Articles | Volume 14, issue 6
Geosci. Model Dev., 14, 3553–3575, 2021
Geosci. Model Dev., 14, 3553–3575, 2021

Model description paper 11 Jun 2021

Model description paper | 11 Jun 2021

A discontinuous Galerkin finite-element model for fast channelized lava flows v1.0

A discontinuous Galerkin finite-element model for fast channelized lava flows v1.0
Colton J. Conroy1,2 and Einat Lev1 Colton J. Conroy and Einat Lev
  • 1Lamont-Doherty Earth Observatory, Columbia University, New York, NY, USA
  • 2Roy M. Huffington Department of Earth Sciences, Southern Methodist University, Dallas, TX, USA

Correspondence: Colton J. Conroy (


Lava flows present a significant natural hazard to communities around volcanoes and are typically slow-moving (<1 to 5 cm s−1) and laminar. Recent lava flows during the 2018 eruption of Kīlauea volcano, Hawai'i, however, reached speeds as high as 11 m s−1 and were transitional to turbulent. The Kīlauea flows formed a complex network of braided channels departing from the classic rectangular channel geometry often employed by lava flow models. To investigate these extreme dynamics we develop a new lava flow model that incorporates nonlinear advection and a nonlinear expression for the fluid viscosity. The model makes use of novel discontinuous Galerkin (DG) finite-element methods and resolves complex channel geometry through the use of unstructured triangular meshes. We verify the model against an analytic test case and demonstrate convergence rates of P+1/2 for polynomials of degree 𝒫. Direct observations recorded by unoccupied aerial systems (UASs) during the Kīlauea eruption provide inlet conditions, constrain input parameters, and serve as a benchmark for model evaluation.

1 Introduction

On 3 May 2018, Kīlauea volcano on the island of Hawai'i began to erupt from new fissures in the lower East Rift Zone at the center of the Leilani Estates subdivision. Before ceasing in early August 2018, the lava flows destroyed over 650 structures and caused significant damage to infrastructure and essential facilities. During the second half of the eruption, the flow field established a complex braided channel system (which is common to many basaltic flows), originating from fissure number 8 (see Fig. 1). The “Fissure 8” flows were unique in the fact that they produced channelized flows reaching speeds as high as 15 m s−1 (Patrick et al.2019). These high speeds, coupled with channel geometry (e.g., constrictions) produced Reynolds numbers (Re>3000) that were significantly higher than typical lava flows. To investigate these extreme dynamics we develop a new channelized lava flow computer model named a discontinuous Galerkin finite-element model for fast channelized lava flows version 1.0.

This paper is organized as follows. In Sect. 1.11.3, we present the motivation for this work, as well as background on the mathematical tools we employ. In Sect. 2, we present the mathematical model and the bottom stress calculation and detail its nuances. We present the discontinuous Galerkin (DG) numerical discretization of the mathematical model in Sect. 3 and verify the model in Sect. 4. In Sect. 5, we evaluate the model against observations of lava flows from the 2018 eruption of Kīlauea volcano. We present misfit errors and root-mean-square errors (RMSEs) for the velocity field from a braided channel section of Fissure 8 and provide quantitative insight into physical quantities of the lava flow field in this area, including its thickness and viscosity. We close the paper in Sect. 6 with some discussion and conclusions.

Figure 1A satellite image (colored, in the background, by DigitalGlobe) overlaid with a thermal aerial orthomosaic (grayscale) where the white and light gray areas reveal the path of the Fissure 8 flow channel as it was on 21 June 2018. Data and map by USGS. The orange rectangle depicts the area of UAS site 8 from where the video we analyzed was captured on 22 June 2018. The flat gray areas south of the active flow channel demarcate the areas inundated by lava during the early stages of the eruption. The top of the image is north. PGV is the Puna Geothermal Ventures power plant that was heavily impacted by the lava.

1.1 Motivation

Typical “operational” lava flow models simulate unconfined lava flow in a 2D plan view (e.g., SCIARA, Crisci et al.2004; MAGFLOW, Vicari et al.2007; LavaPL, Connor et al.2012; VOLCFLOW, Kelfoun and Vargas2015) using either cellular automata or depth-averaged equations in an effort to forecast the area of land inundated by the lava. It is often difficult, however, for these models to accurately reproduce the complicated braided channel network such as those created by Fissure 8. These braided channel networks are common in natural flows (e.g., Dietterich and Cashman2014), and understanding the evolution of the velocity, rheology, and temperature fields (e.g. in response to pulsating effusion) within these channels is critical to hazard mitigation (Patrick et al.2019). Direct measurements of lava properties in situ is usually extremely difficult and dangerous. Modeling lava dynamics within the bounds of an established channel can help to better understand material properties of the flowing lava and inform models and decisions.

Previous attempts to model channelized lava flows have made use of simple heuristic formulas such as Jeffreys equation for laminar flows (Harris and Rowland2015) or Chezy approximations for higher-speed flows (Baloga et al.1995). While convenient, the use of these equations has largely been dictated by the fact that it has been difficult to obtain the physical data necessary for advanced modeling efforts (e.g., channel domain boundaries, inlet boundary conditions, topography). However, with the advent of unoccupied aerial systems (UASs, or “drones”) and their ability to survey active lava fields, we now have access to the data required by sophisticated numerical methods.

1.2 Shallow-water equations for fast lava flows

Commensurate with this development in observational capabilities, we introduce a numerical method for modeling fast-moving lava flows in complex channels. The high Reynolds number associated with these lava flows, coupled with the fact that the total length of the flows (on the order of kilometers) is much greater than the flow depth (on the order of meters), means that the dynamics can be well approximated by two-dimensional depth-integrated equations for mass, momentum, and energy. In particular, we utilize a system of dynamical equations known as the “shallow water equations” (Saint-Venant1851), which quantify average horizontal velocities and the depth of flow. These equations are traditionally used to model free surface flows in coastal oceanic regions, estuaries, and rivers (Dawson and Mirabito2008), although they have been used to model debris flows (George and Iverson2014) and lava flows (Costa and Macedonio2005) as well. The main assumption in the shallow water theory is that the fluid pressure is hydrostatic; gravitational acceleration dominates vertical accelerations in the fluid, and the pressure is calculated via the vertical momentum equation. The formulation of the shallow water equations that we utilize is designed specifically for advection-dominated flows (Kubatko et al.2006), and the pressure gradient term is formulated so that the dynamical equations are well balanced; steady states are preserved, and no artificial motion is induced by numerical artifacts (see Conroy2014, for a full derivation of the dynamical equations from conservation principles).

Lava flows are distinct from hydrological free surface flows in the sense that lava transfers heat to its surroundings; as lava effuses from a vent it cools along lateral flow boundaries and can form solid walls (“levees”) that prevent the lava from spreading to nearby regions. If lava effusion extends for several days, long channels may form that efficiently transport lava from the vent to the flow front. The speed at which the lava flows through the channel system depends on the viscosity of the lava, which in turn is highly dependent on the temperature and chemical composition of the lava (e.g., Griffiths2000). The presence of crystals and/or bubbles in the lava can make its viscosity non-Newtonian (Manga et al.1998; Mader et al.2013) and thus strongly dependent on stress gradients and the thermal properties of the lava. To reflect this strong dependence on temperature, we solve a depth-integrated energy equation that quantifies the thermal evolution of the lava as it interacts with its environment. The depth-integrated energy equation is coupled to the shallow water equations through a thermally dependent nonlinear stress term that reflects the rheology of the lava and can account for the presence of crystals and/or bubbles in the lava flow.

The logistical key to using shallow water equations to model lava flow dynamics rests on the development of the non-Newtonian bottom stress term. Typical friction drag laws do not take into account the viscosity of the fluid (due to the assumption that the fluids inertial acceleration is much greater than its internal resistance). However, in our particular case the flow is not fully turbulent; internal resistance needs to be taken into account in some fashion. Thus, we express the stress at the bottom boundary as a function of the temperature and the vertical stress gradient (which is a function on the vorticity). We solve a thermal boundary layer problem to calculate the temperature at the bottom boundary and utilize the vorticity to determine a virtual length scale over which the interior velocity goes from the depth-averaged value to zero. This results in a bottom stress approximation that is void of a friction factor (e.g., Manning's n) and allows scientists to study physical properties of the lava that are difficult to measure directly (e.g., viscosity). For example, application of the model to the Kīlauea 2018 Fissure 8 lava flows reveals that the lava behaved as a shear thickening fluid due to the high bubble content ( 50 %) with a large capillary number, which agrees well with the recent lab experiments of Lev et al. (2020).

1.3 The discontinuous Galerkin finite-element method

Because closed-form analytic solutions do not exist to the nonlinear shallow water equations and energy equation, we construct approximate solutions to these equations using discontinuous Galerkin (DG) finite-element methods (Cockburn and Shu2001), which have been used to successfully model other geophysical fluid flows including coastal ocean circulation (Kärnä et al.2018), hurricanes (Dawson et al.2011), avalanches (Patra et al.2006), and debris flows (Conroy and George2021). The DG finite-element method differs from continuous Galerkin (CG) finite-element methods in that the DG method solves an integral (or weak) form of the mathematical equations over individual elements and utilizes a solution space that is discontinuous across element boundaries. This allows the DG method to resolve steep gradients that form in the numerical solution, such as the thermodynamic gradients that form at lava channel wall boundaries. Even though the DG method is discontinuous it still conserves mass both locally and globally by utilizing a numerical flux function (introduced by finite-volume methods; see LeVeque (2002) for instance) that takes the discontinuous state of physical properties at element boundaries and creates a consistent flow of information from element to element (Cockburn and Shu2001).

Further, the DG method has been shown to be highly parallelizable using high-performance computing (e.g., Kubatko et al.2009, and Patra et al.2006) and it is amenable to unstructured numerical meshes. The latter feature is important when resolving geometrically complex boundaries of a given fluid domain, such as the flow fields commonly produced by basaltic flows. For instance, the lava flows that effused from Fissure 8 formed a complicated network of braided channels, with multiple locations of branching and merging. In addition, obstructions such as large lava rafts or preexisting structures caused local disruptions in the flow field, making it difficult to evaluate the dynamics using a simplified one-dimensional channelized model, such as those that use a classical rectangular channel geometry (Harris and Rowland2015). To account for these complexities, our new lava flow model discretizes the lava channel domain with an unstructured triangular mesh. This reduces model error as it pertains to a representation of the lava flow domain and is important when reproducing localized flow features of the lava field, such as the jet visible in Fig. 4.

Model verification consists of solving an analytic test case using forcing functions that we choose to exactly satisfy the equations of motion, and results indicate that for smooth solutions the method converges to the exact solution at a rate of P+1/2 for polynomials of degree 𝒫.

2 Mathematical model

Fluid flow on sloped terrain can be quantified in a Cartesian coordinate (x, y, z) system over a time-dependent domain Ω(t)∈ℝ3 by solving Eulerian conservation equations of mass, linear momentum, and energy,


In Eqs. (1)–(3), ρ is the fluid density, u=(u,v,w) is the fluid velocity vector, =x,y,z is the gradient operator, p is the pressure, and τ is the stress term,

(4) τ = τ x x τ x y τ x z τ y x τ y y τ y z τ z x τ z y τ z z .

We denote the x component of gravitational acceleration that is tangential to the sloped surface by gx; gy is the y component of gravitational acceleration that is tangential to the sloped surface, and gz is gravitational acceleration that is normal to the sloped surface. Collectively, these terms form the body force vector fb= (gx, gy, gz) acting on the fluid. T is the temperature of the fluid, cp is the specific heat capacity of the fluid, q=kTTx,kTTy,kTTz is the heat flux through the fluid (kT is a conduction heat transfer coefficient that measures the spread of heat within the fluid), and q˙ quantifies the generation and dissipation of heat within the fluid.

Equations (1) and (2), taken together, are the Navier–Stokes equations and quantify the force dynamics acting on the fluid (see Conroy2014, for a derivation from first principles). Equation (3) is the thermal energy equation: it quantifies the transport of energy through the fluid due to internal temperature gradients and differences between the fluid temperature and the temperature of the surrounding medium (see Moran et al.2003, for instance). To apply Eqs. (1)–(4) to channelized lava flow we need to supplement them with appropriate boundary conditions and define the stress matrix τ. Here, we assume that the system of equations given by Eqs. (1)–(3) are subject to the following kinematic, dynamic, and thermal boundary conditions:

  • Channel wall boundary condition:
    no normal flow un^=0;
    no pressure gradient, p/n^=0;
    slip velocity, ut^=f(τwall,σwall);
    heat loss via conduction, qn^=kT/hwT-Twall.

  • Inlet boundary condition:
    prescribed velocity, un^=prescribed;
    prescribed pressure, p=prescribed;
    prescribed heat content, ρcpT=prescribed.

  • Outlet boundary condition:
    zero change in normal velocity, u/n^=0;
    zero change in pressure, p/n^=0;
    zero change in heat content, ρcpT/n^=0.

  • Free surface boundary condition at z=ζ:
    no relative normal flow,
    atmospheric pressure, p=patm;
    surface shear stress, τsx=-τxx(ζ/x)-τyx(ζ/y)+τzx, τsy=-τxy(ζ/x)-τyy(ζ/y)+τzy;
    heat loss via radiation and convection,

  • Bottom boundary condition at z=-h:
    no slip velocity, u=0;
    bottom shear stress, τbx=τxx(h/x)+τyx(h/y)+τzx, τby=τxy(h/x)+τyy(h/y)+τzy;
    heat loss via conduction, qn^=(kb/hb)T-Tground.

In the boundary conditions above, n^ is the unit normal vector to the wall boundary, outlet boundary, inlet boundary, moving free surface, and bottom boundary, respectively; t^ is the unit tangential vector to the wall; τwall is the tangential shear stress at the wall; and σwall is the normal stress at the wall. We denote the temperature at the wall by Tw, kw is the thermal conductivity of the wall, and hw is the thickness of the thermal boundary layer through which heat is conducted from the lava flow to the wall. We measure the free surface, ζ, relative to a steady depth of flow, h(x,y), that serves as the zero datum in the z direction (see Fig. 2). The surface forces acting on the free surface consist of the atmospheric pressure, patm, along with a surface shear stress applied by the wind τs=(τsx,τsy). Heat transfer from the lava surface to the surrounding atmosphere is dominantly due to radiation and air convection, where ϵ is the emissivity of the lava, σB is the Stefan–Boltzmann constant, Tatm is the temperature of the surrounding atmosphere, and kc is the convection heat transfer coefficient. The main resisting force in dense shallow mass high-speed flows comes from the bottom stress, τb=(τbx,τby), which is a function of the temperature of the lava at the basal boundary, where Tground is the temperature of the ground and hb is the depth of the thermal boundary layer through which heat is transferred from the lava to the ground.

Figure 2A vertical cross section along the center line of the flow, showing the coordinates (x and z, with y the across-flow direction) and the heat transfer mechanisms considered in the model (conduction to the base, advection by the flow, and heat loss by radiation and convection at the surface).


Theoretically, we could define τ and solve Eqs. (1)–(3) along with the prescribed boundary conditions to model channelized lava flows. In practice, however, we need to simplify Eqs. (1)–(3) to make the solution more tractable. More specifically, we assume the following: (i) the lava flow field is incompressible, (ii) vertical accelerations in the lava are dominated by gravity, (iii) lava flow lengths are much greater than the flow depth, and (iv) horizontal flow speeds are large enough that stress gradients are dominated by the first two columns of the stress matrix τ.

Assumption (i) reduces the conservation of mass equation to u=0, while assumptions (ii) and (iv) reduce the z-momentum equation to


which we can leverage to determine the pressure. Integrating from the free surface ζ down to a given z coordinate yields

(5) p = p atm + ρ g z ζ - z .

We assume that gradients in patm are negligible and the horizontal pressure gradient in Eq. (2) becomes


where =x,y. We further simplify the mathematical model and eliminate the vertical dimension by integrating u=0 and Eqs. (2) and (3) over the depth of the lava flow (from h to ζ). We then apply Leibniz's integral rule, utilize the free-surface and bottom boundary conditions, assume the density is constant, and simplify the resulting expression to arrive at the following depth-integrated equations:


where fx and fy are the x component and y component of the depth-averaged gradient of the shear stresses acting on vertical fluid planes, qs is the heat flux through the free surface, qb is the heat flux at the bottom boundary, qi quantifies depth-averaged internal conduction, and q˙ represents depth-averaged internal heat generation and dissipation. The depth-averaged velocity, u=(u,v), and depth-averaged temperature, T, are defined as


The above system of depth-integrated equations (Eqs. 69) can be further simplified due to the dynamics of high-speed flows. More specifically, fx and fy are negligible in high-speed flows except at no-slip and small-slip velocity boundary conditions where large stress gradients form due to the decay of the velocity field to a value of zero (or near zero). In our quantitative analysis of UAS footage from the 2018 Kīlauea eruption, lava flow velocities at channel wall boundaries were much greater than zero, and therefore we utilize a slip (no-flow) channel wall boundary condition and neglect fx and fy in this initial version of the model (see Rao and Rajagopai1999, for an in depth investigation on channel wall boundary conditions in terms of the slip versus no-slip condition). The surface stress terms (τx,τy) in the depth-integrated equations account for wind stress on the lava flow, which we assume to be negligible due to the ratio of the density of air to the density of lava being much less than one. We include the effect of heat conduction at the channel wall boundaries, but we neglect heat conduction in the interior of the lava due to the high speed of the flow, and we set qi=0. Internal heat generation and dissipation can be significant in lava flows with a high crystal content (Griffiths2000) and in lava flows in closed tubes (Costa and Macedonio2005). The Fissure 8 lava flows were hot with limited crust cover, and samples indicate that the crystal content was low in the channel section that we apply the model to (Gansecki et al.2019), and therefore we neglect q˙ in the current model but plan to include it in future releases.

It can be noted that the pressure gradient terms, (1/ρ)p=Hx(gzζ),Hx(gzζ), in Eqs. (7) and (8) are non-conservative product terms that can lead to entropy-violating numerical fluxes if care is not taken in evaluating them numerically (see LeVeque2002, for instance). To circumvent this issue, we make use of the fact H(x,y,t)=ζ(x,y,t)+h(x,y) and rewrite Eqs. (7) and (8) in the conservative form (Kubatko et al.2006),


where P=12gzH2-h2 is the pressure flux and h/x and h/y quantify the gradient in the steady reference depth of flow that ζ is measured relative to. We supplement the system of equations given by Eqs. (13)–(16) with initial conditions along with the channel wall, inlet, and outlet boundary conditions. It can be noted that the depth-integrated mass and momentum equations given by Eqs. (13)–(15) are well studied in the literature and are commonly used to model shallow mass flows such as coastal ocean circulation and hurricane storm surge; see, for example, Dawson et al. (2011) and Kubatko et al. (2006). The addition of the energy equation complicates the solution of Eqs. (14) and (15) due to the fact that the bottom stress term, τb=(τbx,τby), is now a function of both nonlinear velocity gradients and temperature.

2.1 Quantifying the bottom stress term

In the equations of motion (Eqs. 14 and 15), we define the bottom stress term using a Herschel–Bulkley model (Herschel and Bulkley1926),

(17) τ b = τ z x = μ u z + τ yield sgn u z ,

where τb=τzx=(τzx,τzy), τyield is the yield strength of the fluid, sgn denotes the sign of the argument, and μ is the nonlinear viscosity, defined as

(18) μ = K u z n - 1 .

The symbol 𝒦 in Eq. (18) represents the consistency of the lava and can be modeled solely as a function of temperature, while the power law exponent, n, is typically a function of the particle content (crystals and/or bubbles) of the lava (which in turn can be a function of temperature); see Castruccio et al. (2010) and Castruccio et al. (2014) for example. We quantify the temperature dependency of the lava consistency (𝒦) in a fashion similar to Sonder et al. (2006) and use the Vogel–Fulcher–Tammann (VFT) silicate melt model of Giordano et al. (2008),

(19) log K K 0 = A + B T ( K ) - C ,

where A is the value of logK/K0 at infinite temperature, 𝒦0 is a constant set to 1 sn−1, and B and C are parameters that depend on the composition of the lava. The model assumes that A is a constant for all silicate melts regardless of composition, and thus it represents the high temperature limit for silicate melt viscosity. Once the parameter A is fixed then the parameters B and C are determined via a linear ensemble of combinations of oxide components and a subordinate number of multiplicative oxide cross terms; see Giordano et al. (2008) for the full details of the model.

The power law exponent of the viscosity term quantifies the effect that stress gradients have on the material properties of the fluid. A value of n=1 corresponds to a Newtonian fluid, while n<1 or n>1 corresponds to a non-Newtonian fluid. If n>1 then the fluid viscosity increases with increasing shear rate (shear thickening), while if n<1 the fluid viscosity decreases with increasing shear rate (known as shear thinning). Typically, if the lava is sufficiently hot and degassed, then the lava stress can be modeled with a Newtonian approximation and n=1. However, if bubbles and/or crystals are present in the lava (depending on the lava source and the amount of degassing that has occurred), then these structures will deform and realign under an applied shear stress. This consequently causes the viscosity of the lava to become thinner in some situations and thicker in others depending on how the structures rearrange. Lava flows with a high crystal content are typically pseudo-plastic and shear thinning; the crystal structure of the lava resists the flow of the lava and the lava will not flow unless a yield strength is surpassed. In this case, the lava will continue to flow more readily as the shear stress increases. The opposite tends to occur when a lava flow has a high bubble content at higher capillary numbers; large stress gradients in the flow cause the bubbles to rearrange in a fashion that increases the viscosity and the lava behaves as a shear thickening fluid.

The bottom stress term is a function of the velocity gradient evaluated at the bottom boundary, which we do not have access to in the depth-averaged equations, and therefore we define the bottom stress in terms of the depth-averaged velocity as

(20) τ b x μ u δ z + τ yield sgn u δ z μ u z + τ yield sgn u z ,


(21) μ = K u δ z n - 1 .

It can be noted that in the expression above, δz=(δzx,δzy), is a measure of a virtual length over which the shear stress is applied. We determine δz by taking into account the vorticity of the lava flow field, which is defined as

(22) ω = w y - v z i ^ + w x - u z j ^ + u y - v x k ^ ,

where i^, j^, and k^ are unit vectors in the x, y, and z directions, respectively. We solve an auxiliary problem over a pseudo-depth of the lava that consists of an upper mixed layer where u(x,y,z)u(x,y) and a lower layer where (u/z,v/z)(w/x,w/y). We assume that the vorticity in the upper layer is equal to the vorticity in the bottom layer (in terms of magnitude) at the coordinate point where u/z0 (i.e., at the interface between the two layers). This allows us to calculate the vorticity in the upper layer and then use this value to determine δz. In other words, we answer the following question: given a measure of the vorticity associated with the depth-averaged velocity field, what is the associated length scale over which the depth-averaged velocity must decay to a value of zero (the bottom boundary condition) to ensure that the internal vorticity of the flow is conserved? It can be noted that an implicit assumption in depth-integrated models is that internal friction in the vertical is null compared to the friction at flow boundaries. This, along with assumption (i) and a constant density, implies conservation of vorticity about the i^ and j^ directions except at flow boundaries). The key to this approach relies on calculating a measure for the vertical velocity in the upper layer, which we achieve by making use of the kinematic boundary condition,

(23) ζ t + u ζ x z = ζ + v ζ y z = ζ = w z = ζ ,

coupled with the depth-integrated continuity Eq. (13) to obtain a measure of the vertical velocity w. More specifically, expanding derivatives in Eq. (13), solving for ζ/t, while substituting this result into Eq. (23) and noting that in the upper layer (u,v,w)(u(ζ),v(ζ),w(ζ)) yields

(24) w = ζ u x + ζ v y .

The relevant vorticity terms in Eq. (22) include the i^ and j^ components. By definition, u/z=v/z=0 over the upper layer so that the vorticity component about the x axis is w/y and the vorticity component about the y axis is w/x. Because the bottom boundary condition is modeled as a rigid wall where u=0 and the fluid is incompressible, a vorticity layer forms in the fluid near the solid boundary that resists the local rotation of the fluid (this is the reason why the rigid boundary does not deform). The vorticity created at the boundary resists the rotation of the interior and is equal to v/z about the x axis and u/z about the y axis (see Schlichting et al.1968). Now, if we assume that each vorticity component over the bulk of the flow is equal to each vorticity component in the boundary layer at the coordinate point where (u/z,v/z) is no longer equal to zero, then the virtual length over which the shear stress is applied is given by

(25) δ z x = u w / x and δ z y = v w / y .

It can be noted that as (w/x,w/y) goes to 0, the vertical stress in the fluid goes to 0. We can rewrite expression (Eq. 25) solely in terms of the depth-averaged variables using Eq. (24):

(26) δ z x = u x ζ u x + ζ v y - 1 and δ z y = v y ζ u x + ζ v y - 1 .

It can be noted that even though we do not explicitly include horizontal shear stresses in the Kīlauea simulations presented in Sect. 5 (due to the high Re number), the virtual length used to quantify the bottom stress as defined in Eq. (26) is a function of horizontal shear within the fluid. Further, we wish to emphasize that the virtual length, δz, is non-physical and is not necessarily less than the lava flow thickness (H); it merely is a measure to ensure that (u-0)/δzu/z at the bottom boundary in a fashion that conserves internal vorticity, i.e., it ensures that the interior of the flow field is irrotational about the i^ and j^ coordinate axis.

2.2 Heat transfer

As soon as lava effuses from an active vent it begins to degas and transfer heat to its surroundings. Lava cools through the mechanisms of radiation, conduction, and convection in the air above it (we neglect heating from viscous dissipation, which is small compared to heat loss through radiation and conduction for the low-viscosity flows we are considering here, e.g., Harris and Rowland2001). We quantify heat loss due to radiation via Stefan's law (Griffiths2000),

(27) q s = ϵ σ B ρ c p T 4 - T atm 4 ,

where ϵ is the emissivity of the lava, σB is the Stefan–Boltzmann constant, and Tatm is the temperature of the surrounding atmosphere in degrees Kelvin. When lava temperatures fall below the solidus (e.g.,  950 C for Kīlauea lavas), buoyancy-driven convection in the air above the lava becomes the dominant mode of heat transfer at the lava surface instead of radiation (due to crust formation) (Griffiths2000). In this case we set qs in the energy equation to

(28) q s = k c T - T atm 4 / 3 ,

where kc is a heat transfer coefficient (e.g., Patrick et al.2004). We quantify heat transfer from the lava to the ground via conduction (Patrick et al.2004),

(29) q b = k b T - T ground ,

where Tground=f(x) is the temperature of the ground in contact with the lava flow field and kb measures the thermal conductivity of the ground. We utilize Eq. (29) to determine the temperature near the bottom boundary of lava flow field which we use to evaluate the nonlinear viscosity in the bottom stress term in the equations of motion. More specifically, we can rewrite Eq. (29) in terms of a depth-dependent thermal boundary layer temperature, T(z),

(30) T z = k ̃ b h b T - T ground ,

where k̃b is the thermal boundary layer conductivity constant and hb(x,y) is the thickness of the thermal boundary layer (see Fig. 2). We solve Eq. (30) over the thermal boundary layer defined in the z direction from z=zb(x,y) to z=-h(x,y) by setting zb(x,y) to a relative zero. We then integrate Eq. (30) over a thermal boundary coordinate defined from z̃=0 to z̃=-hb(x,y). It can be noted that the relationship between z̃ and z is given by z̃=z-zb so that hb=h-zb. Equation (30) is a non-homogeneous, constant coefficient ordinary differential equation that has the following solution:

(31) T ( z ̃ ) = T int - T ground exp k ̃ b h b z ̃ + T ground ,

which gives an expression for the temperature profile over the thermal boundary layer of the lava (z̃[0,-hb]). Evaluating Eq. (31) at z̃=-hb and setting the interior temperature (Tint) to the depth-integrated value (T), we have

(32) T | z ̃ = - h b = T - T ground exp - k ̃ b + T ground .

We use this temperature to evaluate the consistency in the bottom stress approximation,

(33) K | z = - h = K T | z ̃ = - h b .

The greater the thermal conductivity of the boundary layer, the closer the boundary temperature is to the ground temperature; however, in general, there is usually a steep gradient in the temperature at the interface between the boundary of the flowing lava and the ground that the lava is conducting heat to. It can be noted that we also use an analogous approach to calculate the temperature of the lava at channel wall boundaries.

2.3 Steady reference depth of flow h

We have two options to calculate the steady reference depth of flow (h) of the lava that we use as a zero datum to measure the free surface from. Our particular choice depends on the inflow data available to the model. For instance, if a full set of temporally varying inflow data is available, we set h equal to the time average thickness associated with the data, i.e.,

(34) h = 1 t f - t i t i t f Q in n ^ w in u in n ^ d t ,

where Qinn^ is the inflow flux normal to the boundary and win is the width of the inflow boundary normal to the flow. If, however, the only inflow data available to the model is a set of time-averaged data, then we set h to the solution of the steady, linear system of equations associated with the full nonlinear system of equations given by Eqs. (13)–(16).

3 Numerical discretization

To develop our numerical methods, we rewrite the system of Eqs. (1) in the compact form,

(35) U ( i ) t + F ( i ) ( U ) = S ( i ) ( U ) , i = 1 , 2 , 3 , 4 ,

where U(i), F(i), and S(i) are the i-th row entries of the vectors U, S, and the flux function matrix F, defined as


where P=12gzH2-h2.

3.1 Finite-element partition

To apply a DG spatial discretization to our mathematical model (Eq. 35) over a lava flow channel (see Fig. 5, for example), we begin by introducing a partition of the two-dimensional domain Ω. The complexities of the domain boundary, ∂Ω, are such that an unstructured finite-element partition (or mesh) is necessary to properly capture its intricacies. More specifically, we obtain unstructured triangulations (that we denote by 𝒯h) of the channel domain Ω via an automatic mesh generator known as ADMESH+ (Conroy et al.2012). ADMESH+ solves a number of differential equations to calculate a mesh size function that determines local element sizes based on the curvature of the boundary, channel width, and changes in the topography and domain slope to create a high-quality unstructured simplex mesh (the elements are close to equilateral triangles). The only input required by the program is a list of points defining the boundary and the topography of the domain.

3.2 A weak form and the semi-discrete equations

Given the finite-element partition, 𝒯h, of the domain Ω, we obtain a weak form of Eq. (35) if we first multiply Eq. (35) by a sufficiently smooth test function ψ(x,y)V, integrate over each element Ωj∈𝒯h, and then integrate the flux term by parts,

(36) Ω j U ( i ) t ψ d A - Ω j F ( i ) ψ d A + Ω j F ( i ) n ^ ψ d S = Ω j S ( i ) ψ d A , U ( i ) , ψ V ,

for i=1,2,3,4 and j=1,,N, where 𝒩 is the total number of elements of the triangulation 𝒯h. In the equation above, n^ is the outward unit normal to the element boundary ∂Ωj. Rather than seek solutions to Eq. (36), we search for solutions in the finite dimensional subspace of functions defined as

(37) V h p = ψ : ψ Ω j P Ω j , Ω j ,

where 𝒫 demarcates the space of polynomials of at most degree that is not necessarily continuous across element boundaries. In other words, given a set of basis functions ϕ=ϕ0,ϕ1,,ϕ, we express the trial solution (Uh(i)Vhp) and test function (ψh∈𝒱hp) as

(38) U h ( i ) | Ω j = l = 0 U l ( i ) ( t ) ϕ l ( x ) ,


(39) ψ h | Ω j = l = 0 ψ l ( t ) ϕ l ( x ) ,

where U0(i),U1(i),,U(i) are the time-dependent degrees of freedom of the finite-element solution and i=1,2,3,4. We use products of Jacobi polynomials of degree , Pll=0, as the basis for 𝒱hp. The orthogonal triangular basis is defined in terms of a “collapsed coordinate” system that results in a matrix-free implementation of the method; see Kubatko et al. (2006) for more details. Substituting Uh(i) and ψh into Eq. (36), we arrive at the discrete weak form of the problem: find Uh(i)Vhp such that for all test functions ψh∈𝒱hp for i=1,2,3,4, the expression,

(40) Ω j U h ( i ) t ψ h d A - Ω j F ( i ) ( U h ) ψ h d A + Ω j F ^ ( i ) n ^ ψ h d S = Ω j S ( i ) ( U h ) ψ h d A ,

holds over each element Ωj∈𝒯h, where S(i)(Uh) is the source term evaluated in 𝒱hp and F^(i) is a suitably chosen numerical flux.

3.2.1 Numerical flux

The space of functions defined by Eq. (37) is not necessarily continuous across element boundaries, and thus can be dual-valued (see Fig. 3, for example). To remedy this inconsistency, we replace the dual-valued flux in Eq. (36) with a so-called numerical flux (F^) that makes use of the left and right limits of the trial solution to produce a single-valued flux across a given element's boundary.

Figure 3Jump in numerical solution Uh at an element edge ∂Ωj.


More specifically, given an arbitrary function wh∈𝒱hp at an element boundary point xi, we set the left and right limits of the function to wh-wh(xi-) and wh+wh(xi+), respectively. In this work we utilize the local Lax–Friedrichs (LLF) flux, which defines the numerical flux operator as

(41) F ^ ( i ) n ^ = 1 2 F ( i , + ) + F ( i , - ) n ^ - 1 2 λ max U h ( i , + ) - U h ( i , - ) , for i = 1 , 2 , 3 , 4 ,

where λmax is the maximum eigenvalue of the normal (to the element edges) Jacobian matrix. When solutions to Eq. (35) are sufficiently smooth, we can rewrite Eq. (35) in the quasi-linear form,

(42) U t + J x U x + J y U y = S ,

where the Jacobian matrices Jij=Fixj are




The so-called “normal Jacobian matrix” is then defined by

(43) J n = J x n x + J y n y ,

where nx and ny are the x and y components of the normal edge vector n^. In general, if J is a square (m×m) matrix with m real eigenvalues, then it can be decomposed into its eigensystem,

(44) J x = R ( x ) Λ ( x ) R ( x ) - 1 , and J y = R ( y ) Λ ( y ) R ( y ) - 1 ,

where (⋅) is the matrix of right eigenvectors, Λ(⋅) is the diagonal matrix of eigenvalues, and R()-1 is the matrix of left eigenvectors (LeVeque2002). To determine Λ(x) and Λ(y) we solve for the roots of detJ()-λI=0, which gives the following eigenvalues,

(45) λ 1 , 2 = u n x + v n y , λ 3 = u + g z H n x + v + g z H n y , λ 4 = u - g z H n x + v - g z H n y .

Each eigenvector (ri()) can be determined by solving (J()-λiI)ri=0, where I is the identity matrix, 0 is a vector of zeros, and R()=r1(),r2(),r3(),r4(). Solving for the eigenvectors we have




We use the full eigensystem in the slope-limiting process that stabilizes the method for polynomials of degree greater than or equal to one (Cockburn and Shu2001), and we set λmax in the LLF flux to the maximum value of (λ1,λ2,λ3,λ4).

It can be noted that to mathematically close the solution method of the discrete DG system of equations we need to numerically evaluate expression (Eq. 26) to determine δz. More specifically, we discretize Eq. (26) using a local discontinuous Galerkin (LDG) method (Cockburn and Chi-Wang1998) analogous to the method used in Conroy and Kubatko (2016) to evaluate second-order derivative terms (see Conroy2014, for a detailed discussion on application of the LDG method to shallow mass geophysical fluid flows).

3.3 Strong stability-preserving Runge–Kutta time discretizations

Application of the DG spatial operator to Eq. (40) results in a system of ordinary differential equations (ODEs) for each element,

(46) M ̃ j ( i ) d U ̃ j ( i ) d t = b j ( i ) , i = 1 , 2 , 3 , 4 and j = 1 , , N ,

where 𝒩 is the number of elements in 𝒯h, Ũj(i)=Uj,0(i),Uj,1(i),,Uj,(i) are vectors of the degrees of freedom (i.e., the polynomial basis coefficients), and bj(i)=Rj,0(i),Rj,1(i),,Rj,(i) with,

(47) R j , l ( i ) = Ω j F j ( i ) ϕ l d A - Ω j F ^ j ( i ) n ϕ l d S + Ω j S j ( i ) ϕ l d A .

M̃j(i) is the mass matrix,


which is diagonal due to the choice of basis. Left multiplying (Eq. 46) by the inverse of the mass matrix, we have

(48) d U ̃ ( i ) d t = M ̃ j ( i ) - 1 b j ( i ) = L h p ( U ̃ ) , with i = 1 , 2 , 3 , 4 , and j = 1 , , N ,

where hp is the DG spatial operator. We evaluate the integrals in Eq. (47) using numerical integration rules of sufficiently high degree (Kubatko et al.2006) and discretize Eq. (48) with so-called strong stability-preserving (SSP) Runge–Kutta (RK) methods (Kubatko et al.2014). The unknown polynomial basis coefficients that define the solution over a given element, Ωj, are advanced in time from tm to tm+1 via

  1. Set Ũ0(i)Ũm(i), for i=1,2,3,4.

  2. For each stage r=1,2,,S, set

  3. Finally, set Ũm+1(i)ŨS(i).

It can be noted that Πh is a slope limiter that dampens overshoots and undershoots at solution discontinuities when polynomial approximations greater than 0 are used for the basis (Cockburn and Shu2001), 𝒮 is the number of stages of the RK method, δsΔt is a sub-time step of the time step Δt, and the αrs and βrs are coefficients that define the RK method. In particular, αrs and βrs conform to the following constraints:

  1. αrs=0 if and only if βrs=0;

  2. αrs≥0 and βrs≥0;

  3. s=1rαrs=1.

Because we use explicit RK methods, the time step of the model is limited by a CFL condition; see Kubatko et al. (2014) for more details.

Table 1L2 errors using 𝒫0 for (ζh,uh,vh,Th).

Download Print Version | Download XLSX

Table 2L2 errors using 𝒫1 for (ζh,uh,vh,Th).

Download Print Version | Download XLSX

Table 3L2 errors using 𝒫2 for (ζh,uh,vh,Th).

Download Print Version | Download XLSX

4 Verification

Verification of the DG solution of the mass and momentum equations in the depth-averaged and full three-dimensional case is well documented and can be found in Conroy and Kubatko (2016), Dawson and Aizinger (2005), and Kubatko et al. (2006). To verify our DG solution method for the fully coupled mass, momentum, and energy (depth-averaged) equations, we solve a test problem that is designed to model a free surface wave propagating through a lava channel using the method of manufactured solutions (see Griffiths2000, for a discussion on free surface waves in lava channels at high Re number and Le Moigne et al. (2020) for a detailed investigation on standing waves in lava flow channels). We choose u, v, and H so that the depth-averaged mass equation is satisfied exactly. Specifically, we define

(49) u = u ^ exp - x , v = v ^ exp ( - k x ) , H = h + ζ ^ exp ( i ω t ) ,

with h= constant, u^= constant, v^=u^ky, and ζ^=exp(-iωx^) where x^=exp(kx)/(ku^). The dynamical solution consists of a wave propagating in a direction perpendicular to the y axis with wave number k=5.0×10-3 m−1 and frequency ω=2.0×10-3 s−1. These values were chosen based on velocity and lava flow thickness data recorded by Patrick et al. (2019) during the pulsing effusion regime associated with Fissure 8 during the 2018 Kīlauea event. We then set

(50) T = T ^ exp ( - k T x ) ,

with T^=(y2+Twall) and substitute Eqs. (49) and (50) into the mathematical model (Eqs. 1316) and evaluate the derivative terms using MATLAB's symbolic package. The remainder terms associated with the x momentum equation and the energy equation are then set as artificial source terms that force the numerical solution to be Eqs. (49) and (50). The numerical domain consists of a rectangular channel defined by the Cartesian coordinates x0=0.0 m, xL=200.0 m, y0=0.0 m, yL=30.0 m. We assume symmetry about the centerline (at y=15 m) and only solve the equations over the half-width of the channel. It can be noted that even though the solutions are guaranteed to remain smooth for all time t (because of the forcing functions), the numerical solution is by no means trivial due to the coupling of the equations through the viscosity. We use four different triangular meshes for the verification of the model. The element size of each mesh is 7.50 m (the so-called dx0 mesh), 3.75 m (dx1), 2.50 m (dx2), and 0.625 m (dx3), respectively. Results are displayed in Tables 13, where it can be noted that the method converges to the analytic solution at a rate of approximately +1/2. Further, using 𝒫2 polynomials on the coarsest mesh gives lower errors than 𝒫0 polynomials on the finest mesh. It can be noted that for a given computational mesh, a higher-order polynomial approximation will result in a greater computational expense. However, the goal of using high-order polynomial approximations is to use coarser meshes, which results in better computational efficiency in terms of the number of degrees of freedom necessary to achieve a certain level of accuracy (high-order local polynomials produce more accurate results more efficiently than low-order methods). This is explicitly shown in the works of Kubatko et al. (2006), Kubatko et al. (2009), and Conroy et al. (2018).

Figure 4(Left) Zoomed in view of the section of the braided channel system established by Fissure 8 modeled in this work, near UAS site 8 (See Fig. 1). UAS photo by Ryan Perroy, University of Hawai’i-Hilo. (Right) Map view of the lava surface velocity measured using Optical Flow from videos captured by UAS on 22 June 2018. Colors represent magnitude in m/s. Also shown is the finite-element mesh used to evaluate the model.

5 Evaluation: recent eruption of Kīlauea volcano

We evaluate our model using data captured during the 2018 eruption of Kīlauea volcano, Hawai'i. The East Rift Zone of Kīlauea has erupted repeatedly in historical times and continuously since 1983 (Heliker and Mattox2003; Wolfe1988). A new eruption of unusually large magnitude began on 3 May 2018 in the lower part of the East Rift Zone, with fissures opening in the middle of a residential area (Neal et al.2019). More than 20 fissures opened during the first 12 d of the eruption, erupting slow-moving, unusually high-viscosity lava at low effusion rates. The behavior changed on 18 May, when much hotter and less viscous lava reached the surface. Advance rates and flow lengths increased, widely impacting property and infrastructure. Complete evacuation orders followed within days. Starting on 28 May, activity focused at Fissure 8, located in the heart of the Leilani Estate subdivision. Fissure 8 remained the source of lava for the remainder of the eruption until its abrupt stop on 4 August. The lava that erupted from Fissure 8 soon established a channel that flowed north and east of the vent, forming a moderately branched channel network 4 km from the vent. The flow field exhibited transitions between flow types: a clear transition from pahoehoe to `a`a surface texture occurred downslope and is apparent on the thermal map (see Fig. 1). Overall, the Fissure 8 lava covered an area of 25 km2 and supplied at least 1 3 of lava (out of at least 1.2 total) over 70 d.

5.1 Observational data

During the 2018 Kīlauea eruption, unoccupied aerial systems (UASs) captured a comprehensive time series of overhead videos of channelized lava (the Fissure 8 flow). The videography campaign was purposefully designed to collect data for “remote rheometry” by hovering above specific sites spaced 200–1300 m apart along the length of the open channel and revisiting them throughout the duration of the eruption.

The proximal (near-vent) sites record pahoehoe lava with little crust cover, while the distal sites capture behavior entirely in the `a`a flow regime. Sites within the braided section of the flow recorded video over parallel channels. Over 500 hover videos at the channel sites were acquired over the course of the Fissure 8 eruption between 30 May and 5 August. In this paper, we focus on videos collected at UAS site 8, capturing a junction point where the main channel split into two branches; see Fig. 4.

5.1.1 Velocity field measurements

We analyze the UAS hover videos using the optical flow technique (Horn and Schunck1981; Sun et al.2010) Optical flow is a well-known computer vision technique used to measure velocities of imaged objects based on the motion of brightness within an image sequence or between frames of a video. Lev et al. (2012) used Optical flow to measure the two-dimensional surface velocities of laboratory-scale basaltic lava flows. We follow the same technique as in Lev et al. (2012), tuning parameters to the specifics of the Kīlauea 2018 UAS footage; see Fig. 4. Length scales for video analysis and channel geometry data are provided from camera lens information and the recorded UAS flight altitude and refined using co-registered digital elevation and orthomosaic images produced from additional UAS data collected at the same or very close time.

Figure 5Finite-element partition of the modeled section of the braided channel system. Colors corresponds to (a) topography elevation in meters and (b) the magnitude of the gradient of the topography.


Figure 6Map view of modeled speed using a value of (a) n=1 and (b) n=2 in the viscosity model. Colors represent velocity magnitude in meters per second.


5.2 Model input

We provide our model with a channel geometry, assumed material properties, inlet velocity, and observed temperature. We use topography data from a pre-eruption digital elevation maps (data from the USGS National Elevation Dataset, with a spatial resolution of 10 m/pixel, USGS2002) to calculate the gradient of topography (see Fig. 5). We set the inlet velocity equal to values measured from the UAS video analyzed by optical flow (Fig. 4). Channel edge geometry is obtained from the velocity field (|u|>=0) combined with visible identification of channel boundaries in the UAS image. Figure 5 shows the meshed model domain, with colors depicting the elevation and ground slope used to set up the model.

We set the lava density to ρ=1350 kg m−3, which, with a nominal gas-free density of Hawaiian basalts of 2700 kg m−3 translates to 50 % vesicularity. We set the channel inlet temperature to T=1152C and wall and basal temperatures to Twall=1010C and Tground=477C, respectively. The rheological constants in Eq. (19) are set to A=-4.550, B=5805.30, and C=607.80. These values were calculated using the calculator by Giordano et al. (2008) and are specific for the composition of the basalt that erupted during June 2018 from Fissure 8 as measured by XRF analysis (Gansecki et al.2019). See Table 4 for the coefficient values used in the heat transfer module. It can be noted that because the lava temperature never falls below 950 C in the Kīlauea simulations, surface heat loss is solely due to radiation.

Figure 7Map view of the difference between modeled lava speed and surface speed obtained from UAS video capture for various power law exponents. Colors represent difference in meters per second.


Table 4Value of thermal coefficients.

Download Print Version | Download XLSX

5.3 Model results

All model simulations were executed on a Macbook Pro using Intel's Fortran compiler and on average took 56 min to execute using a time step dt=0.05 s. The finite-element partition of the braided channel system (shown in Fig. 5) consists of 6908 elements with a maximum element size of 8 m and a minimum element size of 1 m. The model's initial conditions were set to the solution of a linearized form of Eqs. (13)–(16) with ζ=0 and h=10 m. Inlet conditions are steady in time, and we set τyield=0 in the Herschel–Bulkley model and time step the nonlinear system of Eqs. (13)–(16) to steady state.

Figure 6 shows the lava velocities calculated for the entire domain by our model using exponent values in the viscosity model (Eq. 18) of n=1 (Newtonian) and n=2 (shear-thickening). Key features of the observed flow field, such as the increase in speed after the constriction in the northern branch and the stagnation at the channel split point, are present in both figures. The overall magnitude of the velocity – up to 9 m s−1 – is also in agreement with the observations.

Figure 8(a) Maximum error in model speed versus power law exponent. (b) Root-mean-square error in model speed versus power law exponent.


Figure 9Map view of the modeled lava thickness for various power law exponents. Colors represent lava thickness in meters. Notice that the lava becomes thinner as the value of n increases.


5.4 Discussion

We evaluate the quality of the fit between the model and the observations by comparing the modeled speed with the observed speed for different power law exponent (n) values shown in Fig. 7. Two areas of relatively large error are clear in the southern branch. We attribute these mostly to uncertainties in the underlying topography data. We use a coarse pre-eruption digital elevation map (DEM) for an area where the overall slopes are very gradual (2–3). The mesh and model resolution is very high compared to the coarseness of the DEM (only 10 DEM grid points across the model), which can lead to inaccuracies in slope estimates. In addition, the DEM is from before the eruption, while the velocity data were captured a few weeks after the channel was established. It is possible that by that time some lava had already deposited on the bottom of the channel and modified the topography.

Figure 10Map view of the effective modeled lava viscosity for various power law exponents. Colors represent viscosity (in Pascal seconds).


An additional source of misfit could be due to the bottom stress calculation. We calculate the thickness of the virtual layer over which the velocity transitions from the depth-integrated velocity to a value of zero (at the bottom boundary) via a two-layer model of vorticity. The two layers correspond to a mixed upper layer and a non-mixed bottom layer where the lava is losing heat due to conduction. While this approach seems to be valid in terms or reproducing lava flow thicknesses observed by USGS surveys (USGS Hawaii Volcano Observatory2019), which place lava thicknesses between 5 and 15 m, (alternative methods produce flows that are 3 times too thick), the two-layer model is still a simplification of reality that most likely introduces some errors.

Further, because we are limited to surface speed observations our error metrics will unavoidably have a misfit in them due to the fact that modeled speeds are depth-averaged. This effect will be small in regions of the channel where the Reynolds number is high because the lava speed will be more uniform over its thickness. However, in areas where the Reynolds number is low(er), model speeds will be lower than observations. This is because the section of channel we modeled has minimal crust cover, and therefore there is zero stress at the top boundary of the lava flow, meaning that lava speed should reach a maximum at the surface. This effect is evident in the northern portion of the northern channel and the portion of the southern channel between −154.886 and −154.8855 W.

An interesting aspect of the model worth drawing attention to is how a change in the value of n in the viscosity model affects numerical results. Figures 7 and 8 reveal that the overall fit improves with increasing n values, which corresponds to shear-thickening in our model. Similarly, the thicknesses predicted by our model for larger power law exponent values (Fig. 9) are closer to the range of thicknesses (5–15 m) measured by USGS survey (USGS Hawaii Volcano Observatory2019), and the apparent viscosity calculated by our model (shown in Fig. 10) for large n values is similar to rough estimates by the USGS (Weston Thelen, personal communication, 2018).

The lower error measures produced by our model for shear-thickening behavior is a departure from other studies that find lava to behave as a shear-thinning fluid (e.g., Castruccio et al.2010; Costa et al.2009; Pinkerton1995). The disparity can most likely be attributed to differences in the particle content of the lava; in the studies of Castruccio et al. (2010), Costa et al. (2009), and Pinkerton (1995), the crystal and bubble content of the lava studied was low (< 20 %), whereas samples taken from the Fissure 8 flow during and after the eruption show a wide range of crystallinity and vesicularity, often with very high vesicle fraction of over 50 % (Halverson et al.2020). The improved overall fit of our model for higher n values is most likely due to this high vesicle fraction and is consistent with the lab experiments of Smith (1997) and Lev et al. (2020) that show high bubble content can produce shear-thickening behavior. In fact, in the experiments of Sayag and Worster (2013), shear-thickening behavior for a constant volume of fluid (as in our investigation) produce flow thicknesses that are less than the shear-thinning case, which is visible in our results in Fig. 9. We conjecture that the increase in the viscosity is due to the fact that as the strain rate increases, the bubbles re-arrange in a fashion that makes it harder for the lava to flow. This effect should be especially pronounced in areas where the strain rate is high (think of the stress associated with solid boundaries in turbulence and how the vorticity created at these boundaries could cause bubbles to run into each other, impeding the flow). This effect is apparent in Fig. 10 for the case n=2, where the effective viscosity is low except in regions of high strain and high(er) Reynolds number, such as near the constriction in the northern channel, at the channel walls where the slope is high in the southern channel, and at the bend area in the southern channel. In the future we will explore mathematical relationships that allow n to be a function of the bubble and crystal content of the lava and examine the sensitivity of the model to non-Newtonian rheological parameters. We plan to infer the best-fitting values for these parameters for a range of locations and times for the Fissure 8 lava flow.

6 Conclusions

We present a novel numerical model for quantifying high-speed lava flows in complex channels. The mathematical model consists of depth-averaged equations of mass, momentum, and energy, and the equations are closed via a nonlinear viscosity model. Because we use discontinuous Galerkin methods to discretize the mathematical model, we are able to capture non-smooth transitions that can occur in lava flows, e.g., jumps in temperature, shear, and viscosity; see Fig. 10. We overcome a major limitation to many depth-integrated models in terms of the need to use an adjustable friction coefficient by solving a heat transfer boundary layer problem coupled with a calculation for the thickness of a virtual layer over which the velocity transitions from the depth-integrated velocity to a value of zero (at the bottom boundary) via a two-layer model of vorticity. This novel approach results in lava flow thicknesses that are in the range of observed values as compared to simple linear schemes that produce flow thicknesses that are 3 times too thick. Further, the use of unstructured triangular meshes allows the model to accurately resolve complex braided channel systems that are commonly produced by basaltic lava flows. This was demonstrated on a section of the complex braided channel system that was created by the Fissure 8 flow from the 2018 Kīlauea lower East Rift Zone eruption, with model results matching observational results quantitatively well. Future work will include using the new versatile model as a tool to infer lava properties and flux during volcanic crises.

Appendix A: List of symbols

Symbol Name Definition
A consistency constant see Sect. 2.1
Aw wetted area Aw=wH
αrs Runge–Kutta coefficient see Sect. 3.3
B consistency constant see Sect. 2.1
b discrete DG spatial matrix see Sect. 3.3
βrs Runge–Kutta coefficient see Sect. 3.3
C consistency constant see Sect. 2.1
dA differential area
dx0,,dx3 mesh size in verification test case
Δt numerical time step
δsΔt sub-time-step in Runge–Kutta method see Sect. 3.3
δz length of virtual bottom boundary layer δz=(δzx,δzy)
ϵ emissivity of lava see Table 4
F flux function matrix see Eq. (35)
F^ numerical flux function see Eq. (41)
fb body force vector fb=gx,gy,gz
fx x comp. of depth-integrated horizontal shear stress fx=-hζτxxx+τyxydz
fy y comp. of depth-integrated horizontal shear stress fy=-hζτxyx+τyyydz
gx, gy, gz x, y, z components of gravitational acceleration
H total depth of flow H=ζ+h
Hilbert space infinite dimensional space with finite energy
h steady reference depth of flow see Sect. 2.3
hb depth of bottom thermal boundary layer see Sect. 2.2
hw thickness of thermal boundary layer
I identity matrix 1000010000100001
i index counter for equations in Eqs. (13)–(16) i=1,,4
i^ unit vector in x coordinate direction
i imaginary number i=-1
J Jacobian matrix see Sect. 3.2.1
Jx Jacobian matrix in x direction see Sect. 3.2.1
Jy Jacobian matrix in y direction see Sect. 3.2.1
j index counter for the set of elements in 𝒯h j=1,,N
j^ unit vector in y coordinate direction
𝒦 lava consistency Pa s
𝒦0 lava consistency constant sn−1
k free surface wave number m−1
k^ unit vector in z coordinate direction
kb thermal conductivity of ground see Table 4
k̃b thermal conductivity of bottom boundary layer see Table 4
kc convection constant see Table 4
kT heat transfer coefficient
L2 L2 error norm ||||2=ΩjU-Uh2dA1/2
hp right hand side of ODE see Sect. 3.3
Λ( ) diagonal matrix of eigenvalues see Sect. 3.2.1
λ an eigenvalue of J see Sect. 3.2.1
λmax maximum eigenvalue of J see Sect. 3.2.1

Symbol Name Definition
l degree of freedom index for polynomial basis see Sect. 3
degree of polynomial basis see Sect. 3
M̃ finite-element mass matrix see Sect. 3.3
m number of eigenvalues of J see Sect. 3.2.1
μ viscosity Pa s
𝒩 number of elements in 𝒯h
n viscosity power law exponent see Sect. 2.1
n^ normal vector vector perpendicular to a plane
Ωj element j in 𝒯h
ω vorticity see Eq. (23)
ω free surface wave frequency s−1
P pressure flux P=12gzH2-h2
𝒫 polynomial space for basis functions see Sect. 3
Pw wetted perimeter Pw=w+2H
p pressure Pa
∂Ω channel domain boundary
∂Ωj boundary of element j in 𝒯h
Πh slope limiter see Cockburn and Shu (2001)
ϕ basis functions combinations of Jacobi polynomials
ψ test function ψ∈𝒱
ψh finite-element approximation of ψ see Eq. (39)
ψl polynomial coefficients of ψh see Eq. (39)
Qin channel inflow flux m3 s−1
q heat flux q=kTTx,kTTy,kTTz
q˙ internal heat generation and dissipation
q˙ internal heat generation and dissipation
qb bottom heat flux see Sect. 2.2
qi depth-integrated internal heat flux
qs surface heat flux see Sect. 2.2
R vector of DG spatial operator see Sect. 3.3
( ) matrix of right eigenvectors of J see Sect. 3.2.1
R()-1 matrix of left eigenvectors of J see Sect. 3.2.1
Re Reynold's number Re=ρ||u||AwμPw
r eigenvector of J see Sect.  3.2.1
ρ density kg m−3
r, s index counters in Runge–Kutta time-stepper
S source vector see Eq. (35)
𝒮 number of Runge–Kutta stages
sgn() sign of argument
σb Stefan–Boltzmann constant see Table 4
σwall normal stress vector at wall (-n^τ)
T temperature K
T depth-averaged temperature T=1H-hζTdz
T^ amplitude of T in verification test see Sect. 4
𝒯h finite-element triangulation of Ω Th=Ωjj=1N
Tatm atmospheric temperature K
Tground ground temperature K
t time s
t^ tangential vector vector parallel to a plane
ti initial time in averaging window s
tf final time in averaging window s

Symbol Name Definition
τ stress tensor see Sect. 2
τb bottom stress vector τb=τbx,τby
τs surface stress vector τs=τsx,τsy
σwall normal stress vector at wall (-n^τ)
(  ) transpose operator (u,v,w)=uvw
U solution vector see Eq. (35)
Ũ vector of polynomial coefficients (degrees of freedom) see Sect. 3.3
Uh finite-element approximation of U see Eq. (38)
Ul polynomial coefficients of Uh see Eq. (38)
u depth dependent velocity u=(u,v,w)
||u|| magnitude of velocity ||u||=(u2+v2+w2)1/2
u depth-averaged velocity vector u=1H-hζudz,1H-hζvdz
u x component of velocity u=dx/dt
u^ amplitude of u in verification test see Sect. 4
𝒱 admissible space of functions 𝒱∈ℋ
𝒱hp finite dimensional subspace 𝒱h∈𝒱
v y component of velocity v=dy/dt
v^ amplitude of v in verification test see Sect. 4
w z component of velocity w=dz/dt
w measure of depth-averaged vertical velocity see Eq. (24)
w̃ right-hand-side function in RK method see Sect. 3.3
win channel width at inlet
x Cartesian coordinate x=(x,y,z)
x0, xL x coordinate of inlet and outlet boundaries
y0, yL y coordinate of channel walls
z̃ thermal boundary layer coordinate z̃0,-hb
zb z coordinate where thermal boundary layer begins see Sect. 2.2
ζ free surface elevation ζ=H-h
ζ^ amplitude of surface elevation in verification test see Sect. 4
τyield yield strength of lava
gradient operator =x,y,z
horizontal gradient operator =x,y
0 vector of zeros 0,,0

Code availability

The computer code and data used in this investigation can be found at (Conroy2014).


The supplement related to this article is available online at:

Author contributions

CJC developed the model code and performed all model simulations. EL assisted in the unoccupied aerial system (UAS) data collection and utilized the optical flow technique to quantify lava velocities from the UAS data. CJC prepared the manuscript with contributions from EL.

Competing interests

The authors declare that they have no conflict of interest.

Financial support

This research has been supported by the National Science Foundation (grant no. EAR-1654588).

Review statement

This paper was edited by Andrew Wickert and reviewed by two anonymous referees.


Baloga, S., Spudis, P. D., and Guest, J. E.: The dynamics of rapidly emplaced terrestrial lava flows and implications for planetary volcanism, J. Geophys. Res.-Sol. Ea., 100, 24509–24519, 1995. a

Castruccio, A., Rust, A. C., and Sparks, R. S. J.: Rheology and flow of crystal-bearing lavas: Insights from analogue gravity currents, Earth Planet. Sc. Lett., 297, 471–480,, 2010. a, b, c

Castruccio, A., Rust, A. C., and Sparks, R. S. J.: Assessing lava flow evolution from post-eruption field data using Herschel-Bulkley rheology, J. Volcanol. Geo. Res., 275, 71–84,, 2014. a

Cockburn, B. and Chi-Wang, S.: The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM J. Num. Anal., 35, 2440–2463, 1998. a

Cockburn, B. and Shu, C.-W.: Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput., 16, 173–261, 2001. a, b, c, d, e

Connor, L. J., Connor, C. B., Meliksetian, K., and Savov, I.: Probabilistic approach to modeling lava flow inundation: a lava flow hazard assessment for a nuclear facility in Armenia, J. Appl. Volcanol., 1, 1–19, 2012. a

Conroy, C. J.: hp discontinuous Galerkin methods for coastal ocean circulation and transport, PhD thesis, The Ohio State University, 2014. a, b, c

Conroy, C. J.: coltonjconroy/DG_2d_lava_flows: DG_LAVA_2D (Version v1.0), Zenodo,, 2020. a

Conroy, C. J. and George, D. L.: Effective normal stress and excess pressure calculations in debris flows with dilatancy, Adv. Water Resour., submitted,, 2021. a

Conroy, C. J. and Kubatko, E. J.: hp discontinuous Galerkin methods for the vertical extent of the water column in coastal settings part I: Barotropic forcing, J. Comput. Phys., 305, 1147–1171, 2016. a, b

Conroy, C. J., Kubatko, E. J., and West, D. W.: ADMESH: An advanced, automatic unstructured mesh generator for shallow water models, Ocean Dynam., 62, 1503–1517, 2012. a

Conroy, C. J., Kubatko, E. J., Nappi, A., Sebian, R., West, D., and Mandli, K. T.: hp discontinuous Galerkin methods for parametric, wind-driven water wave models, Adv. Water Resour., 119, 70–83, 2018. a

Costa, A. and Macedonio, G.: Numerical simulation of lava flows based on depth-averaged equations, Geophys. Res. Lett., 32, L05304,, 2005. a, b

Costa, A., Caricchi, L., and Bagdassarov, N.: A model for the rheology of particle-bearing suspensions and partially molten rocks, Geochem. Geophy. Geosy., 10, Q03010,, 2009. a, b

Crisci, G. M., Rongo, R., Di Gregorio, S., and Spataro, W.: The simulation model SCIARA: the 1991 and 2001 lava flows at Mount Etna, J. Volcanol. Geotherm. Res., 132, 253–267,, 2004. a

Dawson, C. and Aizinger, V.: A discontinuous Galerkin method for three-dimensional shallow water equations, J. Sci. Comput., 22, 245–267, 2005. a

Dawson, C. and Mirabito, C. M.: The Shallow Water Equations, available at: (last access: 21 December 2020), 2008. a

Dawson, C., Kubatko, E. J., Westerink, J. J., Trahan, C., Mirabito, C., Michoski, C., and Panda, N.: Discontinuous Galerkin methods for modeling hurricane storm surge, Adv. Water Resour., 34, 1165–1176, 2011. a, b

Dietterich, H. R. and Cashman, K. V.: Channel networks within lava flows: Formation, evolution, and implications for flow behavior, J. Geophys. Res.-Earth Surf., 119, 1704–1724,, 2014. a

Gansecki, C., Lee, R. L., Shea, T., Lundblad, S. P., Hon, K., and Parcheta, C.: The tangled tale of Kīlauea’s 2018 eruption as told by geochemical monitoring, Science, 366, eaaz0147, 2019. a, b

George, D. L. and Iverson, R. M.: A depth averaged debris flow model that includes the effects of dilatancy. II. Numerical predictions and experimental tests, P. Roy. Soc. A-Math. Phy., 470, 1–31, 2014. a

Giordano, D., Russell, J., and Dingwell, D.: Viscosity of magmatic liquids: A model, Earth Planet. Sc. Lett., 271, 123–134, 2008. a, b, c

Griffiths, R. W.: The dynamics of lava flows, Ann. Rev. Fluid Mech., 32, 477–518, 2000. a, b, c, d, e

Halverson, B., Whittington, A., Hammer, J., deGraffenried, R., Lev, E., Dietterich, H., Patrick, M., Parcheta, C., Carr, B., Zoeller, M., Trusdell, F., and Llewellin, E.: Vesicularity and Rheology of the Kīlauea 2018 Lava Flows, Goldschmidt conference, 21–26 June 2020. a

Harris, A. J. and Rowland, S.: FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel, Bull. Volcanol., 63, 20–44, 2001. a

Harris, A. J. and Rowland, S. K.: FLOWGO 2012: An Updated Framework for Thermorheological Simulations of Channel-Contained Lava, Hawaiian Volcanoes: From Source to Surface, 208, 457–481, 2015. a, b

Heliker, C. and Mattox, T. N.: The First Two Decades of the Pu ‘u ‘O'ö-Küpaianaha, US Geological Survey Professional Paper, 1676, 2003. a

Herschel, W. H. and Bulkley, R.: Konsistenzmessungen von gummi-benzollösungen, Colloid Polym. Sci., 39, 291–300, 1926. a

Horn, B. K. and Schunck, B. G.: Determining optical flow, Artif. Intell., 17, 185–203, 1981. a

Kärnä, T., Kramer, S. C., Mitchell, L., Ham, D. A., Piggott, M. D., and Baptista, A. M.: Thetis coastal ocean model: discontinuous Galerkin discretization for the three-dimensional hydrostatic equations, Geosci. Model Dev., 11, 4359–4382,, 2018. a

Kelfoun, K. and Vargas, S. V.: VolcFlow capabilities and potential development for the simulation of lava flows, Geological Society, London, Special Publications, 426, SP426–8, 2015. a

Kubatko, E. J., Westerink, J. J., and Dawson, C.: hp discontinuous Galerkin methods for advection dominated problems in shallow water flow, Comput. Method. Appl. M., 196, 437–451, 2006. a, b, c, d, e, f, g

Kubatko, E. J., Bunya, S., Westerink, J. J., Dawson, C., and Mirabito, C.: A Performance Comparison of Continuous and Discontinuous Finite Element Shallow Water Models, J. Sci. Comput., 40, 315–339, 2009. a, b

Kubatko, E. J., Yeager, B. A., and Ketcheson, D. I.: Optimal strong-stability-preserving Runge–Kutta time discretizations for discontinuous Galerkin methods, J. Sci. Comput., 60, 313–344, 2014. a, b

Le Moigne, Y., Zurek, J. M., Williams-Jones, G., Lev, E., Calahorrano Di Patre,J., and Anzieta, J.: Standing waves in high speed lava channels: A tool forconstraining lava dynamics and eruptive parameters, J. Volcanol. Geoth. Res., 401, 1–15, 2020. a

Lev, E., Spiegelman, M., Wysocki, R. J., and Karson, J. A.: Investigating lava flow rheology using video analysis and numerical flow models, J. Volcanol. Geoth. Res., 247, 62–73, 2012. a, b

Lev, E., Birnbaum, J., Conroy, C. J., Whittington, A., Halverson, B., Hammer, J., and Llewellin, E.: The Rheology of Three-Phase Lavas and Magmas, Goldschmidt conference, 21–26 June 2020. a, b

LeVeque, R. J.: Finite volume methods for hyperbolic problems, vol. 31, Cambridge university press, 2002. a, b, c

Mader, H., Llewellin, E., and Mueller, S.: The rheology of two-phase magmas: A review and analysis, J. Volcanol. Geoth. Res., 257, 135–158, 2013. a

Manga, M., Castro, J., Cashman, K., and Loewenberg, M.: Rheology of bubble-bearing magmas, J. Volcanol. Geotherm. Res., 87, 15–28, 1998. a

Moran, M. J., Shapiro, H. N., Munson, B. R., and DeWitt, D. P.: Introduction to thermal systems engineering; Thermodynamics, fluid mechanics, and heat transfer, John Wiley & Sons, Inc., 2003. a

Neal, C., Brantley, S., Antolik, L., Babb, J., Burgess, M., Calles, K., Cappos, M., Chang, J., Conway, S., Desmither, L., Dotray, P., Elias, T., Fukunaga, P., Fuke, S., Johanson, I. A., Kamibayashi, K., Kauahikaua, J., Lee, R. L., Pekalib, S., Miklius, A., Million, W., Moniz, C. J., Nadeau, P. A., Okubo, P., Parcheta, C., Patrick, M. R., Shiro, B., Swanson, D. A., Tollett, W., Trusdell, F., Younger, E. F., Zoeller, M. H., Montgomery-Brown, E. K., Anderson, K. R., Poland, M. P., Ball, J. L., Bard, J., Coombs, M., Dietterich, H. R., Kern, C., Thelen, W. A., Cervelli, P. F., Orr, T., Houghton, B. F., Gansecki, C., Hazlett, R., Lundgren, P., Diefenbach, A. K., Lerner, A. H., Waite, G., Kelly, P., Clor, L., Werner, C., Mulliken, K., Fisher, G., and Damby, D.: The 2018 rift eruption and summit collapse of Kīlauea Volcano, Science, 363, 367–374, 2019. a

Patra, A. K., Nichitac, C. C., Bauera, A. C., Pitmanc, E. B., Bursikb, M., and Sheridan, M. F.: Parallel adaptive discontinuous Galerkin approximation for thin layer avalanche modeling, Comput. Geosci., 32, 912–926, 2006. a, b

Patrick, M., Dehn, J., and Dean, K.: Numerical modeling of lava flow cooling applied to the 1997 Okmok eruption: Approach and analysis, J. Geophys. Res.-Sol. Ea., 109, B03202,, 2004. a, b

Patrick, M. R., Dietterich, H. R., Lyons, J. J., Diefenbach, A. K., Parcheta, C., Anderson, K. R., Namiki, A., Sumita, I., Shiro, B., and Kauahikaua, J. P.: Cyclic lava effusion during the 2018 eruption of Kīlauea Volcano, Science, 366, eaay9070,, 2019. a, b, c

Pinkerton, H.: Rheological properties of basaltic lavas at sub-liquidus temperatures: laboratory and field measurements on lavas from Mount Etna, J. Volcanol. Geotherm. Res., 68, 307–323,, 1995. a, b

Rao, I. J. and Rajagopai, K. R.: The effect of the slip boundary condition on the flow of fluids in a channel, Acta Mech., 135, 113–126, 1999. a

De Saint-Venant, A.: Formules et tables nouvelles pour la solution des problèmes relatifs aux eaux courantes. Paris, France, Carilian-Goeury et Vor. Dalmont, 1851. a

Sayag, R. and Worster, M. G.: Axisymmetric gravity currents of power-law fluids over a rigid horizontal surface, Fluid Mech., 716, R5,, 2013. a

Schlichting, H., Gersten, K., Krause, E., Oertel, H., and Mayes, K.: Boundary-layer theory, vol. 7, Springer, 1968. a

Smith, J. V.: Shear thickening dilatancy in crystal-rich flows, J. Volcanol. Geotherm. Res., 79, 1–8, 1997. a

Sonder, I., Zimanowski, B., and Büttner, R.: Non-Newtonian viscosity of basaltic magma, Geophys. Res. Lett., 33, L02303,, 2006. a

Sun, D., Roth, S., and Black, M. J.: Secrets of optical flow estimation and their principles, in: 2010 IEEE computer society conference on computer vision and pattern recognition, IEEE, 2432–2439, 2010. a

USGS: National Elevation Dataset, available at: (last access: 10 April 2019), 2002. a

USGS Hawaii Volcano Observatory: Kīlauea 2018 lower East Rift Zone lava flow thicknesses: a Preliminary Map, available at: (last access: 4 June 2020), 2019. a, b

Vicari, A., Alexis, H., Del Negro, C., Coltelli, M., Marsella, M., and Proietti, C.: Modeling of the 2001 lava flow at Etna volcano by a Cellular Automata approach, Environ. Model. Softw., 22, 1465–1471, 2007. a

Wolfe, E. W.: The Puu Oo eruption of Kīlauea Volcano, Hawaii: Episodes 1–20, January 3, 1983–June 8, 1984., US Geol. Serv. Prof. Paper, 1463, 251 pp., 1988. a

Short summary
Lava flows present a natural hazard to communities around volcanoes and are usually slow-moving (< 1-5 cm/s). Lava flows during the 2018 eruption of Kilauea volcano, Hawai’i, however, reached speeds as high as 11 m/s. To investigate these dynamics we develop a new lava flow computer model that incorporates a nonlinear expression for the fluid viscosity. Model results indicate that the lava flows at Site 8 of the eruption displayed shear thickening behavior due to the flow's high bubble content.