The influence of temperaturedependent viscosity on

0 downloads 0 Views 888KB Size Report
Aug 19, 2013 - refer to as preferred pathways) can affect lava emplacement, and allow flows to reach ...... begin diverging strongly due to the development of trans- ...... Anderson, S. W., S. M. McColley, J. H. Fink, R. K. Hudson (2005), The.
JOURNAL OF GEOPHYSICAL RESEARCH: EARTH SURFACE, VOL. 118, 1516–1532, doi:10.1002/jgrf.20111, 2013

The inf luence of temperature-dependent viscosity on lava f low dynamics S. Diniega,1 S. E. Smrekar,1 S. Anderson,2 and E. R. Stofan 3 Received 30 November 2012; revised 9 July 2013; accepted 11 July 2013; published 19 August 2013.

[1] Many studies have investigated how channels and tubes form within basaltic flows. The

partitioning of a broad lava flow into narrow regions of enhanced flow velocity (which we refer to as preferred pathways) can affect lava emplacement, and allow flows to reach great lengths on Earth and other planets. In this study, we investigate the role that a dynamic instability, driven by large changes in viscosity due to small changes in temperature, can play in the formation and propagation of low-viscosity preferred pathways. We use a fluid dynamics-based model of laminar lava flow emplacement that examines this process exclusively. Analysis of this model shows that preferred pathways will initiate when the temperature dependence of the lava’s viscosity is sufficiently strong. If a preferred pathway develops, it will form and stabilize over a distance that reflects the competition between the cooling rate through the flow’s surfaces and the flow velocity within the preferred pathway. Using velocity and cooling rates of typical basaltic flows on Earth, this model is able to produce preferred pathways within lava flows with reasonable dimensions and timescales. The model also fails to produce such pathways in situations where tubes and channels are not found in nature (e.g., in basalt sheet flows and silicic flows < 1 km in length). These results are sufficiently positive to demonstrate that this thermodynamic instability may play a role in the formation and propagation of preferred pathways within some flows, and thus merit further investigation.

Citation: Diniega, S., S. E. Smrekar, S. Anderson, and E. R. Stofan (2013), The influence of temperature-dependent viscosity on lava flow dynamics, J. Geophys. Res. Earth Surf., 118, 1516–1532, doi:10.1002/jgrf.20111.

1.

Background

[2] Many effusive lava flows begin as a broad, simple flow containing a molten core that extends nearly across the width of the flow [e.g., Greeley, 1971; Kauahikaua et al., 1998]. Some of these flows remain undivided and during emplacement will spread as a laterally coherent sheet. Within other flows, the molten core becomes divided within a single flow unit, with regions of faster flow separated by slower or stagnant regions. For example, lava channels and tubes are obvious and distinct regions of enhanced flow. In addition to channels and tubes, lava flows can also contain less obvious “preferred pathways” where enhanced flow is differentiated from the surrounding flow by only slight differences in velocity, temperature, and viscosity. [3] Understanding the formation of these types of enhanced flow regions within a lava flow is important because the 1 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, California, USA. 2 MAST Institute and Earth Sciences Department, University of Northern Colorado, Greeley, Colorado, USA. 3 Proxemy Research, Rectortown, Virginia, USA.

Corresponding author: S. Diniega, Jet Propulsion Laboratory, California Institute of Technology, M/S 301-250D, 4800 Oak Grove Dr., Pasadena, CA 91109, USA. ([email protected]) ©2013. American Geophysical Union. All Rights Reserved. 2169-9003/13/10.1002/jgrf.20111

formation of such regions within a flow greatly affects lava emplacement. In particular, enhanced flow regions can make possible the great lengths of some lava flows because they serve as vital mechanisms for efficient lava transport [Duraiswami et al., 2004; Hon et al., 1994; Keszthelyi and Self, 1998]. Lava flows containing tubes can remain hot and travel long distances without requiring extremely high effusion rates [Greeley, 1987; Hon et al., 1994; Peterson et al., 1994]. Within many lava flows in Hawaii, for example, lava tubes regularly facilitate transport of basaltic lava over distances of 10–20 km with temperature drops of ≤1°C/km [Kauahikaua et al., 2003]. [4] In this study, we generally refer to all narrow regions of enhanced flow—whether distinctly or subtly divided from the rest of the flow—as “preferred pathways.” We do this because we hypothesize that all of these regions of enhanced flow can be initiated due to the same flow-focusing dynamic (although we recognize that the subsequent evolution processes and end-morphology differs between preferred pathways, channels, and tubes). [5] Here, we investigate how a broad lava flow can become focused into narrow preferred pathways. Several hypotheses have been proposed for the formation of preferred pathways. Portions of a flow that are thicker or that go through topographic lows are more likely to contain regions of sustained flow compared to thinner portions that cool and stagnate [Glaze et al., 2005; Greeley and Hyde, 1972; Hon et al., 1994]. It has also been proposed that within anastomosing

1516

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Figure 1. Correlated measurements of temperature and viscosity for natural melt lavas, from Blatt et al. [2006]. We assume that viscosity has an exponential dependence on temperature: μ ~ exp(β(1  T)), where T is the dimensionless temperature within the temperature range of interest, TC to TH (i.e., T = (Tdim  TC)/(TH  TC)). The curves plotted yield β-values of approximately 3.5, 2.3, 1.4, 1.4, and 0.4 (starting with the upper/rhyolite curve and descending in order) when considering viscosity changes with TH  TC = 100°C. These laboratory measurements were generated at 1 bar air pressure and with volatile-free samples. As discussed within the text, β values are likely much higher within natural lava flows. sheet-flow lobes, cooling from the sides of individual lobes will focus flow into central regions that then connect and coalesce into tubes [Greeley, 1987; Hon et al., 1994; Peterson et al., 1994; Rowland and Walker, 1990]. A dynamic-instability hypothesis has also been proposed [Anderson et al., 1999, 2005, 2012] that is based on the well-established fact that lava viscosity varies exponentially with temperature (Figure 1). Our study examines the viability of this third hypothesis within natural lava flows. [6] Experiments and models with two immiscible fluids of different viscosities, where the less viscous fluid is injected into the more viscous fluid, show formation of stable low-viscosity “fingers” rather than a smooth penetration front [e.g., Saffman and Taylor, 1958]. A similar dynamic occurs within a flow composed of a single fluid with a temperature-dependent viscosity. Laboratory experiments using liquid paraffin [Whitehead and Helfrich, 1991], karo syrup [Wylie et al., 1999], glycerine [Holloway and de Bruyn, 2005], and polyethylene glycol wax [Anderson et al., 2005] have shown that small temperature variations within an initially uniform-velocity flow can cause flow to focus into one or a few preferred pathways. A model of this thermo-viscous instability has been used to explain how fissure eruptions evolve from a curtain of lava to a few localized eruption centers [Costa and Macedonio, 2002; Helfrich, 1995; Wylie and Lister, 1995; Wylie et al., 1999]. [7] In this study, we investigate the way in which small thermal variations within the source flow can lead to spatial variations in viscosity, and the effect that these viscosity

variations can have on the formation of preferred pathways within a lava flow. To evaluate the effect of the thermo-viscous instability in isolation from other pathway formation mechanisms, we model the effect of this type of instability on lava flows moving over featureless terrain. We use a simple two-dimensional fluid dynamics-based model that includes an exponential relationship between temperature and viscosity. Although topography will certainly influence the specific path of a preferred pathway, our model describes the conditions leading to a sufficient viscosity contrast which allows preferred pathways to form and stabilize. Furthermore, using observed ranges of flow velocity and cooling rate from basalt flows that contain channels or tubes, we calculate both the along-flow distance over which a preferred pathway forms and widths of stable preferred pathways. By comparing these model-predicted preferred-pathway dimensions with reported dimensions of typical flow channels, we test the viability of our hypothesis. Furthermore, we investigate the values of model parameters that inhibit the thermo-viscous effect and prevent the generation of preferred pathways within a flow, and use this information to propose explanations for why some flows do not contain channels and continue to spread like a sheet. [8] Our study focuses mainly on basalt flows because tubes and channels form most commonly in these flows. We hypothesize that channels and tubes can initialize as a preferred pathway and become defined as continued flow focusing stabilizes the pathway. Observations have shown that this stabilization can occur very quickly—channels can appear within a few tens of minutes within basaltic flow [Peterson and Swanson, 1974]. Sustained flow over hours to days may allow some of the preferred pathways to evolve into lava tubes [Byrnes and Crown, 2001; Peterson et al., 1994], with large tubes forming over a few weeks [Calvari and Pinkerton, 1998; Hon et al., 1994; Peterson et al., 1994] to months [Harris et al., 1997]. Additionally, stabilized preferred pathways can create an uneven lava surface as a result of differential pressures within the flow. This provides one mechanism for the generation of hummocky topography [Anderson et al., 1999; Hon et al., 1994; Self et al., 1998; Swanson, 1973]. In this study, we also briefly discuss silicic flows because these flows sometimes contain channels. [9] Our model predicts that interior lava flow dynamics induced by temperature-dependent viscosity can lead to the development of preferred pathways. As discussed below, all basalt flows will cool to a point at which they have sufficiently large viscosity contrasts for the thermoviscous instability to influence preferred pathway formation within the main core of the flow. But other effects, such as the influence of topography or changes in effusion rate, can have a stronger influence on flow dynamics. Our results are sufficiently positive, however, to suggest that this dynamic instability may be important in some flows, such as when effusion rates are high or when flows are thick and wide relative to topographic variations. 1.1. Types of Flows Investigated [10] We recognize that all lava flows will have their flow direction and thickness influenced by the presence of pre-flow topography. Flows that are of the scale of the preexisting topography will experience flow focusing due to the topography, and they are likely to contain channels or tubes that

1517

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

follow the local gradient. However, preexisting topography does not explain the formation of preferred pathways in flows with high effusion rates or of sufficient size to be unconstrained by local topography. Within these types of flows, we hypothesize that while topography will influence flow thickness and direction, the lava’s temperature-dependent viscosity also will play a role in determining when and if preferred pathways form within the flow interior. Thermo-viscous instabilities may set the size and number of preferred pathways. [11] We investigate observed and derived flow velocities, cooling rates, and dimensions from basalt flows in Hawaii and on Mt. Etna to provide physically reasonable flow parameter inputs to our model. Etnean and Hawaiian flows contain different amounts of two types of basalt lava flows, pahoehoe and a`a. These two types of basaltic flow are classified separately due to their different flow textures and morphologies. The morphological differences result from differences in how the flow is able to fracture or tear apart its surface during advance [Kilburn, 1993, 2004; Rowland and Walker, 1990] and from differences in the lava rheology [Cashman et al., 1999]. However, the inner cores of these two lava types seem to be of sufficiently similar composition [Cashman et al., 1999; Kilburn, 2004] and temperature to produce the needed viscosity contrast for preferred pathway initiation. Additionally, reported flow velocities and cooling rates are similar for the two types of flow. Because these flow parameters are the relevant quantities for our model, we consider both lava types when making reference to “basalt flows.” Channels are commonly found in both a`a and pahoehoe flows, and tubes are also found in both types of flow although they are more common in pahoehoe flows [Calvari and Pinkerton, 1998; Hon et al., 1994; Kauahikaua et al., 2003]. [12] Our model applies to simple flows, which are flows that cannot be divided into flow units with respect to emplacement or cooling dynamics [Walker, 1971]. Our model applies to only the portion of a flow that is experiencing steady laminar flow. Thus, we focus on the region of flow where the flow is generally of near-constant width [Kilburn and Lopes, 1991] and is likely to be thermally mixed [Crisp and Baloga, 1990; Kilburn and Lopes, 1991]. Although the dynamics that we consider are likely to also occur within near-source and frontal zones of flows, it is unclear how large an effect the temperature-dependent viscosity could exert on the flow relative to flow spreading or turbulence. Investigation of thermo-viscous instabilities within these more complicated flow regions is beyond the scope of this study. [13] In this paper, we use the term “sheet flow” to describe a coherent flow with an advancing undivided molten core that extends across flow width and whose thickness is smaller than its width or length. This is consistent with our simple flow and laminar flow constraints, and is analogous to the flow type examined in laboratory experiments of thermo-viscous instabilities [e.g., Fink and Griffiths, 1992; Wylie et al., 1999]. In the literature, the term “sheet flow” is sometimes used to describe surface morphology. For example, Hon et al. [1994] refer to sheet flows as broad, inflated flows that lack a hummocky surface. Soule et al. [2004] discuss “pahoehoe sheets” that are broad inflating regions within the proximal portion of Hawaiian flows that are either channel-free or between channels. Although

surface morphology may be linked to interior dynamics, our model of thermo-viscous instabilities describes only the behavior of flow within the molten core and does not explicitly predict surface morphology. Thus, our use of the term “sheet flow” refers only to the behavior of the flow within the interior flow core. We make no explicit assumptions about the cooling dynamics along the periphery of the flow or the structure or morphology of the flow’s crust.

2.

Effect of Temperature-Dependent Viscosity

[14] A flow model coupled with a temperature-dependent viscosity developed by Whitehead and Helfrich [1991] was used in the first investigations on the role of thermo-viscous instabilities in magma dynamics. Helfrich [1995] and Wylie and Lister [1995] adapted that model to examine basaltic magma flow through dikes and fissures. These studies more rigorously investigated the flow dynamics predicted with this model. In particular, steady state analysis of the model established the criteria for the development of preferred pathways. The key parameters were the sensitivity of viscosity to temperature and the initial flow rate. These studies demonstrated the plausibility of stable low-viscosity, fast-flowing pathways forming spontaneously within a vertical flow. Additionally, numerical simulation of a vertical flow driven by a specified change in pressure (simulating a fissure eruption) showed that a small local variation within the temperature profile at the influx boundary is sufficient to initiate the thermo-viscous instability. Furthermore, periodic variations within the influx temperature profile yielded a set of widely spaced pathways. This dynamic was shown to be sufficient to cause a uniform curtain of lava within a fissure eruption to evolve into isolated and stable eruption vents [Helfrich, 1995; Wylie and Lister, 1995]. 2.1. Physical Assumptions [15] The equations and assumptions used in this study are very similar to those considered by Helfrich [1995] and Wylie and Lister [1995]. Primary differences in model structure arise from different assumptions about flow orientation and heat loss through the bounding flow surfaces. The prior studies considered vertical flow through a fissure driven by buoyancy forces and a specified pressure difference between the magma chamber and free surface with little heat loss through the fissure walls. In our study, we focus on horizontal rather than vertical flow and thus neglect buoyancy forces. Additionally, we do not consider pressure when defining our influx condition. Instead, we specify an influx velocity and temperature profile. We also include heat loss through the flow’s upper and lower surfaces. This study is the first investigation of how lava’s temperature-dependent viscosity may influence flow emplacement. [16] As with most fluid dynamics and geological problems, we make a number of simplifying assumptions. These assumptions limit the applicability of our model to specific types of flows and the level of complexity we can consider, but they are needed to keep the equations and numerical algorithm tractable. Our assumptions include the following: [17] 1. Newtonian laminar flow. [18] 2. Uniform, shallow flow thickness relative to flow width and length.

1518

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Figure 2. A planform view of a flow with a central hot and low-viscosity preferred pathway. Flow emanates from a linear source (of the same width as the flow) at the left edge (red), with proscribed temperature and velocity profiles. The sides of the flow (black) are held at constant cold temperature. Heat loss is assumed to occur through the flow’s upper and lower surfaces (out of the page). The arrow lengths and directions illustrate, schematically, the flow velocity field. [19] 3. Flow occurs in two-dimensions (planform) and vertical variations in flow temperature, viscosity, and velocity are neglected. [20] 4. Flow is of uniform width and over a featureless, horizontal surface. [21] 5. Effusion rate is constant and from a linear source. [22] 6. Heat lost from the flow through its surfaces is integrated into a single cooling rate. [23] 7. The relationship between temperature and viscosity in the governing equations is a simple exponential function. [24] These assumptions are described in more detail in the following paragraphs. [25] Within the lava core, the mobile material between lower and upper crusts, we assume that flow can be approximated by steady, Newtonian laminar flow. This assumption is supported for pahoehoe flows by the smooth surfaces of some flows [Hon et al., 1994], by laboratory measurements of a small yield stress in these lavas [Shaw, 1969], and by studies which compare the results from this type of simple model with more complex lava flow models [e.g., Glaze and Baloga, 1998; Castruccio et al., 2013]. Additionally, the Reynolds number of basaltic lava flows is typically small enough to suggest laminar flow [Booth and Self, 1973; Wylie and Lister, 1995], and this is a common model assumption [e. g., Moore, 1987; Wylie and Lister, 1995]. The Reynolds number is even smaller for silicic flows, due to their higher viscosities and lower flow rates. Pahoehoe flows can transition into a`a flows that have a non-negligible yield stress as crystallization increases beyond a rheological threshold [Dragoni and Tallarico, 1994; Castruccio et al., 2013]. This does not negatively impact application of our model. If lava were to behave as a Bingham fluid, then flow focusing should be enhanced as cooled lava would be slowed even more than predicted by the Newtonian flow model. Thus, with Bingham flows, our model still yields a reasonable prediction for the onset of pathway formation. [26] We assume that the uniform flow depth is much smaller than the other two flow dimensions (also known as a Hele-Shaw cell geometry or a shallow-flow assumption). Such a geometry is generally found in terrestrial [Anderson et al., 1999; Byrnes and Crown, 2001; Hon et al., 1994; Kauahikaua et al., 1998] and planetary basaltic flows

[Zimbelman, 1998]. Under this assumption, flow can be modeled using a modified Darcy’s law [Bonn et al., 1995; Saffman and Taylor, 1958] where the pressure gradient is proportional to the fluid velocity. [27] On the basis of shallow-flow geometry, we simulate a two-dimensional system instead of a full three-dimensional system. Wylie and Lister [1995] compared results of a twodimensional numerical simulation with a full three-dimensional system and found that the two-dimensional model was sufficient for investigating steady state velocity solutions as long as suitable “average” parameters are considered. [28] For simplicity, we consider flow within an area of uniform width, which we acknowledge may not be realistic as a flow can be significantly wider than its source or vary in width. We also assume flow is over a featureless surface. For the low slopes on which basalt (especially pahoehoe) flows are typically found [Hon et al., 1994], the influence of slope can be indirectly accounted for within our featureless-terrain simulation by varying the influx velocity appropriately (e.g., using Jeffrey’s equation where velocity ~ sin (slope angle) [Jeffreys, 1925]). However, many basaltic regions are dominated by hummocky preeruption topography that is not accounted for in our model. [29] We assume a constant effusion rate from a linear source from which flow emanates in a uniform sheet. In reality, long- and short-term variations in effusion rate (and flow velocity) are common due to changes in supply, upstream blockage (or blockage release), or changes in vesicularity [Harris et al., 2007]. A change in flow velocity could disrupt the growth of an instability (or induce an instability in a previously stable flow), and thus, our results can be applied only to periods of relatively uniform effusion. [30] In addition to diffusion of heat within the flow and from the flow to the side walls, heat is lost through the flow’s upper and lower surfaces (i.e., into or out of the page in Figures 2 and 3). Rather than considering heat loss in different directions and through different mechanisms individually, we integrate all heat loss into a single cooling rate (δ(x,t)). This restricts application of our model to simple flows (as defined by Walker [1971]). In some simulations, we assume a temporally and spatially constant cooling rate. In others, we allow it to smoothly decrease in space and/or time (e.g., as might be the effect of a thickening surface crust as the flow matures). [31] The Arrhenius Model (μ ~ exp(k/T)) is commonly used to relate temperature and viscosity of lavas [Blatt et al., 2006]. In our model, we use a first-order approximation of the Arrhenius Model to describe the relationship between viscosity (μ) and the dimensionless temperature of the lava (T): μ(T) = μH exp(β(1  T)). The dimensionless temperature is defined as: T = (Tdim  TC)/(TH  TC) where Tdim is the temperature at a point (x,y) in the flow, TH is the hottest temperature of the lava at the source, and TC is the steady “coldest” temperature of the side walls. T ranges from 0, when Tdim = TC, to 1, when Tdim = TH. When T = 1, the flow will be at minimum viscosity, μH. When T = 0, we achieve maximum viscosity (μC = μH exp(β)). Thus, the value exp(β) corresponds to the largest viscosity contrast ratio that will exist within the system (μC/μH). For simplicity and consistency in estimating reasonable β values from laboratory and field viscosity measurements, we assume that this ratio (exp(β)) occurs over a 100°C temperature difference, which is

1519

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Figure 3. Numerical simulations illustrating differences between flow developed with small (β = 2) and large (β = 4) dependence of viscosity on temperature. Both simulations use δ = 0.1, uinflux = 1, and T = 1 at (x = 0, y = 0) with a decay to T = 0.9 along the edges of the influx boundary (x = 0, y → ± 10). (left) Temperature contour planform maps (see color bar) within a steady state flow. (right) Corresponding steady state temperature profiles (solid lines) along y = 0, plotted with exponential trend lines predicted by the maximum streamwise velocity (T = exp(x δ/u*); dashed). When β = 4, the steady state solution overpredicts temperatures close to the source as heat and lava are being focused into the low-viscosity preferred pathway (which increases in velocity to 2.3 m s1). Flow focusing continues over the characteristic length scale (x = u*/δ ~ 23 m in this simulation, vertical dashed line; see section 3.1). Beyond that distance, transverse flow velocities go to zero and the model matches with the predicted trend. All distances (x,y) are presented in meters. comparable to the difference found between basaltic flow within lava tubes and the lava tube walls [Kauahikaua et al., 1998]. We emphasize, however, that the important quantity in the generation of preferred pathways is the viscosity contrast ratio, not the specific associated temperatures or temperature difference. In most simulations, we assume a constant β, but we also examined the effect of a β that varied linearly with respect to the temperature. As we will discuss, this had a very small effect on flow dynamics. 2.2. Equations [32] On the basis of these assumptions, we develop the following physical equations: ∇· U ¼ 0 ðconservation of mass=incompressibilityÞ

(1)

∇ P ¼ 12μ U=d ðconservation of momentumÞ

(2)

2

μ ¼ μH expðβ ð1  TÞÞ ðtemperature-dependence of viscosityÞ (4)

where the velocity (U = (u,v)), pressure (P), viscosity (μ), and dimensionless temperature (T) are variables defined for each location (x,y) and evolving in time (t). The coefficients of thermal diffusivity and the cooling rate (δ) may be dependent on location (x,y) and time (t), but are otherwise generally independent of flow conditions. Tt is the time derivative of the dimensionless temperature. The flow thickness (d) is also specified in these equations, but we will show that it is not used in the final model equations. A summary of the typical model input parameters for basalt flows that contain channels and tubes is given in Table 1.

[33] To simplify the system, we convert our velocity vector field (U) into a stream function (ψ, a scalar field where U = ∇ x ψ) because this causes the system to automatically satisfy the incompressibility condition. The energy equation remains similar in structure, becoming

and the equation for the temperature-dependent viscosity is unchanged. After taking the curl of both sides of the momentum equation (Equation 2), we get a modified Laplace’s equation: ∇· ðμ∇ψÞ ¼ 0

(6)

[34] This equation form is simpler to work with and it also allows us to specify a more natural boundary condition for horizontal flow: we need only to specify the influx velocity rather than the pressure at the boundaries of the flow. Additionally, this equation form yields an intuitive idea of how the viscosity relates to velocity. Let us first consider the case where viscosity (μ) is constant; this yields ∇2ψ = 0 (Laplace’s equation) which, to first order, means that the stream function at any point is the average of the values at that point’s nearest neighbors. When viscosity is not a constant, the stream function at any point is instead a weighted average of its nearest-neighbors’ values, where the weighting is inversely related to the nearest-neighbor’s viscosity. This reflects the physical idea that mass is conserved, so if the fluid around a point has average nonzero velocity, then fluid must also be moving at that point. However, if one side of the point has a higher viscosity, then the flow will be retarded on that side, thus skewing the motion.

1520

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS Table 1. The First Set of Model Parameters We Used to Investigate Basaltic Flows That Form Preferred Pathways, for Both Analysis and Numerical Simulation of Model Equationsa

Variables

Parameters

Physical Meaning, Units

Specified Values

U = (u.v)

Velocity, m/s

μ

Viscosity, Pa s

T = (Tdim  TC) / (TH  TC)

Dimensionless temperature

Uinflux = (uinflux, 0) is the influx velocity (constant along influx boundary) μ(1) is the viscosity at highest temperature (i.e., lowest viscosity) TH and TC are the highest (i.e., at source) and lowest (i.e., along wall) temperatures, and thus will yield the range of viscosity values the flow can have; we assume TH - TC = 100 °C

w

flow half-width, m 2 thermal diffusivity, m /s

δ = δmaxf(t)g(x)

β

is heat conduction through the lava and between the lava and walls δmax is the initial/maximal cooling rate, before decay due to increase of surface crust thickness with distance from source f(t) is zero or a monotonically decreasing function g(x) = exp(c x) is change in the cooling rate due to crustal increase in distance from the source

fractional cooling rate (through surface crust), 1 s = absolute cooling rate 1 °C s normalized by the maximum temperature difference within the flow (100 °C) exponent in the equation relating temperature and viscosity, dimensionless

Typical Values 1ms

1

100 Pa s TH = 1150°C (eruption temperatures in Hawaii), TC = 1050 °C

100 m 6 2 1 = 10 m s b 2

1 1

δmax = 10  10 s (as radiative cooling rates from open channels -1 are 1 – 10 °C s ) c See section 3.3 2

c = 10

m

1 c

0–6 once the lava temperature is below the liquidus (~1150°C) (0 would mean no temperature dependence)

a

The ranges of values that were considered in the full course of this study are given in Table 2. Common order-of-magnitude estimate [Wylie and Lister, 1995]. From measured cooling rates of °C/m, after normalizing by the allowed temperature difference (100°C). Example values of δmax and c are derived based on measurements of cooling rate of 8.3°C/m at the vent mouth and 5.8°C/m 40 m downstream and an assumed exponential decay of the cooling rate with distance [Pinkerton et al., 2002]. All other values are explained in the text. b c

[35] Finally, we assume the cooling rate has the form δ = δmax f(t) g(x), where f and g are ≤ 1 and constant or decreasing functions. As we will discuss, we generally focus on just the influence of δmax. Physically, however, the full equation form could correspond to the cooling rate decreasing as insulating crustal thickness increases with distance from the source and with time. 2.3. The Natural Length Scale [36] As the primary rate of temperature change within our system is due to the cooling of the lava, the most natural timescale of the model is 1/δ. The natural length scale is the ratio between the streamwise velocity and the cooling rate (u/δ), which physically reflects the thermal balance that the system achieves between the heat loss through its surfaces and the influx of hotter fluid. This length scale can also be thought of as the distance over which the full viscosity change (μH → μC) would occur within a steady flow if δ was a constant and absolute cooling rate rather than a fractional cooling rate (i.e., if the heat loss term was δ, and not δT). [37] A more intuitive explanation of the length scale can be found by considering the system where viscosity is constant and independent of temperature (i.e., β = 0 so μ = μH exp (0) = μH everywhere). The influx velocity (U = (uinflux,0)) and cooling rate (δ) are also assumed to be constant and uniform. Under these conditions, the equations described in section 2.2 have a trivial steady state solution (i.e., Tt = 0)

in which the dimensionless temperature decays exponentially with distance from the source T ¼ T0 expðδ=uinflux xÞ

(7)

where T0 = T0(y) is the dimensionless temperature at the influx boundary. A similar steady state solution is found when β ≠ 0, but in that case transverse flow will occur near the source due to viscosity contrasts. As a result, the steady state velocity of the system will differ between streamwise profiles. The average streamwise velocity of the flow will be higher than uinflux within a hot, low-viscosity preferred pathway and lower within the cooler flow (this will be discussed further in section 3.1). [38] Thus, our model predicts that the flow will always trend toward T = T0 exp(x/L), where L is the system’s “natural” (or characteristic) length scale (u/δ). These results provide us with a simple way to compare the dynamics observed in different lava fields or as lava flow evolves in time. If we wish to compare two flows that have different flow parameters, we can do this after seeing how their characteristic length scale differs. For example, if Flow 1 has a flow velocity that is higher than Flow 2 by a factor of α, without changing δ, then their dynamics should appear the same after stretching the spatial dimensions of Flow 2 by that factor (e.g., doubling the velocity means that the high temperatures within a pathway will extend twice as far downstream). However, if δ in Flow 1 is also higher than δ in Flow 2 by the same factor α (so that L = u/δ remains constant), then the expected distances over which pathways will form or

1521

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Figure 4. Maximum and minimum streamwise velocity (u) values within the steady state flow (neglecting flow right next to the wall) as a function of viscosity’s dependence on temperature (β). The different lines are for different amounts of preferred pathways that were formed within the flow (solid is with one preferred pathway; dashed is with two; dotted is with three). As β increases, the flow experiences more flow focusing and the steady state streamwise velocity increases within the preferred pathway. Pathways form (i.e., velocities begin diverging strongly due to the development of transverse flow) when β > ~3. extend remain unchanged. (The timescale over which the preferred pathway forms in Flow 2, however, increases by 1/α as the model’s “natural” timescale is 1/δ.) 2.4. The Inf luence of β [39] Helfrich [1995] and Wylie and Lister [1995] both explored the influence of the exponent (β) in the temperature relating temperature and viscosity. Those studies showed that for small β, viscosity changes only slightly with temperature and so little flow focusing develops (i.e., the lava core will retain its “sheet flow” structure). A large β-value yields significant flow focusing into regions of higher temperature (Figure 3) as viscosity changes greatly due to a small difference in temperature, and thus induces large changes in flow velocity (Figure 4). This allows low-viscosity preferred pathways to form and stabilize. As shown in Helfrich [1995] (and reproduced briefly here for completeness), the threshold between these two states can be found by checking for two flow velocities that yield the same flow pressure (for a system with a fixed β and cooling rate (δ)). If two velocities can be found, then it should be possible for neighboring flow regions moving at the different velocities to be stable. As the pressure within each region will be the same, there is no pressure gradient in the transverse direction and thus no deformation of the flow. [40] Phrasing this mathematically: we seek two fixed streamwise velocities u1 and u2 (such that U(x, y1) = (u1,0) and U(x,y2) = (u2,0)) at which the flow pressure at some distance (x) from the source is the same (i.e., (P(x,y1) = P(x,y2)). The total pressure generated by the flow at a distance x from the source (i.e., the integration of equation (2) along a streamwise profile (y = y1) from 0 to distance x) is given by x

x

Pðx; y1 Þ ¼ ∫0 ∇PðsÞds∼∫0 u1  expðβð1  expðsδ=u1 ÞÞds

(8)

where we have used the steady state temperature profile derived in the previous subsection (Equation 7, with uinflux = u1

and T0 = 1), s is a dummy integration variable, and we have dropped all constant coefficients. For a specific β, we can now look for two velocities u1 and u2 where P(x,y1) = P(x,y2). When β < 3, we cannot find two velocity values that will satisfy this equation for any x. For β > ~3 (equivalent to the viscosity contrast ratio μC/μH = exp(β) > ~20), two different velocities can be found that yield the same pressure value for any distance x. Thus, for β values above this threshold value (~3), it is possible for neighboring regions of different velocities to maintain stable flow once a preferred pathway has formed. Thus, the restriction that β > ~3 corresponds to a threshold viscosity contrast ratio within the flow for the stabilization of a preferred pathway. [41] However, this is the case only when both δ and u are not (relatively) extreme values. For all β, when the flow velocity becomes very high and the cooling rate drops to zero (Case I), or when the cooling rate becomes very high and the flow velocity drops to zero (Case II), then it is not possible to find a second velocity that will yield the same pressures. In either case, this is because it is difficult to generate the requisite temperature differences to create a viscosity contrast ratio and initiate pathway formation. In Case I, the influx of hot fluid raises the temperature of the flow at a rate faster than heat can be lost, so the flow will be kept at a near-uniform hot temperature. In Case II, the flow is cooled at a relatively fast rate, so it will be kept at a near-uniform cold temperature.

3.

Analytic and Numerical Investigations

[42] We use our model to predict interior flow dynamics of basaltic flows using two methods. Steady state analysis of model equations allows us to examine conditions under which preferred pathways can form and stabilize due to lava’s thermo-viscous instability, and to identify expected trends and scaling with parameter values. Numerical simulation allows us to examine the development of pathways and the evolving velocities and other system variables. Details of the numerical scheme are given in Appendix A and typical model parameter values are given in Table 1. [43] Whenever possible, model results are compared with observations of terrestrial lava flows (Table 2). However, only a small number of reported measurements are of a type or resolution that we can use. Measurements need to relate to early flow dynamics; thus, only well-monitored, consistently active flows were examined. Also, correlated estimates of velocity, cooling rate, and channel or tube widths are needed from at least a few locations within a single flow. 3.1. Model Study 1: Constant Cooling Rate [44] We first consider the simple case, where the cooling rate is constant in space and time (i.e., δ = δmax). In general, simulations were run with influx temperature profiles having a specified number of spatially periodic temperature highs with variable magnitudes. We assume a maximum difference of 10°C within the influx temperature profile (or 0.1 in our dimensionless temperature), comparable to measured temperature variations in Hawaiian lava pond samples reported by Helz et al. [2003], in Mauna Loa vent flow reported by Lipman and Banks [1987], and in the fluid lava portion (> 4 m depth in 1963) of the Alae lava lake reported by Peck, 1978). However, to be conservative in our modeling, we include local differences of only a few degrees.

1522

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS Table 2. Typical (Order of Magnitude) Values for Different Flow Rheologies and Morphologiesa Silicic Flows Influx velocity (uinflux) Cooling rate (δ) Flow width (w) Exponent in the equation between temperature and viscosity (β) Typical flow length Predicted distance over which pathways develop (L) Predicted timescale over which pathways develop

5

Basaltic Flows with Tubes/Channels

1

1

10 m s 4 1 - 10 °C s

0.01–1 m s 3 1 10 °C s (cooling rate of a well-insulated tube) 1 –10°C s (open channel) 10 m–1 km 3–5+ (once temperature is sufficiently below liquidus) 100 m–10 km 1–100 m, consistent with observations

9

10

300 m 4–9

Basaltic Sheet Flows 1

0.1–5 m s To keep pathways from forming: 3 1 1 10 –10 °C s is sufficient 10 m–10 km 1–5+

1 km 100 m–1 km, up to 100 km 1 km+, consistent with To keep pathways from forming observations when β > 3: > 100 m–1+ km 7 9 3 10–10 s To keep pathways from forming 10 –10 s 4 2 3 (months to years), consistent (seconds to minutes), up to 10 s (few hours), when β > 3: 10 –10 + s (minutes+) with observations consistent with observations -4 1 1–3 m s , consistent with observations N/A Predicted flow velocity within 10 m/s, consistent with observations a stable pathway (u*) Predicted pathway width (f) 30–60 m, smaller than observed 100 m 1–100 m, consistent with observations N/A a Many basaltic flows and some silicic flows contain channels and tubes. Basaltic sheet flows do not. The upper five parameters are input observation values, and the lower four variables are predicted output values, calculated via our flow dynamics model. For silicic and basaltic flows, predicted output variables are generally consistent with observations (except for silicic pathway width), as discussed in the text. Italicized text indicates a prediction where observations were not found in the literature for direct comparison: for basaltic sheet flows, we estimate the cooling rate needed for pathways to not form. All other values are explained in the text.

[45] As discussed in section 2.3, when the cooling rate is constant, the analytical steady state solution along a streamwise profile (y is constant) consists of a constant streamwise velocity (U = (u*,0)) and a dimensionless temperature that decays exponentially with distance from the source (T(x,y) = T0(y) exp( x δmax/u*)). There is of course some small variation between the numerical simulation and this steady state, as the analytical steady state solution assumes no transverse flow. In the numerical simulation, transverse flow occurs near the influx boundary and the streamwise velocity will increase above uinflux within a preferred pathway. This causes the flow’s dimensionless temperature to initially decay faster than predicted by the analytical solution. Beyond the distance x ~ u*/δmax (where we use the maximum velocity within the streamwise profile of the flow for u*), the system is wellapproximated by the analytical steady state solution (Figure 3). [46] The exact value of the streamwise velocity at which the low-viscosity pathway becomes stable (u*) depends on the viscosity contrast ratio achieved by the flow, as shown in Figure 4. This dependence occurs because larger β values will cause larger viscosity contrast ratios for a small temperature difference (the fully developed viscosity contrast ratio is exp (β)) which, in turn, drives higher levels of flow focusing toward regions of higher temperature. Within a pathway, we find that the maximum flow velocity increases nonlinearly with β and varies weakly with respect to the number of thermal peaks included along the temperature profile at the source (Figure 4). [47] We consider the system consisting of a single preferred pathway of half-width f and velocity umax within a flow of half-width w and velocity umin (f and w shown in Figure 2). Using our numerically derived relationship between β and velocities that our simulation generates inside and outside of the preferred pathway (umax and umin, respectively) (Figure 4), we can predict the pathway width relative to the flow width: f umax þ ðw  f Þ umin ¼ wuinflux

(9)

as mass must be conserved and our assumed geometry requires that no breakouts occur within the flow and that

the flow cannot thicken. We are now able to solve for f/w as a function of β (Figure 5): f =w ¼ ðuinflux  umin Þ=ðumax  umin Þ

(10)

[48] As the velocities vary little with the number of pathways (Figure 4), the actual pathway widths will scale with the density of pathways. We note that this does not impose a natural length scale for individual pathway widths, but instead yields an integrated (across all pathways) length scale and a minimum relative spacing between pathways. 3.2. Comparison of Model Study 1 With Observations [49] To compare our model results with observations, we sought out velocity and cooling rate values for typical basalt flows that did not yet appear to have been significantly affected by topographic barriers or by cooling. Generally, we use order-of-magnitude values reported within the early stages of flow and within the middle or channelized portions of flows. We also use only values that were reported for several flows or for “typical” flows. Values reported for flows that were noted as being unusual or greatly varied in morphology, effusion rate, flow velocity, viscosity, etc. were avoided. Values reported for very narrow flows (< 10 m) were also avoided. Finally, when possible, we emphasize direct measurements rather than model-dependent derived values. [50] Basaltic flows that have not yet evolved to form tubes or channels typically have flow velocities on the order of 0.01–1 m s1 [e.g., Guest et al., 1987; Kauahikaua et al., 1998; 2003; Lipman and Banks, 1987; Swanson, 1973]. (Higher velocities have been reported [e.g., Baloga et al., 1995], but such flows are possibly turbulent and thus are outside the scope of our model.) Cooling rate observations within these basaltic flows were difficult to find in the literature, but a loose upper bound can be estimated on the basis of near-vent open-channel radiative surface cooling rates, which are > 1°C m1 [Pinkerton et al., 2002]. A loose lower bound for cooling rates can be estimated from lava

1523

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Lipman and Banks, 1987; Soule et al., 2004] and may sometimes approach 10 m s1 [Greeley, 1971; Kauahikaua et al., 2003]. This suggests that the moving lava has generally been focused into a pathway about 0.1–0.2 times the flow width. Thus, it seems reasonable that this rough estimate for f/w may generally hold over basalt flows. Figure 5 shows that this ratio corresponds to expected β values for basalt flows that develop channels and tubes (β ~ 3–5, as discussed in the next subsection).

Figure 5. Predicted (lower bound) fractional pathway width (f/w) as a function of viscosity’s dependence on temperature (β). The vertical dashed line illustrates a threshold value β ~ 3 (see text for discussion). When β < 2, the approximation used to calculate f breaks down, preventing numerical estimation of f/w. When β = 0, there is no temperature-dependent flow focusing and the preferred pathway is equivalent to the entire flow, so f/w = 1. temperature changes measured within well-insulated tube flows, which cool ~103 °C m1 [Cashman et al., 1999; Harris et al., 2005; Helz et al., 2003; Kauahikaua et al., 1998; Peterson and Swanson, 1974; Riker et al., 2009]. Cooling rates estimated from glass geothermometry and other analyses of the interior of relict channelized flows generally yields values of ~ 102 °C s1 [Cashman et al., 1999]. [51] Using these broad estimates of flow velocities and cooling rates, we find that u/δ ranges from a few to tens of meters and that 1/δ ranges from minutes to a few hours (Table 2). These space and timescales are consistent with observations that channels can begin forming very rapidly after a flow has begun and very close to the source [e.g., Calvari and Pinkerton, 1998; Peterson and Swanson, 1974]. Thus, preferential pathway formation is likely within basaltic flow because typical flow lengths (100 m–10 km) [Hon et al., 1994; Self et al., 1998; Swanson, 1973] and durations (hours to weeks) [Crisp and Baloga, 1990] are much greater than the predicted formation scales. [52] If we assume the existence of only one or a few pathways within a flow, we can compare computed f/w ratios with those observed within basaltic lava flows in Hawaii and at Mt. Etna, Italy. Individual flows are often a few hundred meters wide [Hon et al., 1994; Lipman and Banks, 1987], and channels and tubes can be as large as 10–30 m wide [Greeley, 1971; Harris et al., 1997; Hon et al., 1994; Kauahikaua et al., 1998; Lipman and Banks, 1987]. These dimensions yield f/w ratios of around 0.1–0.2. Some studies report flows containing smaller channels (3–5 m wide), but those studies also note that either the channels are “numerous” [Guest et al., 1987; Peterson and Swanson, 1974] or that the flow base is narrow (i.e., few tens of meters wide) [Mazzarini et al., 2005; Swanson, 1973]. Additionally, estimates of maximum flow velocities within channels and tubes tend generally to be a few times to an order of magnitude higher than the reported flow velocities near the source. Velocities measured within channels and tubes are typically 1–3 m s1 [Cashman et al., 1999; Gregg and Keszthelyi, 2004; Hon et al., 1994; Kauahikaua et al., 1998; 2003;

3.3. Model Study 2: A Variable Exponent (β) in the Temperature-Viscosity Equation [53] As β = ln(μC / μH), the value β can be thought of as the strength of the viscosity’s temperature dependence. This parameter can be found by calculating the “slope” of a plot of temperature versus the natural logarithm of viscosity (Figure 1), over the temperature range of interest (TC to TH). Generally, laboratory experiments involving low viscosity lavas yield small β values; for example, correlated temperature and viscosity laboratory measurements shown in Figure 1 yield β = 1.4 for basalt. Similarly low β values (1–2) are found in other laboratory rheological studies of basalt [e.g., Shaw, 1969]. As these β values are below the threshold we found in section 2.4 for the formation and stabilization of preferred pathways (β ~3), we would not expect preferred pathways to form in basalt. However, β is not independent of the temperature range considered: the strength of the temperature dependence of viscosity increases (i.e., the plotted relationship steepens) as the temperature drops below the liquidus (~1200°C for basalt) and crystallization begins. As the temperature of basalt lava passes below 1150°C, increased crystallinity greatly enhances the viscosity of the lava [Hoover et al., 2001; Pinkerton and Sparks, 1978; Shaw, 1969; Sato, 2005]. For example, laboratory measurements of basalt viscosity [Shaw, 1969, Table 2; Sato 2005] yield β = 1.2 when TH = 1300°C, versus 3.6 for TH = 1200°C. Observations of Hawaiian lava flows yield temperatures of eruptions [Ault et al., 1961; Crisp and Baloga, 1990; Cashman et al., 1999; Fink and Zimbelman, 1990; Pinkerton et al., 2002] and flow in lava tubes and channels [Helz et al., 2003; Pinkerton et al., 2002] of ~1100–1150°C, suggesting that β values for these eruptions are likely at the higher end, and above the threshold β-value of 3 (Table 1). Helfrich [1995] extrapolated from Shaw’s results to estimate β ~ 5 for basalts having temperatures near the solidus ~ 1000°C. Furthermore, there is also some morphological evidence from natural lava flows that very high β values are plausible for basalt flows. For example, Crisp and Baloga [1990] estimated viscosity increases within a 1984 Mauna Loa flow of 4–5 orders of magnitude due to an ~30°C temperature decrease (from ~1130°C) and 20%–30% crystal content increase. Extrapolating this change in viscosity out over a 100°C temperature difference yields β ~ 9! Given this evidence that high β-values are likely for natural basalt flows that have cooled sufficiently below the liquidus, we expect most basalt flows to develop preferred pathways as the lava cools. [54] To evaluate the effect a varying β-value might have on the final pathway configuration, we numerically simulate flows in which β is a linear function of dimensionless temperature and vary its value from 1 to 4. We find that the end-state flow dynamics using this range of β values are very similar to

1524

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

Figure 6. Plot showing how typical ranges in cooling rate (δ) and flow velocity (u) for a lava type can be used to predict the characteristic length scale over which preferred pathways could form within a flow. Boxes correspond to the value ranges considered for the lava types discussed in this study (see Table 2). Dashed diagonals are typical basalt flow lengths of 100 and 1000 m. Our model predicts formation of preferred pathways if the flow velocity and cooling rate yield a characteristic scale of less than the flow length, or if the parameters fall below the relevant diagonal line. Typical parameter values for different types of flows are represented as boxes or lines. The parameter box given for silicic flows is unconstrained vertically because only one flow velocity was found in the literature. those we find when using just the maximum β value of 4. Therefore, β is left as a constant value for most simulation runs. However, if erupted basalt lava is above the liquidus temperature, then this perhaps explains why some basalt flows do not immediately develop preferred pathways (as will be discussed in section 4.2). 3.4. Model Study 3: An Evolving Cooling Rate (δ) [55] Lastly, we consider a system where the cooling rate is allowed to decrease smoothly with distance and with time. As the cooling rate is the parameter least constrained by observations, and because it may vary greatly between types of eruptions and environments, we investigate a range of values and a few functional forms. Generally, we use a constant cooling rate (δ = δmax), but we investigate a cooling rate that decays with distance from the source and with time: δ = δmax g(x) f(t). [56] Some measurements of surface lava temperatures at a variety of distances from the vent suggest that an exponential decrease of the cooling rate with distance is a reasonable approximation [Pinkerton et al., 2002], so we assume g(x) = exp(c x), where c ≥ 0. A lack of reported measurements of cooling rates along channels or tubes as a function of time made it difficult to approximate a time-dependent function of the cooling rate based on observations. Therefore, we assume the cooling rate will have an inverse proportionality to t: f(t) = 1/(1 + n ∙ t), where n ≥ 0. We chose this functional form as the crust thickens as t1/2 [Hon et al., 1994] and conduction of heat through a crust yields Tt ~ Tzz (second derivative of T with respect to crustal thickness) ~ 1/thickness2. The additional “1 +” in the denominator is included to keep the function from blowing up at t = 0.

[57] The use of nonconstant cooling rates did not significantly alter the simulation results. The dynamics of flow remained unchanged and the preferred pathways initialize over the same length and timescales as for δ = δmax. However, as δ changes, the preferred pathways extend farther downstream within the flow, with higher temperatures preserved for longer distances than within the analytical steady state. This behavior is due to the fact that the system’s streamwise characteristic length scale is increasing (L ~ δ -1 ~ exp(c x)(1 + n ∙ t)). [58] This temporal evolution in the characteristic length scale could present a possible additional stabilizing process for preferred pathways that will enhance the formation of mature channels and tubes, as hotter lava is able to extend farther downstream within a preferred pathway. To evaluate this effect, we consider the length (~ c1) and time (~ n1) scales over which the decrease in cooling rate and related pathway extension becomes an important influence. If these scales are generally longer than the scales of pathway initialization, which is the primary focus of this study, or of typical flow durations, then they can be neglected at this time. Estimates of crustal growth over lava lakes [Hon et al., 1994] suggest that n should be very small (perhaps as small as 108/s), which would mean f(t) ~ 1 until a few months or years have passed. Unless n is much larger, the time-dependent component of the cooling rate function can be ignored. Similarly, reasonable c values imply a distance of approximately 70 m before the crustal growth as a function of distance begins to influence the dynamics. As this length is greater than the distance over which pathways are generally expected to initialize in basalt, we can also ignore this component. Thus, it appears sufficient to approximate the cooling rate as a constant (δ = δmax) for basalt flows.

4.

Discussion

[59] Our model predicts that small temperature differences within a hot lava core can influence interior flow dynamics. This effect on interior flow dynamics due to small-scale temperature variations occurs because of lava’s natural thermo-viscous instability. That instability can result in the formation and stabilization of preferred pathways within broad flows. Our model predicts that certain conditions will be favorable for preferred pathway formation. These conditions include (but are not restricted to) the following: [60] 1. A flow thick enough to contain a liquid core, but sufficiently thin for the shallow-flow-geometry approximation to hold. [61] 2. The flow moves over sufficiently featureless ground, so that no large-scale topography will greatly deflect the flow or alter its thickness. The flow moves over a sufficiently constant and shallow slope to avoid slopeinduced changes in pressure [Guest et al., 2012] or thickness [Gregg and Fink, 2000]. [62] 3. Effusion at the source is sufficiently constant to avoid variations in flow conditions that will destabilize preferred pathway development. [63] 4. The viscosity contrast that develops as the flow cools is sufficiently large for focused flow to initialize. For a basalt flow, a sufficient viscosity contrast will arise naturally from typical temperature and crystal-content variations as the flow cools.

1525

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

[64] 5. The flow extends a sufficient length for preferred pathways to develop. The length (u/δ) and time (1/δ) scales for preferred pathway development both depend on the lava’s cooling rate, which depends primarily on the thickness and fracturing of a surface crust. [65] As discussed in section 3, we know that these flow conditions are often met within basalt flows. In particular, point 5 is shown graphically in Figure 6, where a parameter box identifies the flow velocity (u) and cooling rate (δ) values reported for typical basalt flows. For these values, the ratio of u/δ generally falls below the basalt flow length scale of 100 or 1000 m. Thus, our model predicts that many basalt flows are likely to form preferred pathways, and it presents a plausible explanation for why channels and tubes are common features of basaltic flows. [66] In the following discussion, we examine why other flow models do not produce a consistent correlation between flow characteristics and the presence of preferred pathways. We then use our model to investigate lava flows lacking tubes and determine the conditions that will inhibit preferred pathway formation. We also discuss model limitations and possible future extensions. 4.1. Comparison With Another Basalt Flow Model [67] As discussed in section 3, the length scale over which preferred pathways are likely to form depends on both the cooling rate and the flow velocity. This dependence is consistent with an idea that has been investigated by many volcanologists: that competition between development of a lava flow’s crust and molten core is a primary influence on lava flow dynamics. In particular, the “Kilburn model” focused on the competition of deformation forces induced by the crust and core of the flow during emplacement [Kilburn and Lopes, 1988, 1991; Kilburn, 1993, Kilburn, 2004]. In the Kilburn model, stress imparted by the advancing core tries to move the flow front forward whereas crustal strength resists this motion. Therefore, the timescale ratio between flow deformation and crust formation affects flow dynamics. When advection of the flow proceeds faster than crust formation (core-dominated flow), the flow will shear and tear the crust faster than the crust can heal. As a result, flow-front advance will be steady and the flow will lengthen and widen. Conversely, when the crust heals faster than the flow can spread (crust-dominated), the flow motion will be retarded and it will advance in spatial- and temporal-pulses. Core-dominated flow is commonly associated with a`a flows, whereas crust-dominated flow is commonly associated with pahoehoe flows [Kilburn, 1993, Kilburn, 2004]. [68] Although categorization of flow based on a timescale ratio of flow advection to crust formation yields a correlated categorization of flow morphology, the Kilburn model does not differentiate between flows that should or should not contain preferred pathways. We argue that this is because it considers flow velocity and cooling rate as individual flow parameters, rather than considering a coupling between them—as we do in our model. For example, we can explore what might happen to a flow if it experiences an increase in flow velocity and becomes more core-dominated. It is reasonable to assume that an increase in velocity would also increase crustal disruption, which would increase exposure of the molten core and radiative cooling. Thus, we expect both flow velocity (u) and cooling rate (δ) to increase relative

to what these flow parameters would have remained if the flow had not undergone an increase in flow velocity [Cashman et al., 1999]. This leaves it unclear how the length scale of pathway formation (u/δ) will change and we cannot directly predict pathway formation. Similarly, if flow velocity (u) were to decrease and cause the flow to become crustdominated, the crust would be more stable and thus the cooling rate (δ) of the core would also decrease. Thus, we do not expect to find a general correlation between overall flow dimensions, flow advance rates, or flow type (core- or crust-dominated; a`a or pahoehoe) with pathway formation or dimensions. [69] Other studies have also considered competition between a lava flow’s crust and molten core. These studies have proposed models that consider the competition in the thermal budget of long basaltic lava flows to relate the flow’s cooling and effusion rates to flow length [Keszthelyi and Self, 1998; Castruccio et al., 2013] or morphology [Fink and Griffiths, 1990, 1992; Gregg and Fink, 2000]. Most of these models are similar to the Kilburn model in that their predictions about correlations between flow parameters and resultant morphology rely on an estimate of either flow velocity or the cooling rate, but not the ratio of these two parameters. Thus, as with the Kilburn model, we do not expect to find a general correlation between the development (or lack) of preferred pathways and predictions from these other models that consider individual flow parameters. To our knowledge, the occurrence of lava channels or tubes has not been shown to consistently correlate with lava composition, rheology, flow length, effusion rate, or any other individual flow parameters. 4.2. Inhibition of Preferred Pathway Development [70] Although channels and tubes are common in basaltic flows, there are many large basaltic flows that lack a hummocky surface, apparent channels or tubes, or tube-related inflation landforms (such as tumuli [e.g., Duraiswami et al., 2004; Hon et al., 1994; Walker, 1991]). The formation of broad sheet flows implies that they were emplaced with a uniformly distributed liquid core, and observations of active flows support this hypothesis [e.g., Hon et al., 1994; Kauahikaua et al., 1998]. Generally, this type of flow is assumed to occur over near-flat topography; sheet flows are typically emplaced on terrain with slopes < 3 [Hon et al., 1994; Self et al., 1998] and topographic relief < 1 m [Byrnes and Crown, 2001; Swanson, 1973]. This has led to hypotheses suggesting that flow over a very shallow slope and subdued terrain will inhibit the formation of preferred pathways. However, the mechanism by which featureless topography may inhibit pathway formation is not clear. Some studies state that the lack of local relief and spreading of the flow promote low flow velocities that allow sheet flows to spread uniformly [Byrnes and Crown, 2001; Hon et al., 1994]. Other studies state that higher velocities are possible due to lack of flow disruption by topography, and these higher velocities cause rapid emplacement of the lava that stabilizes sheet-like spreading [Crown et al., 1999; Duraiswami et al., 2004; Self et al., 1998]. [71] We try a new approach to evaluate why sheet flows form and sometimes remain stable by considering eruption and rheological factors that could inhibit preferred pathway formation. One possibility for maintaining sheet flow is that the required viscosity contrasts do not occur (i.e., β is too small).

1526

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

For basalt flows, it is possible for a flow to be too hot to contain the required viscosity contrasts. For example, Riker et al. [2009] recount an observation of one flow in the 1859 Mauna Loa eruption emitted at a white heat, implying a temperature well above 1150°C [Williams and McBirney, 1979]. That flow did not separate into a network of lava streams until it had flowed for several kilometers and cooled to blood-red, suggesting a temperature closer to 1000°C [Williams and McBirney, 1979]. If a sheet flow stagnated before cooling sufficiently, then perhaps it could solidify without producing preferred pathways. For example, a sheet flow in the 1969 Mauna Ulu eruption extended only a few hundred meters before the feeding lava column withdrew and the flow stagnated [Swanson, 1973]. [72] Sheet flow stability may also be maintained if the length or timescale for the formation and stabilization of preferred pathways is longer than the flow length. We first consider length scale. Because typical basalt flows are a few hundred meters to a few kilometers in length, which is much longer than the distance over which we generally predict that preferred pathways will form in basalt (~10 m), we seek flow conditions that will increase the characteristic length scale (u/δ) by a few orders of magnitude and thus explain sheet flow observations. [73] Flow moving over a steeper slope may have a velocity higher than is common for basalt flows. For small slope increases, the velocity of a Newtonian fluid increases roughly linearly with the slope angle (i.e., the Jefferys equation states that velocity ~ sin(slope angle) ~ slope angle for small angles [Jeffreys, 1925]) under a constant effusion rate. This will increase the flow velocity, but by only a very small factor as generally we are interested only in low slopes. Steep slopes will increase the velocity by a greater amount, but in a sublinear way as the small-angle approximation will no longer hold and flows may thicken when flowing over higher slopes. Furthermore, gravitational shearing yields a readily fractured surface crust [Gregg and Fink, 2000; Guest et al., 2012], so the cooling rate will also likely increase. As both velocity and cooling rate may increase in flow on a high slope, relative to a typical flow over a horizontal surface, it is unclear what effect higher slope will have on the characteristic length scale of the flow, and thus on the formation of preferred pathways. Thus, we do not expect to find a simple correlation between slope and distance over which preferred pathways may form within a flow. This negative result is supported by observations: some portions of basalt flows on Mt. Etna flow contain channelized a`a despite very steep slopes (~27° [Guest et al., 2012]) and analog wax experiments have also form channelized flow structures over high slopes (up to ~30°) [Gregg and Fink, 2000]. [74] Alternatively, increasing the effusion rate is another way to increase the flow velocity, if the increase in flow volume is not accompanied by an increase in flow thickness or width. If the cooling rate is unchanged, then the flow can have a high u/δ-ratio and preferred pathway development will be slowed. Or, if the flow is moving fast enough, a sufficient viscosity contrast can be prevented from developing (Case I in section 2.4). Effusion rates at least an order of magnitude higher than maximum typical basalt flow rates have been proposed for terrestrial large igneous provinces (104 m3 s1 [Keszthelyi et al., 2006; Self et al., 1998] versus 103 m3 s1 [e.g., Calvari and Pinkerton, 1998; Kauahikaua

et al., 2003]) and peak effusion rates up to 1000 times larger possibly occurred in some of Earth’s largest basaltic fields [Keszthelyi et al., 2006; Swanson et al., 1975]. These rates, however, are very uncommon, especially within the types of flows that we are investigating within this study. Furthermore, high effusion rates are not likely to be sustained over long time periods. After the flow slows, that portion of the flow will then be likely overprinted by later tube-fed flows. Thus, like an increase in slope, it seems implausible that an increase in effusion rate will generate a sufficient increase in length scale to inhibit the formation of preferred pathways in most flows. [75] We now consider changes to the cooling rate, because decreasing it will increase the characteristic length scale. The cooling rate of the lava can be decreased if insulation of the flow interior is increased, such as through formation of a thicker or more uniform and unbroken crust. Formation of a thick, more uniform crust can decrease the cooling rate of a flow without affecting the flow velocity. The cooling rate of lava flow declines rapidly as the flow [Keszthelyi et al., 2006] or the surface crust [Crisp and Baloga, 1990; Hon et al., 1994] thickens, because more heat transfer occurs through conduction relative to radiative cooling or convection. Cooling rates at least 10–100 times less than the rates measured within fresh, uninsulated basalt flows are needed for preferred pathways to form at distances longer than the typical sheet flow length (Table 2; Figure 6). Observations of sheet flows show that a millimeter-thick flexible skin will cover the flow surface within seconds [Hon et al., 1994] and thicken into a centimeter-thick crust within minutes [Hon et al., 1994; Keszthelyi and Denlinger, 1996; Swanson, 1973]. If this crust is uniform and unbroken, it can quickly reduce the cooling rate [Anderson and Fink, 1992; Crisp and Baloga, 1990]. Once a stable crust has formed, cooling rates become very small. An extreme example of a wellinsulated flow is of course that found within a tube, which can cool at ~0.001°C/m [Cashman et al., 1999; Harris et al., 2005; Kauahikaua et al., 1998; Peterson and Swanson, 1974; Riker et al., 2009]. (Unfortunately, we found no direct reports of cooling rates for basalt sheet flows.) Furthermore, as the cooling rate decreases, the timescale for pathway development increases. If the cooling rate drops sufficiently, the development of a viscosity contrast can also be prevented (Case II in section 2.4). [76] However, if the crust is not uniform or is fractured, then the cooling rate will remain high in proportion to the exposed area of hot core [Crisp and Baloga, 1990], which can vary over a few orders of magnitude [Crisp and Baloga, 1994]. Keszthelyi and Denlinger [1996] also show that lava vesicularity and crystal content can affect the insulation provided by the early forming crust. Therefore, it is plausible that in some flows, the insulating crust could quickly decrease the cooling rate so as to inhibit preferred pathway development. In some flows, though, the cooling rate may be kept high, due to crust disruption or the specifics of that flow’s vesicularity and crystal content. [77] In summary, if the time or length scale of preferred pathway formation is larger than the emplacement time or length of a flow, then that flow will retain its sheet flow structure and preferred pathways will not form. This idea can be seen graphically in Figure 6, where a parameter box that is drawn for basalt flows with a slightly higher velocity range

1527

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

and a lower cooling rate range now corresponds to flows whose characteristic length scale will tend to lie above the 100 m diagonal line. We propose that typical sheet flows can remain stable over short distances if they slow (e.g., due to a decrease in influx, decrease in slope, or encountering an obstacle) before cooling sufficiently for the required viscosity contrasts to initiate pathway formation (volumelimited flow). Sheet flows that extend over longer distances could have been stabilized by an atypically rapid decrease in the cooling rate caused by the formation of a thick, stable crust and a relatively uniform emplacement rate. 4.3. Possible Connection Between Preferred Pathway Formation and Surface Morphology [78] As was discussed in section 4.1, no individual flow conditions (such as flow velocity or lava rheology) have been consistently correlated with observations of preferred pathways. However, as discussed in section 4.2, for typical basalt flow parameters, it appears plausible for the cooling rate to be a very influential parameter on pathway formation because it may vary more widely under common eruption conditions. There are many eruption and flow processes and morphologies which can influence the cooling rate of a flow. We propose that within most typical flows, the flow’s cooling rate will depend primarily on growth and disruption of the surficial crust. [79] As described by Crisp and Baloga [1990], the process of crust development and its influence on the thermal evolution of a lava flow are highly complicated. This perhaps explains why it has been difficult to phenomenologically or analytically describe the link between surface morphology and cooling of the molten lava core. However, if crustal metrics such as the mean and variance in thickness, vesicularity, or the fractional area of exposed core material can be better linked to the cooling rate of the molten lava core, then perhaps a clearer connection can be found between crustal morphology and the formation of preferred pathways. [80] Additionally, we propose that some confusing morphology may be the result of a very thick, strong crust forming over a sheet flow before flow conditions evolve to allow preferred pathways to form in the flow interior. In this case, it may take additional time for surface morphology to reflect the existence of subsurface pathways. Within flows where the strength of the surface crust is sufficiently high when preferred pathways form, a relatively smooth surface morphology may be retained despite the nonuniformity of the subsurface flow pattern, leading to contradictory observations. For example, Anderson et al. [1999] and Schaefer and Kattenhorn [2004] both found evidence of a pulsing stress distribution in sheet-flow-margin fractures which they attributed to the formation of low-viscosity pathways in the flow interior, an interpretation at odds with the flat sheet flow surface [Self et al., 1998]. We hypothesize that if the supply of lava to a sheet flow is cut off before newly formed pathways can modify the flow surface, then a flow that is morphologically similar to a sheet may not have been emplaced entirely by sheet-like spreading. [81] In summary, further studies of crust formation over a range of timescales are needed to improve our understanding of how this lava flow process can influence the cooling of the flow and enhance or hinder the formation of preferred pathways. We encourage further investigation into how the

early formation of a thick crust may inhibit the surface expression of subsurface preferred pathways, thus sometimes leading to a disconnect between surface morphology and flow emplacement dynamics. 4.4. A Brief Look at Silicic Lavas [82] In addition to basaltic lavas, we apply our model to silicic lavas. The viscosity of these lavas has a stronger temperature dependence than basalt, and silicic flows can contain channels (e.g., Big Glass Mountain [Fink, 1980]; Caliente dome [Harris et al., 2004]) although such features are rare. Lava tube systems have not yet been found within silicic flows [Williams and McBirney, 1979]. As with basalt sheet flows, we examine why preferred pathway formation may be inhibited within silicic flows. [83] Laboratory studies show that silicic lavas have higher β values at standard eruption temperatures than basalt lavas (Figure 1). For example, Hess and Dingwell [1996] measured β ~ 4.6 for granitic melts with very low water content (TH = 1200°C) and Fink [1980] estimated β between 3.5 and 9 on the basis of viscosity differences between surface and core rhyolitic lavas. These high values of β suggest that a thermo-viscous instability can form preferred pathways within silicic flows. [84] However, the length and timescales (~ u/δ and ~ 1/δ, respectively) of pathway formation in silicic lavas are greatly increased relative to basalt flows. This is due to the much lower cooling rate of these flows, as silicic flows commonly have a thick insulating crust. Harris et al. [2004] estimated a cooling rate of 0.05°C/h (δ ~ 107 s1) for a 10–50 m thick active silicic flow at Caliente dome, and Manley [1992] estimated a cooling rate of 0.15°C yr1 (δ ~1010 s1) on the basis of numerical simulation of heat flow from a 100–300 m thick relict flow in Idaho. Eruption velocities of silicic flows are also much less than those of basaltic flows (u = 1–5 m day1 ~ 105–104 m s1 [Harris et al., 2004; Manley, 1992]), but the difference in the cooling rate dominates the change in the characteristic length and timescales. Silicic flows have preferred pathway formation length scales at least 100X longer (~1 km) and timescales several orders of magnitude longer (~months) than those estimated for basalt flows. These calculated scales are consistent with observations of where and when channels become visible within silicic flows. For example, two >3 km-long silicic flows have channels which appear ~1 km from the vent [Fink, 1980; Harris et al., 2004]. One of these long flows was monitored as it erupted and its channel developed over several months [Harris et al., 2004]. Few silicic flows travel more than 1–2 km from their vent [Williams and McBirney, 1979], which may explain why flow morphologies associated with preferred pathway formation in silicic flows are generally not observed. [85] Channels within silicic lavas are not as narrow as one would expect given the probable β value of these flows. For β > 4–5, we expect f/w < 0.15, but silicic flows can have widths of a few hundred meters with a channel ~ 100 m wide (i.e., f/w ~ 0.3) [Harris et al., 2004]. This discrepancy between predicted and observed widths could reflect that we have computed β incorrectly. It could also reflect that the much thicker insulating crust and longer formation timescales limit the transverse velocities generated by flow focusing, or at least the surface expression of these dynamics (as discussed in section 4.3). Alternatively, it could reflect the

1528

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

invalidity of the model proposed in this study for silicic flows. There are important differences between basaltic and silicic lava rheology and flow geometry that can affect the model dynamics presented in this study. In particular, silicic flows have a much higher crystal and volatile content, which may invalidate some of the physical assumptions used to construct our model. Additionally, the geometry of a silicic flow is not well-approximated by a shallow flow as the much higher viscosity of the lava, the high slopes, and the restriction of surrounding topography create flow thicknesses to tens of meters. As a result, the momentum equation used in this model may be invalid. [86] Further observation and study of active silicic flows are needed to fully determine why channels are rare and tubes are nonexistent within silicic flows. Our model suggests that the expected characteristic scales of pathway formation within silicic flows due to a thermo-viscous instability could provide an explanation (Table 2; Figure 6). Owing to the very low cooling rates of silicic flows, our model predicts that very long flows (> 1 km) are generally needed for preferred pathways to develop. 4.5. Model Limitations and Extensions [87] A key component of any interpretative model analysis is an understanding of the limitations of that model. Here we outline a few physical influences which we have not included in our model and which should be considered when comparing model results with observations of natural flows. [88] 1. By using a simple exponential function with a constant β to relate dimensionless temperature and viscosity, the model assumes a smooth, continuous change in viscosity with temperature. However, this is a broad simplification. For example, as discussed in section 3.3, crystal content increases linearly with a temperature decrease, but its effect on viscosity is nonlinear. Crystal content of a lava has only a small effect on viscosity until crystallization exceeds 40%. At that level, crystals begin locking together and rapidly increase the viscosity [Cashman et al., 1999], resulting in a nonconstant β. High levels of crystallization can also cause other rheological responses, such as the inducement of a yield stress and the transition of basalt flows from pahoehoe to a`a [Cashman et al., 1999]. It can also affect the cooling rate [Dragoni and Tallarico, 1994]. These observations suggest that flow dynamics may be strongly dependent on initial eruption temperature and crystal content [Riker et al., 2009]. [89] 2. In the model, we have considered only coolinginduced changes in viscosity related to increasing crystal content within the flow. However, other conditions, such as volatile outgassing, can also affect lava viscosity significantly. Depending on the rate of deformation and the starting viscosity and temperature of a bubble-free mixture, exsolution of volatiles from the lava can lubricate flow or increase viscosity if the bubbles impede flow [Lejeune et al., 1999] or induce crystallization [Applegarth et al., 2013; Métrich and Rutherford, 1998]. Furthermore, degassing or the latent heat of crystallization can generate viscosity changes without the expected accompanying decrease in temperature [e.g., Crisp and Baloga, 1994; Wylie et al., 1999]. These processes generally do not have dominant effects within typical terrestrial subaerial lava flows (although they can greatly affect the dynamics of magma rise within a vertical conduit and

eruption style [e.g., Szramek et al., 2006; Wylie et al., 1999]). If they are important, it may be possible to use the model proposed here, but with a temperature-viscosity equation that includes both effects of cooling and the effects of temperature-independent processes, such as vesicularity changes. However, to get a reasonable connection between estimates of interior flow dynamics and preferred pathway dimensions, such an equation would need to be carefully fit to the specific rheology, composition, and eruption temperature of a flow. Additionally, if temperature-independent processes dominate the viscosity changes within a lava flow, the measured cooling rate of the natural lava flow (δ) will not be the only parameter affecting the flow’s characteristic length and timescales. [90] 3. Our model assumes a straight channel. Meanders and changes in flow geometry can disrupt the formation of crust and cause wall spallation [Riker et al., 2009], which in turn can affect flow dynamics. Preferred pathways can bifurcate when the lava flow is able to spread laterally, such as when the slope decreases or when the flow front enters a broad and flat region. The effects of a more complicated geometry should be examined in greater detail. [91] 4. Our model assumes flow over a featureless surface. Topography clearly plays a role in lava flow direction and behavior, so we expect that the effect of the temperaturedependent viscosity will be important only within flows with high effusion rates or of sufficient size to not be constrained entirely by local topography. We suggest that our featureless-terrain-model threshold requirements can perhaps be used as a starting point for examining a flow that is influenced by both topography and the thermo-viscous instability. For example, a flow that is only mildly affected by topography can perhaps be investigated using this model with an increased “effective” viscosity contrast ratio (β) or a decreased “effective” length scale of preferred pathway formation (u/δ). For flows that are narrow or thin relative to the surrounding topography, we expect flow dynamics to be dominated by the effect of topography. [92] 5. Our model produces preferred pathways which may become channels or tubes. Once tubes or channels have formed, they can continue to evolve in dimension and shape due to other processes, such as the formation and failure of blockages [Harris et al., 2005; Kauahikaua et al., 1998], changes in effusion rate [Bailey et al., 2006], thermal erosion [Greeley and Hyde, 1972; Kauahikaua et al., 1998; Swanson, 1973], or slope changes [Guest et al., 2012; Mazzarini et al., 2005]. Additionally, channels and tubes may form via other mechanisms. Therefore, comparison of model results with observed active or relict features must be carefully considered. In this study, we have attempted to compare model results with features found within stable and young portions of actively monitored large flows. [93] 6. Our study used estimates for cooling rate and its evolution in time and space appropriate for conduction through a stable and growing surface crust. In a flow where radiation dominates heat loss (such as within an uncrusted channel) or where crust is continually formed and disrupted causing rapidly changing cooling rates, then much higher or more stochastic cooling rates should be considered. As heat may be preferentially lost at high rates from uncrusted core material, core temperatures may become more uniform in the transverse direction which will depress viscosity

1529

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

variations occurring within a flow due to temperature variations. Generally, this type of flow occurs only transiently, especially for pahoehoe flows [Crisp and Baloga, 1990; Hon et al., 1994]. After this initial phase, flow usually settles into a more thermally protected state [Harris et al., 2005]. Therefore, more complex models of cooling can probably be neglected.

5.

Conclusion

[94] In this study, we investigate one mechanism that could contribute toward the formation of preferred pathways within surface lava flows. We hypothesize that a known temperature dependence of lava viscosity can influence large-scale flow dynamics. We test this hypothesis by modeling this effect in isolation, and we compare our model results to observations of channels in natural flows. Our model results are generally consistent with observations from a range of flows, and we therefore suggest that a thermo-viscous instability can influence evolution of a lava flow when certain conditions are met. In particular, we find that [95] 1. preferred pathways within lava flows can form spontaneously due to the natural temperature dependence of lava’s viscosity if a threshold viscosity contrast is attained; and [96] 2. preferred pathways will develop and stabilize over a distance that is related to a balance between flow velocity within the pathway and the cooling rate through the surface crust. [97] The needed viscosity contrast is generally attained within natural basalt flows and the length scale of preferred pathway development is short relative to common flow lengths. Thus, our model predicts that preferred pathways are likely to form spontaneously and rapidly within typical basalt flows. That basalt flows appear to often have the right conditions for the thermo-viscous instability to initiate preferred pathway formation may explain why tubes and channels are so common within basalt flows. Our model also provides plausible explanations for why channels and tubes are not found in some flows, namely basalt sheet flows and silicic lavas. For these kinds of flows, the necessary viscosity contrast may not to be attained due to high eruption temperatures for basalt sheet flows and rheology for silicic flows. Alternatively, even if the necessary viscosity contrast is attained, it is plausible that preferred pathways may not form as the expected time and length scales for pathway development are greater than the duration and length of the flow under reasonable flow conditions. This situation can occur if effusion rates are high, cooling rates decrease quickly, or flow duration is very short. [98] Although we have not shown conclusively that lava’s thermo-viscous instability plays a significant role in lava channel and tube initialization, and we acknowledge that topography plays a role in all flows (and likely the dominant role in flows where flow dimensions are comparable to the surrounding topography), our results are sufficiently positive to suggest that this dynamic should be considered and investigated further. Our model of preferred pathway formation has the potential to identify and estimate some salient conditions for lava channel and tube formation as a function of lava temperature, effusion rate, and crust formation.

Appendix A: Numerical Scheme and Assumptions [99] Within each time step, we directly evolve three scalar parameters within the flow: the stream function (ψ), temperature (T), and viscosity (μ) fields. The temperature field is first updated, using the stream function from the previous time step. The viscosity field is then calculated based on the new temperature field, and finally, the stream function is updated using the new viscosity field. The partial differential equations are evolved via finite-difference approximations. [100] To decrease time and memory requirements, the full domain is divided into subdomains which span the full width, but only a portion of the length of the channel. These subdomains were solved sequentially from the influx condition. As the temperature evolution step relies entirely on previous values, it can be solved in parallel. Due to the elliptic nature of the stream function evolution, that solution step must be done in series. [101] We note that extra care had to be taken in the integration of the energy equation, as spurious oscillations often developed and grew due to overshooting values in regions with steep gradients, yielding physically unrealistic overheating or overcooling. (Numerical instabilities resulting from an inability to properly resolve steep-gradients were also noted in Helfrich, [1995].) Thus, a flux-limiting scheme (which forces all flow to remain pseudolaminar) was necessary for the first-order upwinding advection scheme used in both the x- and y-directions: in the x-direction, we choose the slope value that is the minimum (in magnitude) between the upstream difference (Tx = (Ti,j  Ti,j1)/dx) and the central difference (Tx = (Ti,j + 1  Ti,j  1)/2dx). Within transverse flow, pure upwinding is used and flow is prevented (set = 0) if a local extrema is found (i.e., if the column profile T*,j is not monotonic). These flux-limiting algorithms were developed by van Leer [1979] for stable computations of compressible flow that are total variation diminishing and monotonicity preserving. Flux limiter schemes are commonly used in studies where strong gradients (i.e., schocks) can arise (e.g., steep viscosity gradients developing within incompressible flow in planetary interiors [Choblet et al., 2007] or sharp fronts in models of breaking waves [Bradford, 2012]). [102] A second-order spatial, first-order temporal finitedifference scheme was used for thermal diffusion in the y-direction. In calculations of the stream function, a fivegridpoint scheme was used to solve the modified Laplace’s equation, yielding a second-order solution. The solution was done explicitly, via the Jacobi method. [103] The influx boundary stream function was chosen so as to yield a constant influx velocity (ψinflux(y) =  uinflux y). Along the influx boundary, the temperature was defined as T = (1  perturbation(y)), where the perturbation function was a superposition of periodic functions with a whole number of periods (1–6 so as to have adequate resolution across each pathway) across the channel and with maximum magnitude (up to 0.1) as it approaches the side walls. Along the side walls, temperature was set to 0 and the stream function values were fixed, reflecting a free-slip condition. (Note that with the free-slip condition, the cold side walls influence the flow only through heat diffusion and we are neglecting the effects of along-wall turbulence in the model.) Within the channel, temperature was initialized to exponentially decay (sometimes

1530

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS

steeply) with distance from the source—this was found to be necessary to prevent numerical issues that resulted with an initial stepwise change in temperature near the source. [104] With all components of the algorithm, care was taken to insure that all Courant-Friedrichs-Lewy (CFL) conditions were satisfied in setting the time steps (relative to the grid spacing). We also checked that mass was conserved. Additionally, the resulting simulation was tested for solution convergence through comparison of model results (specifically, the shape and size of pathways) for a range of grid sizes. Generally, 40–50 elements in the transverse direction were found to be sufficient in detail and convergence, yielding similar results as simulations run with 2–4 times as many elements. Results also were comparable to those reported by Helfrich [1995]. We monitored numerical simulations to ensure that fronts and pathways did not focus down to regions smaller than a few elements; for example, β was limited to < 6 when 50 elements were used in the transverse direction. [105] Acknowledgments. Diniega was supported by an appointment to the NASA Postdoctoral Program, administered by Oak Ridge Associated Universities, at the California Institute of Technology Jet Propulsion Laboratory under a contract with NASA. Anderson acknowledges support from grant NNG05GL55G from the NASA Mars Fundamental Research Program. We thank two reviewers and two editors for their extensive comments which have greatly improved the manuscript. We also thank Christophe Sotin for suggesting the flux-limiting method.

References Anderson, S. W., and J. H. Fink (1992), Crease structures as indicators of emplacement rates and surface stress regimes of lava flows, Geol. Soc. Am. Bull., 104, 615–626. Anderson, S. W., E. R. Stofan, S. E. Smrekar, J. E. Guest, and B. Wood (1999), Pulsed inflation of pahoehoe lava flows: Implications for flood basalts, Earth Planet. Sci. Lett., 168, 7–18. Anderson, S. W., S. M. McColley, J. H. Fink, R. K. Hudson (2005), The development of fluid instabilities and preferred pathways in lava flow interiors: Insights from analog experiments and fractal analysis, in Kinematics and dynamics of lava flows, Geo. Soc. America Special Paper 396, edited by M. Manga and G. Ventura, 147–161. Anderson, S. W., S. E. Smrekar, and E. R. Stofan (2012), Tumulus development on lava flows: Insights from observations of active tumuli and analysis of formation models, Bull. Volcanol., doi:10.1007/s00445-012-0576-2. Applegarth, L. J., H. Tuffen, M. R. James, H. Pinkerton, and K. V. Cashman (2013), Direct observations of degassing-induced crystallization in basalts, Geology, 41(2), 243–246, doi:10.1130/G33641.1. Ault, W. U., J. P. Eaton, and D. H. Richter (1961), Lava temperatures in the 1959 Kilaeau eruption and cooling lake, Geol. Soc. Am. Bull., 72, 791–794. Bailey, J. E., A. J. L. Harris, J. Dehn, S. Calvari, and S. K. Rowland (2006), The changing morphology of an open lava channel on Mt. Etna, Bull. Volcanol., 68, 497–515, doi:10.1007/s00445-005-0025-6. Baloga, S., P. D. Spudis, and J. E. Guest (1995), The dynamics of rapidly emplaced terrestrial lava flows and implications for planetary volcanism, JGR, 100(12), 24,509–24,519. Blatt, H., R. Tracy, and B. Owens (2006), Petrology: Igneous, Sedimentary and Metamorphic, 3rd ed., Freeman, New York. Bonn, D., H. Kellay, M. Bräunlich, M. Ben Amar, and J. Meunier (1995), Viscous fingering in complex fluids, Physica A, 220, 60–73. Booth, B., and S. Self (1973), Rheological features of the 1971 Mount Etna lavas, Phil. Trans. R. Soc. London, A, 274, 99–106. Bradford, S. F. (2012), Improving the efficiency and accuracy of a nonhydrostatic surf zone model, Coastal Engineering, 65, 1–10, doi:10.1016/j.coastaleng.2012.02.004. Byrnes, J. M., and D. A. Crown (2001), Relationships between pahoehoe surface units, topography, and lava tubes at Mauna Ulu, Kilauea Volcano, Hawaii, J. Geophys. Res., 106(B2), 2139–2151. Calvari, S., and H. Pinkerton (1998), Formation of lava tubes and extensive flow field during the 1991–1993 eruption of Mount Etna, J. Geophys. Res., 103(B11), 27,291–27,301. Cashman, K. V., C. Thornber, and J. P. Kauahikaua (1999), Cooling and crystallization of lava in open channels, and the transition of pahoehoe lava to `a`a, Bull. Volcanol., 61, 306–323.

Castruccio, A., A. C. Rust, and R. S. J. Sparks (2013), Evolution of crustand core-dominated lava flows using scaling analysis, Bull. Volcanol., 75(681), 1–15, doi:10.1007/s00445-012-0681-2. Choblet, G., O. Čadek, F. Couturier, and C. Dumoulin (2007), Oedipus: A new tool to study the dynamics of planetary interiors, Geophys. J. Int., 170, 9–30. Costa, A., and G. Macedonio (2002), Nonlinear phenomena in fluids with temperature-dependent viscosity: An hysteresis model for magma flow in conduits, Geophys. Res. Lett., 29(10), 1402, doi:10.1029/2001GL014493. Crisp, J., and S. Baloga (1990), A model for lava flows with two thermal components, J. Geophys. Res., 95(B2), 1255–1270. Crisp, J., and S. Baloga (1994), Influence of crystallization and entrainment of cooler material on the emplacement of basaltic aa lava flows, J. Geophys. Res., 99(B6), 11,819–11,831. Crown, D. A., and S. M. Baloga (1999), Pahoehoe toe dimensions, morphology, and branching relationships at Mauna Ulu, Kilauea Volcano, Hawai’i, Bull. Volcanol., 61, 288–305. Dragoni, M., and A. Tallarico (1994), The effect of crystallization on the rheology and dynamics of lava flows, J. Volcanol Geoth. Res, 59, 241–252, doi:10.1016/0377-0273(94)90098-1. Duraiswami, R. A., N. R. Bondre, and G. Dole (2004), Possible lava tube system in a hummocky lava flow at Daund, western Deccan Volcanic Province, India, J. Earth System Sci., 113(4), 819–829, doi:10.1007/ BF02704040. Fink, J. (1980), Surface folding and viscosity of rhyolite flows, Geology, 8, 250–254. Fink, J. H., and R. W. Griffiths (1990), Radial spreading of viscous-gravity currents with solidifying crust, J. Fluid Mech., 221, 485–509. Fink, J. H., and R. W. Griffiths (1992), A laboratory analogue study of the surface morphology of lava flows extruded from point and line sources, J. Volcanol. Geotherm. Res., 54, 19–32. Fink, J. H., and J. Zimbelman (1990), Longitudinal variations in rheological properties of lavas: Puu Oo basalt flows, Kileauea Volcano, Hawaii, in Lava Flows and Domes, Emplacement Mechanisms and Hazard Implications, pp. 157–173, Springer-Verlag, New York. Glaze, L. S., and S. M. Baloga (1998), Dimensions of Pu’u O’o lava flows on Mars, J. Geophys. Res., 103(E6), 13,659–13,666, doi:10.1029/98JE01427. Glaze, L. S., S. W. Anderson, E. R. Stofan, S. Baloga, and S. E. Smrekar (2005), Statistical distribution of tumuli on pahoehoe flow surfaces: Analysis of examples in Hawaii and Iceland and potential applications to lava flows on Mars, J. Geophys. Res., 110, B08202, doi:10.1029/ 2004JB003564. Greeley, R. (1971), Observations of actively forming lava tubes and associated structures, Hawaii, Modern Geo., 2, 207–223. Greeley, R. (1987), The role of lava tubes in Hawaiian volcanoes,USGS Prof. Paper, 1350, 1590–1602. Greeley, R., and J. H. Hyde (1972), Lava tubes of the Cave Basalt, Mount St. Helens, Washington, Geol. Soc. Am. Bull., 83, 2397–2418. Gregg, T. K. P., and J. H. Fink (2000), A laboratory investigation into the effects of slope on lava flow morphology, J. Volcanol. Geotherm. Res., 96, 145–159. Gregg, T. K. P., and L. P. Keszthelyi (2004), The emplacement of pahoehoe toes: Field observations and comparison to laboratory simulation, Bull. Volcanol., 66, 381–391. Guest, J., A. Duncan, E. R. Stofan, and S. W. Anderson (2012), Effect of slope on development of pahoehoe flow fields: Evidence from Mount Etna, J. Volcanol. Geotherm. Res., 219–220, 52–62, doi:10.1016/j. jvolgeores.2012.01.006. Guest, J. E., C. R. J. Kilburn, H. Pinkerton, and A. M. Duncan (1987), The evolution of lava flow-fields: Observations of the 1981 and 1983 eruptions of Mount Etna, Sicily, Bull. Volcanol., 49, 527–540. Harris, A., S. Blake, and D. Rothery (1997), A chronology of the 1991 to 1993 Mount Etna eruption using advanced very high resolution radiometer data: Implications for real-time thermal volcano monitoring, J. Geophys. Res., 102(B4), 7985–8003. Harris, A. J. L., L. P. Flynn, O. Matias, W. I. Rose, and J. Cornejo (2004), The evolution of an active silicic lava flow field: An ETM + perspective, J. Volcanol. Geotherm. Res., 135, 147–168, doi:10.1016/j.jvolgeores.2003.12.011. Harris, A., J. Bailey, S. Calvari, and J. Dehn (2005), Heat loss measured at a lava channel and its implications for down-channel cooling and rheology,. in Kinematics and dynamics of lava flows, Geo. Soc. America Special Paper 396, edited by M. Manga and G. Ventura, 125–146. Harris, A., J. Dehn, and S. Calvari (2007), Lava effusion rate definition and measurement: A review, Bull. Volcanol., 70, 1–22, doi:10.1007/s00445007-0120-y. Helfrich, K. R. (1995), Thermo-viscous fingering of flow in a thin gap: A model of magma flow in dikes and fissures, J. Fluid Mech., 305, 219–238, doi:10.1017/S0022112095004605. Helz, R. T., C. Heliker, K. Hon, and M. Mangan (2003), Thermal efficiency of lava tubes in the Pu’u ‘O’o-Kupaianaha eruption, USGS Prof. Paper 1676, 105–120.

1531

DINIEGA ET AL.: THERMO-VISCOUS LAVA FLOW DYNAMICS Hess, K. U., and D. B. Dingwell (1996), Viscosities of hydrous leucogranitic melts: A non-arrhenian model, Am. Mineral., 81, 1297–1300. Holloway, K. E., and J. R. de Bruyn (2005), Viscous fingering with a single fluid, Can. J. Phys., 83, 551–564, doi:10.1139/P05-024. Hon, K., J. Kauahikaua, R. Denlinger, and K. Mackay (1994), Emplacement and inflation of pahoehoe sheet flows: Observations and measurements of active lava flows on Kilauea Volcano, Hawaii, Geol. Soc. Am. Bull., 106, 351–370. Hoover, S. R., K. V. Cashman, and M. Manga (2001), The yield strength of subliquidus basalts—Experimental results, J. Volcanol. Geotherm. Res., 107, 1–18. Jeffreys, H. (1925), The flow of water in an inclined channel of rectangular section, Phil. Mag., 49, 793–807. Kauahikaua, J. P., K. V. Cashman, T. N. Mattox, K. Hon, C. C. Heliker, M. T. Mangan, and C. R. Thornber (1998), Observations on basaltic lava streams in tubes from Kilauea Volcano, Hawai’i, J. Geophys. Res., 103, 27,303–27,324. Kauahikaua, J., D. R. Sherrod, K. V. Cashman, C. Heliker, K. Hon, T. N. Mattox, and J. A. Johnson (2003), Hawaiian lava flow dynamics during the Pu’u ‘O’o—Kupaianaha eruption: A tale of two decades, USGS Professional Paper 1676, p. 63–88. Keszthelyi, L., and R. Denlinger (1996), The initial cooling of pahoehoe flow lobes, Bull. Volcanol., 58, 5–18. Keszthelyi, L., and S. Self (1998), Some physical requirements for the emplacement of long basaltic lava flows, J. Geophys. Res., 103(B11), 27,447–27,464. Keszthelyi, L., S. Self, and T. Thordarson (2006), Flood lavas on Earth, Io, and Mars, J. Geol. Soc. London, 163, 253–264. Kilburn, C. R. J. (1993), Lava crusts, aa flow lengthening and the pahoehoe-aa transition, in Active Lavas, edited by C. R. J. Kilburn and G. Luongo, UCL Press Limited, London, 263–280. Kilburn, C. R. J. (2004), Fracturing as a quantitative indicator of lava flow dynamics, J. Volcanol. Geotherm. Res., 132, 209–224, doi:10.1016/ S0377-0273(03)00346-9. Kilburn, C. R. J., and R. M. C. Lopes (1988), The growth of aa lava flow fields on Mount Etna, Sicily, J. Geophys. Res., 93(B12), 14,759–14,772. Kilburn, C. R. J., and R. M. C. Lopes (1991), General patterns of flow field growth: Aa and blocky lavas, J. Geophys. Res., 96(B12), 19,721–19,732. van Leer, B. (1979), Towards the ultimate conservative difference scheme, J. Comput. Phys., 32(1), 101–136. Lejeune, A. M., Y. Bottinga, T. W. Trull, and P. Richet (1999), Rheology of bubble-bearing magmas, Earth Planet. Sci. Lett., 166, 71–84. Lipman, P. W., and N. G. Banks (1987), A’a flow dynamics, Mauna Loa 1984, USGS Prof. Paper 1350, 1527–1567. Manley, C. R. (1992), Extended cooling and viscous flow of large, hot rhyolite lavas: Implications of numerical modeling results, J. Volcanol. Geotherm. Res., 53, 27–46. Mazzarini, F., M. T. Pareschi, M. Favalli, I. Isola, S. Tarquini, and E. Boschi (2005), Morphology of basaltic lava channels during the Mt. Etna September 2004 eruption from airborne laser altimeter data, Geophys. Res. Lett., 32, L04305, doi:10.1029/2004GL021815. Métrich, N., and M. J. Rutherford (1998), Low pressure crystallization paths of H2O-saturated basaltic-hawaiitic melts from Mt Etna: Implications for open-system degassing of basaltic volcanoes, Geochim. Cosmochim. Acta, 62(7), 1195–1205. Moore, H. J. (1987), Preliminary estimates of the rheological properties of 1984 Mauna Loa lava, USGS Prof. Paper 1350, 1569–1588. Peck, D. L. (1978), Cooling and vesiculation of Alae Lava Lake, Hawaii, USGS Prof. Paper 935-B.

Peterson, D. W., and D. A. Swanson (1974), Observed formation of lava tubes during 1970–71 at Kilauea Volcano, Hawaii, Stud. Speleol., 2(6), 209–222. Peterson, D. W., R. T. Holcomb, R. I. Tilling, and R. L. Christiansen (1994), Development of lava tubes in the light of observations at Mauna Ulu, Kilauea Volcano, Hawaii, Bull. Volcanol., 56, 343–360. Pinkerton, H., and R. S. J. Sparks (1978), Field measurements of the rheology of lava, Nature, 276, 383–385. Pinkerton, H., M. James, and A. Jones (2002), Surface temperature measurements of active lava flows on Kilauea volcano, Hawai’i, J. Volcanol. Geotherm. Res., 113, 159–176, doi:10.1016/S0377-0273(01)00257-8. Riker, J. M., K. V. Cashman, J. P. Kauahikaua, and C. M. Montierth (2009), The length of channelized lava flows: Insight from the 1859 eruption of Mauna Loa Volcano, Hawai’i, J. Volcanol. Geotherm. Res., 183, 139–156. Rowland, S. K., and G. P. L. Walker (1990), Pahoehoe and aa in Hawaii: Volumetric flow rate controls the lava structure, Bull. Volcanol., 52, 615–628. Saffman, P. G., and G. Taylor (1958), The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous liquid, Proc. Roy. Soc. Lond Math Phys. Sci., 245(1242), 312–329. Sato, H. (2005), Viscosity measurement of subliquidus magmas: 1707 basalt of Fuji volcano, J Miner Petrol Sci, 100, 133–142. Schaefer, C. J., and S. A. Kattenhorn (2004), Characterization and evolution of fractures in low-volume pahoehoe flows, eastern Snake River Plain, Idaho, Geol. Soc. Am. Bull., 116, 322–336. Self, S., L. Keszthelyi, and T. Thordarson (1998), The importance of pahoehoe, Annu. Rev. Earth Planet. Sci., 26, 81–110. Shaw, H. R. (1969), Rheology of basalt in the melting range, J. Petrol., 10(3), 510–535. Soule, S. A., K. V. Cashman, and J. P. Kauahikaua (2004), Examining flow emplacement through the surface morphology of three rapidly emplaced, solidified lava flows, Kilauea Volcano, Hawaii, Bull. Volcanol., 66, 1–14, doi:10.1007/s00445-003-0291-0. Swanson, D. A. (1973), Pahoehoe flows from the 1969–1971 Mauna Ulu eruption, Kilauea, Hawaii, Geol. Soc. Am. Bull., 84, 615–626. Swanson, D. A., T. L. Wright, and R. T. Helz (1975), Linear vent systems and estimated rates of magma production and eruption for the Yakima Basalt on the Columbia Plateau, Am. J. Sci., 275, 877–905. Szramek, L., J. E. Gardner, and J. Larsen (2006), Degassing and microlite crystallization of basaltic andesite magma erupting at Arenal Volcano, Costa Rica, J. Volcanol. Geotherm. Res., 157, 182–201. Walker, G. P. L. (1971), Compound and simple lava flows and flood basalts, Bull. Volcanologique, 35(3), 579–590. Walker, G. P. L. (1991), Structure, and origin by injection of lava under surface crust, of tumuli, “lava rises,” “lava-rise pits,” and “lava-inflation clefts” in Hawaii, Bull. Volcanol., 53, 546–558. Whitehead, J. A., and H. R. Helfrich (1991), Instability of flow with temperature-dependent viscosity: A model of magma dynamics, J. Geophys. Res., 96, 4145–4155. Williams, H. and A. R. McBirney (1979), Volcanology, Freeman, Cooper, and Co., 391 pp. Wylie, J. J., and J. R. Lister (1995), The effects of temperature-dependent viscosity on flow in a cooled channel with application to basaltic fissure eruptions, J. Fluid Mech., 305, 239–261. Wylie, J. J., K. R. Helfrich, B. Dade, J. R. Lister, and J. F. Salzig (1999), Flow localization in fissure eruptions, Bull. Volcanol., 60, 432–440. Zimbelman, J. R. (1998), Emplacement of long lava flows on planetary surfaces, J. Geophys. Res., 103, 27,503–27,516, doi:10.1029/98JB01123.

1532