Lava flow in a channel with a bifurcation

Pbys. Chem. Earth (A), Vol. 24, No. 1 l-12, pp. 953-956,

1999 0 1999 Elsevier Science Ltd All rights reserved 1464-1895/99/$ - see front matter

Pergamon

PII:

S1464-1895(99)00141-6

Lava Flow in a Channel with a Bifurcation G. Macedonia’ and A. Longo’ ‘Osservatorio Vesuviano, Via Manzoni, 249, I-80123 Napoli, Italy *Dipartimento Scienze della Terra, Universiti degli Studi di Pisa, Via S. Maria, 53 I-56126 Pisa, Italy Received 21 May 1997; accepted 05 July 1999

Abstract. Etna poses risks to inhabited areas with its frequent effusive eruptions. During the 1989 and 1991 eruptions the Italian Department of Civil Protection diverted the lava from its natural path into an artificial channel, reducing the risk of lava inundation. The intervention resulted in the creation of a bifurcation between the natural and the artificial channels. In this paper magma dynamics in the bifurcation is investigated by solving the equations of mass and momentum balance with a simplified two-dimensional geometry, describing magma as an homogeneous, isothermal incompressible fluid with Newtonian rheology. Results show the important role played by the slope of the artificial channel and the effect of the width of the artificial mouth on the efficiency of the diversion. 0 1999 Elsevier Science Ltd. All rights reserved 1

Introduction

Effusive eruptions on Etna produce lava flows which can pose hazard to several towns. During the 1989 and 1991 eruptions a diversion of lava from its natural path was needed in order to reduce the hazard of lava invasion. The method of diversion adopted by the Italian Department of Civil Protection agency consisted in the creation of an artificial channel lateral to the natural one and in the breaking of the wall separating the two channels using explosives, as reported in Barberi et al. (1992). The resulting shape of the system was a channel with a bifurcation as shown in Fig. 1. The breaking of the separation wall allowed magma to flow into the artificial channel with a decrease in the mass flow rate in the natural channel. In order to study the effectiveness of such an intervention, a numerical model was developed to investigate the dynamics of the flow through the natural and the artiCorrespondence

to:

viano, Via Manzoni,

G. Macedonia,

Osservatorio 249, I-89123 Napoli, Italy

Vesu-

ficial branches of the bifurcation. The fluid motion is described by the mass and momentum balance equations, which are solved by the Finite Element Method (FEM) adopting a Streamline-Upwind/Petrov-Galerkin (SU/PG) scheme. Lava is treated as an incompressible, isothermal, homogeneous, and constant viscosity fluid with Newtonian rheology. This last assuption on magma rheology, strongly simplifies our problem by the numerical point of view. Although more realistic nonNewtonian reologies were adopted in the past by many authors, to model lava flows (eg. Dragoni et al. (1986): Crisp and Baloga (1990); etc.), previous works considered only very simple geometries such as an inclined plane. In this work, on the other hand, we focused our attention on the effects on lava flows produced by more complex channel shapes. Simulations were performed to study the efficiency of the diversion as a function of the geometry of the system.

2

The physical

model

The dynamics of lava flow in a channel with a bifurcation can be described by a set of partial differential equations describing the mass and momentum balance: v.v = 0

(1)

p~+p(v-v)v+vp-v~-~~=.0

(2)

where t is the time, p the pressure, p is the magma density, 1-1the kinematic viscosity, g is gravity force, V is the gradient operator. The viscous drag is expressed by the stress tensor 7 defined as: 7=/J(V7J+KW) The fluid was assumed to be incompressible, and with Newtonian rheology.

953

(3) isothermal,

954

G. Macedonio and A. Longo: Lava Flow in a Channel with a Bifurcation

The simplified geometry of the system is shown in Fig. 1. The width of the natural and the artificial channels are w and W respectively, whereas the width of the artificial mouth is L. The slopes of the natural and the artificial channels are LYand 0, and the components of the gravity acceleration along the two channels are gsincY and gsing. In this work, the efficiency of the diversion was investigated as a function of the shape of the system, defined by the above parameters. The computational domain is two-dimensional and a Cartesian coordinate system was used. We assumed that magma completely fills both the channels and has an infinite extension along the vertical direction; this is equivalent to neglecting the drag at the floor of the channels and to setting free-slip boundary conditions at the bottom of the channels. The vertical walls of the channels are taken as rigid and impermeable with no-slip conditions on the velocity (V = 0). An appropriate gravitational field (gsincr and gsing) is set in each channel to reproduce its slope. Conditions of free mass flow are imposed at the inlet of the natural channel, and at the outlet of both the natural and the artificial channels. The velocity field at the initial instant (t = 0) is everywhere null, and initial null pressure field are is imposed throughout the domain. Equations (1) and (2) are solved with a Fortran code based on the Finite Element Method (FEM) and adopting the Streamline-Upwind/Petrov-Galerkin (SU/PG) formulation as suggested by Brooks and Hughes (1982). The two-dimensional computational domain is discretized in rectangular elements (four nodes) and an isoparametric description of the velocity field is adopted using bilinear interpolating functions. The interpolating functions are continuous across element boundaries. The pressure field, in contrast, is assumed to be constant inside each element, resulting in a discontinuous function across elements. This scheme is particularly suitable for the description of fully incompressible fluids. Spatial numerical integration of mass and momentum balance equations (1,2) was performed using a Gauss 2x2 quadrature, yielding to a system of algebraic equations:

Ma+Cv+iV(v) Gtv - D = 0

-Gp-F=O

(4) (5)

In the above equations, v and a are vectors that represent the nodal values of the velocity V; and of their time derivatives ai = iti. The vector p contains the value of the pressure at element centers, M is the consistent generalized mass matrix, and C accounts for the diffusive term. N(v) represents the convective term and, since it is non linear, the elements of the matrix are functions of the unknown velocities. G represents the gradient matrix, whereas Gt is the transposed divergence matrix. F is a generalized force vector and D accounts

W

W

/

/

g sinp

gsina

T

-

f

‘/ I

L

/

W I /

Fig.

1. Sketch of the bifurcating lava channel

for the prescribed boundary velocity on the continuity equation. The conservation equations (4) and (5) are solved using the Predictor/Multicorrector transient algorithm as suggested by Brooks and Hughes (1982). The integrations in time starts from the velocity and pressure fields prescribed at t = 0. The solution method is explicit in the velocity and implicit in the pressure, this last being computed through the solution of a Poisson equation. Convergence of the algorithm is assured by a control criterion based on the mass conservation. A uniform grid was used, with a resolution of 0.5 m, the total number of elements being about 500. The time step varied from 0.01 s to 0.005-s and the maximum number of iterations necessary for the convergence of the algorithm ranged from 2 to 6.

3

Results

The FEM code is used to perform a parametric study on the width L of the artificial mouth and the slope p of the artificial channel. The width w and the slope a of the natural channel, and the width W of the artificial channel are kept fixed. The other input quantities that are maintained constant are:

955

G. Macedonio and A. Longo: Lava Flow in a Channel with a Bifurcation 0.50

1H 0.45 0 $ k$ ;

0.40

I 0.35

3.0

4.0

5.0

6.0

7.0

8.0

9.0

-I

10.0

Width of the artificial mouth (m)

Fig. 3. Variations of Qart./Qtot. as a function of the mouth of the artificial channel.

Fig. 2. from

Velocity

the

field distribution

beginning

P

=

P w

= zz

of the

2700 100 5 8”

the bifurcation

m

The investigated = =

intervals of L and p are:

2+8 m 50~120

The aim of the parametric study is to analyze the effects of the width L of the artificial mouth and the slope /3 of the artificial channel on the efficiency of the diversion Q, that is defined as the ratio between the mass flow rate per unit vertical length through the artificial channel and the total mass flow rate per unit vertical length at the entrance of the natural channel: Q

=

Qartificiat Qtotu1

=

I,,

at 1100 s

simulation

kg/m3 Pas m

;I,

L p

inside

of the width,

Sk” V(x)dx Sd” v(xPx

where V(z) and v(z) are the velocities along the outlet of the artificial channel and along the inlet of the natural channel respectively. The velocity field in the bifurcation resulting at t=llOO seconds from the beginning of the calculation is shown in Fig. 2. The development of a parabolic velocity profile across the channels 40 s from the beginning of the simulation marked the steady state flow. The values obtained for the maximum velocities range in the interval of 1+5 m/s, whereas the velocity observed in the artificial channel during the 1991 eruption was about 1 m/s. It should be pointed out that the distortion of the lava streamlines during the diversion of the fluid from the natural into the artificial channel causes much viscous

dissipation, which can be minimized by increasing the width, L, of the artificial mouth. As can be seen in Fig. 2, the formation of the vortex implies a narrowing of the effective width of the mouth and lava is forced to flow through a reduced section of the channel, decreasing the efficiency of the diversion. Variations of Q as a function of L are depicted in Fig. 3. It should be observed that the efficiency of the diversion can be increased by enlarging the width of the artificial mouth. Results shown in Fig. 3 are obtained in a system with cy = fi = 10” and w = W = 5 m. The dependence of Q on the width of the artificial mouth is not strong, however we observe an increase in the efficiency as the width of the mouth, L, increases. The maximum efficiency of the diversion for channels of equal width is l/2 and, as shown in Fig. 3, it is approached when L is about 10 m, that is twice the width of the channels (w = W = 5 m). Further increases of L will result in little change in Q. Results of the parametric study of Q as a function of the slope /I of the artificial channel are shown in Fig. 4. In this case the width of the channels and the artificial mouth are set equal (w = W = L = 5 m), with cy = 8”. It can be observed that the dependence of Q on /I is strong and variations of only a few degrees in the slope of the artificial channel result in a great change in the efficiency of the diversion: as an example, doubling the angle of inclination p from 6’ up to 12’ produces a doubling in the diversion efficiency.

4

Conclusions

The fluid dynamics of lava flow in a channel with a bifurcation was simplified and numerically simulated. This problem has practical interest for the planning of lava diversions such as those performed at Etna by the Italian Civil Protection Department during the 1989 and 1991

G. Macedonio and A. Longo: Lava Flow in a Channel with a Bifurcation

956 0.6

diversion (Fig. 4), and that the viscous drag at the artificial channel mouth results in the narrowing of the effective section of the channel (Fig. 2). The dissipation can be minimized by increasing the width L of the artificial mouth about twice the width of the artificial channel (Fig. 3). The adopted model for describing the lava flow is based on simplifying assumptions regarding magma rheology. For a more realistic description, Non-Newtonian reology models, and magma cooling need to be considered in the future.

k

5.0

6.0

7.0

6.0

9.0

10.0

Slope of the artificialchannel

Fig. 4. Variations of Q.?t./QtOt. the artificial channel.

11.0

12.0

(deg)

as a function of the

slope, /3, of

Acknowledgments. The authors wish to thank Gruppo Nazionale per la Vulcanologia, Consiglio Nazionale delle Ricerche, Italy, for financial support.

References eruptions. The model is based on the solution of the mass and momentum balance equations for an incompressible, homogeneous, constant viscosity and isothermal magma with Newtonian rheology. The conservation equations are solved with a code based on the Finite Element Method, and adopting a SU/PG formulation. The efficiency Q of the lava diversion is analyzed by a parametric study performed varying the width L of the arti&ial mouth and the slope p of the artificial channel. Results showed that the slope of the artificial channel plays a major role in the efficiency of the

Barberi, F., Carapezza, M., Valenza, M., and Villari, L., L’erwione 1991-1992 dell%tna e gli interuenti per fernore o titardare I’ouanzato della lava, GNV-CNR, Giardini, Pisa, 1992. Brooks, A. and Hughes, T., Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Meths. Appl. Mech. Engrg., 32, 199-259, 1982. Crisp, J. and Baloga, S., A model for lava flows with two thermal components, J. Gwphys. Res., 95, 1255-1270, 1990. Dragoni, M., Bonafede, M., and Boschi, E., Downslope flow models of a bingham liquid: implications for lava flows, J. Vo’olcanol. Gwtherm. Res., 30, 305-325, 1986.