A discontinuous Galerkin finite-element model for fast channelized lava flows v1.0
- 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 (email@example.com)
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 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.
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.1–1.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.
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 Vargas, 2015) 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 Cashman, 2014), 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 Rowland, 2015) 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-Venant, 1851), 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 Mirabito, 2008), although they have been used to model debris flows (George and Iverson, 2014) and lava flows (Costa and Macedonio, 2005) 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 Conroy, 2014, 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., Griffiths, 2000). 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 Shu, 2001), 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 George, 2021). 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 Shu, 2001).
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 Rowland, 2015). 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 for polynomials of degree 𝒫.
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,
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, is the heat flux through the fluid (kT is a conduction heat transfer coefficient that measures the spread of heat within the fluid), and 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 Conroy, 2014, 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 ;
no pressure gradient, ;
slip velocity, ;
heat loss via conduction, .
Inlet boundary condition:
prescribed velocity, ;
prescribed pressure, p=prescribed;
prescribed heat content, ρcpT=prescribed.
Outlet boundary condition:
zero change in normal velocity, ;
zero change in pressure, ;
zero change in heat content, .
Free surface boundary condition at z=ζ:
no relative normal flow,
atmospheric pressure, p=patm;
surface shear stress, , ;
heat loss via radiation and convection,
Bottom boundary condition at :
no slip velocity, u=0;
bottom shear stress, , ;
heat loss via conduction, .
In the boundary conditions above, is the unit normal vector to the wall boundary, outlet boundary, inlet boundary, moving free surface, and bottom boundary, respectively; 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 . 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, , 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.
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 , 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
We assume that gradients in patm are negligible and the horizontal pressure gradient in Eq. (2) becomes
where . We further simplify the mathematical model and eliminate the vertical dimension by integrating 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, quantifies depth-averaged internal conduction, and represents depth-averaged internal heat generation and dissipation. The depth-averaged velocity, , and depth-averaged temperature, , are defined as
The above system of depth-integrated equations (Eqs. 6–9) 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 Rajagopai, 1999, for an in depth investigation on channel wall boundary conditions in terms of the slip versus no-slip condition). The surface stress terms 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 . Internal heat generation and dissipation can be significant in lava flows with a high crystal content (Griffiths, 2000) and in lava flows in closed tubes (Costa and Macedonio, 2005). 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 in the current model but plan to include it in future releases.
It can be noted that the pressure gradient terms, , 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 LeVeque, 2002, for instance). To circumvent this issue, we make use of the fact and rewrite Eqs. (7) and (8) in the conservative form (Kubatko et al., 2006),
where is the pressure flux and and 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, , is now a function of both nonlinear velocity gradients and temperature.
2.1 Quantifying the bottom stress term
where , τyield is the yield strength of the fluid, sgn denotes the sign of the argument, and μ is the nonlinear viscosity, defined as
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),
where A is the value of 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
It can be noted that in the expression above, , 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
where , , and 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 and a lower layer where . 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 (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 and 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,
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 , while substituting this result into Eq. (23) and noting that in the upper layer yields
The relevant vorticity terms in Eq. (22) include the and components. By definition, over the upper layer so that the vorticity component about the x axis is and the vorticity component about the y axis is . 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 about the x axis and 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 is no longer equal to zero, then the virtual length over which the shear stress is applied is given by
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 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 and 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 Rowland, 2001). We quantify heat loss due to radiation via Stefan's law (Griffiths, 2000),
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) (Griffiths, 2000). In this case we set qs in the energy equation to
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),
where 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 to by setting zb(x,y) to a relative zero. We then integrate Eq. (30) over a thermal boundary coordinate defined from to . It can be noted that the relationship between and z is given by so that . Equation (30) is a non-homogeneous, constant coefficient ordinary differential equation that has the following solution:
which gives an expression for the temperature profile over the thermal boundary layer of the lava (). Evaluating Eq. (31) at and setting the interior temperature (Tint) to the depth-integrated value (), we have
We use this temperature to evaluate the consistency in the bottom stress approximation,
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.,
where 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).
To develop our numerical methods, we rewrite the system of Eqs. (1) in the compact form,
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
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 , integrate over each element Ωj∈𝒯h, and then integrate the flux term by parts,
for and , where 𝒩 is the total number of elements of the triangulation 𝒯h. In the equation above, 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
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 , we express the trial solution () and test function (ψh∈𝒱hp) as
where are the time-dependent degrees of freedom of the finite-element solution and . We use products of Jacobi polynomials of degree ℓ, , 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 and ψh into Eq. (36), we arrive at the discrete weak form of the problem: find such that for all test functions ψh∈𝒱hp for , the expression,
holds over each element Ωj∈𝒯h, where S(i)(Uh) is the source term evaluated in 𝒱hp and 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 () 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.
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 and , respectively. In this work we utilize the local Lax–Friedrichs (LLF) flux, which defines the numerical flux operator as
where the Jacobian matrices are
The so-called “normal Jacobian matrix” is then defined by
where nx and ny are the x and y components of the normal edge vector . In general, if J is a square (m×m) matrix with m real eigenvalues, then it can be decomposed into its eigensystem,
where ℛ(⋅) is the matrix of right eigenvectors, Λ(⋅) is the diagonal matrix of eigenvalues, and is the matrix of left eigenvectors (LeVeque, 2002). To determine Λ(x) and Λ(y) we solve for the roots of det, which gives the following eigenvalues,
Each eigenvector () can be determined by solving , where I is the identity matrix, 0 is a vector of zeros, and . 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 Shu, 2001), and we set λmax in the LLF flux to the maximum value of .
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-Wang, 1998) analogous to the method used in Conroy and Kubatko (2016) to evaluate second-order derivative terms (see Conroy, 2014, 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,
where 𝒩 is the number of elements in 𝒯h, are vectors of the degrees of freedom (i.e., the polynomial basis coefficients), and with,
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
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
Set , for .
For each stage , set
Finally, set .
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 Shu, 2001), 𝒮 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:
αrs=0 if and only if βrs=0;
αrs≥0 and βrs≥0;
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 Griffiths, 2000, 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 , , and H so that the depth-averaged mass equation is satisfied exactly. Specifically, we define
with h= constant, constant, , and where . The dynamical solution consists of a wave propagating in a direction perpendicular to the y axis with wave number m−1 and frequency 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
with and substitute Eqs. (49) and (50) into the mathematical model (Eqs. 13–16) 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 1–3, where it can be noted that the method converges to the analytic solution at a rate of approximately . 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).
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 Mattox, 2003; Wolfe, 1988). 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 Schunck, 1981; 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.
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, USGS, 2002) 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 ∘C and wall and basal temperatures to Twall=1010 ∘C and Tground=477 ∘C, respectively. The rheological constants in Eq. (19) are set to , 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.
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.
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.
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 Observatory, 2019), 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 Observatory, 2019), 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; Pinkerton, 1995). 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.
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.
|A||consistency constant||see Sect. 2.1|
|α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|
|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|
|ϵ||emissivity of lava||see Table 4|
|F||flux function matrix||see Eq. (35)|
|numerical flux function||see Eq. (41)|
|fb||body force vector|
|fx||x comp. of depth-integrated horizontal shear stress|
|fy||y comp. of depth-integrated horizontal shear stress|
|gx, gy, gz||x, y, z components of gravitational acceleration||–|
|H||total depth of flow|
|ℋ||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||index counter for equations in Eqs. (13)–(16)|
|unit vector in x coordinate direction||–|
|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|
|unit vector in y coordinate direction||–|
|𝒦||lava consistency||Pa s|
|𝒦0||lava consistency constant||sn−1|
|k||free surface wave number||m−1|
|unit vector in z coordinate direction||–|
|kb||thermal conductivity of ground||see Table 4|
|thermal conductivity of bottom boundary layer||see Table 4|
|kc||convection constant||see Table 4|
|kT||heat transfer coefficient|
|L2||L2 error norm|
|ℒ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|
|l||degree of freedom index for polynomial basis||see Sect. 3|
|ℓ||degree of polynomial basis||see Sect. 3|
|finite-element mass matrix||see Sect. 3.3|
|m||number of eigenvalues of J||see Sect. 3.2.1|
|𝒩||number of elements in 𝒯h||–|
|n||viscosity power law exponent||see Sect. 2.1|
|normal vector||vector perpendicular to a plane|
|Ωj||element j in 𝒯h||–|
|ω||vorticity||see Eq. (23)|
|ω||free surface wave frequency||s−1|
|𝒫||polynomial space for basis functions||see Sect. 3|
|∂Ω||channel domain boundary||–|
|∂Ωj||boundary of element j in 𝒯h||–|
|Πh||slope limiter||see Cockburn and Shu (2001)|
|ϕ||basis functions||combinations of Jacobi polynomials|
|ψh||finite-element approximation of ψ||see Eq. (39)|
|ψl||polynomial coefficients of ψh||see Eq. (39)|
|Qin||channel inflow flux||m3 s−1|
|internal heat generation and dissipation||–|
|internal heat generation and dissipation||–|
|qb||bottom heat flux||see Sect. 2.2|
|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|
|matrix of left eigenvectors of J||see Sect. 3.2.1|
|r||eigenvector of J||see Sect. 3.2.1|
|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|
|amplitude of in verification test||see Sect. 4|
|𝒯h||finite-element triangulation of Ω|
|tangential vector||vector parallel to a plane|
|ti||initial time in averaging window||s|
|tf||final time in averaging window||s|
|τ||stress tensor||see Sect. 2|
|τb||bottom stress vector|
|τs||surface stress vector|
|σwall||normal stress vector at wall|
|( )′||transpose operator|
|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|
|magnitude of velocity|
|depth-averaged velocity vector|
|u||x component of velocity|
|amplitude of in verification test||see Sect. 4|
|𝒱||admissible space of functions||𝒱∈ℋ|
|𝒱hp||finite dimensional subspace||𝒱h∈𝒱|
|v||y component of velocity|
|amplitude of in verification test||see Sect. 4|
|w||z component of velocity|
|measure of depth-averaged vertical velocity||see Eq. (24)|
|right-hand-side function in RK method||see Sect. 3.3|
|win||channel width at inlet||–|
|x0, xL||x coordinate of inlet and outlet boundaries||–|
|y0, yL||y coordinate of channel walls||–|
|thermal boundary layer coordinate|
|zb||z coordinate where thermal boundary layer begins||see Sect. 2.2|
|ζ||free surface elevation|
|amplitude of surface elevation in verification test||see Sect. 4|
|τyield||yield strength of lava||–|
|horizontal gradient operator|
|0||vector of zeros|
The supplement related to this article is available online at: https://doi.org/10.5194/gmd-14-3553-2021-supplement.
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.
The authors declare that they have no conflict of interest.
This research has been supported by the National Science Foundation (grant no. EAR-1654588).
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, https://doi.org/10.1016/j.epsl.2010.06.051, 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, https://doi.org/10.1016/j.jvolgeores.2014.02.004, 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
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. and George, D. L.: Effective normal stress and excess pressure calculations in debris flows with dilatancy, Adv. Water Resour., submitted, https://doi.org/10.13140/RG.2.2.34282.67527, 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., Caricchi, L., and Bagdassarov, N.: A model for the rheology of particle-bearing suspensions and partially molten rocks, Geochem. Geophy. Geosy., 10, Q03010, https://doi.org/10.1029/2008GC002138, 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, https://doi.org/10.1016/S0377-0273(03)00349-4, 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: https://users.oden.utexas.edu/~arbogast/cam397/dawson_v2.pdf (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, https://doi.org/10.1002/2014JF003103, 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 https://doi.org/10.1126/science.aaz0147, 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
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
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, https://doi.org/10.5194/gmd-11-4359-2018, 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
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
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, https://doi.org/10.1029/2003JB002537, 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, https://doi.org/10.1126/science.aay9070, 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, https://doi.org/10.1016/0377-0273(95)00018-7, 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
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
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 Hawaii Volcano Observatory: Kīlauea 2018 lower East Rift Zone lava flow thicknesses: a Preliminary Map, available at: https://www.usgs.gov/maps/k-lauea-2018-lower-east-rift-zone-lava-flow-thicknesses-a-preliminary-map (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