HAL Id: hal-01059334
https://hal.archives-ouvertes.fr/hal-01059334v1
Preprint submitted on 29 Aug 2014 (v1), last revised 2 Sep 2013 (v4)
HAL
is a multi-disciplinary open access archive for the deposit and dissemination of sci- entific research documents, whether they are pub- lished or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire
HAL, estdestinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.
A new shallow approximation for tridimensional non-isothermal viscoplastic lava flows
Noé Bernabeu, Pierre Saramito, Claude Smutek
To cite this version:
Noé Bernabeu, Pierre Saramito, Claude Smutek. A new shallow approximation for tridimensional
non-isothermal viscoplastic lava flows. 2014. �hal-01059334v1�
A new shallow approximation for
tridimensional non-isothermal viscoplastic lava flows
Noé Bernabeu
†,‡, Pierre Saramito
†,∗and Claude Smutek
‡Abstract – A new shallow reduced model for the non-isothermal tridimen- sional viscoplastic fluid flowing on a general topography is presented in this paper. Both the consistency and yield stress are supposed temperature- dependent. An asymptotic analysis leads to reduce the 3D problem to a 2D surface one with depth-averaged equations. These equations are numer- ical approximated by an autoadaptive finite element method, based on the Rheolef C++ library, allowing to track accurately the front position. The proposed approach is first evaluated by comparing numerical prediction with non-isothermal experimental measurements for a silicone oil dome. Next, the December 2010 eruption of Piton de la Fournaise (La Réunion island) is numerically reproduced and compared with available data.
The risk assessments for volcanic lava flow pose a difficult challenge to numerical methods.
Indeed, these flows with free surface couple the fluid dynamic equations (Navier-Stokes like) with thermal effects, as diffusion-convection in lava, radiation and convection in air, diffusion in the substrate while the fluid rheology is complex and temperature-dependent.
In order to overcome these complexities, some authors proposed a probabilistic approach based on the topography: the flow path is determined by the maximum slope direction.
See the DOWNFLOW [28] or the ELFM [9] codes. Such codes require few computational time and few input data, which is a definitive advantage. However, these approaches do not permit to predict the evolution of the lava, its thickness and its arrested state. The deterministic approaches adopt more complete descriptions of the phenomena. Harris et al. [14] proposed a kinematic one-dimensional thermo-rheological model of a lava flow in a confined channel (the FLOWGO code). In [16], a three-dimensional model for Newtonian fluids, including the energy conservation, solidification and free surface is presented. A variant of the deterministic approach is the cellular automata method for 2D computation systems, based on a spacial partition into cells whose state evolves according to that of their neighbours and a transition function which describes the exchanges of lava and heat:
see the MAGFLOW code [29] that evaluates the lava exchange from the steady state solution of the Bingham fluid equations on an inclined plane, or the SCIARA code [27]
that bases on a minimization rule of the difference in height between neighbouring cells while the lava cooling takes into account the radiation at the free surface and the heat exchanges due to lava mixing.
The lava flow extends over several kilometres while its thickness remains relatively small.
An asymptotic analysis based on this aspect ratio allows to reduce rigorously the three- dimensional set of conservation laws and constitutive equations to a surface bidimensional
†Laboratoire Jean Kuntzmann, UMR 5524 Univ. J. Fourier - Grenoble I and CNRS, BP 53, F-38041 Grenoble cedex, France
∗Corresponding author (e-mail: [email protected])
‡Lab. géosciences – IPGP and La Réunion university, 15, av. René Cassin, CS 92003, 97744 Saint- Denis cedex 09, France
system of equations. For inertia dominated flows, we obtain the classical Saint-Venant shallow water model [3], where viscous terms are neglected. Gerbeau and Perthame [13]
proposed a modified Saint-Venant system where viscous effects are maintained thought a wall friction term. Costa and Macedomio [8] extended this viscous variant by including non-isothermal and viscous heating effects for volcanic lava flows. Notice that free sur- face viscoplastic lava flows are highly viscous and these inertia dominated models could be improved. Huppert [17] investigated Newtonian viscous dominated flows. For more complex non-Newtonian rheologies, shallow approximations of the horizontal dam break problem were first studied for a viscoplastic fluid by Lui and Mei [18] and revisited by Balmforth and Craster [1]. Computations for some specific tridimensional topographies was next performed by using a specific axisymmetric coordinate systems: a curved chan- nel [20] and a conical surface [31]. Recently, an extension for an arbitrarily topography is presented by the present authors in [5]. Tacking into account thermal effects is more complex with the asymptotic analysis approach: the problem does not reduce to a sur- face bidimensional one, it remains fully three-dimensional. In order to overcome this difficulty, different approaches based on a depth-average version of the heat equation was investigated by Bercovici and Lin [4] for a Newtonian fluid to model the cooling of mantle plume heads with temperature-dependent buoyancy and viscosity. Balmforth and Craster applied this approach to the evolution of a lava dome [2] for viscoplastic fluid with temperature-dependent consistency and yield stress.
The aim of this paper is to bring a new robust and efficient numerical method based on the asymptotic analysis approach for the resolution of the shallow approximation of tridimensional non-isothermal viscoplastic flow problem with temperature-dependent consistency and yield stress on a general topography. The first section presents the tridi- mensional problem, its surface bidimensional reduction and the numerical resolution of the set of equations. The second section validates our approach based on comparisons of numerical simulations with experimental measurements performed on a non-isothermal silicone oil dome laboratory experiment. The last section develops in details some com- parison of numerical predictions with data available for a real volcanic lava flow, the December 2010 Piton de la Fournaise lava one.
1 Problem statement and numerical method
1.1 Tridimensional formulation
The viscoplastic Herschel-Bulkley constitutive equation [15] expresses the deviatoric part τ of the stress tensor versus the rate of deformation tensorγ˙ =∇u+∇uT as:
τ=K(θ)|γ|˙ n−1γ˙+τy(θ)|˙γγ|˙ when γ˙ 6= 0,
|τ| ≤τy. otherwise. (1a)
whereuis the velocity field,θis the temperature,K(θ)>0the temperature-dependent consistency,n >0 is the power-law index andτy(θ)is the temperature-dependent yield stress. The consistency and the yield stress are in general decreasing with the temper- ature. Here |τ| =
(1/2)P3
i,j=1τij21/2
denotes the conventional norm of a symmetric tensor τ in mechanics. The total Cauchy stress tensor is σ =−p.I +τ where p is the
pressure andI the identity tensor. Whenτy = 0andn= 1, the fluid is Newtonian and Kis the viscosity. Whenτy= 0andn >0the fluid is a quasi-Newtonian power-law one.
Whenτy >0 and n= 1, the model reduces to the Bingham one [6]. The constitutive equation (1a) should be completed by the mass, momentum and energy conservation laws:
div u= 0, (1b)
ρ(∂tu+u.∇u)−div(−p.I+τ) =ρg, (1c) ρCp(∂tθ+u.∇θ)−div(k∇θ)−τ : ˙γ
2 = 0. (1d)
where ρ is the constant density, Cp the constant specific heat, k the thermal conduc- tivity andgthe gravity vector. Notice that there are four equations (1a)-(1d) and four unknownsτ,u= (ux, uy, uz),pandθ. The corresponding problem is closed by defining the boundary and initial conditions.
Figure 1: Eruption flow on a variable topography with cooling.
The flow over a variable topography is considered with cooling and mass eruption from a vent (see Fig.1). For any timet >0, the flow domain is represented by:
Q(t) ={(x, y, z)∈Ω×R; f(x, y)< z < f(x, y) +h(t, x, y)}
whereΩis an open and bounded subsetR2. The functionf denotes the topography and his the flow height. Notice that, sincehis a function,his mono-valued, which excludes the small front cusp often observed on real flows. The eruption zone is described by an open subset Ωe of Ω (see Fig. 1) and Ωs = Ω\Ωe denotes its complementary. The boundary ∂Q(t)of the flowing lava volume Q(t) splits in four parts: the bottom relief in the eruption zone Γe, the bottom relief out of the eruption zone Γs and the top free
surface Γf(t):
Γe={(x, y, z)∈Ωe×R; z=f(x, y)}, Γs={(x, y, z)∈Ωs×R; z=f(x, y)},
Γf(t) ={(x, y, z)∈Ω×R; z=f(x, y) +h(t, x, y)}.
Also, S={(x, y, z)∈Ωs×R; z < f(x, y)} denotes the substratum. For anyt >0, the boundary conditions are a non-slip (Dirichlet) condition on the bottom for the velocity field and natural (Neumann) one on the free surface:
ux=uy= 0 anduz=weonΓe∪Γs, (1e)
σ.νf =0onΓf(t). (1f)
whereνfdenotes the unit outward vector of∂Q(t)onΓf(t). Here,weis the lava eruption velocity, defined on]0,+∞[×Γe, and satisfyingwe6= 0onΓe andwe= 0onΓs. Notice that, from the Dirichlet condition, we have u= 0 on Γs and a vertical profile on the eruption zone Γe. For the temperature, a Dirichlet condition in the eruption zone, a conduction flux on the bottom with the substratum (out of eruption zone) and radiative and convection fluxes with air on the free surface are considered:
θ=θeonΓe, (1g) kνs.∇θ=ksνs.∇θsonΓs, (1h) kνf.∇θ+ǫσSB θ4−θ4a
+λ(θ−θa) = 0onΓf(t). (1i) where θe is the initial eruption temperature, θs is the temperature in the substratum, θa is the temperature in the air,ks is the thermal conductivity of the substratum,νsis the unit outward vector of∂Q(t)onΓs,ǫis the emissivity,σSB is the Stefan-Boltzmann constant and λ is the convective heat transfer coefficient. Here, the term ksνs.∇θs
represents the heat flux from the substratum: it will be estimated in the forthcoming paragraph 1.2. It remains to describe the evolution of the free surface. It is convenient to introduce the level set functionϕthat is defined for allt >0and(x, y, z)∈Ω×Rby:
ϕ(t, x, y, z) =z−f(x, y)−h(t, x, y).
Notice that the zero level, where ϕ(t, x, y, z) = 0, coincides with the free surface Γf(t). Let us write that the level set function is transported by the u velocity field:
∂tϕ+u.∇ϕ= 0. OnΓf(t), wherez=f +h, this relation becomes:
∂th+ux∂x(f +h) +∂y(f+h) =uz in ]0,+∞[×Ω (1j) This is a first-order transport equation for the height hthat should be completed by an initial condition:
h(t= 0, x, y) =hinit, ∀(x, y)∈Ω, (1k) where hinit is given. It should also be completed by an inflow boundary condition forh:
h=hext on]0,+∞[×∂Ω− where∂Ω−={(x, y)∈∂Ω| uxνw,x+uyνw,y<0}.
We assume that the domainΩis sufficiently large that the flow never reaches the bound- aries of Ω. Then, on ∂Ω, ux = uy = 0 for any time and Ω− = ∅ and thus this last
boundary condition is here empty. The set of equation is finally completed by an initial condition for the velocityuand for the temperatureθ:
u(t= 0) =uinit and θ(t= 0) =θinit inQ(0). (1l) Finally, the tridimensional problem expresses: findh,τ, u,pandθ satisfying (1a)-(1l).
For applications to lava flows, the finite element of finite difference discretisations of this problem lead to huge nonlinear and time-dependent set of equations. The next paragraph present a reduction to a bidimensional problem.
1.2 Reduction to a bidimensional problem
The asymptotic analysis approach developped here was initiated by Lui and Meil [18] and then revisited by Balmforth and Craster [1] for an isothermal bidimensional viscoplastic flow on a constant slope. Balmforth and Craster then extended this analysis to the non-isothermal case for an axisymmetric geometry and presented applications to a lava dome [2]. In [5], the present authors extended this asymptotic analysis to the case of an isothermal tridimensional viscoplastic flow on an arbitrarily topography. In the present paper, this analysis is extended to the case of a non-isothermal flow: the consistency and the yield stress are supposed temperature-dependent.
For any function g of]0,+∞[×Q(t), let us denoteg its vertical-averaged value, defined for all(t, x, y)∈]0,+∞[×Ωby
g(t, x, y) =
1 h
Z f+h f
g(t, x, y, z)dz whenh6= 0,
0 otherwise.
For simplicity, assume that the consistency and the yield stress depend only of the vertical-averaged temperature: K(θ) =K(θ)andτy(θ) =τy(θ). It means that the effects of vertical variation of temperature on the consistency and the yield stress are here ne- glected, while their effects in the horizontal directionsxandyis still taken into account.
The problem is reformulated with dimensionless quantities and unknowns, denoted with tildes. The temperature is expressed asθ=θa+ (θe−θa)˜θ. The temperature-dependent consistency and yield stress are rescaled asK(θ) =KeK(˜˜ θ)andτy(θ) =τy,eτ˜y(˜θ)where Ke=K(θe)andτy,e=τy(θe). Remark thatK(1) = 1˜ and˜τy(1) = 1. LetH be a charac- teristic flow height andLbe a characteristic horizontal length of the bidimensional flow domainΩ. Let us introduce the dimensionless aspect ratio:
ε=H/L.
Let U = ρgH3/(ηL) a characteristic flow velocity in the horizontal plane, where η=K0(U/H)n−1 is a characteristic viscosity and g =|g| denotes the norm of gravity vector. The characteristic velocity expands as
U =
ρgH2 K0L
n1
H.
LetW =εU be a characteristic velocity in the vertical direction,T =L/U be a charac- teristic time and P =ρgH a characteristic pressure. We consider the following change
of variables:
x=Lx, y˜ =L˜y, z=Hz, t˜ =T˜t, p=Pp, h˜ =H˜h, ux=Uu˜x, uy =Uu˜y, uz=Wu˜z.
Remark the non-isotropic scaling procedure for the vertical coordinatezand the vertical vector componentuzof the velocity vector. In the rest of this paper, only the dimension- less problem is considered: for simplicity and since there is no ambiguity, the tilde are omitted on the dimensionless variables. After the change of variable and an asymptotic analysis (see [5] details for isothermal case) the problem reduces to:
(P): find handθ satisfying:
∂th+div|| hu||
=we in ]0,+∞[×Ω (2a) h(t= 0) =hinit in Q(0), (2b)
∂(f+h)
∂n = 0 in ]0,+∞[×∂Ω, (2c)
∂tθ+ux∂xθ+uy∂yθ+uz∂zθ− 1
P e∂zzθ= 0 in ]0,+∞[×Q(t), (2d) θ(t= 0) =θinit in Q(0), (2e)
∂θ
∂z +R pµ(θ)θ+Nuθ= 0 on]0,+∞[×Γf(t), (2f)
−∂θ
∂z+ks
k P es
πt 12
θ= 0 on]0,+∞[×Γs, (2g) θ= 1 on]0,+∞[×Γe. (2h) The problem involves six dimensionless numbers:
Bi = τy,eH
ηU the Bingham,
P e = ρCpU L
k the Péclet number,
P es = ρsCp,sU L ks
the Péclet number in the substrate,
Nu = λH
k the Nüsselt number,
R = HǫσSB(θe−θa)3
k a radiation number,
µ = θa
θe−θa
a temperature ratio number, for radiation.
Also,pµ(θ) =θ3+ 4µθ2+ 6µ2θ+ 4µ3 for anyθ≥0. Remark that problem(P)involves onlyhandθ: the velocity componentsu||= (ux, uy)anduzadmits then explicit expres- sions involvinghandθ:
u||=
n
n+1|∇||(f+h)|n1dir(∇||(f+h))K−1(θ)h
(f+hc(θ)−z)n+1n −h
n+1
cn
i whenz∈[f, f+hc(θ)],
−n+1n |∇||(f +h)|n1dir(∇||(f+h))K−1(θ)h
n+1 n
c
whenz∈]f +hc(θ), f+h],
and
uz(t, x, y, z) =we− Z z
f
div||(u||)dz where hc(t, x, y) = max
0, h− Bi τy(θ)
|∇||(f+h)|
. For convenience, we also denote as dir(v||) = v||/|v||| the direction of any nonzero plane vector. In (2g), the heat flux from the substratum has been estimated as in [11, p. 6], thanks to an autosimilar so- lution for the temperature in the substratumθs(t, z)given in [7, p. 53]. This approach leads to an explicitly expression of θs versus θ and then we obtain a Robin boundary condition in terms ofθonly.
1.3 Reduction of the heat equation
Notice that the reduced problem still remain tridimensional, as (2d) is still defined on the tridimensional flow domain Q(t). It remains to integrate this equation in the ver- tical direction and express the problem in terms of the averaged temperature θ. For any t > 0 and (x, y) ∈ Ω, the heat equation (2d) is integrated from z=f(t, x, y) to z=f(t, x, y) +h(t, x, y). Then, using (2a) anddivu= 0leads to:
h∂tθ+div h θu||
−div hu||
θ−we 1−θ
− 1
P e[∂zθ]f+hf = 0. (3) In order to obtain a well-posed problem in term ofθ instead ofθ, an additional relation should be introduced: the so-called closure relation, that expresses θ versusθ. A verti- cal profile θ is chosen according to θ as θ(t, x, y, z) =ϕ(t, x, y, z)θ(t, x, y) where ϕis a unknown function satisfyingϕ= 1and the boundary conditions in the vertical direction:
∂zϕ+Rpµ(θϕ)ϕ+Nu ϕ= 0 onΓf(t),
−∂zϕ+ks
k P es
πt 12
ϕ= 0 onΓsandθϕ= 1 onΓe. With these notations, (3) becomes:
h ∂tθ+ϕu||.∇θ
+div h(ϕu||−u||)
θ−we 1−θ
− 1
P e[∂zϕ]f+hf θ= 0. (4) Notice thatϕu||interprets as a weighted averaged velocity. The closure equation (4) was first introduced by Bercovici and Lin [4] in the context of the cooling mantle plume heads:
these authors showed that, restrictingϕ(t, x, y, z)to be a second degree polynomial inz for any fixed(t, x, y), leads to a well-posed reduced problem when replacing (2d) by (4) in problem(P). Moreover, the computation of the ϕpolynomial coefficient at any(t, x, y) was easy and explicit. However, this second order polynomial approximation in z do not permit to eliminate the term div h(ϕu||−u||)
θ from equation (3). There is also no evidence that its always positive: it can thus generate an exponential growth of the averaging temperature and this possible behavior is difficult to handle in numerical simulation. Finally, the physical meaning of this therm is unclear and Bercovici and Lin neglected div h(ϕu||−u||)
θ in the previous equation. These authors showed that a real numerical error is committed with this hypothesis: its order of error is of about 10%
for their specific problem [4, p. 3307]. Notice that 10% of errors is acceptable for some applications ; nevertheless, from a mathematical point of view, all convergence properties versusεare definitively lost.
In the present paper, we propose a variant of this approach that conserve the convergence properties versusε: we chooseϕas a third degree polynomial inz. This choice allows one additional degree of freedom at each(t, x, y)and we choose to impose exactlyϕu||=u||at any(t, x, y)as an additional constraint. Let(t, x, y)be fixed andϕ(z) =az3+bz2+cz+d:
the four unknown coefficients a, b, cand d are the only solution of the following four equations;
ϕ = 1, (5a)
ϕu|| − u|| = 0, (5b)
∂ϕ
∂z + (R pµ(θϕ) +Nu)ϕ = 0 onz=f+h, (5c)
−∂ϕ
∂z +ks
k P es
πt 12
ϕ= 0onΓsandθϕ = 1 onΓe, onz=f (5d) The reduced problem is then obtained from (P)by replacing (2d) by
h ∂tθ+u||.∇θ
−we(1−θ).−
3ah2+ 2bh P e
θ= 0. (6a)
and the initial and boundary conditions (2e)-(2h) by
θ(t= 0) = θinit onΩ, (6b)
∂θ
∂n = 0 on]0,+∞[×∂Ω. (6c)
1.4 Numerical resolution
The nonlinear reduced problem inhandθ is first discretized versus time by an implicit second order variable step finite difference scheme (see [5]). It leads to a sequence of nonlinear subproblems that depends only of the horizontal coordinates (x, y)∈ Ω. An under-relaxed fixed point algorithm is used for solving these nonlinear subproblem: it allows to decouple the parabolic evolution equation evolution, in terms of h, and the averaged heat equation, in term of θ. Then, these equations are discretized by an au- toadaptive finite element method: the practical implementation bases on the Rheolef finite element library [24]. Such autoadaptive mesh method was first introduced in [25]
for viscoplastic flows and then extended in [22], and we refer to theses articles for the implementation details. The procedure bases on a mesh adaptation loop at each time step: your goal is to catch accurately the front evolution, wereh= 0(see Fig.2). At the front, both thehandθgradients are sharp.
Figure 2: Uniform (a) and dynamic (b) auto-adaptive meshes. Zoom near the front (c).
This mesh was automatically generated for the volcanic lava flow simulation presented in the forthcoming section 3.
2 Comparison with a silicone oil dome experiment
2.1 Experimental setup and physical parameters
Physical quantities symbol value unit
euption flow rate Q 2.2×10−8 m−3.s−1
initial temperature of the eruption fluid θe 315.15 K initial temperature of air and polystyrene θa 293.151 K
fluid density ρ 954 kg.m−3
air density ρa 1.2 kg.m−3
fluid viscosity atθetemperature Ke 3.4 Pa.s
fluid emissivity ǫ 0.96 -
thermal conductivity of the fluid k 0.15 W.m−1.K−1
fluid specific heat Cp 1500 J.m−1.K−1
convective heat transfer coefficient with air λ 1 − 3 W.m−2.K−1 thermal conductivity of the polystyrene ks 0.03 W.m−1.K−1 thermal diffusivity of the polystyrene κs 6×10−7 m2.s−1
vent radius re 2 − 4 mm
constant in Arrhenius law α 0.00808044 K−1
constant in Huppert [17] formula a 0.715 dimensionless Table 1: Physical parameters for the laboratory experiment [10,12].
1 10
15 20 25 30 35 40 45 50 55 60 θin ˚C
K(θ)in Pa.s
Figure 3: Evolution of the silicon oil viscosity (in Pa.s) vs temperature (in◦C).
The evolution of lava dome was numerical studied by Balmforth and Craster using yield stress fluids [2]. In 2012, Garel presented [12] a laboratory experiment that reproduce analogously the lava dome growth with silicon oil (Newtonian fluid) and allowed accurate temperature measurements. The fluid is injected through a vertical vent in an horizontal polystyrene plane at constant flow rate Q. The fluid is initially heated at θe above the ambient temperature θa and a system of optical and infrared cameras are used to
follow the dome growth and the surface temperature. Table1groups the relevant physical parameters for this experiment. The fluid viscosity was measured for the full experimental temperature range (see Fig. 3). Data was first presented in [10, pp. 60–62]. Observe that the viscosity can be approximated by an Arrhenius lawK(θ) =Keeα(θe−θ) where α= 0.00808044K−1 was obtained by a linear regression. The flow thought the vent is supposed to be a steady Poiseuille flow. As the vent is a circular hole with radius re, the vertical velocitywe is a second order polynomial versus the radius. Then we(r) = cmax(0, r2e−r2)wherec >0 is such thatRre
0 we(r)rdr=Qi.e. c= 2Q/(πre4).
2.2 Comparison between experiments and simulations
20 25 30 35 40
2 4 6 8 10
(a)t= 7480s,λ= 2
r(cm) surfacetemperatureθ(z=f+h)in◦C P2simulation
P3simulation K(θ)-P3simulation measurement
20 25 30 35 40
2 4 6 8 10
(b)t= 7480s
surfacetemperatureθ(z=f+h)in◦C
r(cm) λ= 1,K(θ)-P3simulation λ= 3,K(θ)-P3simulation measurement
Figure 4: Comparison between the present simulations and a laboratory experiment from [12]: surface temperature vs radius. (a) influence of some model variants (with λ= 2); (b) influence of theλconvection parameter.
Fig.4.a shows the computed surface temperature as obtained by numerical simulations for different variants of our shallow model. These simulations are also compared with experimental measurements, from [12]. Observe the good correspondence between all simulations and the experimental data. The labelP2is associated to a constant viscosity model similar to those of by Bercovici and Lin [4]: the temperature profile θ =ϕ θ is approximated with a second order polynomial ϕin the vertical direction and the term div h(ϕu||−u||)
θis neglected. The labelP3denotes a constant viscosity model where ϕis approximated by a third order polynomial that satisfiesϕu||=u||. Finally, the label K(θ)-P3 denotes a non-constant viscosity model based on the Arrhenius’s law. For this silicone oil dome problem, the results of the simulations are very close for all these model variants ; the comparison is presented at time t = 7480 s but this observation is still valid for other times. The change from aP2 to aP3 model has few effects on numerical computations, despite this change is more satisfactory from a mathematical point of view.
The change from aP3 to a K(θ)-P3 model has also few effects: the difference between the two curves is imperceptible on Fig. 4.a. Remark that in the present experiment,
20 25 30 35 40
2 4 6 8 10
t= 160s
r(cm) surfacetemperatureθ(z=f+h)in◦C K(θ)-P3simulation
measurement
20 25 30 35 40
2 4 6 8 10
t= 600s
r(cm) surfacetemperatureθ(z=f+h)in◦C K(θ)-P3simulation
measurement
20 25 30 35 40
2 4 6 8 10
t= 1200s
r(cm) surfacetemperatureθ(z=f+h)in◦C K(θ)-P3simulation
measurements
20 25 30 35 40
2 4 6 8 10
t= 3000s
r(cm) surfacetemperatureθ(z=f+h)in◦C K(θ)-P3simulation
measurement
20 25 30 35 40
2 4 6 8 10
t= 7480s
r(cm) surfacetemperatureθ(z=f+h)in◦C K(θ)-P3simulation
measurement
Figure 5: Comparison between the present simulations and a laboratory experiment from [12]: surface temperature vs radius (λ= 2 and K(θ)-P3 model). Time evolution fort= 160s,600s,1200s,3000sand7480s.
the viscosity varies about a factor two in the temperature range from 20◦C to 60◦C (see Fig.3). Based on an autosimilar solution for the lava dome problem, Huppert [17]
proposed an explicit formula of the reference heighthref and dome front radiusrf(t)for the isothermal case:
href=a
(ρ−ρa)gQ3 3K
1/4
and rf(t) =a2/3
3KQ (ρ−ρa)g
1/8
t1/2, (7) whereais a constant andρais the air density. Observe that the viscosityKappears with some 1/4 and 1/8 exponents: as a consequence, the dependence upon the temperature is finally negligible for the present silicone oil dome evolution experience. This situation will dramatically change for real lava flows, as presented in the next section.
The experimental estimate of the coefficient of convective heat transfer λranges from 1 to 3 W.m−2.K−1. Fig.4.b plots the computed surface temperature for the two extreme values ofλ. Observe that the corresponding numerical results are close. Moreover, the two curves wrap completely the experimental data. In the rest of this paragraph, the median valueλ= 2is used for all the simulations.
Fig.5, present the surface temperature evolution versus time for theK(θ)-P3model. All these simulations appear in good concordance with experimental data for times larger than t = 600s. At t = 160 s, we observe some discrepancy between the experimental data and the theoretical model. Following [12, p. 7], "we interpret the faster experimental cooling as due to lateral heat conduction in the substrate that initially widens the size of the thermal anomaly. This effect decreases as the temperature decreases at the front of the current, and becomes negligible fort >320s, when the experimental observations are well reproduced by theory."
3 Comparison with a real volcanic lava flow
3.1 The December 2010 Piton de la Fournaise lava flow
Figure 6: The Piton de la Fournaise 2010 lava flow (credit N. Bernabeu, May 2014).
Piton de la Fournaise, a volcano from La Réunion island, is among the most active volcanoes in the world. A lava flow occurred in December 2010 on the north flank of
Physical quantities symbol value unit
average eruption flow rate Q 9.7 m−3.s−1
eruption duration de 15 : 15 h:min
initial temperature of the fluid θe 1423 K
initial temperature of air and substrate θa 303 K
lava density ρ 2200 kg.m−3
air density ρa 1.2 kg.m−3
lava viscosity at θetemperature Ke 104 Pa.s lava yield stress atθe temperature τy,e 102 Pa.s
lava emissivity ǫ 0.95 -
lava thermal conductivity k 2 W.m−1.K−1
lava specific heat Cp 1225 J.m−1.K−1
convective heat transfer coefficient with air λ 80 W.m−2.K−1
flow characteristic thickness H 1 m
flow characteristic length L 1000 m
constant in viscosity Arrhenius law α 0.016447 K−1 constant in yield stress Arrhenius law β 0.016447 K−1
Table 2: Physical parameters for the lava flow. See [19,21,26,30] and [23] for the 2010 eruption data.
volcano. It was a basalt flow lasted few hours which outcome of a fissure on the volcano surface. Fig. 6 shows the flow at arrested state in 2014. See [23] for more information about the 2010 lava flow. Table 2 groups all the relevant physical quantities suitable for the simulation. The fluid is represented by a Bingham viscoplastic rheological model (the power law index is n= 1). The substrate and the lava present similar properties, such as density and thermal conductivity, since the substrate is composed by old cooled lava. The consistency and the yield stress depend strongly upon the temperature and vary on a large range of magnitude. More precisely, we suppose that the viscosity and the yield stress follow some Arrhenius laws: K(θ) =Keeα(θe−θ)andτy(θ) =τy,eeβ(θe−θ) where α=β = 0.016447K−1.
3.2 Digital elevation model and flow conditions
Volcano observatory from La Réunion and the laboratory of geosciences at La Réunion provided us a tridimensional digital elevation model (DEM) of the Piton de la Fournaise with a 5 meter horizontal resolution (see Fig. 7), covering the zone of deposit. This topographic survey being older than 2010, it can be used to define the relief functionf in our model. The volcano observatory also provided us the detailed contours of the final lava flow at arrested state (see Fig. 7 right): this data will be used for comparison with numerical simulations. The average lava flow rate Qduring the eruption is also known (see table 2). An expedition was necessary to complete some missing data and localize accurately the different vents (see Fig. 8.a). Observe on Fig.8.b that the 2010 volcanic flow is splited in two disjoint zones, denoted asAandB, where vents are numbered from 1 to 6. The flow rates of these two disjoint flows are estimated as proportional to their relative areas, i.e. QA =κ Q and QB = (1−κ)Q, whereκ =|A|/(|A|+|B|) and |A|
and|B|denote the areas of theAandB zones, respectively. Notice that these areas are
Figure 7: Digital elevation model of the Piton de la Fournaise volcano in 2008, with a 5 meter horizontal resolution (left). Zoom with the observed arrested flow contours (right).
known from the detailed contours of the lava flow at arrested state. On theAzone, we suppose that the vent1is a circular cone with a radius of 20 meters and is active during the full eruption time. On theB zone, we suppose that the vents from 2to6 are active one after the other and that they are circular. Let us denote bydethe eruption duration.
We suppose that the vent 2 is active from t= 0 to de/5 with a radius of 10 meters and the vent i, 3 ≤i ≤6, is active from t= (i−2)de/5 to (i−1)de/5 with a radius of 20 meters.
Figure 8: (a) View of the aligned cones along the fault line associated to different vents at arrested state (credit N. Bernabeu, 2014). (b) Schematic view of the different vent positions and the two separated zones of the observed arrested state.
3.3 Numerical simulations
Numerical simulations are performed for this volcanic lava flow: the arrested state is reached att= 25h and is represented on Fig.9over the DEM. Observe that the predicted arrested state is relatively closed to that of the observed one, represented by a thin white
Figure 9: Simulation of the volcanic lava flow: prediceted arrested state represented over the DEM with a colormap showing the flow height h. The contour of the observed arrested state is represented by a thin white line.
line. Some discrepancies remain between the simulation and the observation. The final deposit is overestimated in different places. Also, the widths of the different flow branches are not always well predicted and some flow bifurcations do not appear in the simulation.
There are many explanations to understand these discrepancies. The main problem for the simulation is the lack of data concerning the 2010 lava flow, especially concerning the vent positions and their corresponding flow rates during the time. These missing data was here estimated and we observed that different choices slightly influence the final result.
A second source of error is the accuracy of the present DEM: its 5 meter horizontal resolution should be reported to the characteristic flow height that is of about 1 meter.
At this scale, a lot of ground details are lost and a single block of rock or a little relief with a size of 2 or 3 meters is able to generate a bifurcation while it is invisible on the present DEM. In 2012, a new DEM of the Piton de la Fournaise has been realized with a one meter spatial resolution: it can very interesting to test it for future eruptions. A third source of errors could be due to the model itself: the vertical variation of the temperature is here neglected in the viscosity and the yield stress temperature dependent functions, since only a vertical averaged value of the temperature is used. A future research direction will be to reconsider this vertical variation. Nevertheless, considering theses lack of physical data and our simplified shallow viscoplastic flow model, the numerical simulations are in relatively good agreement with observations. Finally, Figs.10and11shows the evolution of the volcanic lava flow simulation for various times until the arrested state.
Figure 10: Simulation of the volcanic lava flow att= 5h and10h (left) heighth; (right) vertical-averaged temperatureθ. The contour of the real covered zone is represented by a thin black line.
Figure 11: Simulation of the volcanic lava flow at 15h and25h: (left) heighth; (right) vertical-averaged temperature θ. The contour of the real covered zone is represented by a thin black line.
Conclusion
A new reduced model for the non-isothermal free-surface tridimensional viscoplastic fluid for shallow flows on a general topography has been presented in this paper. Both the consistency and the yield stress are supposed temperature-dependent. The governing equations was numerical approximated by an autoadaptive finite element method, al- lowing to track accurately the front position. The proposed approach was evaluated by comparing numerical predictions with available data for both a laboratory experiment with a non-isothermal silicon oil flow and a real volcanic lava flow. These two tests showed that the numerical simulations are relatively in good agreement with observa- tions. Future works will consider a new 1 meter resolution DEM of the Piton de la Fournaise volcano for more accurate simulations and also to explore lava flow simulation for others volcanoes. Another research direction is to improve our model, by tacking into account the vertical dependence of the temperature in the viscosity and yield stress functions.
References
[1] N. J. Balmforth, R. V. Craster, A. C. Rust, and R. Sassi. Viscoplastic flow over an inclined surface. J. Non-Newt. Fluid Mech., 139(1):103–127, 2006.
[2] N. J. Balmforth, R. V. Craster, and R. Sassi. Dynamics of cooling viscoplastic domes. J. Fluid Mech., 499:149–182, 2004.
[3] A. J. C. Barré de Saint-Venant. Théorie et équations générales du mouvement non permanent des eaux courantes. Comptes Rendus des séances de l’Académie des Sciences, Paris, France, Séance 17, 73:147–154, 1871.
[4] D. Bercovici and J. Lin. A gravity current model of cooling mantle plume heads with temperature-dependent buoyancy and viscosity. J. Geophys. Res., 101(B2):3291–
3309, 1996.
[5] N. Bernabeu, P. Saramito, and C. Smutek. Numerical modeling of non-newtonian viscoplastic flows: part II. Viscoplastic fluids and general tridimensional topogra- phies. Int. J. Numer. Anal. Model., 11(1):213–228, 2013.
[6] E. C. Bingham. Fluidity and plasticity. Mc Graw-Hill, New-York, USA, 1922.
[7] H. S. Carslaw and J. C. Jaeger. Conduction of heat in solids. Oxford University Press, UK, second edition, 1946.
[8] A. Costa and G. Macedonio. Numerical simulation of lava flows based on depth- averaged equations. Geophys. Res. Lett., 32(5), 2005.
[9] M. L. Damiani, G. Groppelli, G. Norini, E. Bertino, A. Gigliuto, and A. Nucita. A lava flow simulation model for the development of volcanic hazard maps for mount Etna (Italy). Comput. Geosci., 32(4):512–526, 2006.
[10] F. Garel. Modélisation de la dynamique et du refroidissement des coulées de lave : utilisation de la télédétection thermique dans la gestion d’une éruption effusive. PhD thesis, Institut de physique du globe, Paris, 2012.
[11] F. Garel, E. Kaminski, S. Tait, and A. Limare. An experimental study of the surface thermal signature of hot subaerial isoviscous gravity currents: Implications for thermal monitoring of lava flows and domes.J. Geophys. Res., 117:B02205, 2012.
[12] F Garel, E Kaminski, S Tait, and A Limare. An experimental study of the sur- face thermal signature of hot subaerial isoviscous gravity currents: implications for thermal monitoring of lava flows and domes. J. Geophys. Res., 117:B02205, 2012.
[13] J.-F. Gerbeau and B. Perthame. Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete and continuous dynamical systems - serie B, 1(1):89–102, 2001.
[14] A. J. Harris and S. Rowland. FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel. Bull. Volcanol., 63(1):20–44, 2001.
[15] W. H. Herschel and T. Bulkley. Measurement of consistency as applied to rubber- benzene solutions. volume 26, pages 621–633, 1926.
[16] M. Hidaka, A. Goto, S. Umino, and E. Fujita. VTFS project: development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification. Geochem. Geophys. Geosyst., 6:Q07008, 2005.
[17] H. E. Huppert. The propagation of two-dimensional and axisymmetric viscous grav- ity currents over a rigid horizontal surface. J. Fluid Mech., 121:43–58, 1982.
[18] K. F. Liu and C. C. Mei. Slow spreading of a sheet of Bingham fluid on an inclined plane. J. Fluid Mech., 207:505–529, 1989.
[19] A. R. McBirney and T. Murase. Rheological properties of magmas.Ann. Rev. Earth Planet. Sci., 12:337–357, 1984.
[20] C. C. Mei and M. Yuhi. Slow flow of a Bingham fluid in a shallow channel of finite width. J. Fluid Mech., 431:135–159, 2001.
[21] H. Pinkerton and G. Norton. Rheological properties of basaltic lavas at sub-liquidus temperatures: laboratory and field measurements on lavas from mount Etna. J.
Volcanol. Geotherm. Res., 68(4):307–323, 1995.
[22] N. Roquet and P. Saramito. An adaptive finite element method for Bingham fluid flows around a cylinder. Comput. Meth. Appl. Mech. Engrg., 192(31):3317–3341, 2003.
[23] G. Roult, A. Peltier, B. Taisne, T. Staudacher, V. Ferrazzini, and A. di Muro. A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985–2010 period. Search and analysis of short-term precursors from a broad-band seismological station. J. Volcanol. Geotherm. Res., 241:78–104, 2012.
[24] P. Saramito. Efficient C++ finite element computing with Rheolef. CNRS and LJK, 2013. http://cel.archives-ouvertes.fr/cel-00573970.
[25] P. Saramito and N. Roquet. An adaptive finite element method for viscoplastic fluid flows in pipes. Comput. Meth. Appl. Mech. Engrg., 190(40):5391–5412, 2001.
[26] H. R. Shaw. Rheology of basalt in the melting range. J. Petrology, 10(3):510–535, 1969.
[27] W. Spataro, M. V. Avolio, V. Lupiano, G. A. Trunfio, R. Rongo, and D. d’Ambrosio.
The latest release of the lava flows simulation model SCIARA: first application to mt Etna (Italy) and solution of the anisotropic flow direction problem on an ideal surface. Procedia Computer Science, 1(1):17–26, 2010.
[28] S. Tarquini and M. Favalli. Mapping and DOWNFLOW simulation of recent lava flow fields at mount Etna. J. Volcanol. Geotherm. Res., 204(1):27–39, 2011.
[29] A. Vicari, H. Alexis, C. del Negro, M. Coltelli, M. Marsella, and C. Proietti. Model- ing of the 2001 lava flow at Etna volcano by a cellular automata approach.Environ.
Model. Softw., 22(10):1465–1471, 2007.
[30] N. Villeneuve, D. R. Neuville, P. Boivin, P. Bachelery, and P. Richet. Magma crystallization and viscosity: a study of molten basalts from the Piton de la Fournaise volcano (La Réunion island). Chemical Geology, 256(3):242–251, 2008.
[31] M. Yuhi and C. C. Mei. Slow spreading of fluid mud over a conical surface. J. Fluid Mech., 519:337–358, 2004.
