Layout 6


Cooling of a channeled lava flow with non-Newtonian rheology:
crust formation and surface radiance

Andrea Tallarico1,2,*, Michele Dragoni3, Marilena Filippucci1,4, Antonello Piombo3,
Stefano Santini5, Antonella Valerio3

1 Università di Bari, Centro Interdipartimentale di Ricerca per il Rischio Sismico e Vulcanico,Bari, Italy
2 Università di Bari, Dipartimento di Geologia e Geofisica, Bari, Italy
3 Università di Bologna, Dipartimento di Fisica, Bologna, Italy
4 Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania, Osservatorio Etneo, Catania, Italy
5 Università di Urbino, Dipartimento di Scienze di Base e Fondamenti, Urbino, Italy

ANNALS OF GEOPHYSICS, 54, 5, 2011; doi: 10.4401/ag-5335

ABSTRACT

We present here the results from dynamical and thermal models that
describe a channeled lava flow as it cools by radiation. In particular, the
effects of  power-law rheology and of  the presence of  bends in the flow are
considered, as well as the formation of  surface crust and lava tubes. On
the basis of  the thermal models, we analyze the assumptions implicit in
the currently used formulae for evaluation of  lava flow rates from satellite
thermal imagery. Assuming a steady flow down an inclined rectangular
channel, we solve numerically the equation of  motion by the finite-volume
method and a classical iterative solution. Our results show that the use of
power-law rheology results in relevant differences in the average velocity
and volume flow rate with respect to Newtonian rheology. Crust formation
is strongly influenced by power-law rheology; in particular, the growth
rate and the velocity profile inside the channel are strongly modified. In
addition, channel curvature affects the flow dynamics and surface
morphology. The size and shape of  surface solid plates are controlled by
competition between the shear stress and the crust yield strength: the
degree of  crust cover of  the channel is studied as a function of  the
curvature. Simple formulae are currently used to relate the lava flow rate
to the energy radiated by the lava flow as inferred from satellite thermal
imagery. Such formulae are based on a specific model, and consequently,
their validity is subject to the model assumptions. An analysis of  these
assumptions reveals that the current use of  such formulae is not consistent
with the model.

Introduction
Laboratory studies on basaltic melts show that lava

rheology can show nonNewtonian behavior under certain
conditions, which include vesicularity [Spera et al. 1988, Stein
and Spera 1992, Badgassarov and Pinkerton 2004], crystal
concentration [Pinkerton and Stevenson 1992, Smith 2000,

Sonder et al. 2006, Champallier et al. 2008], and temperature
and shear rates [Shaw et al. 1968]. The literature indicates
that magmas have nonNewtonian, pseudoplastic behavior,
with the exception of Smith [2000], who attributes a dilatant
rheology to lava at high crystal concentrations. Pseudoplasticity
and dilatancy can be described by power-law rheology. 

When solving the problem of gravity-driven lava flow,
power-law rheology introduces nonlinearity into the
diffusion term of the momentum equation, and an analytical
solution of the governing differential equations is not
possible. This has given rise to various approximate solutions
to the problem. The fully developed laminar flow of power-
law fluids has been studied numerically using the finite
element method [Syrjala 1995] and the finite volume method
[Capobianchi 2008] for a pressure driven flow in a horizontal
rectangular duct. 

The equation of motion for gravity-driven lava flow
with a power-law rheology was solved by Filippucci et al.
[2010] using the finite volume method for the cases in which
the rheology is temperature independent. Several models
describing lava tube formation have been proposed [Greeley
1971, Peterson and Swanson 1974, Rowland and Walker
1990, Hon et al. 1994, Dragoni et al. 1995, Valerio et al. 2008,
Valerio et al. 2010]. In the present study, the cooling of a
channel of lava with power-law rheology is considered. A
two-dimensional thermal model with heat flux assigned at
the upper surface is introduced to describe lava cooling due
to radiation and convection into the atmosphere. The lava
crust is considered as a plastic body, and its rheology is
described through the introduction of a yield strength as a
function of temperature. 

Article history
Received November 11, 2010; accepted May 31, 2011.
Subject classification:
Rheology, Magmas, Thermodynamics, Lava flow, Effusion rate.

510

Special Issue: V3-LAVA PROJECT



TALLARICO ET AL.

Bends in lava flow are often observed, and these can
strongly affect the flow dynamics and formation of lava tubes
[Greely 1971, Peterson et al. 1994]. In the present study, we
consider the effects of the curvature of a channel on the lava
velocity and shear stress, and on the formation of the solid
crust at the flow surface. To allow an analytical solution of
the Navier-Stokes equation in the presence of bends, the lava
is regarded as a Newtonian, homogeneous, isotropic and
incompressible fluid [Valerio et al. 2011]. Heat radiation and
convection into the atmosphere are considered as the main
cooling mechanisms. The model analyses the effects of
curvature on the development of surface solid plates.

The availability of high-resolution thermal imagery of
active lava flows has stimulated the use of radiance maps for
evaluation of lava effusion rates [Mouginis-Mark et al. 1994,
Harris et al. 1995, Calvari et al. 2010]. This has been made
possible by the use of simple formulae that relate the lava
flow rate to the energy radiated per unit time from the
planimetric surface of the flow [Pieri and Baloga 1986, Harris
et al. 1997, Harris et al. 2005a]. Such formulae are based on
a specific flow model, and consequently, the results depend

on the flow-model assumptions. An analysis of these
assumptions was performed by Dragoni and Tallarico [2009],
which revealed that the current use of these formulae is not
consistent with the model.

Crust formation in a lava channel
with power-law rheology

The constitutive equation for a power-law fluid is:

(1)

where vij is the stress tensor, ėij is the strain rate tensor, k is
the fluid consistency, n is the power-law exponent, which is
a measure of the nonlinearity, and:

(2)

where I2 is the second invariant of the strain rate tensor. The
apparent viscosity of the fluid is given by:

(3)

If n < 1, the fluid is pseudoplastic and it thins (i.e. it
deforms more rapidly) with an increase in stress. If n=1, the
fluid is Newtonian. If n>1, the fluid is dilatant and it thickens
with an increase in stress.

An illustration of the model with the coordinate system
is shown in Figure 1. The equation of motion for a gravity-
driven flow down an inclined rectangular channel is:

(4)

where vx is the x component of velocity, g is the gravity
acceleration, and a is the slope angle. The apparent viscosity is:

(5)

511

Figure 1. Coordinate system and geometrical parameters.

Figure 2. a) Average velocity with respect to a for a fixed geometry (a = 10 m, h = 4 m, t �= 2650 kg m−3, k = 104 Pa sn); b) Average velocity with respect to
a for a fixed q (a = 10 m, q = 500 m3s−1, t = 2650 kg m−3, k = 104 Pa sn); c): Thickness h with respect to q (a = 10 m, a = 20˚, t = 2650 kg m−3, k = 104 Pa sn).

2ke e1ij
n

ij=v
-o o

2e I2=o

ke 1app
n=h -o

sin 0g x y
v

y
v
zapp

x
app

x+ + =
2
2

2
2

2
2

2
2

t a h hc `m j

k y
v

z
v2 2

app
x x

2
1n

= +
2
2

2
2

h

-

c `m j; E



512

The boundary conditions are the nonslip at the walls,
and free surface at the top of the flow:

(6)

(7)

The dynamic problem has been solved using the finite
volume method with an iterative solver [Patankar 1980]. The
numerical tests carried out by Filippucci et al. [2010]
indicated that the numerical solution, with respect to the
analytical solution, that is available for the Newtonian
rheology produces an error of about 0.03% with a grid of
51 × 51 control volumes (CV). The error due to the mesh
refinement tested with the power-law rheology, for which
the analytical solution is not available, is about 0.03% with a
grid of 51 × 51 CVs, so the grid of 51 × 51 CVs was adopted
for all of the following computations, as a good compromise

between precision and time for calculation. 
The effects of the slope a on the average velocity and

the flow rate can be remarkably different for Newtonian and
nonNewtonian fluids. Figure 2a shows the effects of a on the
average velocity, with varying n, and fixed k, channel width
a, and thickness h. Differences between the linear and the
nonlinear rheology are enhanced by increases in a. The
effects of the slope have been evaluated also for a constant
volume flow rate q (Figure 2b). The effects of nonlinearity
on vx vary with a and k. Finally, we evaluated the dependence
of the volume flow rate q on h for fixed values of a and a
(Figure 2c).

We simulate two possible scenarios of nonlinearity
(pseudoplastic, with n = 0.5, and dilatant, with n = 1.5) and
compare these with the Newtonian case (n = 1) for the
Laghetto lava channel (Figure 3a, b) and for another active
lava channel, both of which were formed during the 2001
eruption of Mount Etna (Figure 4). Our results show that the

COOLING OF CHANNELED LAVA FLOW

Figure 3. Velocity profiles of the Laghetto channel (Mount Etna, 2001). a) Horizontal profiles at z = 0. b) Vertical profiles, at y = 0. The velocity profiles
are plotted for three different values of the power-law exponent n (0.5, 1, 1.5). The parameters of the Etna channel are: a = 12 m, h = 6 m, a = 5˚ (data
from Ferlito and Sievert [2006]); the parameters of the lava are: t = 2030 kg m−3 and k = 8625 Pa sn (data from Pinkerton and Sparks [1978]).

Figure 4. Velocity profiles of the channel flow of the 2001 Mount Etna eruption. a) Horizontal profiles at z = 0. b) Vertical profiles at y = 0. The velocity
profiles are plotted for three different values of the power-law index n (0.5, 1, 1.5). The parameters of the Etna channel are: a = 3 m, h = 1.2–1.5 m, a =
7.5˚; the parameters of the lava are: t = 2030 kg m−3 and k = 1.1 × 104 Pa sn (data from Harris et al. [2005a]).

( , ) 0; 2 , 0v y h v
a zx x= =!- ` j

0, 0; , 0 0y
v

z z
v

yx x= =
2
2

2
2^ ^h h



velocity profiles and average velocities can vary significantly,
depending on the assumption of linear or nonlinear
rheology. 

We have also considered the May 1997 flow from Pu'u
'O'O (Kilauea, Hawai'i). The apparent viscosity varies from 9
to 879 Pa s, and decreases with an increase in the shear rate,
so the system is pseudoplastic [Weed et al. 1986]. The
channel geometry was described by Harris and Rowland
[2001]. We simulated the flow with pseudoplastic rheology,
setting n = 0.9, 0.7 and 0.5. Given the large variations in
channel width and the possible variations in channel
thickness, we considered two possible geometries for the
channel section. In this way, we can observe how the average
velocities are modified by a change in the aspect ratio b =
h/a of the channel. As a result, with an increase in b the
average velocity increases, and this effect is greater as n
decreases (Figure 5a). Moreover, the volume flow rate q,
which we hypothesized as constant with respect to b,
increases very rapidly as n decreases (Figure 5b). 

Finally, we considered the 1984 Mauna Loa lava
channels. We used power law rheology to fit the flow
velocity data measured from videotapes by Sakimoto and
Gregg [2001]. As a result, we found n in the dilatant range
(Figure 6a, b).

Cooling is modeled for a lava flow exiting a vent with
an effusion temperature T0 = 1,273 K. We assume that the
flow is laterally isothermal and that the cooling process
occurs from the upper surface with a fixed value of the heat
flux q0. This describes the average heat loss from the surface
that is covered by a solid crust and from the crust-free shear
regions [Valerio et al. 2008]. The energy equation governing
the heat transfer is: 

(8)

where T is the temperature and | is the thermal diffusivity:

(9)

TALLARICO ET AL.

513

Figure 5. Pu'u 'O'O (Kilauea, Hawai'i) case study. a) Average velocity va vs n for two values of b. b) Volume flow rate q vs n, with q equal for the two values
of b. The parameters of the channel are: a = 5 m, h = 3 m, a = 10˚; the parameters of the lava are: t = 2600 kg m−3 and k = 895 Pa sn.

Figure 6. Matches between the measurements of the flow velocity (crosses) and the velocity profiles (lines) computed with a dilatant model (n > 1) for
the channel flow of April 2, 1984, Mauna Loa eruption. a) Horizontal profiles at z = 0. b) Vertical profiles at y = 0. The parameters of the channel are: a =
14.6 m, h = 4.1 m, a = 8˚; the parameters of the lava are: t = 1200 kg m−3 and k = 850–1400 Pa sn (data from Sakimoto and Gregg [2001]; viscosity from
Harris and Allen III [2008]).

t
T

y
T

z
T

2

2

2

2
= +

2
2

2
2

2
2

|c m

c
K

p
=| t



514

where cp is the specific heat, and K is the thermal conductivity.
The boundary conditions are the adiabatic conditions at

the walls and the constant heat flux at the upper surface.
Initially, the fluid has a uniform temperature T0. At the time
t>0, a constant heat flux q0 is set and the thermal boundary
conditions become:

(10)

(11)

In addition to the spatial integration, the discretization
of the energy equation involves integration of Equation (8)
over the time interval from t to t + Dt. This integration is
operated by using a fully implicit scheme, which ensures
simplicity and physical consistency [Patankar 1980]. Details
of the numerical test of this solution are given below.

We substituted the time dependence of the solution in
T with the space dependence by the relation:

(12)

where:
(13)

We first considered the conductive cooling of a steady
lava flow, and the consequent surface crust formation.
Indeed, many basaltic lava channels have a surface crust that
appears dark due to its lower radiative temperature, and that
is carried along by the flow. We evaluated the effects of the
nonlinearity of the lava rheology on crust formation, by
modeling the plastic behavior of the crust using a
temperature-dependent yield strength x, given by:

(14)

where Ts is the solidus temperature, x0 is the maximum yield
strength, and b is a constant.

We call d the width of the flow surface covered by crust.
For temperatures below the solidus, crust formation and
growth are possible in the region |y|< d/2, where the
horizontal shear stress vxy < x. Shearing movements occur
and crust formation is inhibited where vxy>x. As the surface
temperature reaches the solidus, x starts to increase up to
the threshold. The distance from the vent where this happens
strongly depends on the assumption on the rheology (Figure 7a).
Simultaneously, the crust width d/a also grows, and once it
attains its maximum width, it stops widening. In particular,
for the dilatant fluid, the crust stops widening earlier in space
with respect to the Newtonian and the pseudoplastic fluid,
which will cover a much longer distance before the crust
reaches its maximum width (Figure 7b).

Models of formation of lava tubes as a consequence of
the crust welding to the channel levees have been proposed
that consider lava as a Bingham liquid [Dragoni et al. 1995]
or a Newtonian liquid [Cashman et al. 2006, Valerio et al.
2008]. We also explored how variations in the channel width,
ground slope and volume flow rate can affect the formation
of a lava tube, assuming a power-law rheology. As a result,
we find that for topographic and morphological variations
of the lava channel, the nonlinearity of the rheology has only
a slight influence on the formation of a lava tube, and the
results for the Newtonian fluid appear to also be valid for
power-law fluids.

Numerical test on the solution of the heat equation
The solution of the heat equation was tested to study

the accuracy, by comparing it with the analytical solution
with constant surface flux q0 [Turcotte and Shubert 1982]:

(15)

COOLING OF CHANNELED LAVA FLOW

Figure 7. a) Yield strength x as a function of x´ for different values of n. b) Crust coverage on the lava surface d/a as function of x´ for different values of
n (a = 10 m, h = 3 m, a = 0.2 rad, t = 2800 kg m−3, k = 104 Pa sn, cp =837 J kg

−1 K−1, K =3 W K−1 m−1, q0 =10
4 W m−2, x0 =8100 Pa, b =170). x´ is the

point on the flux direction where the surface temperature reaches the solidus, and it varies with varying n.

( , ) 0; 2 , 0z
T y h y

T a z= =!
2
2

2
2

- ` j
0, 0; , 0y

T z z
T y K

q0= =
2
2

2
2

-^ ^h h

( )t v z
x

x
= r

1 ,v a v y z dy2

2

x x
a

a

=
-

r r ^ h#

( 1)b-
1 ,

0,
e T T

T T
0 s

s

T
T
s

=
#

$
x

x -6 @) ( , ) 2 2 erfc 2T z t T K
q t

e z
t

z
0

4
an

r t
z2

= +r
|

|
-

-|
- c m



We evaluated the difference DTu between the upper
temperatures of the two solutions and the difference DT
between the average temperatures [Syrjala 1996,
Capobianchi and Wagner 2010], as given by: 

(16)

As we expected, both and DTu always decrease as
the number NCV of the control volumes increases, i.e. by
thickening of the computational mesh, which means that
the numerical temperature gradually approaches the
analytical value (Figures 8, 9, respectively). In this test,      
is always about two orders of magnitude less than DTu, and
this difference is a consequence of the cooling of the lava
channel, which at the time considered, involves only the
upper layers of the flow; the lower part remains at a
constant temperature. It can be seen that is very low also
for a coarse grid, as it reaches 0.01% for a grid 61 × 61 CVs
(Figure 8), while DTu needs a grid of 201 × 201 CVs to drop
under 0.05% (Figure 9). In Figures 8 and 9, the errors on 
and Tu due to the mesh refinement are also plotted. It can
be seen that due to mesh refinement, is a decreasing
function of NCV, starting from the grid of 51 × 51 CVs. For
Figure 10, we evaluated the trend of with NCV at
different times. This test indicates that decreases as the
time increases, which indicates that the discretization error
does not propagate with time and the solution is
convergent.

Crust formation in the presence of a bend in the channel
Bends in lava flows arise from local variations in the

topography, and they are commonly observed in volcanic
fields. The curvature of a channel affects the flow dynamics
and morphology [Greeley 1971].

An analytical model was proposed by Valerio et al.
[2011] to explain the effects of curvature on the velocity and
viscous stress. In a bending channel, a cylindrical coordinate
system is assumed (Figure 11). A flow segment with a
constant curvature is limited by arcs of concentric
circumferences, with their centers in the origin of the
coordinate system. We studied the equation of motion for a
viscous, homogenous, isotropic, Newtonian fluid in the
gravity field. In the near-vent portion of the flow, neglecting
nonNewtonian effects is a reasonable assumption, as a
consequence of the high temperature. The constant density
t and viscosity h are used. Lava flows in a rectangular
channel with a constant slope. The steady-state motion is
assumed to be unidirectional in the azimuthal direction: the
flow is parallel to the levees. This implies that no radial
component of the gravity force acts on it. This condition
often applies when a channel changes direction as a
consequence of local topography. Further assumptions are
introduced: the velocity only depends on the radial

TALLARICO ET AL.

515

Figure 8. Numerical tests on the heat conduction equation. Percentage
errors     as function of the number NCV of control volumes of the
computational mesh at 301 s of cooling. Solid line,    between the
numerical and analytical solutions. Dashed line, between the actual
grid and the coarser grid computed with a grid step of 10 × 10 CVs.

Figure 9. Numerical tests on the heat conduction equation. Percentage
errors DTu as function of the number NCV of control volumes of the
computational mesh at 301 s of cooling. Solid line, DTu between the
numerical and analytical solutions. Dashed line, DTu between the actual
grid and the coarser grid computed with a grid step of 10 × 10 CVs.

Figure 10. Numerical tests on the heat conduction equation. Percentage
errors     as a function of the number NCV of control volumes of the
computational mesh at different times t = 101, 201, 301 s of cooling.

1 ( )T h T z dz0

h

=
-

r #
TDr

TDr

TDr

TDr

TDr
TDr

Tr

TDr
TDr

TDr

TDr



516

coordinate, and the pressure gradient is negligible. This
obtains equation of motion as:

(17)

where u is the velocity, and:

(18)

By solving the equation with no slip boundary
conditions at the lateral walls of the channel, we obtain the
velocity as a function of the radius at the surface of the flow:

(19)

where:
(20)

The flow is not symmetric with respect to rc, the center
of the channel. The gradient of the velocity is greater close
to the levee with the higher curvature. This behavior is
particularly noticeable in comparison with the velocity in a
rectilinear channel (Figure 12), which is calculated
considering analogous hypotheses: the unidirectional flow
of a Newtonian fluid in a gravity field that depends on the
horizontal coordinate, while dependence on depth and
direction of flow is neglected. In this case, the velocity is
symmetric with respect to the center of the channel, where
it reaches its maximum. The maximum velocity in the
curved channel is shifted toward the internal levee. The
extent of the displacement from rc is a decreasing function
of the radius of curvature, which shows that the asymmetry
of the velocity is higher for sharp, rather than smooth, bends,
and it increases with channel width (Figure 13).

We calculate the viscous stress as a function of the
radius from the constitutive equation of a Newtonian fluid,
obtaining:

(21)

In analogy with the velocity gradient, the shear stress
reaches greater values in the region close to the internal
levee. In particular, it reaches its maximum in the internal
levee. We compared the viscous stress calculated in a
curvilinear channel with the stress in a rectilinear channel
(Figure 14). Here, the curvature produces a reduction in
stress close to the external levee, and an increase close to the
internal levee. The maximum shear stress is calculated as a
function of rc: it is a decreasing function of the radius of

COOLING OF CHANNELED LAVA FLOW

Figure 11. Sketch of the model for bends and the coordinate system.

Figure 12. Comparison between the velocity as a function of the
horizontal coordinate in a rectilinear channel (dashed line), and in a
curvilinear channel (solid line). The vertical line indicates the center of
the channel rc.

Figure 13. Distance between the point where the velocity is maximum
and the center of the channel rc, as a function of the radius of curvature
of the channel. Each curve is calculated for a different value of channel
width a.

1 0dr
d u

r dr
du

2

2
+ + =c

sin
g

=c h
t

a

( ) 4 ( ) lnu r r r c r
r2

1
2

1
= +

c
- -

4 ln
c

r r
r r

1

2
1
2

2

2=
c -

( ) 4 1 lnr r r
r

r r
rc1 2

1
2

1
= + +v h

c
- -c `m j; E



curvature, and reaches higher values with sharp bends. 
We consider radiation and convection into the

atmosphere to be the dominant mechanisms in the cooling
process. We calculate iteratively the surface heat-flow density
as a function of temperature along the down-flow direction. 

(22)

We use the model of cooling of a conductive half-space,
described by the time-dependent heat equation for the heat
flux [Turcotte and Schubert 1982]. 

Substituting the Fourier law in the solution, and on the
assumption that the temperature is initially uniform in the
medium, we obtain the temperature. Substituting the time
t according to Equation (12), we obtain the flow
temperature as a function of the coordinates z and L, which
represent the vertical coordinate and the distance covered,
respectively:

(23)

where T0 is the initial uniform temperature, and ū is the
average velocity.

As long as the lava cools, crust platforms form at the
flow surface. These platforms of solidified lava are driven by
the underlying fluid, and laterally they are limited by the
action of the shear stress. Greeley [1971] described some
effects of the curvature of a channel on the lava surface
morphology in Hawaiian volcanic fields. The crust formed
at a distance of a few meters from the eruptive vent. At the
first bend in the channel, this crust breaks into plates, which
float on the flow surface. The increased velocity gradient
provokes fragmentation of the plates and their inclusion in
the flow through remelting. Peterson et al. [1994] also noted

that the crust plates break in the lateral flow region when
they run across a curve. 

A common phenomenon in bends is the presence of a
different level of deposit on the two levees of the channel,
which can be observed when the eruption ends and the
lava drains. The level of deposits is considerably higher on
the external levee rather than on the internal levee, as
observed on Kilauea volcano, Hawaii [Heslop et al. 1989]
and on Pico Partido volcano, Lanzarote [Woodcock and
Harris 2006]. The centripetal force provokes only a slight
difference in height between the two lateral walls. Its
effects are evaluated from the radial component of the
Navier-Stokes equation, which considers the variation in
pressure due to the variation in height in the flow surface.
The calculated superelevation explains an increase in the
height lower than 10% of the superelevation inferred from
the data. The present model applies to cases where the
flow is parallel to the levees, whereas the measured
different amounts of deposits occur in flows where the
radial component of the gravity force acts in bends, due to
the local topography.

Surface radiance
and implications for satellite thermal imagery

The effusion rate of the lava from an eruption vent is
the primary quantity that controls the evolution of the lava
flow that ensues. For this reason, a lot of effort has been
devoted to the evaluation of this quantity, which involves
measurement of the flow velocity and the cross-sectional area
[e.g. Calvari et al. 2002]. The knowledge of effusion rates also
has a major role in real-time simulations of lava flow paths
that are carried out when lava flows threatens inhabited areas
[e.g. Ishihara et al. 1990, Miyamoto and Sasaki 1997, Vicari et
al. 2007, Rongo et al. 2008, Hérault et al. 2009]. 

Direct measurements of effusion rates in the field are
difficult, and calculations from other flow parameters would
require the measurement of such parameters [Tallarico et al.
2006], which can be equally problematic. In recent years, the
availability of thermal images of volcanoes during eruptions
has stimulated the evaluation of lava effusion rates from the
measurement of their heat radiation. The thermal data
obtained from high-resolution radiometers mounted in
Earth satellites were first used to monitor active lavas by
Mouginis-Mark et al. [1994] and Harris et al. [1995]. Thermal
images from hand-held radiometers on the ground have also
been used [Harris et al. 2005b]. 

The use of radiance maps for the evaluation of lava
effusion rates is made possible by simple formulae that relate
the lava flow rate to the energy radiated per unit time from
the planimetric surfaces of the flow. Such formulae are based
on a specific flow model, and consequently, their validity is
subject to the model assumptions. The evaluation of lava
effusion rates is based on a formula originally proposed by

TALLARICO ET AL.

517

Figure 14. Comparison between the viscous stress as a function of the
horizontal coordinate in a rectilinear channel (dashed line), and in a
curvilinear channel (solid line). The vertical line indicates the center of
the channel rc.

( ) ( ) ( )q T E T T c T T0
4 4

a conv a= +R - -

( , )
2

2 erfc 2T z L T K
q

u
L

e z z L
u

0
0 4 L

z u2

=
r

|
|- -

|
-

r
rr c m; E



518

Pieri and Baloga [1986] that relates the flow rate to the
planimetric area of a flow:

(24)

where q is the volume flow rate, W is the heat radiated per
unit time from the flow surface, and DT is the difference
between the lava temperatures at the vent and at the front.
If the flow rate is assumed to coincide with the effusion rate,
this formula states that the effusion rate is proportional to
the heat radiated per unit time by the surface of the flow.
Harris et al. [1997] modified the formula, by including the
thermal contribution W1 of convection in the air, and the
degree {0 of crystallization of the lava. The contribution W2
of heat conduction to the ground was included later,
obtaining [Harris et al. 2005a]:

(25)

where cL is the latent heat of crystallization.
The radiance of a hot body is related to the surface

temperature at a given instant of time. If the body is a
flowing liquid, the connection between the radiated heat and
the flow rate is not straightforward. The singularity of
Equation (24) was noted by Wright et al. [2001]. They
suggested that the apparent success of the formula in giving
reasonable values of effusion rates is due to numerical
coincidence between a constant by which Harris et al. [1997]
multiply space-based estimates of lava flow area, and the
empirical ratio between effusion rates and flow areas
reported by Pieri and Baloga [1986] for a suite of Hawaiian
flows. Harris et al. [2007] asserted that the technique gives
time-averaged discharge rates.

Dragoni and Tallarico [2009] proposed a different
explanation. Indeed, the formula of Pieri and Baloga [1986]
can only be obtained in the framework of a lava flow model
based on several simplifying assumptions. They explicitly
addressed such assumptions, to ascertain which of them are
acceptable approximations of real lava flows, and which
instead impose strong limits to the applicability of the
model. They also considered whether current use of the
formula is consistent with the model itself. The analysis of
the thermal model from which the formulae by Pieri and
Baloga [1986] and Harris et al. [1997] are obtained shows that
the measurement of only the total instantaneous heat flow
from the flow surface is not sufficient to calculate the lava
flow rate. It is thus necessary to provide further information,
i.e., the difference between the lava temperatures at the vent
and at the front. If this difference is assumed to be the same
for all lava flows, as is currently assumed, then the formula
states that the effusion rate is an increasing function of flow
length and a decreasing function of flow thickness, which is
unrealistic. Apparently reasonable results are obtained as the
assumption of a uniform crust temperature is compensated

for by the inconsistent use of the formulae. In the real world,
the average surface temperature is lower for longer flows,
and is higher for thicker flows. The measured heat flow
incorporates these effects, which happen to counterbalance
the use of a constant temperature difference between the
vent and the front of the flow.

Conclusions
The models presented in the present study indicate the

relevance of thermal and rheological processes on the
dynamics of channeled lava flows. In particular, the effects
of nonlinearity on the velocity and the flow rate of a gravity-
driven lava flow are not negligible, especially for steep slopes.
Moreover, the effects of nonlinearity on the average velocity
are enhanced by an increase in the fluid consistency
parameter. Generally, the average velocity tends to increase
with the aspect ratio of the channel, and this effect grows by
decreasing the exponent n of the power-law rheology.

The rapidity of crust formation due to the constant heat
flow at the channel surface strongly depends on the degree
of non-linearity of the rheology. The crust grows, attains its
maximum width, and stops widening, depending on the
rheology. While for the dilatant fluid, the crust stops
widening very close to the point where the temperature
reaches solidus, a pseudoplastic fluid will cover a much
longer distance before the crust reaches its maximum width.
This result is a consequence because with our assumptions,
the pseudoplastic fluid flows with a greater velocity than the
Newtonian fluid, which in turn flows with a greater velocity
than the dilatant fluid. 

Numerical tests were carried out with the purpose of
testing the stability and accuracy of the solution of the
thermal equation. We find that the numerical error decreases
as the number of cells of the mesh increases, and it does not
propagate with time.

The presence of bends in the lava path significantly
affects the flow dynamics. A simplified form of the equation
of motion for a Newtonian fluid has been solved, which
analytically obtained the velocity and stress fields in a
cylindrical coordinate system as a function of the radial
coordinate. Fluid velocity and viscous stress show
asymmetric behavior with respect to the center of the
channel. The maximum of the velocity is shifted towards the
internal levee, in proximity of which the shear stress also
reaches higher values. Although we use strong assumptions
to simplify the dynamic, rheological and thermal aspects of
lava flows, the model provides a possible explanation of
some field observations. An analysis was carried out of the
formulae currently used to evaluate the effusion rate of lava
flows on the basis of the measurement of the instantaneous
heat flow from the flow surface. The analysis of the model
from which the formulae are derived shows that further
information should be supplied, i.e. the difference between

COOLING OF CHANNELED LAVA FLOW

q c T
W
p

=
t D

( )q c T c
W W W1

p L 0

2= +
+ +

t {D



the lava temperatures at the vent and at the front. If this
difference is assumed to be the same for all lava flows, as is
currently assumed, the formula states that the effusion rate
is an increasing function of the flow length and a decreasing
function of flow thickness, which is not realistic.

Consistent application of the models presented here to
actual eruptions needs reliable data from both field and
laboratory measurements of viscosity, which deserve much
more effort in the future.

Acknowledgements. This study benefited from funding provided
by the Italian Presidenza del Consiglio dei Ministri Dipartimento della
Protezione Civile (DPC), Scientific papers funded by the DPC do not
represent its official opinion and policies. The authors wish to thank Ciro
Del Negro, Pierfrancesco Dellino, and an anonymous reviewer for their
very useful comments on the manuscript.

References
Bagdassarov, N. and H. Pinkerton (2004). Transient phe-

nomena in vesicular lava flows based on laboratory ex-
periments with analogue materials, J. Volcanol. Geoth.
Res., 132, 115-136.

Calvari, S., M. Neri and H. Pinkerton (2002). Effusion rate
estimations during the 1999 summit eruption on Mount
Etna, and growth of two distinct lava flow fields, J. Vol-
canol. Geoth. Res., 119, 107-123.

Calvari, S., L. Lodato, A. Steffke, A. Cristaldi, A.J.L. Harris, L.
Spampinato and E. Boschi (2010). The 2007 Stromboli
eruption: event chronology and effusion rates using ther-
mal infrared data, J. Geophys. Res., 115, B04201; doi:
10.1029/2009JB006478.

Capobianchi, M. (2008). Pressure drop predictions for lami-
nar flows of extended modified power law fluids in rec-
tangular ducts, Int. J. Heat Mass Tran., 51, 1393-1401.

Capobianchi, M. and D. Wagner (2010). Heat transfer in lam-
inar flows of extended modified power law fluids in rec-
tangular ducts, Int. J. Heat Mass Transfer, 53, 558-563.

Cashman, K.V., R.C. Kerr and R.W. Griffiths (2006). A labo-
ratory model of surface crust formation and disruption
on channelized lava flows, B. Volcanol., 68, 753-770.

Champallier R., M. Bystricky and L. Arbaret (2008). Experi-
mental investigation of magma rheology at 300 MPa:
from pure hydrous melt to 75 vol. % of crystals, Earth
Planet. Sc. Lett., 267, 571-583.

Dragoni, M., A. Piombo and A. Tallarico (1995). A model for
the formation of lava tubes by roofing over a channel, J.
Geophys. Res., 100, 8435-8447.

Dragoni, M. and A. Tallarico (2009). Assumptions in the eval-
uation of lava effusion rates from heat radiation, Geo-
phys. Res. Lett., 36, L08302; doi: 10.1029/2009GL037411.

Ferlito, C. and J. Siewert (2006). Lava channel formation dur-
ing the 2001 eruption on Mount Etna: evidence for me-
chanical erosion, Phys. Rev. Lett., 96, 028501; doi:
10.1103/Phys Rev Lett. 96.028501.

Filippucci, M., A. Tallarico and M. Dragoni (2010). A three-
dimensional dynamical model for channeled lava flow
with non-linear rheology, J. Geophys. Res., 115, B05202.

Greeley, R. (1971). Observations of actively forming lava
tubes and associated structures, Hawaii, Modern Geol-
ogy, 2, 207-223.

Harris, A.J.L., R.A. Vaughan and D.A. Rothery (1995). Vol-
cano detection and monitoring using AVHRR data: the
Krafla eruption, 1984, Int. J. Remote Sens., 16, 1001-1020.

Harris, A.J.L., S. Blake and D.A. Rothery (1997). A chronol-
ogy of the 1991 to 1993 Mount Etna eruption using ad-
vanced very high resolution radiometer data:
implications for real-time thermal volcano monitoring, J.
Geophys. Res., 102, 7985-8003.

Harris, A.J.L. and S.K. Rowland (2001). FLOWGO: a kine-
matic thermo-rheological model for lava flowing in a
channel, B. Volcanol., 63, 20-44.

Harris, A.J.L., J. Bailey, S. Calvari and J. Dehn (2005a). Heat
loss measured at a lava channel and its implications for
down-channel cooling and rheology, In: M. Manga and
G. Ventura (eds.), Kinematics and Dynamics of Lava
Flows, Geological Society of America Special Paper, 396,
125-146; doi: 10.1130/2005.2396(09).

Harris, A.J.L., J. Dehn, M. Patrick, S. Calvari, M. Ripepe and
L. Lodato (2005b). Lava effusion rates from hand-held
thermal infrared imagery: an example from the June 2003
effusive activity at Stromboli, B. Volcanol., 68, 107-117.

Harris, A.J.L., J. Dehn and S. Calvari (2007). Lava effusion
rate definition and measurement: a review, B. Volcanol.,
70, 1-22.

Harris, A.J.L. and J.S. Allen III (2008). One-, two- and three-
phase viscosity treatments for basaltic lava flows, J. Geo-
phys. Res., 113, B09212; doi: 10.1029/2007JB005035.

Helsop, S.E., L. Wilson, H. Pinkerton and J.W. Head III
(1989). Dynamics of a confined lava flow on Kilauea vol-
cano, Hawaii, B. Volcanol., 51, 415-432.

Hérault, A., A. Vicari, A. Ciraudo and C. Del Negro (2009).
Forecasting lava flow hazards during the 2006 Etna erup-
tion: Using the MAGFLOW cellular automata model, Com-
put. Geosci.-UK, 35 (5), 1050-1060; doi: 10.1016/j.cageo.
2007.10.008.

Hon, K., J. Kauahikaua, R. Denlinger and K. McKay (1994).
Emplacement and inflation of pahoehoe sheet flows: ob-
servations and measurements of active flows an Kilauea
Volcano, Hawaii, Geol. Soc. Am. Bull., 106, 351-370.

Ishihara, K., M. Iguchi and K. Kamo (1990). Numerical sim-
ulation of lava flows on some volcanoes in Japan, In: J.H.
Fink (ed.), Lava Flows and Domes: Emplacement Mech-
anisms and Hazard Implications (IAVCEI Proceedings in
Volcanology), Springer-Verlag, 174-207.

Miyamoto H. and S. Sasaki (1997). Simulating lava flows by
an improved cellular automata method, Comput. Geosci.-
UK, 23 (3), 283-292.

TALLARICO ET AL.

519



520

Mouginis-Mark, P.J., H. Garbeil and P. Flament (1994). Effects
of viewing geometry on AVHRR observations of volcanic
thermal anomalies, Remote Sens. Environ., 48, 51-60.

Patankar, S.V. (1980). Numerical Heat Transfer and Fluid
Flow, Hemisphere Series on Computational Methods in
Mechanics and Thermal Science, Taylor & Francis, 214 pp.

Peterson, D.W. and D.A. Swanson (1974). Observed forma-
tion of lava tubes during 1970-1971 at Kilauea volcano,
Hawaii, Studies in Speleology, 2, 209-222.

Peterson, D.W., R.T. Holcomb, R.I. Tilling and R.L. Chris-
tiansen (1994). Development of lava tubes in the light of
observations at Mauna Ulu, Kilauea Volcano, Hawaii, B.
Volcanol., 56, 343-360.

Pieri, D.C. and S.M. Baloga (1986). Eruption rate, area and
length relationships for some Hawaiian lava flows, J. Vol-
canol. Geoth. Res., 30, 29-45.

Pinkerton, H. and R.S.J. Sparks (1978). Field measurements
of the rheology of lava, Nature, 276, 383-385.

Pinkerton, H. and R. Stevenson (1992). Methods of deter-
mining the rheological properties of magmas at sub-
solidus temperatures, J. Volcanol. Geoth. Res., 53, 47-66.

Rongo, R., W. Spataro, D. D'Ambrosio, M.V. Avolio, G.A. Trun-
fio and S. Di Gregorio (2008). Lava flow hazard evaluation
through cellular automata and genetic algorithms: an ap-
plication to Mt Etna volcano, Fund. Inform., 87, 247-268.

Rowland, S.K. and G.P.L. Walker (1990). Pahoehoe and aa in
Hawaii: volumetric flow rate controls the lava structure,
B. Volcanol., 52 (8), 631-641.

Sakimoto, S.E.H. and T.K.P. Gregg (2001). Channeled flow:
analytic solutions, laboratory experiments, and applica-
tions to lava flows, J. Geophys. Res., 106, 8629-8644.

Shaw, H.R., T.L. Wright, D.L. Peck and R. Okamura (1968).
The viscosity of basaltic magma: an analysis of field
measurements in Makaopuhi lava lake, Hawaii, Am. J.
Sci., 266, 255-264.

Smith, J.V. (2000). Textural evidence for dilatant (shear thick-
ening) rheology of magma at high crystal concentrations,
J. Volcanol. Geoth. Res., 99, 1-7.

Sonder I., B. Zimanowski and R. Büttner (2006). Non-New-
tonian viscosity of basaltic magma, Geophys. Res. Lett.,
33, L02303; doi: 10.1029/2005GL024240.

Spera, F., A. Borgia, J. Strimple, and M. Feigenson (1988).
Rheology of Melts and Magmatic Suspensions 1. Design
and Calibration of Concentric Cylinder Viscometer with
Application to Rhyolitic Magma, J. Geophys. Res., 93 (B9),
10273-10294.

Stein, D.J. and F.J. Spera (1992). Rheology and microstruc-
ture of magmatic emulsions: theory and experiments, J.
Volcanol. Geoth. Res., 49, 157-174.

Syrjala, S. (1995). Finite-element analysis of fully developed
laminar flow of power-law non- Newtonian fluid in a rec-
tangular duct, Int. Commun. Heat Mass, 22, 549-557.

Syrjala, S. (1996). Further finite element analysis of fully de-

veloped laminar flow of power-law non-Newtonian fluid
in rectangular ducts: heat transfer predictions, Int. Com-
mun. Heat Mass, 23, 6, 799-807.

Tallarico, A., M. Dragoni and G. Zito (2006). Evaluation of lava
effusion rate and viscosity from other flow parameters,
J. Geophys. Res., 111, B11205; doi: 10.1029/2005JB003762.

Turcotte, D.L. and G. Schubert (1982). Geodynamics: appli-
cations of continuum physics to geological problems,
New York, 450 pp.

Valerio, A., A. Tallarico and M. Dragoni (2008). Mechanisms
of formation of lava tubes, J. Geophys. Res., 113, B08209;
doi: 10.1029/2007JB005435.

Valerio, A., A. Tallarico, and M. Dragoni (2010). A model for
the formation of lava tubes by the growth of the crust
from the levees, J. Geophys. Res., 115, B09208; doi:
10.1029/2009JB006598.

Valerio, A., A. Tallarico and M. Dragoni (2011). Effects of
curvature of a lava channel on velocity and stress fields,
Geophys. J. Int., 187 (2), 825-832; doi: 10.1111/j.1365-246X.
2011.05166.x.

Vicari, A., A. Hérault, C. Del Negro, M. Coltelli, M. Marsella
and C. Proietti (2007). Modeling of the 2001 lava flow at
Etna volcano by a cellular automata approach, Environ.
Modell. Softw., 22, 1465-1471.

Weed, H.C., F.J. Ryerson and A.J. Piwinskii (1986). Rheolog-
ical Properties of Molten Kilauea Iki Basalt Containing
Suspended Crystals, In: K.S. Vorres (ed.), Mineral Matter
and Ash in Coal, ACS Symposium Series, 301, 223-233.

Woodcock, D. and A.J.L. Harris (2006). The dynamics of a
channel-fed lava flow on Pico Partido volcano, Lanzarote,
B. Volcanol., 69, 207-215.

Wright, R., S. Blake, A.J.L. Harris and D.A. Rothery (2001). A
simple explanation for the space-based calculation of lava
eruption rates, Earth Planet. Sci. Lett., 192, 223-233.

*Corresponding author: Andrea Tallarico,
Università di Bari, Centro Interdipartimentale di Ricerca per il Rischio
Sismico e Vulcanico, Bari, Italy; email: tallarico@geo.uniba.it.

© 2011 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights
reserved.

COOLING OF CHANNELED LAVA FLOW