A discontinuous Galerkin ﬁnite-element model for fast channelized lava ﬂows v1.0

. Lava ﬂows present a signiﬁcant natural hazard to communities around volcanoes and are typically slow-moving ( < 1 to 5 cm s − 1 ) and laminar. Recent lava ﬂows 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 ﬂows formed a complex network of braided channels departing from the classic rectangular channel geometry often employed by lava ﬂow models. To investigate these extreme dynamics we develop a new lava ﬂow model that incorporates nonlinear advection and a nonlinear expression for the ﬂuid viscosity. The model makes use of novel discontinuous Galerkin (DG) ﬁnite-element methods and resolves complex channel geometry through the use of unstructured triangular meshes. We verify the model against an analytic test case and demonstrate convergence rates of P + 1 / 2 for polynomials of degree P . 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.


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

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 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 Published by Copernicus Publications on behalf of the European Geosciences Union.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 ac-tive lava fields, we now have access to the data required by sophisticated numerical methods.

Shallow-water 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 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 the-ory 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 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).

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, andPatra 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 P + 1/2 for polynomials of degree P. https://doi.org/10.5194/gmd-14-3553-2021 Geosci.Model Dev., 14, 3553-3575, 2021 Fluid flow on sloped terrain can be quantified in a Cartesian coordinate (x, y, z) system over a time-dependent domain (t) ∈ R 3 by solving Eulerian conservation equations of mass, linear momentum, and energy, (3) In Eqs.(4) 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, 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 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 u • n = 0; no pressure gradient, ∂p/∂ n = 0; slip velocity, u • t = f (τ wall , σ wall ); heat loss via conduction, q • n = (k T /h w ) (T − T wall ).
-Inlet boundary condition: prescribed velocity, u • n = prescribed; prescribed pressure, p = prescribed; prescribed heat content, ρc p T = prescribed. - In the boundary conditions above, n is the unit normal vector to the wall boundary, outlet boundary, inlet boundary, moving free surface, and bottom boundary, respectively; t is the unit tangential vector to the wall; τ wall is the tangential shear stress at the wall; and σ wall is the normal stress at the wall.We denote the temperature at the wall by 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 τ s = (τ s x , τ s y ) .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 high-speed flows comes from the bottom stress, τ b = (τ b x , τ b y ) , 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 τ .(5) We assume that gradients in p atm are negligible and the horizontal pressure gradient in Eq. ( 2) becomes where ∇ = ∂ ∂x , ∂ ∂y .We further simplify the mathematical model and eliminate the vertical dimension by integrating ∇ • u = 0 and Eqs. ( 2) and (3) over the depth of the lava flow (from −h to ζ ).We then apply Leibniz's integral rule, utilize the free-surface and bottom boundary conditions, assume the density is constant, and simplify the resulting expression to arrive at the following depth-integrated equations: where f x and f y are the x component and y component of the depth-averaged 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, ∇ • q i quantifies depth-averaged internal conduction, and q represents depthaveraged internal heat generation and dissipation.The depthaveraged velocity, u = (u, v), and depth-averaged temperature, T , are defined as https://doi.org/10.5194/gmd-14-3553-2021Geosci.Model Dev., 14, 3553-3575, 2021 The above system of depth-integrated equations (Eqs.6-9) can be further simplified due to the dynamics of high-speed flows.More specifically, f x and f y are negligible in highspeed 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 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 no-slip condition).The surface stress terms (τ x , τ y ) 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 ∇ • q i = 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 q 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 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) = ζ (x, y, t) + h(x, y) and rewrite Eqs. ( 7) and ( 8) in the conservative form (Kubatko et al., 2006), where P = 1 2 g z H 2 − h 2 is the pressure flux and ∂h/∂x and ∂h/∂y quantify the gradient in the steady reference depth of flow that ζ is measured relative to.We supplement the system of equations given by Eqs. ( 13)-( 16) with initial conditions along with the channel wall, inlet, and outlet boundary conditions.It can be noted that the depth-integrated mass and momentum equations given by Eqs. ( 13)-( 15) are well studied in the literature and are commonly used to model shallow mass flows such as coastal ocean circulation and hurricane storm surge; see, for example, Dawson et al. (2011) and Kubatko et al. (2006).The addition of the energy equation complicates the solution of Eqs. ( 14) and ( 15) due to the fact that the bottom stress term, τ b = (τ b x , τ b y ) , is now a function of both nonlinear velocity gradients and temperature.

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 τ b = τ zx = (τ zx , τ zy ) , τ yield is the yield strength of the fluid, sgn denotes the sign of the argument, and µ is the nonlinear viscosity, defined as The symbol K 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 (K) 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 log K/K 0 at infinite temperature, K 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 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 where It can be noted that in the expression above, δ z = (δ z x , δ z y ) , 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 k are unit vectors in the x, y, and z directions, respectively.We solve an auxiliary problem over a pseudo-depth of the lava that consists of an upper mixed layer where u(x, y, z) ≡ u(x, y) and a lower layer where (∂u/∂z, ∂v/∂z) (∂w/∂x, ∂w/∂y) .We assume that the vorticity in the upper layer is equal to the vorticity in the bottom layer (in terms of magnitude) at the coordinate point where ∂u/∂z = 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 depth-averaged 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 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 ∂ζ /∂t, while substituting this result into Eq.( 23) and noting that in the upper layer (u, v, w) The relevant vorticity terms in Eq. ( 22) include the î and ĵ components.By definition, ∂u/∂z = ∂v/∂z = 0 over the upper layer so that the vorticity component about the x axis is ∂w/∂y and the vorticity component about the y axis is ∂w/∂x.Because the bottom boundary condition is modeled as a rigid wall where u = 0 and the fluid is incompressible, a vorticity layer forms in the fluid near the solid boundary that resists the local rotation of the fluid (this is the reason why the rigid boundary does not deform).The vorticity created at the boundary resists the rotation of the interior and is equal to ∂v/∂z about the x axis and ∂u/∂z about the y axis (see https://doi.org/10.5194/gmd-14-3553-2021 Geosci.Model Dev., 14, 3553-3575, 2021 Schlichting et al., 1968).Now, if we assume that each vorticity component over the bulk of the flow is equal to each vorticity component in the boundary layer at the coordinate point where (∂u/∂z, ∂v/∂z) is no longer equal to zero, then the virtual length over which the shear stress is applied is given by It can be noted that as (∂w/∂x, ∂w/∂y) goes to 0, the vertical stress in the fluid goes to 0. We can rewrite expression (Eq.25) solely in terms of the depth-averaged variables using Eq. ( 24): and It can be noted that even though we do not explicitly include horizontal shear stresses in the Kīlauea simulations presented in Sect. 5 (due to the high Re number), the virtual length used to quantify the bottom stress as defined in Eq. ( 26) is a function of horizontal shear within the fluid.Further, we wish to emphasize that the virtual length, δ z , is non-physical and is not necessarily less than the lava flow thickness (H ); it merely is a measure to ensure that (u − 0)/δ z ≈ ∂u/∂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 î and ĵ coordinate axis.

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 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), 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 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 depth-dependent thermal boundary layer temperature, T (z), where kb 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 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 z = 0 to z = −h b (x, y).It can be noted that the relationship between z and z is given by z = z − z b so that 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 (z ∈ [0, −h b ]).Evaluating Eq. ( 31) at z = −h b and setting the interior temperature (T int ) to the depth-integrated value (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.

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 Q in • 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 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).

Numerical discretization
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 ,

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 twodimensional 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 T 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.

A weak form and the semi-discrete equations
Given the finite-element partition, T h , of the domain , we obtain a weak form of Eq. ( 35) if we first multiply Eq. ( 35) by a sufficiently smooth test function ψ(x, y) ∈ V, integrate over each element j ∈ T h , and then integrate the flux term by parts, for i = 1, 2, 3, 4 and j = 1, . .., N , where N is the total number of elements of the triangulation T h .In the equation above, n is the outward unit normal to the element boundary ∂ j .
Rather than seek solutions to Eq. ( 36), we search for solutions in the finite dimensional subspace of functions defined as where P demarcates the space of polynomials of at most degree that is not necessarily continuous across element boundaries.In other words, given a set of basis functions φ = (φ 0 , φ 1 , . .., φ ) , we express the trial solution (U and where i) are the time-dependent degrees of freedom of the finite-element solution and i = 1, 2, 3, 4. We use products of Jacobi polynomials of degree , {P l } l=0 , as the basis for V 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 U (i) h and ψ h into Eq.( 36), we arrive at the discrete weak form of the problem: https://doi.org/10.5194/gmd-14-3553-2021 Geosci.Model Dev., 14, 3553-3575, 2021 find h ∈ V hp such that for all test functions ψ h ∈ V hp for i = 1, 2, 3, 4, the expression, holds over each element j ∈ T h , where S (i) (U h ) is the source term evaluated in V hp and F(i) is a suitably chosen numerical flux.

Numerical flux
The space of functions defined by Eq. ( 37) is not necessarily continuous across element boundaries, and thus can be dual-valued (see Fig. 3, for example).To remedy this inconsistency, we replace the dual-valued flux in Eq. ( 36) with a so-called numerical flux ( F) that makes use of the left and right limits of the trial solution to produce a single-valued flux across a given element's boundary.More specifically, given an arbitrary function w h ∈ V hp at an element boundary point x i , we set the left and right limits of the function to w − h ≡ w h (x − i ) and w + h ≡ w h (x + i ), 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 J ij = ∂F i ∂x j are , and The so-called "normal Jacobian matrix" is then defined by where n x and n y are the x and y components of the normal edge vector n.In general, if J is a square (m × m) matrix with m real eigenvalues, then it can be decomposed into its eigensystem, where R (•) is the matrix of right eigenvectors, (•) is the diagonal matrix of eigenvalues, and R −1 (•) is the matrix of left eigenvectors (LeVeque, 2002).To determine (x) and (y) we solve for the roots of det J (•) − λI = 0, which gives the following eigenvalues, λ 1,2 = un x + vn y , Each eigenvector (r i ) can be determined by solving (J (•) − λ i I)r i = 0, where I is the identity matrix, 0 is a vector of zeros, and 4 .Solving for the eigenvectors we have and 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 (λ 1 , λ 2 , λ 3 , λ 4 ).
It can be noted that to mathematically close the solution method of the discrete DG system of equations we need to numerically evaluate expression (Eq.26) to determine δ z .More specifically, we discretize Eq. ( 26) using a local discontinuous Galerkin (LDG) method (Cockburn and Chi-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).

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 N is the number of elements in T h , are vectors of the degrees of freedom (i.e., the polynomial basis coefficients), and b M(i) j 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 with i = 1, 2, 3, 4, and j = 1, . .., N , where L 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 t m to t m+1 via 1.Set 2. For each stage r = 1, 2, . .., S, set α rs w rs ,

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), S is the number of stages of the RK method, δ s t is a sub-time step of the time step t, and the α rs and β rs are coefficients that define the RK method.In particular, α rs and β rs conform to the following constraints: 1. α rs = 0 if and only if β rs = 0; 2. α rs ≥ 0 and β rs ≥ 0; 3. r s=1 α rs = 1.
Because we use explicit RK methods, the time step of the model is limited by a CFL condition; see Kubatko et al. (2014) for more details.

Verification
Verification of the DG solution of the mass and momentum equations in the depth-averaged and full three-dimensional case is well documented and can be found in Conroy and Kubatko (2016), Dawson and Aizinger (2005), and Kubatko et al. (2006).To verify our DG solution method for the fully coupled mass, momentum, and energy (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 u, v, and H so that the depth-averaged mass https://doi.org/10.5194/gmd-14-3553-2021 Geosci.Model Dev., 14, 3553-3575, 2021 equation is satisfied exactly.Specifically, we define with h = constant, û = constant, v = ûky, and ζ = exp(−iω x) where x = exp(kx)/(k û).The dynamical solution consists of a wave propagating in a direction perpendicular to the y axis with wave number k = 5.0 × 10 −3 m −1 and frequency ω = 2.0 × 10 −3 s −1 .These values were chosen based on velocity and lava flow thickness data recorded by Patrick et al. (2019) during the pulsing effusion regime associated with Fissure 8 during the 2018 Kīlauea event.We then set with T = (y 2 + T 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.0m, y 0 = 0.0 m, y L = 30.0m.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 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 + 1/2.Further, using P 2 polynomials on the coarsest mesh gives lower errors than P 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), andConroy et al. (2018).

Evaluation: recent eruption of Kīlauea volcano
We evaluate our model using data captured during the 2018 eruption of Kīlauea volcano, Hawai'i.The East Rift Zone of Kīlauea has erupted repeatedly in historical times and continuously since 1983 (Heliker and 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 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.

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 videog-raphy 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.

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. https://doi.org/10.5194/gmd-14-3553-2021 Geosci.Model Dev., 14, 3553-3575, 2021

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 T = 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 = −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.

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 (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.

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

Conclusions
We present a novel numerical model for quantifying highspeed 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 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 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.

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

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

Figure 3 .
Figure 3. Jump in numerical solution U h at an element edge ∂ j .

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

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

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

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

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

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

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

Table 4 .
Value of thermal coefficients.