Model description paper 11 Jun 2021
Model description paper  11 Jun 2021
A discontinuous Galerkin finiteelement model for fast channelized lava flows v1.0
 ^{1}LamontDoherty Earth Observatory, Columbia University, New York, NY, USA
 ^{2}Roy M. Huffington Department of Earth Sciences, Southern Methodist University, Dallas, TX, USA
 ^{1}LamontDoherty Earth Observatory, Columbia University, New York, NY, USA
 ^{2}Roy M. Huffington Department of Earth Sciences, Southern Methodist University, Dallas, TX, USA
Correspondence: Colton J. Conroy (cjc2235@columbia.edu)
Hide author detailsCorrespondence: Colton J. Conroy (cjc2235@columbia.edu)
Lava flows present a significant natural hazard to communities around volcanoes and are typically slowmoving (<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) finiteelement 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 $\mathcal{P}+\mathrm{1}/\mathrm{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.
 Article
(8286 KB) 
Supplement
(292 KB)  BibTeX
 EndNote
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 finiteelement 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 rootmeansquare 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.
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 Vargas, 2015) using either cellular automata or depthaveraged 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 higherspeed 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 Shallowwater equations for fast lava flows
Commensurate with this development in observational capabilities, we introduce a numerical method for modeling fastmoving 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 twodimensional depthintegrated equations for mass, momentum, and energy. In particular, we utilize a system of dynamical equations known as the “shallow water equations” (SaintVenant, 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 advectiondominated 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 nonNewtonian (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 depthintegrated energy equation that quantifies the thermal evolution of the lava as it interacts with its environment. The depthintegrated 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 nonNewtonian 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 depthaveraged 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 finiteelement method
Because closedform 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) finiteelement 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 finiteelement method differs from continuous Galerkin (CG) finiteelement 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 finitevolume 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 highperformance 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 onedimensional 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 $\mathcal{P}+\mathrm{1}/\mathrm{2}$ for polynomials of degree 𝒫.
Fluid flow on sloped terrain can be quantified in a Cartesian coordinate (x, y, z) system over a timedependent domain Ω(t)∈ℝ^{3} by solving Eulerian conservation equations of mass, linear momentum, and energy,
In Eqs. (1)–(3), ρ is the fluid density, $\mathit{u}=(u,v,w{)}^{\prime}$ is the fluid velocity vector, $\mathrm{\nabla}={\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\right)}^{\prime}$ is the gradient operator, p is the pressure, and τ is the stress term,
We denote the x component of gravitational acceleration that is tangential to the sloped surface by g_{x}; g_{y} is the y component of gravitational acceleration that is tangential to the sloped surface, and g_{z} is gravitational acceleration that is normal to the sloped surface. Collectively, these terms form the body force vector f_{b}= (g_{x}, g_{y}, g_{z})^{′} acting on the fluid. T is the temperature of the fluid, c_{p} is the specific heat capacity of the fluid, $\mathit{q}={\left({k}_{T}\frac{\partial T}{\partial x},{k}_{T}\frac{\partial T}{\partial y},{k}_{T}\frac{\partial T}{\partial z}\right)}^{\prime}$ is the heat flux through the fluid (k_{T} is a conduction heat transfer coefficient that measures the spread of heat within the fluid), and $\dot{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 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 $\mathit{u}\cdot \widehat{\mathit{n}}=\mathrm{0}$;
no pressure gradient, $\partial p/\partial \widehat{\mathit{n}}=\mathrm{0}$;
slip velocity, $\mathit{u}\cdot \widehat{\mathit{t}}=f({\mathit{\tau}}_{\text{wall}},{\mathit{\sigma}}_{\text{wall}})$;
heat loss via conduction, $\mathit{q}\cdot \widehat{\mathit{n}}=\left({k}_{T}/{h}_{\text{w}}\right)\left(T{T}_{\text{wall}}\right)$. 
Inlet boundary condition:
prescribed velocity, $\mathit{u}\cdot \widehat{\mathit{n}}=\text{prescribed}$;
prescribed pressure, p=prescribed;
prescribed heat content, ρc_{p}T=prescribed. 
Outlet boundary condition:
zero change in normal velocity, $\partial \mathit{u}/\partial \widehat{\mathit{n}}=\mathrm{0}$;
zero change in pressure, $\partial p/\partial \widehat{\mathit{n}}=\mathrm{0}$;
zero change in heat content, $\mathit{\rho}{c}_{p}\partial T/\partial \widehat{\mathit{n}}=\mathrm{0}$. 
Free surface boundary condition at z=ζ:
no relative normal flow,
$\partial \mathit{\zeta}/\partial t=u(\partial \mathit{\zeta}/\partial x)v\partial (\mathit{\zeta}/\partial y)+w$;
atmospheric pressure, p=p_{atm};
surface shear stress, ${\mathit{\tau}}_{{\text{s}}_{x}}={\mathit{\tau}}_{xx}(\partial \mathit{\zeta}/\partial x){\mathit{\tau}}_{yx}(\partial \mathit{\zeta}/\partial y)+{\mathit{\tau}}_{zx}$, ${\mathit{\tau}}_{{\text{s}}_{y}}={\mathit{\tau}}_{xy}(\partial \mathit{\zeta}/\partial x){\mathit{\tau}}_{yy}(\partial \mathit{\zeta}/\partial y)+{\mathit{\tau}}_{zy}$;
heat loss via radiation and convection,
$\mathit{q}\cdot \widehat{\mathit{n}}=\mathit{\u03f5}{\mathit{\sigma}}_{B}\left({T}^{\mathrm{4}}{T}_{\text{atm}}^{\mathrm{4}}\right)+{k}_{\text{c}}{\left(T{T}_{\text{atm}}\right)}^{\frac{\mathrm{4}}{\mathrm{3}}}$. 
Bottom boundary condition at $z=h$:
no slip velocity, u=0;
bottom shear stress, ${\mathit{\tau}}_{{\text{b}}_{x}}={\mathit{\tau}}_{xx}(\partial h/\partial x)+{\mathit{\tau}}_{yx}(\partial h/\partial y)+{\mathit{\tau}}_{zx}$, ${\mathit{\tau}}_{{\text{b}}_{y}}={\mathit{\tau}}_{xy}(\partial h/\partial x)+{\mathit{\tau}}_{yy}(\partial h/\partial y)+{\mathit{\tau}}_{zy}$;
heat loss via conduction, $\mathit{q}\cdot \widehat{\mathit{n}}=({k}_{\text{b}}/{h}_{\text{b}})\left(T{T}_{\text{ground}}\right)$.
In the boundary conditions above, $\widehat{\mathit{n}}$ is the unit normal vector to the wall boundary, outlet boundary, inlet boundary, moving free surface, and bottom boundary, respectively; $\widehat{\mathit{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 T_{w}, k_{w} is the thermal conductivity of the wall, and h_{w} 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, p_{atm}, along with a surface shear stress applied by the wind ${\mathit{\tau}}_{\text{s}}=({\mathit{\tau}}_{{\text{s}}_{x}},{\mathit{\tau}}_{{\text{s}}_{y}}{)}^{\prime}$. 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, T_{atm} is the temperature of the surrounding atmosphere, and k_{c} is the convection heat transfer coefficient. The main resisting force in dense shallow mass highspeed flows comes from the bottom stress, ${\mathit{\tau}}_{\text{b}}=({\mathit{\tau}}_{{\text{b}}_{x}},{\mathit{\tau}}_{{\text{b}}_{y}}{)}^{\prime}$, which is a function of the temperature of the lava at the basal boundary, where T_{ground} is the temperature of the ground and h_{b} 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 $\mathrm{\nabla}\cdot \mathit{u}=\mathrm{0}$, while assumptions (ii) and (iv) reduce the zmomentum 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 p_{atm} are negligible and the horizontal pressure gradient in Eq. (2) becomes
where $\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}={\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)}^{\prime}$. We further simplify the mathematical model and eliminate the vertical dimension by integrating $\mathrm{\nabla}\cdot \mathit{u}=\mathrm{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 freesurface and bottom boundary conditions, assume the density is constant, and simplify the resulting expression to arrive at the following depthintegrated equations:
where f_{x} and f_{y} are the x component and y component of the depthaveraged gradient of the shear stresses acting on vertical fluid planes, q_{s} is the heat flux through the free surface, q_{b} is the heat flux at the bottom boundary, $\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}\cdot {\stackrel{\mathrm{\u203e}}{q}}_{i}$ quantifies depthaveraged internal conduction, and $\dot{\stackrel{\mathrm{\u203e}}{q}}$ represents depthaveraged internal heat generation and dissipation. The depthaveraged velocity, $\stackrel{\mathrm{\u203e}}{\mathit{u}}=(\stackrel{\mathrm{\u203e}}{u},\stackrel{\mathrm{\u203e}}{v})$, and depthaveraged temperature, $\stackrel{\mathrm{\u203e}}{T}$, are defined as
The above system of depthintegrated equations (Eqs. 6–9) can be further simplified due to the dynamics of highspeed flows. More specifically, f_{x} and f_{y} are negligible in highspeed flows except at noslip and smallslip 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 (noflow) channel wall boundary condition and neglect f_{x} and f_{y} 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 noslip condition). The surface stress terms $({\mathit{\tau}}_{x},{\mathit{\tau}}_{y}{)}^{\prime}$ in the depthintegrated 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 $\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}\cdot {\stackrel{\mathrm{\u203e}}{q}}_{i}=\mathrm{0}$. 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 $\dot{\stackrel{\mathrm{\u203e}}{q}}$ in the current model but plan to include it in future releases.
It can be noted that the pressure gradient terms, $(\mathrm{1}/\mathit{\rho})\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}p={\left(H\frac{\partial}{\partial x}\left({g}_{z}\mathit{\zeta}\right),H\frac{\partial}{\partial x}\left({g}_{z}\mathit{\zeta}\right)\right)}^{\prime}$, in Eqs. (7) and (8) are nonconservative product terms that can lead to entropyviolating 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 $H(x,y,t)=\mathit{\zeta}(x,y,t)+h(x,y)$ and rewrite Eqs. (7) and (8) in the conservative form (Kubatko et al., 2006),
where $P=\frac{\mathrm{1}}{\mathrm{2}}{g}_{z}\left({H}^{\mathrm{2}}{h}^{\mathrm{2}}\right)$ is the pressure flux and $\partial h/\partial x$ and $\partial h/\partial 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 depthintegrated 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, ${\mathit{\tau}}_{\text{b}}=({\mathit{\tau}}_{{\text{b}}_{x}},{\mathit{\tau}}_{{\text{b}}_{y}}{)}^{\prime}$, 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 Bulkley, 1926),
where ${\mathit{\tau}}_{\text{b}}={\mathit{\tau}}_{zx}=({\mathit{\tau}}_{zx},{\mathit{\tau}}_{zy}{)}^{\prime}$, τ_{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 $\mathrm{log}\mathcal{K}/{\mathcal{K}}_{\mathrm{0}}$ at infinite temperature, 𝒦_{0} is a constant set to 1 s^{n−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 nonNewtonian 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 pseudoplastic 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 depthaveraged equations, and therefore we define the bottom stress in terms of the depthaveraged velocity as
where
It can be noted that in the expression above, ${\mathit{\delta}}_{z}=({\mathit{\delta}}_{{z}_{x}},{\mathit{\delta}}_{{z}_{y}}{)}^{\prime}$, 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 $\widehat{i}$, $\widehat{j}$, and $\widehat{k}$ are unit vectors in the x, y, and z directions, respectively. We solve an auxiliary problem over a pseudodepth of the lava that consists of an upper mixed layer where $\mathit{u}(x,y,z)\equiv \stackrel{\mathrm{\u203e}}{\mathit{u}}(x,y)$ and a lower layer where $(\partial u/\partial z,\partial v/\partial z{)}^{\prime}\gg (\partial w/\partial x,\partial w/\partial y{)}^{\prime}$. 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 $\partial \mathit{u}/\partial z\ne \mathrm{0}$ (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 depthaveraged velocity field, what is the associated length scale over which the depthaveraged 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 depthintegrated 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 $\widehat{i}$ and $\widehat{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,
coupled with the depthintegrated continuity Eq. (13) to obtain a measure of the vertical velocity w. More specifically, expanding derivatives in Eq. (13), solving for $\partial \mathit{\zeta}/\partial t$, while substituting this result into Eq. (23) and noting that in the upper layer $(\stackrel{\mathrm{\u203e}}{u},\stackrel{\mathrm{\u203e}}{v},\stackrel{\mathrm{\u203e}}{w})\equiv \left(u\right(\mathit{\zeta}),v(\mathit{\zeta}),w(\mathit{\zeta}\left)\right)$ yields
The relevant vorticity terms in Eq. (22) include the $\widehat{i}$ and $\widehat{j}$ components. By definition, $\partial \stackrel{\mathrm{\u203e}}{u}/\partial z=\partial \stackrel{\mathrm{\u203e}}{v}/\partial z=\mathrm{0}$ over the upper layer so that the vorticity component about the x axis is $\partial \stackrel{\mathrm{\u203e}}{w}/\partial y$ and the vorticity component about the y axis is $\partial \stackrel{\mathrm{\u203e}}{w}/\partial 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 $\partial v/\partial z$ about the x axis and $\partial u/\partial 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 $(\partial u/\partial z,\partial v/\partial z{)}^{\prime}$ 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 as $(\partial \stackrel{\mathrm{\u203e}}{w}/\partial x,\partial \stackrel{\mathrm{\u203e}}{w}/\partial y)$ goes to 0, the vertical stress in the fluid goes to 0. We can rewrite expression (Eq. 25) solely in terms of the depthaveraged variables using Eq. (24):
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 nonphysical and is not necessarily less than the lava flow thickness (H); it merely is a measure to ensure that $(\stackrel{\mathrm{\u203e}}{\mathit{u}}\mathit{0})/{\mathit{\delta}}_{z}\approx \partial \mathit{u}/\partial 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 $\widehat{i}$ and $\widehat{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 lowviscosity 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 T_{atm} 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), buoyancydriven 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 q_{s} in the energy equation to
where k_{c} 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),
where T_{ground}=f(x) is the temperature of the ground in contact with the lava flow field and k_{b} 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 depthdependent thermal boundary layer temperature, T(z),
where ${\stackrel{\mathrm{\u0303}}{k}}_{\text{b}}$ is the thermal boundary layer conductivity constant and h_{b}(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={z}_{\text{b}}(x,y)$ to $z=h(x,y)$ by setting z_{b}(x,y) to a relative zero. We then integrate Eq. (30) over a thermal boundary coordinate defined from $\stackrel{\mathrm{\u0303}}{z}=\mathrm{0}$ to $\stackrel{\mathrm{\u0303}}{z}={h}_{\text{b}}(x,y)$. It can be noted that the relationship between $\stackrel{\mathrm{\u0303}}{z}$ and z is given by $\stackrel{\mathrm{\u0303}}{z}=z{z}_{\text{b}}$ so that ${h}_{\text{b}}=h{z}_{\text{b}}$. Equation (30) is a nonhomogeneous, 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 ($\stackrel{\mathrm{\u0303}}{z}\in [\mathrm{0},{h}_{\text{b}}]$). Evaluating Eq. (31) at $\stackrel{\mathrm{\u0303}}{z}={h}_{\text{b}}$ and setting the interior temperature (T_{int}) to the depthintegrated value ($\stackrel{\mathrm{\u203e}}{T}$), 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 ${\mathbf{Q}}_{\text{in}}\cdot \widehat{\mathit{n}}$ is the inflow flux normal to the boundary and w_{in} 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 timeaveraged 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 ith row entries of the vectors U, S, and the flux function matrix F, defined as
where $P=\frac{\mathrm{1}}{\mathrm{2}}{g}_{z}\left({H}^{\mathrm{2}}{h}^{\mathrm{2}}\right)$.
3.1 Finiteelement 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 twodimensional domain Ω. The complexities of the domain boundary, ∂Ω, are such that an unstructured finiteelement 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 highquality 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 semidiscrete equations
Given the finiteelement 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 $\mathit{\psi}(x,y)\in \mathcal{V}$, integrate over each element Ω_{j}∈𝒯_{h}, and then integrate the flux term by parts,
for $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$ and $j=\mathrm{1},\mathrm{\dots},\mathcal{N}$, where 𝒩 is the total number of elements of the triangulation 𝒯_{h}. In the equation above, $\widehat{\mathit{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
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 $\mathit{\varphi}={\left({\mathit{\varphi}}_{\mathrm{0}},{\mathit{\varphi}}_{\mathrm{1}},\mathrm{\dots},{\mathit{\varphi}}_{\mathrm{\ell}}\right)}^{\prime}$, we express the trial solution (${U}_{h}^{\left(i\right)}\in {\mathcal{V}}_{hp}$) and test function (ψ_{h}∈𝒱_{hp}) as
and
where ${\left({U}_{\mathrm{0}}^{\left(i\right)},{U}_{\mathrm{1}}^{\left(i\right)},\mathrm{\dots},{U}_{\mathrm{\ell}}^{\left(i\right)}\right)}^{\prime}$ are the timedependent degrees of freedom of the finiteelement solution and $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$. We use products of Jacobi polynomials of degree ℓ, ${\left\{{\mathcal{P}}_{l}\right\}}_{l=\mathrm{0}}^{\mathrm{\ell}}$, as the basis for 𝒱_{hp}. The orthogonal triangular basis is defined in terms of a “collapsed coordinate” system that results in a matrixfree implementation of the method; see Kubatko et al. (2006) for more details. Substituting ${U}_{h}^{\left(i\right)}$ and ψ_{h} into Eq. (36), we arrive at the discrete weak form of the problem: find ${U}_{h}^{\left(i\right)}\in {\mathcal{V}}_{hp}$ such that for all test functions ψ_{h}∈𝒱_{hp} for $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$, the expression,
holds over each element Ω_{j}∈𝒯_{h}, where S^{(i)}(U_{h}) is the source term evaluated in 𝒱_{hp} and ${\widehat{\mathbf{F}}}^{\left(i\right)}$ 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 dualvalued (see Fig. 3, for example). To remedy this inconsistency, we replace the dualvalued flux in Eq. (36) with a socalled numerical flux ($\widehat{\mathbf{F}}$) that makes use of the left and right limits of the trial solution to produce a singlevalued flux across a given element's boundary.
More specifically, given an arbitrary function w_{h}∈𝒱_{hp} at an element boundary point x_{i}, we set the left and right limits of the function to ${w}_{h}^{}\equiv {w}_{h}\left({\mathbf{x}}_{i}^{}\right)$ and ${w}_{h}^{+}\equiv {w}_{h}\left({\mathbf{x}}_{i}^{+}\right)$, respectively. In this work we utilize the local Lax–Friedrichs (LLF) flux, which defines the numerical flux operator as
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 quasilinear form,
where the Jacobian matrices $\left({\mathbf{J}}_{ij}=\frac{\partial {F}_{i}}{\partial {x}_{j}}\right)$ are
and
The socalled “normal Jacobian matrix” is then defined by
where n_{x} and n_{y} are the x and y components of the normal edge vector $\widehat{\mathit{n}}$. 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 ${\mathcal{R}}_{(\cdot )}^{\mathrm{1}}$ is the matrix of left eigenvectors (LeVeque, 2002). To determine Λ_{(x)} and Λ_{(y)} we solve for the roots of det$\left({\mathbf{J}}_{(\cdot )}\mathit{\lambda}\mathcal{I}\right)=\mathrm{0}$, which gives the following eigenvalues,
Each eigenvector (${\mathit{r}}_{i}^{(\cdot )}$) can be determined by solving $({\mathbf{J}}_{(\cdot )}{\mathit{\lambda}}_{i}\mathbf{I}){\mathit{r}}_{i}=\mathbf{0}$, where I is the identity matrix, 0 is a vector of zeros, and ${\mathcal{R}}_{(\cdot )}=\left[{\mathit{r}}_{\mathrm{1}}^{(\cdot )},{\mathit{r}}_{\mathrm{2}}^{(\cdot )},{\mathit{r}}_{\mathrm{3}}^{(\cdot )},{\mathit{r}}_{\mathrm{4}}^{(\cdot )}\right]$. Solving for the eigenvectors we have
and
We use the full eigensystem in the slopelimiting 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 $({\mathit{\lambda}}_{\mathrm{1}},{\mathit{\lambda}}_{\mathrm{2}},{\mathit{\lambda}}_{\mathrm{3}},{\mathit{\lambda}}_{\mathrm{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 ChiWang, 1998) analogous to the method used in Conroy and Kubatko (2016) to evaluate secondorder derivative terms (see Conroy, 2014, for a detailed discussion on application of the LDG method to shallow mass geophysical fluid flows).
3.3 Strong stabilitypreserving 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}, ${\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{j}^{\left(i\right)}={\left(\phantom{\rule{0.33em}{0ex}}{U}_{j,\mathrm{0}}^{\left(i\right)},{U}_{j,\mathrm{1}}^{\left(i\right)},\mathrm{\dots},{U}_{j,\mathrm{\ell}}^{\left(i\right)}\phantom{\rule{0.33em}{0ex}}\right)}^{\prime}$ are vectors of the degrees of freedom (i.e., the polynomial basis coefficients), and ${\mathbf{b}}_{j}^{\left(i\right)}={\left(\phantom{\rule{0.33em}{0ex}}{\mathit{R}}_{j,\mathrm{0}}^{\left(i\right)},{\mathit{R}}_{j,\mathrm{1}}^{\left(i\right)},\mathrm{\dots},{\mathit{R}}_{j,\mathrm{\ell}}^{\left(i\right)}\phantom{\rule{0.33em}{0ex}}\right)}^{\prime}$ with,
${\stackrel{\mathrm{\u0303}}{\mathbf{M}}}_{j}^{\left(i\right)}$ 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 socalled strong stabilitypreserving (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 t_{m} to t_{m+1} via

Set ${\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{\mathrm{0}}^{\left(i\right)}\leftarrow {\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{m}^{\left(i\right)}$, for $i=\mathrm{1},\phantom{\rule{0.25em}{0ex}}\mathrm{2},\phantom{\rule{0.25em}{0ex}}\mathrm{3},\phantom{\rule{0.25em}{0ex}}\mathrm{4}$.

For each stage $r=\mathrm{1},\mathrm{2},\mathrm{\dots},\mathcal{S}$, set
$$\begin{array}{rl}& {\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{r}^{\left(i\right)}\leftarrow {\mathrm{\Pi}}_{h}\left(\sum _{s=\mathrm{1}}^{r}{\mathit{\alpha}}_{rs}{\mathbf{w}}^{rs}\right),\\ & {\mathbf{w}}^{rs}={\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{s\mathrm{1}}^{\left(i\right)}+{\displaystyle \frac{{\mathit{\beta}}_{rs}}{{\mathit{\alpha}}_{rs}}}\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\mathbf{L}}_{h}\phantom{\rule{0.25em}{0ex}}\left({\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{s\mathrm{1}},{t}_{m}+{\mathit{\delta}}_{\text{s}}\mathrm{\Delta}t\right).\end{array}$$ 
Finally, set ${\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{m+\mathrm{1}}^{\left(i\right)}\leftarrow {\stackrel{\mathrm{\u0303}}{\mathit{U}}}_{\mathcal{S}}^{\left(i\right)}$.
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 subtime 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;

${\sum}_{s=\mathrm{1}}^{r}\phantom{\rule{0.25em}{0ex}}{\mathit{\alpha}}_{rs}=\mathrm{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.
Verification of the DG solution of the mass and momentum equations in the depthaveraged and full threedimensional 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 (depthaveraged) 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 $\stackrel{\mathrm{\u203e}}{u}$, $\stackrel{\mathrm{\u203e}}{v}$, and H so that the depthaveraged mass equation is satisfied exactly. Specifically, we define
with h= constant, $\widehat{u}=$ constant, $\widehat{v}=\widehat{u}ky$, and $\widehat{\mathit{\zeta}}=\mathrm{exp}(\text{i}\mathit{\omega}\widehat{x})$ where $\widehat{x}=\mathrm{exp}\left(kx\right)/\left(k\widehat{u}\right)$. The dynamical solution consists of a wave propagating in a direction perpendicular to the y axis with wave number $k=\mathrm{5.0}\times {\mathrm{10}}^{\mathrm{3}}$ m^{−1} and frequency $\mathit{\omega}=\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{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
with $\widehat{T}=({y}^{\mathrm{2}}+{T}_{\text{wall}})$ 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 x_{0}=0.0 m, x_{L}=200.0 m, y_{0}=0.0 m, y_{L}=30.0 m. We assume symmetry about the centerline (at y=15 m) and only solve the equations over the halfwidth 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 socalled dx_{0} mesh), 3.75 m (dx_{1}), 2.50 m (dx_{2}), and 0.625 m (dx_{3}), 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 $\mathrm{\ell}+\mathrm{1}/\mathrm{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 higherorder polynomial approximation will result in a greater computational expense. However, the goal of using highorder 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 (highorder local polynomials produce more accurate results more efficiently than loworder 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 slowmoving, unusually highviscosity 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 km^{2} 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 (nearvent) 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 wellknown 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 twodimensional surface velocities of laboratoryscale 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 coregistered 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 preeruption 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 gasfree density of Hawaiian basalts of 2700 kg m^{−3} translates to 50 % vesicularity. We set the channel inlet temperature to $\stackrel{\mathrm{\u203e}}{T}=\mathrm{1152}$ ^{∘}C and wall and basal temperatures to T_{wall}=1010 ^{∘}C and T_{ground}=477 ^{∘}C, respectively. The rheological constants in Eq. (19) are set to $A=\mathrm{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.
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 finiteelement 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 (shearthickening). 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.
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 preeruption 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 depthintegrated velocity to a value of zero (at the bottom boundary) via a twolayer model of vorticity. The two layers correspond to a mixed upper layer and a nonmixed 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 twolayer 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 depthaveraged. 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 shearthickening 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 shearthickening behavior is a departure from other studies that find lava to behave as a shearthinning 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 shearthickening behavior. In fact, in the experiments of Sayag and Worster (2013), shearthickening behavior for a constant volume of fluid (as in our investigation) produce flow thicknesses that are less than the shearthinning 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 rearrange 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 nonNewtonian rheological parameters. We plan to infer the bestfitting 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 highspeed lava flows in complex channels. The mathematical model consists of depthaveraged 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 nonsmooth 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 depthintegrated 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 depthintegrated velocity to a value of zero (at the bottom boundary) via a twolayer 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.
Symbol  Name  Definition 
A  consistency constant  see Sect. 2.1 
A_{w}  wetted area  A_{w}=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  – 
$\text{d}{x}_{\mathrm{0}},\mathrm{\dots},\text{d}{x}_{\mathrm{3}}$  mesh size in verification test case  – 
Δt  numerical time step  – 
δ_{s}Δt  subtimestep in Runge–Kutta method  see Sect. 3.3 
δ_{z}  length of virtual bottom boundary layer  ${\mathit{\delta}}_{z}=({\mathit{\delta}}_{zx},{\mathit{\delta}}_{zy}{)}^{\prime}$ 
ϵ  emissivity of lava  see Table 4 
F  flux function matrix  see Eq. (35) 
$\widehat{F}$  numerical flux function  see Eq. (41) 
f_{b}  body force vector  ${\mathit{f}}_{\text{b}}={\left({g}_{x},{g}_{y},{g}_{z}\right)}^{\prime}$ 
f_{x}  x comp. of depthintegrated horizontal shear stress  ${f}_{x}={\int}_{h}^{\mathit{\zeta}}\frac{\partial {\mathit{\tau}}_{xx}}{\partial x}+\frac{\partial {\mathit{\tau}}_{yx}}{\partial y}\phantom{\rule{0.25em}{0ex}}\text{d}z$ 
f_{y}  y comp. of depthintegrated horizontal shear stress  ${f}_{y}={\int}_{h}^{\mathit{\zeta}}\frac{\partial {\mathit{\tau}}_{xy}}{\partial x}+\frac{\partial {\mathit{\tau}}_{yy}}{\partial y}\phantom{\rule{0.25em}{0ex}}\text{d}z$ 
g_{x}, g_{y}, g_{z}  x, y, z components of gravitational acceleration  – 
H  total depth of flow  $H=\mathit{\zeta}+h$ 
ℋ  Hilbert space  infinite dimensional space with finite energy 
h  steady reference depth of flow  see Sect. 2.3 
h_{b}  depth of bottom thermal boundary layer  see Sect. 2.2 
h_{w}  thickness of thermal boundary layer  – 
I  identity matrix  $\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]$ 
i  index counter for equations in Eqs. (13)–(16)  $i=\mathrm{1},\mathrm{\dots},\mathrm{4}$ 
$\widehat{i}$  unit vector in x coordinate direction  – 
i  imaginary number  $i=\sqrt{\mathrm{1}}$ 
J  Jacobian matrix  see Sect. 3.2.1 
J_{x}  Jacobian matrix in x direction  see Sect. 3.2.1 
J_{y}  Jacobian matrix in y direction  see Sect. 3.2.1 
j  index counter for the set of elements in 𝒯_{h}  $j=\mathrm{1},\mathrm{\dots},\mathcal{N}$ 
$\widehat{j}$  unit vector in y coordinate direction  – 
𝒦  lava consistency  Pa s 
𝒦_{0}  lava consistency constant  s^{n−1} 
k  free surface wave number  m^{−1} 
$\widehat{k}$  unit vector in z coordinate direction  – 
k_{b}  thermal conductivity of ground  see Table 4 
${\stackrel{\mathrm{\u0303}}{k}}_{\text{b}}$  thermal conductivity of bottom boundary layer  see Table 4 
k_{c}  convection constant  see Table 4 
k_{T}  heat transfer coefficient  
L^{2}  L^{2} error norm  $\left\right\cdot {}_{\mathrm{2}}={\left({\int}_{{\mathrm{\Omega}}_{j}}{\left(\mathit{U}{\mathit{U}}_{h}\right)}^{\mathrm{2}}\phantom{\rule{0.25em}{0ex}}\text{d}A\right)}^{\mathrm{1}/\mathrm{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 
$\stackrel{\mathrm{\u0303}}{\mathbf{M}}$  finiteelement 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 
$\widehat{\mathit{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=\frac{\mathrm{1}}{\mathrm{2}}{g}_{z}\left({H}^{\mathrm{2}}{h}^{\mathrm{2}}\right)$ 
𝒫  polynomial space for basis functions  see Sect. 3 
P_{w}  wetted perimeter  ${P}_{w}=w+\mathrm{2}H$ 
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}  finiteelement approximation of ψ  see Eq. (39) 
ψ_{l}  polynomial coefficients of ψ_{h}  see Eq. (39) 
Q_{in}  channel inflow flux  m^{3} s^{−1} 
q  heat flux  $\mathit{q}={\left({k}_{T}\frac{\partial T}{\partial x},{k}_{T}\frac{\partial T}{\partial y},{k}_{T}\frac{\partial T}{\partial z}\right)}^{\prime}$ 
$\dot{\mathit{q}}$  internal heat generation and dissipation  – 
$\dot{\stackrel{\mathrm{\u203e}}{q}}$  internal heat generation and dissipation  – 
q_{b}  bottom heat flux  see Sect. 2.2 
${\stackrel{\mathrm{\u203e}}{q}}_{i}$  depthintegrated internal heat flux  – 
q_{s}  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 
${\mathcal{R}}_{\left(\phantom{\rule{0.25em}{0ex}}\right)}^{\mathrm{1}}$  matrix of left eigenvectors of J  see Sect. 3.2.1 
Re  Reynold's number  $\mathit{Re}=\frac{\mathit{\rho}\left\right\mathit{u}\left\right{A}_{\text{w}}}{\mathit{\mu}{P}_{w}}$ 
r  eigenvector of J  see Sect. 3.2.1 
ρ  density  kg m^{−3} 
r, s  index counters in Runge–Kutta timestepper  – 
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  $(\widehat{\mathit{n}}\mathit{\tau}{)}^{\prime}$ 
T  temperature  K 
$\stackrel{\mathrm{\u203e}}{T}$  depthaveraged temperature  $\stackrel{\mathrm{\u203e}}{T}=\frac{\mathrm{1}}{H}{\int}_{h}^{\mathit{\zeta}}T\text{d}z$ 
$\widehat{T}$  amplitude of $\stackrel{\mathrm{\u203e}}{T}$ in verification test  see Sect. 4 
𝒯_{h}  finiteelement triangulation of Ω  ${\mathcal{T}}_{h}={\left\{{\mathrm{\Omega}}_{j}\right\}}_{j=\mathrm{1}}^{\mathcal{N}}$ 
T_{atm}  atmospheric temperature  K 
T_{ground}  ground temperature  K 
t  time  s 
$\widehat{\mathit{t}}$  tangential vector  vector parallel to a plane 
t_{i}  initial time in averaging window  s 
t_{f}  final time in averaging window  s 
Symbol  Name  Definition 
τ  stress tensor  see Sect. 2 
τ_{b}  bottom stress vector  ${\mathit{\tau}}_{\text{b}}={\left({\mathit{\tau}}_{{\text{b}}_{x}},{\mathit{\tau}}_{{\text{b}}_{y}}\right)}^{\prime}$ 
τ_{s}  surface stress vector  ${\mathit{\tau}}_{\text{s}}={\left({\mathit{\tau}}_{{\text{s}}_{x}},{\mathit{\tau}}_{{\text{s}}_{y}}\right)}^{\prime}$ 
σ_{wall}  normal stress vector at wall  $(\widehat{\mathit{n}}\mathit{\tau}{)}^{\prime}$ 
( )^{′}  transpose operator  $(u,v,w{)}^{\prime}=\left[\begin{array}{c}u\\ v\\ w\end{array}\right]$ 
U  solution vector  see Eq. (35) 
$\stackrel{\mathrm{\u0303}}{\mathit{U}}$  vector of polynomial coefficients (degrees of freedom)  see Sect. 3.3 
U_{h}  finiteelement approximation of U  see Eq. (38) 
U_{l}  polynomial coefficients of U_{h}  see Eq. (38) 
u  depth dependent velocity  $\mathit{u}=(u,v,w{)}^{\prime}$ 
$\left\right\mathit{u}\left\right$  magnitude of velocity  $\left\right\mathit{u}\left\right=({u}^{\mathrm{2}}+{v}^{\mathrm{2}}+{w}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ 
$\stackrel{\mathrm{\u203e}}{\mathit{u}}$  depthaveraged velocity vector  $\stackrel{\mathrm{\u203e}}{\mathit{u}}={\left(\frac{\mathrm{1}}{H}{\int}_{h}^{\mathit{\zeta}}u\phantom{\rule{0.25em}{0ex}}\text{d}z,\phantom{\rule{0.25em}{0ex}}\frac{\mathrm{1}}{H}{\int}_{h}^{\mathit{\zeta}}v\phantom{\rule{0.25em}{0ex}}\text{d}z\right)}^{\prime}$ 
u  x component of velocity  $u=\text{d}x/\text{d}t$ 
$\widehat{u}$  amplitude of $\stackrel{\mathrm{\u203e}}{u}$ in verification test  see Sect. 4 
𝒱  admissible space of functions  𝒱∈ℋ 
𝒱_{hp}  finite dimensional subspace  𝒱_{h}∈𝒱 
v  y component of velocity  $v=\text{d}y/\text{d}t$ 
$\widehat{v}$  amplitude of $\stackrel{\mathrm{\u203e}}{v}$ in verification test  see Sect. 4 
w  z component of velocity  $w=\text{d}z/\text{d}t$ 
$\stackrel{\mathrm{\u203e}}{w}$  measure of depthaveraged vertical velocity  see Eq. (24) 
$\stackrel{\mathrm{\u0303}}{\mathbf{w}}$  righthandside function in RK method  see Sect. 3.3 
w_{in}  channel width at inlet  – 
x  Cartesian coordinate  $\mathbf{x}=(x,y,z{)}^{\prime}$ 
x_{0}, x_{L}  x coordinate of inlet and outlet boundaries  – 
y_{0}, y_{L}  y coordinate of channel walls  – 
$\stackrel{\mathrm{\u0303}}{z}$  thermal boundary layer coordinate  $\stackrel{\mathrm{\u0303}}{z}\in \left[\mathrm{0},{h}_{\text{b}}\right]$ 
z_{b}  z coordinate where thermal boundary layer begins  see Sect. 2.2 
ζ  free surface elevation  $\mathit{\zeta}=Hh$ 
$\widehat{\mathit{\zeta}}$  amplitude of surface elevation in verification test  see Sect. 4 
τ_{yield}  yield strength of lava  – 
∇  gradient operator  $\mathrm{\nabla}={\left(\frac{\partial \phantom{\rule{0.25em}{0ex}}}{\partial x},\frac{\partial \phantom{\rule{0.25em}{0ex}}}{\partial y},\frac{\partial \phantom{\rule{0.25em}{0ex}}}{\partial z}\right)}^{\prime}$ 
$\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}$  horizontal gradient operator  $\stackrel{\mathrm{\u203e}}{\mathrm{\nabla}}={\left(\frac{\partial \phantom{\rule{0.25em}{0ex}}}{\partial x},\frac{\partial \phantom{\rule{0.25em}{0ex}}}{\partial y}\right)}^{\prime}$ 
0  vector of zeros  ${\left(\mathrm{0},\mathrm{\dots},\mathrm{0}\right)}^{\prime}$ 
The computer code and data used in this investigation can be found at https://doi.org/10.5281/zenodo.3863306 (Conroy, 2014).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1435532021supplement.
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. EAR1654588).
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 crystalbearing 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 posteruption field data using HerschelBulkley rheology, J. Volcanol. Geo. Res., 275, 71–84, https://doi.org/10.1016/j.jvolgeores.2014.02.004, 2014. a
Cockburn, B. and ChiWang, S.: The local discontinuous Galerkin method for timedependent convectiondiffusion systems, SIAM J. Num. Anal., 35, 2440–2463, 1998. a
Cockburn, B. and Shu, C.W.: Runge–Kutta discontinuous Galerkin methods for convectiondominated 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, https://doi.org/10.5281/zenodo.3863306, 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, 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, winddriven water wave models, Adv. Water Resour., 119, 70–83, 2018. a
Costa, A. and Macedonio, G.: Numerical simulation of lava flows based on depthaveraged equations, Geophys. Res. Lett., 32, L05304, https://doi.org/10.1029/2004GL021817, 2005. a, b
Costa, A., Caricchi, L., and Bagdassarov, N.: A model for the rheology of particlebearing 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/S03770273(03)003494, 2004. a
Dawson, C. and Aizinger, V.: A discontinuous Galerkin method for threedimensional 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. AMath. 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 thermorheological 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 ChannelContained 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 gummibenzollö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 threedimensional hydrostatic equations, Geosci. Model Dev., 11, 4359–4382, https://doi.org/10.5194/gmd1143592018, 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 strongstabilitypreserving Runge–Kutta time discretizations for discontinuous Galerkin methods, J. Sci. Comput., 60, 313–344, 2014. a, b
Le Moigne, Y., Zurek, J. M., WilliamsJones, 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 ThreePhase 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 twophase 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 bubblebearing 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., MontgomeryBrown, 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 subliquidus temperatures: laboratory and field measurements on lavas from Mount Etna, J. Volcanol. Geotherm. Res., 68, 307–323, https://doi.org/10.1016/03770273(95)000187, 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 SaintVenant, A.: Formules et tables nouvelles pour la solution des problèmes relatifs aux eaux courantes. Paris, France, CarilianGoeury et Vor. Dalmont, 1851. a
Sayag, R. and Worster, M. G.: Axisymmetric gravity currents of powerlaw fluids over a rigid horizontal surface, Fluid Mech., 716, R5, https://doi.org/10.1017/jfm.2012.545, 2013. a
Schlichting, H., Gersten, K., Krause, E., Oertel, H., and Mayes, K.: Boundarylayer theory, vol. 7, Springer, 1968. a
Smith, J. V.: Shear thickening dilatancy in crystalrich flows, J. Volcanol. Geotherm. Res., 79, 1–8, 1997. a
Sonder, I., Zimanowski, B., and Büttner, R.: NonNewtonian viscosity of basaltic magma, Geophys. Res. Lett., 33, L02303, https://doi.org/10.1029/2005GL024240, 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: http://nationalmap.usgs.gov/ (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: https://www.usgs.gov/maps/klauea2018lowereastriftzonelavaflowthicknessesapreliminarymap (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