Numerical modeling of slope flows entraining bottom material

COLTEC-01993; No of Pages 10 Cold Regions Science and Technology xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Cold Regions Science and Technology journal homepage: www.elsevier.com/locate/coldregions

Numerical modeling of slope flows entraining bottom material M.E. Eglit a,⁎, A.E. Yakubenko b a b

Lomonosov Moscow State University, Moscow, Russia Institute of Mechanics, Lomonosov Moscow State University, Moscow, Russia

a r t i c l e

i n f o

Article history: Received 15 December 2013 Received in revised 24 May 2014 Accepted 4 July 2014 Available online xxxx Keywords: Slope flows Entrainment of bottom material Laminar and turbulent flows Bingham and Hrschel-Bulkley fluids

a b s t r a c t Entrainment of bottom material has a substantial influence on the dynamics of slope flows such as snow avalanches and debris flows, and should therefore be taken into account by mathematical models. An approach based on the assumption that the entrainment takes place when the bed shear stress is equal to the shear strength of the bottom material was proposed by Issler and Pastor Pérez (2011). In this paper, this approach is applied to model laminar flows of Bingham and Herschel–Bulkley fluids as well as turbulent flows with entrainment. Full (not depthaveraged) equations are used to calculate the flow parameters including the shear stress at the bottom, the velocity profile and the flow depth variation with time due to entrainment. The paper focuses on studying of qualitative features of motion with basal entrainment. Modeling of real flows from start to stop is not considered here. It is shown that in entraining flows moving down long constant slopes the flow velocity increases. The velocity profiles in laminar flows described by Bingham and Herschel–Bulkley models at large time become almost linear except in a relatively narrow layer near the upper flow surface and the entrainment rate tends to a constant value. This value depends on the coefficients describing the flow material properties, the shear strength of the bottom material, and the slope angle. Asymptotic formulae for the entrainment rate and the shear rate in the linear part of the velocity profile are derived analytically and confirmed by numerical calculations. © 2014 Elsevier B.V. All rights reserved.

1. Introduction It is well known that slope flows (such as snow avalanches or debris flows) typically entrain slope material during motion so that their mass increases substantially (Sovilla et al., 2001). The dynamics of slope flows is substantially influenced by the entrainment of bed material. Besides, the depth of the layer that could be entrained by flows may be of great interest itself, for example if pipes or other objects crossing the flow path are buried in the ground to protect them. Various mathematical models for flows with entrainment have been proposed depending on different assumptions about the entrainment mechanism and the properties of the flow and bottom materials (see, e.g., Barbolini et al., 2003; Bouchut et al., 2008; Bovet et al., 2010; Briukhanov et al.,1967; Crosta et al., 2009; Eglit, 1968, 1983; Eglit and Demidov, 2005; Fischer et al., 2012; Gauer and Issler, 2004; Issler and Pastor Pérez, 2011; Naaim et al., 2003, 2004; Sovilla, 2004; Sovilla and Bartelt, 2002; Sovilla et al., 2006, 2007; Tai and Kuo, 2008). The majority of these models are based on the depth-averaged flow equations as well as on the assumption that the entrainment rate is a known function of the flow parameters, e.g., the depth-averaged velocity, or the flow depth. One of the first models for motion of snow avalanches with entrainment was presented by Briukhanov et al. (1967) and by Eglit (1968) ⁎ Corresponding author at: Faculty of Mechanics and Mathematics, Lomonosov Moscow State University, Moscow 119991, Russia. E-mail address: [email protected] (M.E. Eglit).

(see also Eglit, 1983; Eglit and Demidov, 2005). It was assumed in this model that entrainment takes place in a narrow zone adjacent to the avalanche front. This assumption was based on observation data confirmed now by measurements (Sovilla, 2004). Since the length of this zone is much less, than the avalanche length, it may be modeled as a “destruction front,” where the snow cover undergoes the transition from solid to fluid while it is set in motion. This type of entrainment is now referred to as frontal or plug entrainment. An important assumption is that destruction of the snow cover and its entrainment into the flow occur when the stress exerted on it is equal to its compressive strength. This condition together with the equations of mass and momentum conservation across the front and the depth-averaged momentum and continuity differential equations allowed calculating the flow dynamics including the front speed and height. To model continuous entrainment from the bottom, empirical formulae for the entrainment rate have been proposed by analogy with flows of different physical nature (see Eglit and Demidov, 2005 for a review). However, these formulae are not derived from physical laws and the relations between the values of their coefficients and the properties of the moving material and substrate remain unknown. Issler and Pastor Pérez (2011) proposed a model based on the assumption that the entrainment from the bottom takes place when the value of the shear stress at the flow bed, τb, is equal to the shear strength of the bottom material, τc. The shear stress at the bed and the entrainment rate in their model should be calculated by solving the full (not depth-averaged) flow equations. Issler and Pastor Pérez (2011)

http://dx.doi.org/10.1016/j.coldregions.2014.07.002 0165-232X/© 2014 Elsevier B.V. All rights reserved.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

2

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

presented numerical solutions for 2-dimensional laminar unsteady motion of linearly viscous and Bagnold fluids down an infinitely long slope with constant incline. In this paper, the entrainment model proposed by Issler and Pastor Pérez (2011) is applied to investigate entraining unsteady laminar flows of Bingham and Herschel–Bulkley fluids as well as turbulent flows. The Prandtl turbulence model and the three-parameter turbulence model developed by Lushchik et al. (1978, 2004), hereinafter referred to as LPY model, are used to describe turbulent flows. It should be emphasized that we do not consider here the full problem of the flow motion from start to stop but focus only on some aspects of continuous entrainment along the bottom of the flow. The outline of the paper is as follows. In Section 2, the general equations for 2-dimensional unsteady motion down an infinitely long slope with constant incline are given with a focus on the Bingham rheology. The results from numerical simulations based on the theory presented in Section 2 are described in Section 3. Section 4 presents analytical asymptotic solution when the entrainment rate becomes constant at larger times. In Section 5 flows of Herschel–Bulkley fluids are considered. Entrainment by turbulent flows is studied in Sections 6 and 7. In the Conclusions, the qualitative predictions of theoretical models are summarized.

entrainment will occur (except in limiting cases where the gravitational force and the inertial pseudo-force combined balance τc and the entrainment rate smoothly tends to 0 while |τb| remains locked at τc). In any case |τb| cannot be larger than τc because at this stress level, the texture of the bed material is destroyed immediately and any excess shear stress from the flow is absorbed by accelerating the eroded material. The boundary between the flow and the intact slope material is a front where the solid bottom material is ground up into a (granular) fluid and set in motion. If entrainment takes place, this destruction front moves down along the normal to the slope. The propagation speed of this front is the entrainment rate. It should be found by solving the full problem. The momentum equation reads ∂vx 1 ∂τ xz ¼ g sin θ þ ρ ∂z ∂t

at 0 ≤ z ≤ hðt Þ

ð1Þ

where τxz is the shear stress. For all flows in our study, we will neglect friction at the upper surface z = 0, τxz ¼ 0

at z ¼ 0;

ð2Þ

assume the no-slip condition at the bottom z = h(t), 2. Statement of the problem

vx ¼ 0

The aim of this paper is to model the process of basal entrainment and to study its influence on the flow dynamic parameters. That is why we consider the simplest statement of the problem. We do not take into account many features of real flows: complex slope topography, existence of the leading front, variation of flow depth and velocity along the flow. Instead, we consider only a part of the flow assuming that longitudinal variation of flow depth, flow velocity, inclination angle of the slope and properties of the bottom material in this part can be neglected. It is reasonable because most natural slope flows are long and the derivatives of the flow depth and velocity over the downstream coordinate are small except possibly in separated narrow zones, which are not considered here. The strict statement of the problem is as follows. An unsteady twodimensional laminar or turbulent flow moving down a slope with constant incline and constant bed material properties is considered. The moving material is incompressible. The flow depth h is constant along the flow but may vary in time due to entrainment of bottom material, h = h(t) . The x-axis is parallel to the flow velocity and is located at the upper surface of the flow; the z-axis is directed downward along the normal to the bottom (Fig. 1). The flow is laminar, vz = 0, vx = v(t, z). In turbulent flows (Sections 6 and 7), similar relations are assumed for Reynolds-averaged velocity components. If the value of the shear stress, |τb|, at the bottom of the flow is less than the shear strength τc of the bottom material, the flow does not entrain the bottom material. If |τb| = τc,

and an additional boundary condition at the bottom for the entraining flow

Fig. 1. Schematic representation of a Bingham flow on an inclined plane: 1 – the solid layer, 2 – the shear layer; 3 – the bottom material; AB – the destruction front; h(t) – the flow depth; θ – the slope angle; r(t) – the solid layer depth.

at z ¼ hðt Þ;

jτxz j ¼ τ c at z ¼ hðt Þ:

ð3Þ

ð4Þ

The condition (4) will allow us to find the rate of change of the flow depth, i.e., the entrainment speed. We use here the classical no-slip condition (3) at the bottom of the flow. Solutions of many problems about motion of viscous fluids obtained with the use of this boundary condition are in good agreement with experimental data. It was used also in modeling natural flows (Dent and Lang, 1980, 1982; Issler and Pastor Pérez, 2011). However, for some flows the no-slip condition is not suitable, e.g., for motion of rarefied gas or motion of concentrated suspensions (Soltani and Yilmazer, 1998). It is a topic for discussion whether the no-slip condition should be replaced by a certain “slip” condition in modeling geophysical flows. The slip condition states that there is a jump of velocity at the bottom: tangential velocity equals zero in the bottom material and equals vslip at the bottom of the flow. Measurements of velocity profiles (Kern et al., 2004; Naaim et al., 2004) show that the flow velocity is not small already at short distance above the bottom. However, these results do not reject the applicability of the no-slip condition. The slip condition was used by Bovet et al. (2010), Kern et al. (2004), Lang and Dent (1980), and Naaim et al. (2004). Unfortunately, the data about vslip are very poor. That is why in this study we use the classical no-slip condition, though calculations could be easily done with a slip-condition if we knew the value of vslip. To close the systems (1)–(4) we need to formulate the constitutive rheological relations for the flow material. The latter are different for flows of different physical nature: dry and wet dense snow avalanches, powder avalanches, mudflows, debris flows, volcanic lava, landslides, water flows etc. The simplest model for dense snow avalanches is the model of incompressible isotropic linear-viscous fluids with high value of the viscosity coefficient (Bovet et al., 2010; Dent and Lang, 1980; Issler and Pastor Pérez, 2011; Lang and Dent, 1980; Nishimura, 1996). The components τij of the tensor of viscous stresses in this model are linear functions of the components of the tensor of strain rate, eij, and vanish at eij = 0: τij = 2µeij, µ being the viscosity coefficient. For the simple shear flow in our study, the latter relations are reduced to the equation τxz ¼ μ

∂vx ∂vx ; ¼ 2exz ; ∂z ∂z

ð5Þ

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

3

where τxz is the shear stress and exz is the shear strain rate tensors components. The disadvantage of this model is that it cannot describe the stopping of an avalanche on a slope. This can be described by introducing additional stresses that do not vanish at v = 0. The simplest such model is the so-called Bingham fluid. The Bingham fluid behaves as a solid if the values of the tangential stresses are less than the so-called yield stress τ0; otherwise it moves as a linear viscous fluid. For the simple shear flow considered in our study the rheological equation for a Bingham fluid takes the form If jτ xz j≤τ0 then




∂v
∂vx ¼ 0; if jτxz jNτ0 then jτxz j ¼ τ 0 þ μ 0

x

: ∂z ∂z

ð6Þ

Here τ0 is the yield stress, μ 0 is the viscosity coefficient. Note that if τ0 = 0, the fluid reduces to a Newtonian linear-viscous fluid. Eq. (6) can be rewritten in the form similar to Eq. (5) with the viscosity µ depending on the strain rate

τxz

8 > > > <∞

∂v ¼ μ x ;μ ¼ τ0 > ∂z > > : μ 0 þ

∂v =∂z

x

if if




∂v
jτxz j b τ 0 ;

x

¼ 0
∂z

∂v
jτ xz j≥τ0 ;

x

≥0 ∂z

ð7Þ

A generalization of the Bingham fluid is the so-called Herschel– Bulkley fluid with the following rheological equation (again for the simple shear flow we consider): If jτ xz j≤τ0 then



α
∂v
∂vx ¼ 0; if jτxz jNτ0 then jτxz j ¼ τ 0 þ K

x

; ∂z ∂z

ð8Þ

where K, α are coefficients called the consistency and the power–law index respectively. Bingham and Herschel–Bulkley models have been considered as possible rheological models for snow avalanches, e.g., by Bovet et al. (2010), Dent and Lang (1982), Issler (2003), and Kern et al. (2004) and for debris flows by Coussot et al. (1998). All attempts to formulate rheological relations for different kinds of flows are based on measurements of velocity profiles mainly in steady flows (Dent and Lang, 1982; Dent et al., 1998; Naaim et al., 2004; Nishimura and Maeno, 1989; Tiefenbacher and Kern, 2004). However, in practice different models may explain the same measured profile. That is why many different rheological models have been proposed for slope flows (see, e.g., Bovet et al., 2010; Coussot et al., 1998; Dent and Lang, 1983; Issler and Pastor Pérez, 2011; Kern et al., 2004; Maeno et al., 1980; Naaim et al., 2003, 2004; Nishimura, 1996; Norem et al., 1989; Rognon et al., 2007, and references therein). Below we study the entraining flow of a Bingham fluid. Since according to condition (3) τ = 0 at the upper surface, there exists a layer 0 ≤ z ≤ r(t) adjacent to the flow surface where |τ| ≤ τ0. This layer moves without deformation as a solid (Fig. 1). We will call it a solid layer or a plug layer. The depth r(t) of the solid layer is not known in advance. It can be calculated from the condition |τ| = τ0 at z = r(t). In our calculations, we approximated the standard constitutive Eqs. (6) and (7) for the Bingham fluid by the equations for a fluid with nonlinear viscosity (Fig. 2). We suppose that the viscosity coefficient μ depends on the second invariant of the strain rate tensor so that in our simple problem

μ¼

8 > > > < τ 0 =ε τ0 −μ 0 ε > > > : μ 0 þ

∂v =∂z

x

Note that jτxz j ¼ τ 0

if if




∂vx


bε


∂z

∂vx


≥ε
∂z





x
at
∂v
¼ ε and jτ xz j≤τ0 ∂z

ð9Þ




x
at
∂v
≤ε. The pa∂z

rameter ε is chosen to be small (0.001 s− 1 in our calculations).

Fig. 2. Dashed line – standard Bingham fluid, Eq. (7); full line – approximation (9).

Henceforth we will refer to the layer where |τxz| ≤ τ0 as a solid layer since the shear rates are small in it. At given values of τ0, τc, ρ, θ and given initial flow depth and velocity, we should find the variations of the velocity profile with time, as well as r(t) and h(t). The function h(t), and therefore the domain where the solution is to be constructed, is unknown in advance. Relation (4) should be used to calculate the variation of the flow depth h(t) due to bottom material entrainment. Note that Eq. (4) does not contain h(t) explicitly. For this reason we reformulated the problem in terms of a new coordinate z ¼ hðzt Þ. Then additional terms containing h(t) and its time derivative appear in Eq. (1) and in the boundary condition (4), see Eqs. (12)–(14) below. All calculations have been done in dimensionless variables, designated by bars over the letters and defined by the following relations vx ρ μ ε t τ h τ z ; ρ ¼ ; μ ¼ ; ε ¼ ; t ¼ ; τ 0 ¼ 0 ; h ¼ ; τ xz ¼ xz ; z ¼ : v τc τc ρ0 μ0 ε t z hðt Þ

vx ¼

ð10Þ Here ρ0 is the fluid density and v ¼

rffiffiffiffiffiffi τc v2 v v v z μ ; z ¼  ; t  ¼  ; ε ¼  ; Re ¼   ; ν 0 ¼ 0 : ρ0 g g z ν0 ρ0

ð11Þ

The equations and boundary conditions (1)–(4), (9) in new variables are (henceforth bars over the letters are omitted)   ∂vx ∂vx z dh 1 ∂ ∂v þ sin θ þ 2 ¼ at 0 ≤ z ≤ 1 μ x ∂t ∂z h dt ∂z ρh Re ∂z

μ¼

8 > > > < τ 0 Re =ε

if

τ 0 Re −μ 0 ε > > > : 1 þ

∂v =∂z



if

x

∂vx ¼ 0 at z ¼ 0; vx ¼ 0 ∂z

μ

∂vx

¼ h at z ¼ 1 Re
∂z





∂vx


b εh


∂z

∂vx


≥εh
∂z


at z ¼ 1;

ð12Þ

ð13Þ

ð14Þ

Relation (14) is the new form of the entrainment condition (4). Introduction of the coordinate z ¼ hðzt Þ has two advantages: first, we now have a problem in a fixed known domain and second, the condition (14) now contains the unknown flow depth explicitly and can be used to calculate it.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

4

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

3. Simulations of laminar Bingham flows Software has been created to solve the Eqs. (12)–(14) and a series of simulations have been run. To make calculations we need to set the values of the viscosity coefficient, µ0, yield stress, τ0, and the shear strength of the bed material, τc. The values of these parameters depend on the physical nature of the flow and on the properties of the bottom material. We tried to make calculations having in mind wet and dry dense snow avalanches. Unfortunately, the available data about possible values of µ0, τ0, and τc are poor and sometimes contradictory. We chose the orders of values of these parameters following the information from Dent and Lang (1980, 1982, 1983), Issler (2003), Issler and Pastor Pérez (2011), Kern et al. (2004), Lang and Dent (1980), Maeno and Nishimura (1979), Maeno et al. (1980), Nishimura (1996), Nishimura and Maeno (1989), Sovilla et al. (2006, 2007), and Tiefenbacher and Kern (2004). It is worth noting that, e.g., the values of kinematic viscosity ν = µ/ρ proposed by Lang and Dent (1980) are about 0.015 m2 s−1 for a small-scale experimental avalanche and 0.5 m2 s−1, i.e., 33 times larger, for a natural avalanche with 1 m depth. The authors even speak about dependence of kinematic viscosity on the flow depth. In our opinion, a more natural explanation of their result is that large-scale flows are turbulent, the turbulent viscosity being much larger than the laminar viscosity. To demonstrate the features of the model we do not limit the depth of the entrainable layer. Of course, in real flows the entrainment ceases when there is no entrainable material left in the slope. Some examples of simulations results are presented below. In these simulations it is assumed that at t = 0 the flow is stationary but its depth is a little larger than the depth of the stationary flow with shear stress τc at the bottom. Then the entrainment starts immediately. The basic system of the other input parameters is 

2 −2

θ ¼ 30 ; τ0 =ρ ¼ 2 m s hð0Þ ¼ 0:5 m:

2 −2

; τc =ρ ¼ 2:3 m s

2 −1

; ν ¼ μ=ρ ¼ 0:005 m s

;

ð15Þ

The value of τ0/ρ is taken following Tiefenbacher and Kern (2004), Kern et al. (2004), who performed experiments with wet snow with τ0 = 740 ± 100 Pa , ρ = 400 ± 50 kg/m3, and Dent and Lang (1980), who found τ0 = 540 Pa , ρ = 250 − 300 kg/m3. The value of the kinematic viscosity is chosen following Maeno et al. (1980), who suggested ν ~ 0.001 m2 s− 1 for fluidised snow, and Dent and Lang (1982), who estimated ν as 0.004 m2 s−1 . The order of τc/ρ is in accordance with the data from Shapiro et al. (1997), Mellor (1975), Sovilla et al. (2007), Voitkovsky (1975). However, it should be emphasized again that in this paper we do not simulate real avalanches, but focus on studying qualitative behavior of the flow if the entrainment is determined by condition (4), and, besides this, on studying the influence of various input parameters. Results obtained using the input parameters (15) are presented in Figs. 3–6. Fig. 3 shows that the velocity at the flow surface indicated as vmax and the depth-averaged velocity vmean increase with time due to entrainment. Velocity profiles at different instants of motion with entrainment are plotted in Fig. 4 . There is a “solid” or “plug” layer adjacent to the upper surface of the flow; the thickness of this layer increases with time, but more slowly than the full flow depth does, so that the ratio of the plug layer depth to the flow depth tends to zero (Figs. 5 and 6). The velocity profile in laminar entraining flow is almost linear near the destruction front. At large time the profile tends to be linear everywhere except in a relatively narrow zone near the upper surface. Fig. 5 shows that asymptotically, as t increases, the flow depth becomes a linear function of t, so the entrainment rate tends to a constant. For a

Fig. 3. Surface velocity vmax (curve 1), and depth-averaged velocity vmean (curve 2) versus time for a Bingham fluid.

Newtonian fluid, these facts were found by Issler and Pastor Pérez (2011) and they were confirmed by our calculations. Fig. 7 shows that the entrainment rate depends on the value of the flow material yield stress. Lower values of τ0/ρ lead to greater entrainment rates. Fig. 8 demonstrates the influence of the kinematic viscosity coefficient value on the value of the entrainment rate. The entrainment is faster if the viscosity is higher. 4. Asymptotic formulae for the entrainment rate and for the velocity profile. Flow of Bingham fluid It is possible to construct an asymptotic solution for large values of t in the domain just behind the destruction front. Assume that the entrainment rate we ¼ dh is constant and the motion is steady in the coordt dinate system linked with the bottom of the flow. Let us introduce a new coordinate z′ = h(t) − z, z′ being the distance from the bottom. In this coordinate system vz' = const ≠ 0. It follows from the boundary condition at the bottom that vz0 ¼ we . The momentum Eq. (1) with viscosity coefficient (7) or (9) takes the form we

dvx d2 vx 0 ¼ g sin θ þ μ 0 dz dz02 x The other boundary conditions at the bottom are vx = 0, τ 0 þ μ 0 dv ¼ dz0

x τc for the standard Bingham fluid and τ0 −μ 0 ε þ μ 0 dv ¼ τc for its apdz0 proximation (9). The solution to the momentum equation that does

θ 0 g sin θ not grow exponentially with z′ is vx ¼ g sin we z , i.e. we ¼ dvx Taking dz0

Fig. 4. Velocity profiles in an entraining flow of a Bingham fluid at different instances.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

5

Fig. 5. Flow depth (curve 1) and solid (plug) layer depth (curve 2) versus time.

Fig. 7. Influence of the yield stress value: 1 – τ0/ρ = 2 m2 s−2; 2 – τ0/ρ = 1.5 m2 s−2.

into account the relations between the shear stress and the shear rate and the boundary conditions, we obtain the formulae ð16Þ

square root of time. The same isptrue ffiffiffiffiffi for the depth of the plug layer r(t): at large values of time, r ðt Þ  c νt where c is a constant, that depends on the physical properties of the flow and bottom material but does not depend on the bottom incline (Figs. 5, 6, and 11). These results follow from the asymptotic solution to the total problem that was constructed by use of dimensional analysis in discussions with A.G. Kulikovsky. We do not include the details of this solution in this paper.

ð17Þ

5. Flow of Herschel–Bulkley fluid

vx ¼

τc −τ0 0 z; μ0

we ¼

dh μ 0 g sin θ ¼ dt τc −τ0

for the standard Bingham fluid and vx ¼

τc −τ0 þ μ 0 ε 0 z; μ0

we ¼

μ 0 g sin θ τc −τ0 þ μ 0 ε

for its approximation, respectively. Calculations confirm this behavior of slope flows (Figs. 4, 5, and 7–10). The value of the entrainment rates at large time are close to the values given by formula (16) which are 0.083 m s−1 for input parameters (15) and 0.019 ms−1 if τ0/ρ = 1 m2 s−2. In fact, our calculations gave slightly overestimated values; this overestimation was smaller for smaller values of τ0 (Figs. 9 and 10). Formulae (16) and (17) show that for large values of time the velocity profile is linear except in the relatively thin layer adjacent to the upper surface, and the value of the shear rate in the linear part of the 0 flow is τc −τ μ0 : So it does not depend on the slope angle, depending only on the physical properties of the flow and bottom material. If the flow and bottom parameters are given by relations (15) then it is approximately equal to 60 m s−1. The asymptotic value of the entrainment rate is proportional to the fluid viscosity coefficient and to the slope angle. These effects are observed also in turbulent flows, see Fig. 16. The acceleration of all particles in the zone of linear velocity profile is equal to g sin θ. While the total flow depth increases with time almost linearly due to entrainment (Figs. 5, 7, and 8), the depth of the upper layer, where the velocity profile differs from linear, is approximately proportional to the

Fig. 6. Ratio of the solid layer depth to the flow depth vs. time.

We considered also the flows obeying the Herschel–Bulkley model (8) with the power–law index α equal to 2 (Kern et al., 2004). The behavior of an entraining laminar flow is qualitatively similar to the behavior of laminar Bingham fluid flows. The fluid behaves as a solid in a layer adjacent to the surface (Fig. 12). At large time the flow depth grows approximately proportionally to time while the velocity profiles tend to be linear in the lower layer (Figs. 12 and 13). If α = 2, an asymptotic formula for the entrainment rate at large time is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K we ¼ g sin θ: τc −τ0

ð18Þ

Comparison of calculation results for Bingham and Herschel–Bulkley fluids is presented in Figs. 12–14. The values of model coefficients for the Bingham fluid are as in (15). For the Herschel–Bulkley fluid (see Eq. (8)) they are chosen following Tiefenbacher and Kern (2004): 

2

2

2

2

−4

θ ¼ 30 ; τ0 =ρ ¼ 2 m =s ; τc =ρ ¼ 2:3 m =s ; K=ρ ¼ 0:889  10 hð0Þ ¼ 0:5 m; α ¼ 2:

2

m ; ð19Þ

Fig. 8. Influence of the value of the coefficient of kinematic viscosity: 1 – ν = 0.005 m2 s−1; 2 – ν = 0.01 m2 s−1.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

6

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

Fig. 9. Entrainment rate vs. time for input parameters (15).

Note that we use the no-slip condition at the lower boundary of the flow in our calculations. Figs. 12–14 show that the Bingham model and the Herschel–Bulkley model with α = 2 predict almost identical behavior of the entraining flows described by the sets of input parameters (15) and (19), respectively. This holds at least for the first 10 s from the start of entrainment. Moreover, the asymptotic values of the entrainment rate for both models are very close to each other: 0.083 m s− 1 for the Bingham model and 0.086 m s−1 for the Herschel–Bulkley model. This is a rather unexpected result because we just took the input parameters (15) an (19) from the literature and did not choose them specially to obtain such closeness.

6. Turbulent flows – Prandtl turbulence theory

Fig. 11. Ratio of the plug layer depth to parameters are as in (15).

pffiffiffiffiffi νt vs. time: 1 – θ = 30°; 2 – θ = 45°. The other

viscosity values found in small-scale experiments often lead to unrealistically high values of the velocity. In this section, we consider a turbulent entraining flow of a linearly viscous fluid. Assume that the averaged flow parameters satisfy the conditions and equations indicated in Section 2 but the shear stress is now the sum of the averaged viscous and turbulent (Reynolds) stresses. Note that the term “averaged” in this context does not mean depth-averaged, but time-averaged or ensemble-averaged, as is usual in turbulence theories. Below this kind of averaging is referred to as Reynolds averaging. Turbulent stresses depend on the intensity of turbulence, which in turn depends not only on properties of the moving medium but also on the external conditions. That is why different semi-empirical theories for turbulent stresses have been proposed for different classes of turbulent flows (Davidson, 2004). Prandtl (1925) proposed one of the first and the simplest theories for flows above a solid surface. Measurements of velocity profiles in the wind above the Earth surface, and in pipe flows are in accordance with the predictions of the Prandtl theory. In this section, we use the Prandtl turbulence theory. The averaged momentum equation, called the Reynolds equation, formally coincides with (1), if the velocity is replaced by the Reynolds-averaged velocity. The shear stress τxz is defined in the Prandtl theory as follows:

It is known that, as a general rule, open channel flows are laminar if the Reynolds number Re ¼ vh ν is less than 500 (French, 1985). Let, e.g., v= 25 m s− 1, h = 2 m; then the flow is laminar if ν N 0.1 m2 s− 1. Most natural flows are turbulent. In particular, the values of the moving snow viscosity evaluated from experimental data are much smaller than 0.1 m2 s− 1 (Kern et al., 2004; Lang and Dent, 1980; Dent and Lang,1983; Maeno et al., 1980; Nishimura, 1996). The turbulent nature of slope flows is confirmed by the fact that many of them are successfully described by models with viscous drag proportional to velocity squared (Voellmy-type models). If they were laminar, this drag would depend on velocity linearly. Moreover, attempts to compute the velocity of a stationary fullscale flow under the assumptions that the flow is laminar and with

Here μ 0 is the fluid viscosity coefficient, τxzðturbÞ is the turbulent stress, l characterizes the distance at which the flow particles jump due to turbulence, and κ is the empirical constant approximately equal to 0.4 for all flows; h − z is the distance from the bottom. So formally, the

Fig. 10. Entrainment rate vs. time for parameters (15) except τ0/ρ = 1m2s−2.

Fig. 12. Velocity profiles at different instances in Bingham and Herschel–Bulkley fluids.

τ xz ¼ μ 0

∂vx þ τxzðturbÞ ; ∂z



2
∂v
∂v τxzðturbÞ ¼ ρl

x

x ; ∂z ∂z

l ¼ κ ðh−zÞ

ð20Þ

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

Fig. 13. Flow depth vs. time in entraining flows of Bingham and Herschel–Bulkley fluids.

equations for turbulent flows in the Prandtl theory coincide with equations for laminar flows with a certain nonlinear viscosity depending not only on the shear rate but also on the distance from the bottom

τxz



∂v 2 2
∂v
¼ μ x ; μ ¼ μ 0 þ μ T ; μ T ¼ ρκ ðh−zÞ

x

; ∂z ∂z

ð21Þ

µT is called the turbulent viscosity coefficient. The boundary conditions at the surface and at the bottom in an entraining flow in consideration are:

z¼0:




∂v
∂vx ¼ 0; z ¼ hðt Þ : vx ¼ 0; μ

x

¼ τc : ∂z ∂z

ð22Þ

As in the previous sections, we reformulated the system (20)–(22), introducing dimensionless variables given by formulae (10), (11). The examples of calculations results shown below are for θ = 30°, τc/ ρ = 2.3 m2/s2 and different values of ν = μ/ρ, h(0). The behavior of the flow after the start of entrainment is shown in Figs. 15–17; time is counted from the beginning of entrainment. Figs. 16 and 17 demonstrate the influence of the initial flow depth and the value of the kinematic viscosity coefficient on the entrainment rate, respectively. At large time the dependence of the flow depth on time tends to be linear, and the entrainment rate tends to a constant not depending on the initial flow depth. The entrainment rate is larger for larger values of the viscosity coefficient. The velocity at the surface and the Reynolds number are plotted as functions of time in Figs. 18 and 19, respectively. Fig. 18 shows that at large time the velocity at the surface, i.e., maximum velocity in a given cross-section, increases proportionally to time.

Fig. 14. Velocity increase due to entrainment in Bingham and Herschel–Bulkley flows.

7

Fig. 15. Velocity profiles in a turbulent flow with entrainment (Prandtl turbulence theory) at different instances; h(0) = 1 m, ν = 0.001 m2 s−1.

7. Turbulent flows – LPY-model One of the basic assumptions of the Prandtl theory is the absence of a longitudinal pressure gradient and of gravity. Though this assumption is not valid for many flows, the Prandtl theory is still widely used. However, many more complicated models appeared. Here we present the results of simulations performed with the use of the three-parameter model referred to as the Lushchik–Paveliev–Yakubenko (LPY) model (Leontiev et al., 2009; Lushchik et al., 1978, 2004). This model was created to describe a variety of fluid flows along impermeable and permeable walls in presence of a pressure gradient, heat and mass transfer and other physical processes. In slope flows the pressure gradient as a driving force is replaced by gravity. We consider an unsteady onedimensional flow as in the previous sections. Three parameters that describe the flow in the LPY model are the kinematic turbulent   shear stress, T ¼ 1ρ τ turb ¼ − v0x v0z , the turbulent energy density, E ¼  0 0   0 0  1 vx vx þ vz vz , and the parameter ω, which is linked to the scale of 2 turbulence L by the relation ω ¼ LE2 . Here angular brackets mean Reynolds averaging and the prime stands for the difference between the true and Reynolds averaged parameter values. The parameter ω was first introduced by Kolmogorov (1942). It characterizes the square of the frequency of turbulence. The model contains empirical coefficients, which were determined by matching calculation results and experimental data. The model coefficients, once determined, were afterwards successfully used to describe different turbulent flows without any additional calibration (Leontiev et al., 2009; Lushchik et al., 1978, 2004). The system of LPY equations for our problem is ρ

∂v ∂τ ∂v ∂v ¼ ρg sin θ þ ; τ ¼ τxz ¼ ρν þ τ turb ¼ ρν þ ρT ∂t ∂z ∂z ∂z

ð23Þ

Fig. 16. Depth of entrained layer vs. time for two values of initial flow depth, h0 = h(0).

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

8

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

Fig. 17. Dependence of the entrainment rate on the fluid viscosity coefficient. Fig. 19. Reynolds number vs. time – Prandtl turbulence theory.

 

pffiffiffi E ∂E ∂ ∂E ∂v ¼ − c E L þ c1 ν 2 þ DE þT ∂t ∂z ∂z ∂z L

ð24Þ

 

pffiffiffi T ∂T ∂ ∂T ∂v ¼ − c5 EL þ c6 ν 2 þ Dτ þ c7 E ∂t ∂z ∂z ∂z L

ð25Þ

   

pffiffiffi ω ∂ω ∂ ∂ω T ∂v ∂v ¼ − 2c EL þ 1:4c1 ηf ω 2 þ Dω þ 2c4 sign þ ω E ∂t ∂z ∂z ∂z ∂z L

ð26Þ pffiffiffi Here v = vx, Dϕ ¼ aϕ EL þ α ϕ ν ðϕ ¼ E; T; ωÞ, f ω ¼ 1− 2c11 The values of the dimensionless coefficients are



L ∂E E ∂z

2 .

c ¼ 0:3; c1 ¼ 5π=4; c4 ¼ 0:04; c5 ¼ 3c; c6 ¼ 9c1 ; c7 ¼ 0:2;

ð27Þ

aE ¼ aω ¼ 0:06; aτ ¼ aE c5 =c; α E ¼ α τ ¼ 1; α ω ¼ 1:4:

ð28Þ

The boundary conditions for an entraining flow are at the upper surface z ¼ 0 : τ ¼ 0;

at the bottom z ¼ h : v ¼ 0;

E ¼ 0;

∂E ¼ 0; ∂z

T ¼ 0;

∂ω ¼ 0; ∂z

τ ¼ τc :

Fig. 18. Velocity at the surface vs. time.

∂v ¼ 0; ∂z ð29Þ

ð30Þ

The problem was rewritten in the dimensionless variables introduced by formulae (10) and (11). Software has been created to solve Eqs. (23)–(30). Variation of all parameters including the flow depth, velocity, turbulence energy and turbulent stresses in unsteady motion with entrainment were investigated numerically. Examples of calculation results in dimensional form and comparison with results obtained by Prandtl theory are presented in Figs. 20–23. The results shown in Figs. 19–21 are obtained for θ = 30°, hð0Þ ¼ 1 m, τc/ρ = 2.3 m2 s−2, ν ¼ μρ ¼ 0:001 m2 s−1 . Entrainment begins when |τ| = τc at the bottom. Time indicated in all figures is counted from the beginning of entrainment. The initial flow velocity is larger than the initial velocity of corresponding flow simulated by the Prandtl model, but the shear stress at the bottom is the same. Fig. 22 shows the dependence of the entrainment rate on the value of the bottom material shear strength. The velocity profiles are nonlinear, but the velocity at the surface, the depth averaged velocity and the flow depth increase with time approximately linearly due to entrainment as in the previous models. Fig. 22 shows that both the Prandtl model and the LPY model give almost the same value of the depth of the entrained layer for the first 10 s of the entrainment process. It should be emphasized that we did not perform calibration of model coefficients; they are supposed to be the same for all flows above a solid bottom. To make calculations we need to know only the rheological properties of the flow and the value of shear strength of the bottom material.

Fig. 20. Velocity profiles at different instances. LPY model.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

Fig. 21. Velocity at the surface vs. time. LPY model.

9

Fig. 23. Influence of the value of bottom material shear strength in the LPY turbulence model.

8. Conclusions This paper presents some results of simulations of laminar flows of Bingham and Herschel–Bulkley fluids and turbulent flows of linear viscous fluid entraining the bottom material. The entrainment hypothesis proposed by Issler and Pastor Pérez (2011) is used. Calculations have been made for idealized situations: homogeneous flows on homogeneous slopes. Besides, the values of the coefficients defining the rheological properties of the flow and of the bottom material do not directly correspond to any real materials. Therefore, the results cannot be directly applied to describe a real event. However, they show some qualitative features of flows entraining the bottom material that can be confirmed (or refuted) by observations. First of them is that the flow velocity increases if the flow entrains the slope material. Entrainment starts if the flow depth h becomes larger than the depth hst of the stationary flow with |τb| = τc, or if the strength of the bottom material τc becomes less than the value of the bed friction τb. At large times the velocity profile in entraining flows of both Bingham and Herschel– Bulkley fluids is almost linear except in a relatively narrow boundary layer near the upper surface. For turbulent flows on the contrary, a linear velocity profile is found only near the bed. Still in all studied models, the entrainment rate tends to a constant while the flow velocity and depth increase. Acknowledgements This work was supported by Russian Foundation of Basic Research (projects 12-01-00960-а and 13-08-00084). The authors would like to thank D. Issler and an anonymous reviewer for their valuable review of this paper, their helpful comments and suggestions. We thank also A.G. Kulikovsky, N.E. Leontiev and J.S. Zaiko for fruitful discussions.

Fig. 22. Flow depth increase due to entrainment. Full line – LPY model; dotted line – Prandtl model.

References Barbolini, M., Cappabianca, F., Issler, D., Gauer, P., Eglit, M.E., 2003. WP5: model development and validation. Erosion and deposition processes in snow avalanche dynamics: Report on the state of the art. SATSIE project Avalanche studies and model validation in Europe. Bouchut, F., Fernández-Nieto, E.D., Mangeney, A., Lagréee, P.-Y., 2008. On new erosion models of Savage–Hutter type for avalanches. Acta Mech. 199, 181–208. Bovet, E., Chiaia, B., Preziosi, L., 2010. A new model for snow avalanche dynamics based on non-Newtonian fluids. Meccanica 45 (6), 753–765. Briukhanov, A.V., Grigorian, S.S., Miagkov, S.M., Plam, M.Ya, Shurova, I.Ya, Eglit, M.E., Yakimov, Yu.L., 1967. On some new approaches to the dynamics of snow avalanches. In: Oura, H. (Ed.), Physics of Snow and Ice. Institute of Low Temperature. Science. Hokkaido University, Sapporo, Japan, pp. 1221–1241. Coussot, P., Laigle, D., Arattano, M., Deganutti, A., Marchi, L., 1998. Direct determination of rheological characteristics of debris flow. J. Hydraul. Eng. 124 (8), 865–868. Crosta, G.B., Imposimato, S., Roddeman, D., 2009. Numerical modelling of entrainment deposition in rock and debris-avalanches. Eng. Geol. 109 (1–2), 135–145. Davidson, P.A., 2004. Turbulence: An Introduction for Scientists and Engineers. Oxford University Press, Oxford, United Kingdom. Dent, J.D., Lang, T.E., 1980. Modeling of snow flow. J. Glaciol. 26 (94), 131–140. Dent, J.D., Lang, T.E., 1982. Experiments on the mechanics of flowing snow. Cold Reg. Sci. Technol. 5 (3), 253–258. Dent, J., Lang, T., 1983. A biviscous modified Bingham model of snow avalanche motion. Ann. Glaciol. 4, 42–46. Dent, J.D., Burrel, K.J., Schmidt, D.S., Louge, M.Y., Adams, E.E., Jazbutis, T.G., 1998. Density, velocity and friction measurements in dry-snow avalanche. Ann. Glaciol. 26, 247–252. Eglit, M.E. 1968. Theoretical approaches to calculations of snow avalanche motion. VINITI. Hydrology and Glaciology. pp. 60–97. (in Russian). English translation in M.E. Eglit. Theoretical approaches to avalanche dynamics.Glaciological Data. Report GD-16, World Data Center A for Glaciology (Snow and Ice), November 1984, pp. 63–116. Eglit, M.E., 1983. Some mathematical models of snow avalanches. In: Shahinpoor, M. (Ed.), Advances in the Mechanics and the Flow of Granular Materials, vol.2. Gulf Publishing Company, Houston, Texas, pp. 577–588. Eglit, M.E., Demidov, K.S., 2005. Mathematical modeling of snow entrainment in avalanche motion. Cold Reg. Sci. Technol. 43 (1–2), 10–23. Fischer, J.-T., Granig, M., Jörg, P., 2012. 2012. Snow entrainment in applied avalanche modeling using SamosAT. Geophys. Res. Abstr. 14, EGU2012–EGU12620. French, R.H., 1985. Open-Channel Hydraulics. McGraw-Hill Book Company, New York. Gauer, P., Issler, D., 2004. Possible erosion mechanism in snow avalanches. Ann. Glaciol. 38, 384–392. Issler, D., 2003. Experimental information on snow avalanches. In: Hutter, K., Kirchner, N. (Eds.), Dynamic Response of granular and porous materials under large and catastrophic deformations. Springer, pp. 109–160. Issler, D., Pastor Pérez, M., 2011. Interplay of entrainment and rheology in snow avalanches; a numerical study. Ann. Glaciol. 52 (58), 143–147. Kern, M.A., Tiefenbacher, F., McElwaine, J.N., 2004. The rheology of snow in large chute flows. Cold Reg. Sci. Technol. 39, 181–192. Kolmogorov, A.N., 1942. Equations of turbulent motion of incompressible fluid. Izv. Akad. Nauk SSSR 6 (1–2), 56–58. Lang, T.E., Dent, J.D., 1980. Scale modeling of snow avalanche impact on structures. J. Glaciol. 26 (94), 189–196. Leontiev, A.I., Lushchik, V.G., Yakubenko, A.E., 2009. A heat-insulated permeable wall with suction in a compressible gas flow. Int. J. Heat Mass Transf. 52 (17–18), 4001–4007. Lushchik, V.G., Paveliev, A.A., Yakubenko, A.E., 1978. Three-parameter model of shear turbulence. Fluid Dyn. 13 (3), 350–362. Lushchik, V.G., Paveliev, A.A., Yakubenko, A.E., 2004. Transition to turbulence in the boundary layer on a plate in the presence of a negative free-stream pressure gradient. Fluid Dyn. 39 (2), 250–259. Maeno, N., Nishimura, K., 1979. Fluidization of snow. Cold Reg. Sci. Technol. 1 (2), 109–120.

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002

10

M.E. Eglit, A.E. Yakubenko / Cold Regions Science and Technology xxx (2014) xxx–xxx

Maeno, N., Nishimura, K., Kaneda, Y., 1980. Viscosity and heat transfer in fluidized snow. J. Glaciol. 26 (94), 263–274. Mellor, M., 1975. A review of basic snow mechanics. IAHS Publ. 114, 251–291 (Symposium at Grindelwald 1974 – Snow Mechanics). Naaim, M., Faug, T., Naaim-Bouvet, F., 2003. Dry granular flow modelling including erosion and deposition. Surv. Geophys. 24, 569–585. Naaim, M., Naaim-Bouvet, F., Faug, T., Bouchet, A., 2004. Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects. Cold Reg. Sci. Technol. 39, 193–204. Nishimura, K., 1996. Viscosity of fluidized snow. Cold Reg. Sci. Technol. 24 (2), 117–127. Nishimura, K., Maeno, N., 1989. Contribution of viscous forces to avalanche dynamics. Ann. Glaciol. 13, 202–206. Norem, H., Irgens, F., Schieldrop, B., 1989. Simulation of snow-avalanche flow in run-out zones. Ann. Glaciol. 13, 218–225. Prandtl, L., 1925. Über die ausgebildete turbulenz. ZAMM 5, 136–139. Rognon, P., Chevoir, F., Bellot, H., Ousset, F., Naaim, M., Coussot, P., 2007. Rheology of dense snow flows. Inferences from steady state chute-flow experiments. J. Rheol. 52 (3), 729–748. Shapiro, L.H., Johnson, J.B., Sturm, M., Blaisdell, G.L., 1997. Snow mechanics: review of the state of knowledge and applications. CRREL Report 97–3. U.S. Army Corps of Engineers Cold Regions Research and Engineering Laboratory, Hanover, New Hampshire.

Soltani, F., Yilmazer, U., 1998. Slip velocity and slip layer thickness in flow of concentrated suspensions. J. Appl. Polym. Sci. 70, 515–522. Sovilla, B., 2004. Field experiments and numerical modelling of mass entrainment. PhD thesis. ETH Dissertation nr. 15462Swiss Federal Institute of Technology, Zurich, Switzerland. Sovilla, B., Bartelt, P., 2002. Observations and modelling of snow avalanche entrainment. Nat. Hazards Earth Sci. 2, 169–179. Sovilla, B., Sommavilla, F., Tomaselli, A., 2001. Measurements of mass balance in dense snow avalanche events. Ann. Glaciol. 32, 230–236. Sovilla, B., Burlando, P., Bartelt, P., 2006. Field experiments and numerical modeling of mass entrainment in snow avalanches. J. Geophys. Res. 111 (F03007), 1–16. Sovilla, B., Margreth, S., Bartelt, P., 2007. On snow entrainment in avalanche dynamics calculations. Cold Reg. Sci. Technol. 47, 69–79. Tai, Y.C., Kuo, C.Y., 2008. A new model of granular flows over general topography with erosion and deposition. Acta Mech. 199, 71–96. Tiefenbacher, F., Kern, M.A., 2004. Experimental devices to determine snow avalanche basal friction and velocity profiles. Cold Reg. Sci. Technol. 38, 17–30. Voitkovsky, K.F., 1975. Mekhanicheskie svoistva snega (Mechanical properties of snow). Nauka, Moscow (126 pp.).

Please cite this article as: Eglit, M.E., Yakubenko, A.E., Numerical modeling of slope flows entraining bottom material, Cold Reg. Sci. Technol. (2014), http://dx.doi.org/10.1016/j.coldregions.2014.07.002