Kseniya Ivanova
email: kseniya.ivanova@slf.ch

Andrin Caviezel

Yves Bühler

Perry Bartelt

Numerical modeling of turbulent geophysical flows using a hyperbolic shear shallow water model: Application to powder snow avalanches

Keywords: shear shallow water equations, powder-snow avalanche, ground-based photogrammetry, air-entrainment April 8, 2021

come

Introduction

Shallow water-type equations are commonly used to describe fluid flow in rivers, lakes, coastal areas and oceans [START_REF] Pinder | Finite element simulation in surface and subsurface hydrology[END_REF][START_REF] Tan | Shallow water hydrodynamics[END_REF]. Because they reduce the spatial dimension of a problem by depth-averaging, they offer many computational advantages over full three-dimensional models. Essentially, the long wave approximation allows the reduction in dimension and, because of the hyperbolicity property, the propagation of discontinuities (shocks) can be accurately modelled [START_REF] Savage | The motion of a finite mass of granular material down April 8, 2021, 2:09pm : X -21 a rough incline[END_REF][START_REF] Mageney-Castlenau | Numerical modelling of avalanches based on saint venant equations using a kinetic scheme[END_REF]. Subsequently, shallow water equations are commonly used to simulate geophysical mass movements such as debris flows [START_REF] Pudasaini | A general two-phase debris flow model[END_REF][START_REF] Iverson | A depth-averaged debris-flow model that includes the effects of evolving dilatancy. i. physical basis[END_REF], snow avalanches [START_REF] Naaim | Dry granular flow modelling including erosion and deposition[END_REF][START_REF] Christen | Ramms: Numerical simulation of dense snow avalanches in three-dimensional terrain[END_REF], pyroclastic flows [START_REF] Sheridan | Evaluating titan2d mass-flow model using the 1963 little tahoma peak avalanches, mount rainier, washington[END_REF] and rock/ice avalanches [START_REF] Hungr | A model for the runout analysis of rapid flow slides, debris flows and avalanches[END_REF][START_REF] Bartelt | Thermomechanical modelling of rock avalanches with debris, ice and snow entrainment[END_REF]. When dispersive effects are negligible (this is the case, for example, for the modelling of hydraulic jumps for large Froude numbers, or tsunami waves), one usually employs the classical Saint-Venant (SV) equations [START_REF] Saint-Venant | Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et a l'introduction de marées dans leurs lits[END_REF]. The Saint-Venant equations have also found a wide range of applications in the geosciences, especially in fields where practical engineering solutions are required, see for example [START_REF] Graf | An overview of a decade of applied debris-flow runout modeling in switzerland: challenges and recommendations[END_REF]. However, the question has arisen if the classical shallow water equations can be equally well-applied to model strongly turbulent phenomena, such as powder avalanches (see Fig. S1). The answer is negative, because the underlying hypothesis in the derivation of the Saint-Venant equations is the flow potentiality. Without significant modification, the classical shallow water equations are not able to provide information hypothesis of weakly sheared flows allows the authors to keep the second-order depthaveraged correlations in the governing equations but neglect the third-order correlations, and thus to close the governing system in the dissipationless limit [START_REF] Davidson | Turbulence: An introduction for scientists and engineers[END_REF]. The turbulent energy is separated into two different scales: the turbulent kinetic energy in the boundary layer, and, the turbulent kinetic energy in the roller near the free surface.

In the one-dimensional case, the corresponding equations of shear shallow water model coincide with the equations describing non-isentropic gas flows with a special equation of state. However, in the multi-dimensional case, the system differs significantly from the gas dynamics model. This is a 2D hyperbolic non-conservative system of equations which is reminiscent of a generic Reynolds averaged model of barotropic turbulent flows. The model has three families of characteristics corresponding to the propagation of surface waves, shear waves and average flow.

One area to test the application of depth-averaged turbulence models is the simulation of powder snow avalanches. Powder snow avalanches are strongly turbulent geophysical flows (Fig. S1) that essentially meet the long-wave approximation (the avalanches flow April 8, 2021, 2:09pm X -4 : several kilometers with heights typically less than 100 m). They are composed of a mixture of ice-dust and air and can reach high propagation velocities. Powder avalanches exhibit turbulent flow features, such as velocity fluctuations [START_REF] Ancey | Powder snow avalanches: approximation as non-boussinesq clouds with a richardson number-dependent entrainment function[END_REF][START_REF] Carroll | Frontal dynamics of powder snow avalanches[END_REF][START_REF] Sovilla | The intermittency regions of powder snow avalanches[END_REF] and billows and clefts [START_REF] Simpson | Gravity currents in the environment and laboratory[END_REF] causing significant streamwise height variations in the flow [START_REF] Dreier | Comparison of simulated powder snow avalanches with photogrammetric measurements[END_REF]. Kinetic energy is concentrated at the avalanche front, creating a wind blast, a shock that can easily exert pressures of over 10 kP a, enough to blow down trees [START_REF] Bartelt | Dynamic magnification factors for tree blow-down by powder snow avalanche air blasts[END_REF] and destroy weak buildings [START_REF] Grigoryan | Experimental investigation of an avalanche air blast[END_REF][START_REF] Sukhanov | The mechanism of avalanche air blast formation as derived from field measurements[END_REF]. The tail of the avalanche can be considered a turbulent wake; that is, a suspended cloud of ice dust containing only turbulent energy. Three-dimensional models exist to simulate these avalanches [START_REF] Dutykh | Mathematical modeling of powdersnow avalanche flows[END_REF][START_REF] Sampl | Avalanche simulation with samos[END_REF]), but they are often costly and unpractical to apply, especially for an engineering analysis where different hazard scenarios need to be quickly investigated. Thus, the new model of shear shallow water equations represent an alternative solution to the free surface Navier-Stokes system and lead to a description of the vertical profile of the horizontal velocity while preserving the computational efficiency of the shallow water approximations.

The purpose of this paper is to test if indeed this newly developed turbulence model can be applied to the powder avalanche problem. For this purpose we make use of highresolution photogrammetric powder cloud measurements obtained at the Swiss Vallèe de la Sionne test site in the early 2000s [START_REF] Dreier | Comparison of simulated powder snow avalanches with photogrammetric measurements[END_REF]. Powder-avalanche clouds are ideal for this purpose because the ice-dust serves as a tracer to highlight the cloud turbulence (Fig. S1). The particle tracers define the location of the upper surface of the cloud, a primary variable (h) in the shallow water formulation. For now we restrict our attention to the one-dimensional problem, comparing longitudinal section profiles of measured powder cloud height to model predictions. These heights are governed by turbulent April 8, 2021, 2:09pm : X -5 air-entrainment, and therefore a direct indication of model performance. We concentrate on the prediction of turbulent energy fluxes, and the separation of translational kinetic energy from the energy associated with the fine and large scale turbulent structures. To concentrate the results on the turbulence modelling, we model only the powder cloud, treating the avalanche core as well-defined source of mass and momentum. We believe this simplification is warranted for now, because of the complexity of the combined mixed flowing/powder avalanche problem.

The structure of the article is organized as follows. The governing equations of shear shallow-water theory are presented in the next section with a brief overview of the numerical solution scheme. The following section compares numerical calculations with photogrammetric measurements for two data sets (the avalanches depicted in Fig. S1).

The paper is rounded-off with concluding remarks and an outlook to future work.

Theory, governing equations and numerical solution

Consider a flow regime in which the vertical length scale is much smaller than the horizontal length scale. When the flow is turbulent there is a mean component of motion plus a random component (fluctuations): ũ = U + u , where U is the mean horizontal velocity averaged over the depth in the vertical z-direction U = u = 1 h h 0 u(t, x, z)dz and u is its fluctuation. By depth averaging of incompressible Euler equations of fluid dynamics without assuming potential flow one obtains in one-dimensional case the following system of shear shallow water equations applied here for the avalanche powder-cloud with the source terms representing the balance between gravity and friction as well as the mass exchange between the core of the avalanche, air and the powder-cloud [START_REF] Teshukov | Gas dynamics analogy for vortex free-boundary flows[END_REF][START_REF] Richard | A new model of roll waves: comparison with Brock's experiments[END_REF][START_REF] Richard | Elaboration d'un modèle d'écoulements turbulents en faible profondeur: application au ressaut hydraulique et aux trains de rouleaux[END_REF][START_REF] Ivanova | Formation and coarsening of roll-waves in shear shallow water flows down an inclined rectangular channel[END_REF][START_REF] Gavrilyuk | Multi-dimensional shear shallow water flows: Problems and solutions[END_REF] April 8, 2021, 2:09pm

X -6 : ?, ?):

∂h ∂t + ∂(hU ) ∂x = Q 0 , (1)
∂ ĥ

∂t + ∂( ĥU ) ∂x = Q 0 + Q(U, ϕ, Φ), (2)
∂(hU ) ∂t + ∂ ∂x hU 2 + p = Q 0 W 0 -C f U |U |, (3)
∂ (hE) ∂t + ∂ (hU E + pU ) ∂x = -C e |U | 3 + Q 0 W 0 U + Q 0 gh - U 2 2 + 3 2 (ϕ + Φ)h 2 . ( 4
)
Here h = h(x, t) is the avalanche height without taking in account the air-entrainment, ĥ = ĥ(x, t) is the powder cloud height by taking into account the air-entrainment, t is the time, g is the gravitational acceleration, C f is the Chézy coefficient (friction coefficient) corresponding to the dissipation in the momentum equation.

Let us note u 2 the variance of the horizontal velocity

u 2 = 1 h h 0 u 2 dz
and let us decompose

u 2 = u 2 s + u 2 L .
This decomposition allows to distinguish the turbulent kinetic energy in the boundary layer, and, the turbulent kinetic energy in the eddy near the free surface, and, the definition of two types of enstrophies: Φ = u 2 L h 2 is the enstrophy (squared vorticity) of large eddies formed in the roller and ϕ = u 2 s h 2 is the enstrophy of small eddies developed near the bottom. The wall enstrophy ϕ is supposed to be constant. Here we take

C e = C f + C r Φ ϕ + Φ
. The coefficient C r is the dissipation coefficient in the large scale eddies near the free surface.

April 8, 2021, 2:09pm : X -7

The equation system (1)-(4) represents a system of non-linear conservation laws of mass, momentum and energy with source terms describing the balance between gravity and friction.

The total specific energy E, internal energy e, and the total "pressure" p (we call it "pressure" by analogy with the Euler equations of compressible fluids even if the dimension of p is not the same as that of a pressure) are defined as

E = 1 2 U 2 + e, e = 1 2 gh + (ϕ + Φ)h 2 , (5)
p = gh 2 2 + (ϕ + Φ)h 3 . ( 6
)
The second term in the expression of internal energy e is the sum of the "turbulent" kinetic energy in the roller 1 2 Φh 2 , and "turbulent" kinetic energy in the boundary layer

1 2 ϕh 2 .
This model can be considered as an extension of the standard shallow water model, involving some additional unknowns, the enstrophy (squared vorticity) Φ of large eddies formed in the roller, and an additional parameter, the enstrophy of small vortexes at the vicinity of the bottom, ϕ.

The system (1)-( 4) admits the following equation for enstrophy Φ:

DΦ Dt = 2 h 3 (C f -C e )|U | 3 < 0. ( 7
)
The system (1)-( 4) is a time-dependent system of non-linear hyperbolic partial differential equations with characteristic speeds given by U and U ± a s , where a s is the speed of surface waves. It plays the role of sound speed in this model:

a s = gh + 3(ϕ + Φ)h 2 . ( 8
)
The equation ( 7) means that the enstrophy is decreasing along the trajectories, if C f < C e . Since the large-scale eddies appears at the jump toe and disappears after some April 8, 2021, 2:09pm the depth averaged horizontal velocity u(t, x) and shear stress tensor P in multi-D case are

u(t, x) = 1 h h 0 ũ(t, x, y, z)dz, P = 1 h h 0 (ũ -u) ⊗ (ũ -u) dz.
The tensor P is the stress tensor which measures the distortion of the instantaneous horizontal velocity profile ũ(t, x, y, z) depending of the vertical coordinate z. P is symmetric and positive definite. The sign ⊗ means the tensor product, and I is the identity tensor.

The positive definiteness of P is a consequence of the Cauchy-Schwarz inequality. The enstrophy (squared vorticity) Φ is a one-dimensional analogue of det(P) h 2 .

For numerical solution of the governing equations we use here standard conservative, finite volume Godunov type scheme on a fixed grid. It requires the solution of the Riemann problem at every cell boundary at each time step [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF], [START_REF] Leveque | Numerical methods for conservation laws[END_REF], [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF], [START_REF] Russo | Central schemes for conservation laws with application to shallow water equations[END_REF]. The MUSCL-Hancock extension of the Godunov method is used with the MinMod limiter for the depth, the velocity and the pressure. A mesh convergence study has been performed in order to guarantee a convergent solution.

Let us consider a fixed grid of size ∆x = x i+1/2 -x i-1/2 , the time increment is defined as ∆t = t n+1 -t n that must respect the Courant-Friederichs-Lewy's (CFL) condition.

We rewrite the system (1) -( 4) in the following conservative form

U t + F(U) x = S(U), (9)
where the vectors of conservative variables U, fluxes F(U), and source term S(U) are

U =     h ĥ hU hE     , F(U) =     hU ĥU h + p hU E + pU     , S(U) =       Q 0 Q Q 0 W 0 -C f |U |U -C e |U | 3 + Q 0 W 0 U + Q 0 gh - U 2 2 + 3 2 (ϕ + Φ)h 2       . April 8, 2021, 2:09pm X -10 :
The discrete values of the vector-function U(x, t) at (x i , t n ) will be denoted by

U n i ≡ U(x i , t n ).
The first step (hyperbolic one) consists in computing the source term-free system

U t + F x = 0, (10)
with the initial condition for the complete problem U(x,

t n ) = U n . Integrating in space and time [x i-1/2 , x i+1/2 ] × [t n , t n+1
] the conservation laws (10), one obtains a conservative finite volume Godunov scheme on a fixed grid

Ūn+1 i = U n i - ∆t ∆x F * ,n i+1/2 -F * ,n i-1/2 , ( 11
)
where F * ,n i+1/2 and F * ,n i-1/2 are numerical fluxes. They are constant across interfaces between cells during the time step. For computing the fluxes F * ,n i+1/2 and F * ,n i-1/2 , we solve the Riemann problems between cells i, i + 1 and i -1, i, respectively. The HLLC Riemann solver is used for this aim [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF].

The last step is to integrate the differential equation

dU dt = S(U), (12)
with the initial condition U| t=0 = Ūn+1 given by (11). To integrate the source term, we use Strang splitting [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] :

U n+1 = S 1 2 ∆t M ∆t S 1 2 ∆t U n ,
where M is the MUSCL-Hancock operator, S is the 4 th order explicit Runge-Kutta numerical operator [START_REF] Press | Numerical recipes in Fortran : the art of scientific computing[END_REF].

Numerical simulations and comparison with powder cloud height profiles

In this section we focus on the numerical simulation of the powder cloud avalanche using the proposed system of differential equations Eqs. 1-4. , 2004, 2005). The two avalanches are depicted if Fig. S1 and are denoted No. 509 (2002) and No. 628 (2003). A complete description of the avalanches is contained in the reports, including release conditions and mass balance. Photogrammetric techniques were applied to track the cloud surface over time [START_REF] Vallet | Observations on powder snow avalanches using videogrammetry[END_REF][START_REF] Wicki | Photogrammetrische erfassung vonfliess-und staublawinen mit digitalen amateur-kameras[END_REF]. Previous attempts to simulate these avalanches with shallow water type approaches [START_REF] Bartelt | Configurational energy and the formation of mixed flowing/powder snow and ice avalanches[END_REF] led to all too smooth cloud surfaces, that accurately represented the mean velocity and height of the cloud, but not the turbulent billow and cleft structures [START_REF] Dreier | Comparison of simulated powder snow avalanches with photogrammetric measurements[END_REF].

The governing set of differential equations is specifically set up to be independent of the avalanche core, which is characterized by the source term Q 0 (x, t) for mass influx and velocity W 0 . These right-hand side terms are set up to be analytical functions that are independent of any particular model for the core. Our goal with these assumptions is to concentrate only on the production of turbulence in the powder cloud by explicitly defining source terms for mass and momentum transfer to the cloud. For the mass flux we choose

Q 0 (x, t) = ρ ρ 0 • W 0 (x, t)
Here ρ is the reference density (corresponding to the height h) and ρ 0 is the density of the ice-dust mixture as it exists in the core (corresponding to h 0 ). That is, as the core air is loaded with ice-dust and injected into the cloud, it mixes with some air, which reduces its density. We take these density values as constants (ρ=7 kg/m 3 and ρ 0 =10 kg/m 3 ).

The mean density of the cloud ρ(x, t) (corresponding to the height ĥ(x, t) is given by the relationship ρh = ρĥ . To test the model we first take the injection velocity W 0 (x, t) to be a Gaussian pulse of half-length L 0 and amplitude A 0 ,

W 0 (x, t) = A 0 exp - (x -x 0 -S 0 t) 2 2 • L 2 0 . ( 13
)
April 8, 2021, 2:09pm

X -12

:

The pulse propagates with the speed of the core S 0 . The position x 0 defines the start position of the core at time t=0. For simplicity we assume the amplitude of the pulse to be constant; however, we later modify the amplitude to include an oscillating frequency

W 0 (x, t) = A 0 sin 2 (ω 0 t) exp - (x -x 0 -S 0 t) 2 2 • L 2 0 . ( 14
)
where ω 0 is the frequency. We take the square of the sinus function to ensure mass injection from the core; mass is never extruded from the cloud back into the core. As we shall show the pulse frequency is introduced to model the billow/cleft structure of the observations.

The source term Q(x, t) models air-entrainment and is a simple function of the mean avalanche speed U and turbulent fluctuation velocity u ,

Q(x, t) =    α U U + α T u 2 ρ ρ = α U U + α T (ϕ + Φ)h 2 ρ ρ , if ρ > 2ρ air &U > 0 0, else. ( 15
)
where ρ air is the density of air, ρ air =1.3 kg/m 3 and α U , α T are the steady and turbulent entrainment coefficients, respectively. This function fulfills two properties. Firstly, the higher the mean avalanche speed U and fluctuation velocity u , the larger the airentrainment and, secondly, the denser the ice-dust mixture, the larger the air-entrainment [START_REF] Davidson | Turbulence: An introduction for scientists and engineers[END_REF]. One of the primary reasons for modelling the cloud with turbulence models is to improve air-entrainment models. We note that an air-entrainment models for depth-averaged powder cloud modelling have been proposed by Russian authors, see [START_REF] Bozhinskiy | An energy-based method to calculate streamwise density variations in snow avalanches[END_REF]. These models did not include turbulence.

The computational domain is Ω = [0, 3000]m and is discretized with an uniform Cartesian mesh composed of 3000 grid cells. We have to impose the initial conditions for x belonging to the interval Ω. The initial conditions for height are:

h(x, t = 0) = h 0 = 0.1 m, ĥ(x, t = 0) = h 0 (16)
April 8, 2021, 2:09pm : X -13

and for velocity and large scale enstrophy,

U (x, t = 0) = 0, Φ(x, t = 0) = 0. ( 17
)
We take Neumann boundary conditions at the beginning and end of the computational domain Ω.

We began our numerical investigations with the use of a constant pulse Eq. 13 to simulate avalanche No. 628 (2004), specifically using A 0 = 0.2 and taking the velocity and length of the core to be constant S 0 = 35 m/s and L 0 = 7 m, respectively. This use of a constant pulse, allows us to model (1) the mean height of the cloud and (2) the turbulent wake of the avalanche (Fig. S3). We define the turbulent wake of the avalanche to the region in the avalanche cloud with mostly turbulent energy but little or no translational It is important to note, that we did not change the speed of the core (S 0 = 35 m/s) or the half-length of the pulse (L 0 = 7 m). The amplitude of the pulse increases from A 0 = 0.2 to A 0 = 0.4 such that the total momentum injection does not change from the first simulation (constant pulse) to the oscillating pulse. The free surface profiles predicted by the model ( 1)-( 4) are found to be in good agreement with the photogrammetric measurements.

Although the exact location of the billows is sometimes offset, the surface undulations April 8, 2021, 2:09pm The model predicts high translational kinetic energies at the avalanche front, near the source of momentum Q 0 (x, t). Turbulent energy is also created at the avalanche front, an indication that the injection of momentum is strongly linked to the production of turbulence (velocity fluctuations). Immediately behind the front there is a strong decrease in kinetic energy, as well as a strong increase in turbulent energy. In fact, at the tail of the powder avalanche, all energy is turbulent. Thus, the model clearly divides the powder avalanche into a fast-moving front followed by a turbulent wake. This large difference in translational velocity between the front and tail of the avalanche is well documented [START_REF] Dreier | Comparison of simulated powder snow avalanches with photogrammetric measurements[END_REF]. Note that the source Q 0 (x, t) can disappear, but the front continues to propagate with speed. How the front velocity decreases is dependent on the drag coefficient C f .

The parameterization of the air-entrainment is also of interest. Here, we assume that the air-entrainment is more strongly dependent on the sum of the small and large scale enstrophies (ϕ + Φ) than the translational velocity U , see Eq. 15, α T >> α U . The model predicts that turbulent energy is concentrated at specific locations in the avalanche tail.

Because more air is entrained at these locations, the height ĥ of the cloud grows. This leads to the strong surface oscillations observed in the experiments.

To simulate both avalanches we made only two changes to the parameters. In comparison to avalanche No. 509, the amplitude of the source Q 0 (x, t) for avalanche No.

628 decreased (from A 0 =0.5S 0 to A 0 =0.4S 0 ) while, at the same time, the small scale wall enstrophy increased (from ϕ=2 to ϕ = 10). The half-length L 0 changed slightly April 8, 2021, 2:09pm : X -15 from L 0 =17.0m (No. 509) to L 0 =14.0m (No. 628). All other model parameters did not change. The unchanged parameters control the air-entrainment (α U , α T ), as well as the production of turbulent energy from shearing (C e ), and its eventual decay (the dissipation coefficient C r ). In our view, these parameters should indeed be constants, not varying from avalanche to avalanche, but are intrinsic properties describing the behavior ice-dust mixture. It appears that the injection of momentum into the cloud will vary from avalanche to avalanche. Both the strength (amplitude A 0 , length L 0 ) of the injection and production of small scale turbulence (ϕ) may change. Because we consider avalanches on the same track, it is not possible to identify how terrain features (steepness, curvature, roughness) influence these parameters, but they could certainly be influenced by avalanche mass, snow entrainment, or the temperature of the snow cover. More important, however, is the fact that a pulse with an oscillating amplitude is required to drive the turbulent powder cloud model. There are, in fact, models that would predict such behavior, treating the avalanche core as a type of parametric oscillator with periodic injections of mass and momentum into the powder cloud, see (Buser & Bartelt, 2015;[START_REF] Bartelt | Configurational energy and the formation of mixed flowing/powder snow and ice avalanches[END_REF].

Understanding that there is a random component to the exact location of the height oscillations, we used a Fast Fourier Transform to convert the spatial data of both the experiments and numerical results into a wavelength spectra describing the billow length.

The key idea of Fourier's theory is that any periodic function, however complex it is along the period, can be exactly (i.e. with no information loss) represented as a weighted sum of simple sinusoids. Spectra are often computed from gridded data to determine the horizontal-scale dependence of quantities such as kinetic and turbulent energies, vertical and horizontal velocity, or flow height. On the Figure S9 one can see the amplitude spectrum as a function of the wavelength and the comparison between the calculated (black) and measured (red) data. In order to get better comparability we normalized the April 8, 2021, 2:09pm amplitude to the highest amplitude spectrum value within the plotting domain. A good agreement between the simulations and the measurements wavelength is apparent.

This result indicates that modelling the mass and momentum exchange between the core and the cloud as an oscillating pulse with known frequency (Eq. 14) is too simple. In reality the exchange is driven by other factors such as terrain features which would induce a different spatial distribution. Of significance, however, is the fact that the distribution of billow lengths is fairly well represented by the model. This fact is a direct consequence of constructing the air-entrainment model to be a function of the turbulence variables (ϕ, Φ). jets (injection) and air entrainment [START_REF] Davidson | Turbulence: An introduction for scientists and engineers[END_REF]. It was clearly not possible to model the structure of the powder cloud with a steady source of momentum. Instead, it was necessary to select a specific injection frequency, with periods varying between T =2s to T =3s. We can only speculate on the physical source of these oscillations, which must arise from the mechanics of the granular core. Our future work will certainly be oriented to the development of a full multi-dimensional model capable to describe multi-dimensional powder-snow avalanche dynamics. The numerical validations show the capability of the shear shallow water model to reproduce the complex non-stationary avalanche dynamics. April 8, 2021, 2:09pm

Conclusion

.

  to the roller length), the large-scale enstrophy created by the shock must decrease thereafter. Since the equations are reminiscent of the Euler equations of compressible flows, the conservation laws imply standard Rankine-Hugoniot relations.At the shock front, the enstrophy is increasing analogously to the entropy increase for the Euler equations of compressible flows. The enstrophy production through shocks corresponds physically to the vortex formation. Then the enstrophy dissipates over the length of the roller according to (7). When one takes ϕ = 0, Φ = 0, and C f = C e , the system is reduced to the classical Saint-Venant (shallow water) equations.The shear shallow water model complemented by friction terms provides a more realistic description of some flow phenomena, such as the turbulent hydraulic jumps, travelling waves down inclined plane, shock instabilities in SWASI experiment or a turbulent powdersnow avalanches.As we mentioned above, the underlying hypothesis in the derivation of the Saint-Venant equations is the flow potentiality. The horizontal vorticity (parallel to the bottom) in shallow water approximation is related with the horizontal velocity shear: ω || ≈ ũz , where ũ is the instantaneous (non-averaged) horizontal velocity, and the index z means the derivative in the vertical direction. The absence of the vorticity means the absence of the horizontal velocity shear. In 2D the vorticity vector in shear shallow water approximation is ω = -For the one-dimensional problem we study here, we have ω = ∂u ∂z .The corresponding multi-dimensional model of shear shallow water flows is a hyperbolic system of equations which is reminiscent of a generic Reynolds-averaged Euler equations for barotropic compressible turbulent flows. The multi-dimensional case is much more challenging because one important subset of evolution equations is nonconservative and the nonconservative products also act across genuinely nonlinear fields. The definitions of April 8, 2021, 2:09pm : X -9

  energy. The constant pulse model does not reproduce the billow/cleft structure observed in the experimental measurements. Subsequently, we simulated both avalanches No. 509 (2003) and No. 628 (2004) using an oscillating pulse, with circular frequency ω 0 = 1.3 rad/s. The pulses inject momentum into the cloud at a natural frequency of f = 0.4 Hz; that is, with a period of approximately T = 2.5 s . The results are depicted in figs. S4-S5 for avalanche No. 509 and figs. S6-S8 for avalanche No. 628. Here we present comparisons of profile height for different times t as well as the distribution of kinetic and turbulent energies from the avalanche front to tail.

  to the measurements. The height of the billows, as well as the deepness of the clefts (e.g. Fig. S8, S5), are well represented. Another interesting model result is the calculated steepness of the cloud front, which is reproduced by the model (Fig. S8, S5).

A

  photogrammetric measurements. These measurements show highly dynamic surface un-

Figure S3 .Figure S4 .

 S3S4 Figure S3. Avalanche No. 628, constant amplitude of pulse function W 0 , time t = 55 s.

Acknowledgement Kseniya Ivanova thank Ilya Peshkov for useful discussion. This work was financially supported by RAMMS.