Layout 6 1 ANNALS OF GEOPHYSICS, 62, 2, VO226, 2019; doi: 10.4401/ag-7745 “CRUST DEVELOPMENT INFERRED FROM NUMERICAL MODELS OF LAVA FLOW AND ITS SURFACE THERMAL MEASUREMENTS„ Igor Tsepelev1, Alik Ismail-Zadeh2,3,*, Yulia Starodubtseva1, Alexander Korotkii1,4, Oleg Melnik5 (1) Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences, Yekaterinburg, Russia (2) Institute of Applied Geosciences, Karlsruhe Institute of Technology, Karlsruhe, Germany (3) Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Moscow, Russia (4) Institute of Natural Sciences and Mathematics, Ural Federal University, Yekaterinburg, Russia (5) Institute of Mechanics, Lomonosov Moscow State University, Moscow, Russia 1. INTRODUCTION A flow of lava is a common volcanic feature on the Earth’s surface. A lava flow can destroy villages or se− riously damage cities as happened in and around Catania, Sicily, during the eruption of Etna volcano in 1669 [Branca et al., 2013]. Although lavas have been used as a construction material for centuries and have been a source of nutrients in agriculture, lava flows remain a threat to human activity [Poland et al., 2016]. The hazard is not negligible as hot lava kills vegeta− tion, destroys infrastructure, and may trigger floods due to melting of snow/ice [e.g., Harris, 2015]. Lava flow fields are widespread, e.g., the Permian−Triassic Article history Receveid April 30, 2018; accepted February 20, 2019. Subject classification: Lava crust; Nonlinear heat flux; Lava rheology; Numerical modelling; Data assimilation. ABSTRACT Propagation of a lava flow is governed by slope topography, magma rheology, heat exchange with the atmosphere and the underlying ter− rain, and the rate of the eruption. Highly viscous crust is formed due to cooling and solidification of the uppermost layer of the flow. We consider here two numerical model problems for lava flows, both based on the fundamental physics of a hot fluid flow: a model problem, where thermal conditions (e.g. temperature and heat flow) at the lava surface are unknown a priori (a direct model problem), and a model problem, where the lava surface conditions are known and determined from observations (an inverse model problem). In both models, the lava viscosity depends on temperature and the volume fraction of crystals. By way of solving the direct model problem, we perform a para− metric study of steady state lava flows to investigate the influence of the heat flux, viscosity, and effusion rate on the lava crust devel− opment. Numerical experiments show that a lava crust becomes thicker in the case of the nonlinear heat transfer compared to the case of a linear heat flow at the interface of lava with the atmosphere. Also, the crust thickens at lower lava effusion rates, while higher rates re− sult in a rapid lava advection, slower cooling, and development of a thinner crust. Moreover, a lava crust becomes thicker with a higher coefficient of conductive heat transfer, or a higher lava viscosity, or the growth of effective emissivity of the lava surface. By way of solv− ing the inverse model problem, we use an assimilation technique (that is, a method for an optimal combination of a numerical model of lava flows with observations) to propagate the temperature and heat flow, inferred from measurements at the interface between lava and the atmosphere, into the lava flow interior and to analyse the evolving lava crust. Results of thermal data assimilation illustrate that the physical parameters of lava flows, including the thickness of it crust, can be recovered from measured surface thermal data well enough at least for slow effusion rates. Siberian and Cretaceous−Tertiary Deccan large ig− neous provinces, the Quaternary Yellowstone rhyolitic lava flows or the Holocene lava flows in Iceland [Christiansen, 2001; Self et al., 2008; Bryan et al., 2010; Pinton et al., 2018]. Several regions in the world are prone to lava flow hazard including in Hawaii/USA [Patrick et al., 2017], Iceland [Pedersen et al., 2018], Italy [Corsaro et al., 2009], Réunion/France [e.g., Sol− dati et al., 2018], Kamchatka/Russia [Belousov and Be− lousova, 2018], and some other regions of the active or historic effusive volcanism [e.g., Chevrelet al., 2016; Dietterich et al., 2018]. Volcanic eruptions produce a variety of lava flows depending on the chemical composition and tempera− ture of the erupting material, and the topography of the surface over which the lava flows [e.g., Griffiths, 2000; Rumpf et al., 2018]. Eruption (effusion) rates control lava flow dynamics: with higher effusion rates a lava flows the more rapidly and longer [Walker, 1973; Row− land and Walker, 1990; Harris et al., 1998; Castruccio and Contreras, 2016]. The rapid development of ground−based thermal cameras, drones and satellite data allows collection of repeated thermal images of the surface of active lava flows during a single lava flow eruption [Calvari et al., 2005; Wright et al., 2010; Kelfoun and Vargas, 2015]. These data require devel− opment of appropriate quantitative methods to link subsurface dynamics with surface observations and measurements. Numerical modelling plays an essential role in un− derstanding lava flow patterns as well as their mor− phology and thermal evolution [e.g., Costa and Macedonio, 2005; Cordonnier et al., 2015; and refer− ences herein]. Numerical modelling of lava flows was advanced for the last few decades moving from one − (1D) to three−dimensional (3D) flows. Cordonnier et al., [2015] provided a recent overview of the volcanologi− cal community’s codes for modelling of lava flow em− placement. We highlight here only those codes which are related to deterministic models (the mathematical models defined by a set of governing equations and boundary/initial conditions) as they are related to the type of modelling we used in our study. The 1D chan− nel lava flow code developed by Harris and Rowland [FLOWGO, 2001] is very fast to run as it does not cal− culate dynamic properties of fluid flow; but its limita− tion is that it does not consider vertical variations within a lava flow. Numerical codes for depth−average models of lava flow (using a shallow−water approach) were developed by Macedonio et al. [2005] and Kelfoun and Druitt [2005]; the codes provide rather fast compu− tations, but their limitation is that the lava viscosity does not vary vertically. Smooth particle hydrodynam− ics codes have been also employed in modelling of lava flows [e.g., Herault et al., 2011]. The codes use meshless numerical methods to calculate flow of individual par− ticles, for which physical characteristics are specified. 3D finite−volume or finite−element fluid−dynamics codes allow for studying lava flows on a complex to− pography, using complex rheology, and various bound− ary conditions, although they are slower compared to those mentioned earlier. Tsepelev et al. [2016] used the OpenFOAM toolbox (ver. 2.3.0, http://www.open− foam.org) to develop a code for a 3D isothermal lava flow with crustal pieces in a channel and on a real to− pography. Isothermal models of viscous flows have demonstrated how the lava generated by effusive vol− canic eruptions would advance in the absence of cool− ing, and these have been supplemented by more sophisticated models including complex rheology, cool− ing, solidification, and dynamic interaction of lava and its crust [e.g., Griffiths, 2000; Costa and Macedonio, 2005; Fujita and Nagai, 2015; Tsepelev et al., 2016]. A mathematical model of lava flow can be described by a set of partial differential equations and boundary and initial conditions defined in a specific domain. A model links the causal characteristics of a lava flow (e.g., lava viscosity, thermal conditions at lava inter− faces with the atmosphere and the terrain) with its ob− served/measured properties (e.g., lava temperature and its flow rate). The aim of the direct model problem is to determine the properties of a lava flow based on the knowledge of its causes, and hence to find a solution to the mathematical problem for a given set of parameters and coefficients. An inverse model problem is consid− ered when there is a lack of information on the causal characteristics but some information on lava flow prop− erties exists. For example, a mathematical problem of determination of the thermal and dynamic characteris− tics of lava flow from measurements of the absolute temperature on its upper surface belongs to a class of inverse problems [e.g., Kirsch, 1996; Kabanikhin, 2011], and hence, a small error in measured temperature at the lava surface can result in a significant deviation from the true solution of the problem (e.g. lava temperature). To solve the inverse problem, special techniques for data assimilation [e.g., a variational data assimilation; Is− mail−Zadeh et al., 2016] should be employed to con− TSEPELV ET AL. 2 strain the unknown thermal condition at the lower boundary of the lava flow from observations at the lava surface [e.g., Korotkii et al., 2016]. The basic principle of this assimilation in the case of lava flows is to consider a thermal condition (temperature or heat flow) at the lava bottom (the interface between moving lava and the underlying terrain over which it flows) as a control variable and to find the optimal condition at the bot− tom by minimizing the misfit between the measure− ments at the lava interface with the atmosphere and the model solution at the same interface. To numerically model a lava flow, boundary condi− tions should be known. Meanwhile the temperature or heat flux at the lava bottom is unknown as it is almost impossible to measure it. However, airborne and space measurements of temperature (heat flux) at the interface of lava flow with the atmosphere allow to determine thermal conditions at the lava bottom using data as− similation. Once the boundary conditions at the lava bottom are determined, the lava flow can be modelled to estimate its extent, temperature and flow rate. In this paper, we develop two−dimensional numer− ical models of steady−state lava flows to analyse the development of the lava crust depending on the heat flux into the atmosphere from the lava surface, the lava viscosity, its discharge rate, and the lava front propagation rate (Sections 2−4). We assimilate syn− thetic thermal data (a model of thermal measurement at the lava interface with the atmosphere) into the lava interior to determine the unknown thermal condition at the lava bottom (Sections 5−8). We develop direct and inverse models of lava flows introducing the fol− lowing complexities: (i) the lava viscosity depending on temperature and the volume fraction of crystals; (ii) nonlinear heat transfer mechanism at the surface of lava with the atmosphere, including the radiant and convective parts of the total heat flux; (iii) varying flow conditions at the lava surface; and (iv) lava flow on a real surface topography. 2. THE DIRECT MODEL PROBLEM OF LAVA FLOW To state mathematically a problem of lava flow, we need to consider the following components of a model: computational domain (e.g. topography on which lava flows or slope of the inclined surface or channel as well as vent geometry); basic equations which govern the lava flow; the boundary condition on flow veloc− ity (e.g., effusion rate, free−slip or no−slip or free sur− face conditions); the boundary conditions on temperature (e.g., effusive temperature as well as ra− diative, convective or conductive heat flow at the in− terface with the atmosphere and with the land); and physical properties of the lava (that is, lava viscosity, density, thermal conductivity, rheology etc.). Here we present a mathematical statement of the direct problem of steady−state lava flow in model domain Ω ⊂ R2 (Figure 1). The model domain covers a part of lava flow at some distance from the volcanic vent and the lava flow front. The boundary of the model domain consists of the following parts: Г1 is a line seg− ment connecting points A and E (xA = (x1 A, x2 A) = (0.0,1.5), xE = (x1 E, x2 E) = (0.0,2.5)); Г2 is a circular arc connecting points A and B (xB=(x1 B, x2 B) = (10.0,0.0)); Г3 is a line segment connecting points B and C (xC = (x1 C, x2 C) = (10.0,0.5)); and Г4 is a circular arc con− 3 LAVA CRUST DEVELOPMENT FIGURE 1. Model geometry. � � � � E D �1 2 4 3 10h 2h h h/2 A B C vent flow front TSEPELV ET AL. 4 necting points C and E. All coordinates are in meters. Although the lava rheology can be more complicated, we assume that the lava behaves as a Newtonian fluid with a temperature− and volume−fraction−of−crystals− dependent viscosity and temperature−dependent density and thermal conductivity. We note that crys− tallisation is shown to be the primary process, which is responsible for the increase of the viscosity during emplacement of mafic lava flows [Chevrel et al., 2013]. The flow is described by the Navier−Stokes, continuity, and heat equations [e.g., Hidaka et al., 2005; Ismail− Zadeh and Tackley, 2010]. Although a lava flow can be non−steady state depending on an effusion rate, for a simplicity of the mathematical problem’s description we assume a steady−state lava flow in the modelling [a justification to this assumption is provided by Korotkii et al., 2016]. In the case of viscous flows with the small Reynolds number (Re<1), the governing equations can be represented in the following form: , (1) , (2) . (3) Here x = (x1, x2) ∈ Ω are the Cartesian coordinates; u = (u1(x),u2(x)) is the vector velocity; p = p(x) is the pressure; T = T(x) is the temperature; η is the viscosity; r is the temperature−dependent density; k is the tem- perature-dependent heat conductivity; g is the accele- ration due to gravity; c* (= 1000 J kg -1K-1) is the spe- cific heat capacity; e2 = (0,1) is the unit vector; ∇, T, and 〈·,·〉 denote the gradient vector, the transposed matrix, and the scalar product of vectors, respectively. Lava viscosity depends on temperature, the volume fraction of crystals, water content and some other phy- sical parameters [e.g., Dragoni, 1989; Giordano and Din- gwell, 2003; Harris and Allen, 2008; Costa et al., 2009]. In this study, we consider two models of viscosity: a tem- perature-dependent viscosity [Dragoni, 1989] , (4) where n = 4 x 10-2 K-1, and Tm (= 1333 K) is the typi− cal lava melting temperature; and a volume-fraction- of-crystals−dependent viscosity [Marsh, 1981; Costa et al., 2009] , (5) where , (6) , (7) , (8) f is the volume fraction of crystals; f* (= 0.384) is the specific volume fraction of crystals; θ = (T-Ts)/(Tm-Ts), Ts (= 1053 K) is the solidus temperature; B (=2.5) is the Einstein coefficient’s theoretical value (it was experi− mentally determined that the Einstein coefficient varies from 1.5 to 5; Jeffrey and Acrivos, 1976); δ = 13-γ, ξ = 2.0x10-4, γ = 7.701 [Lejeune and Richet, 1995; Costa et al., 2009], θ1 = 0.5, b2 = 3/2 [Wright and Okamura, 1977], and b1 = √30 [Marsh, 1981]. Finally, we use in the modelling the following lava viscosity restricting its exponential growth with tem− perature by the prescribed value η0 (Figure 2): (9) where η* (=10 6 Pa s) is the typical lava viscosity, and η0 takes values from 10 3 to 106 in the modelling. In addition to these equations, we consider the fol− lowing relationships for the thermal conductivity and the density: , (10) , , 11) where, T = 1473 K; rm (= 2750 kg m -3) and rc (= 2950 kg m-3) are the typical lava melt density and the typ− ical crystal density corresponding to the crystals composed of 50% olivine and 50% plagioclase, re− spectively. As the effective density of the lava crust with about 20% volume of vesicles is estimated to be about 2200 kg m-3 [Kilburn, 2000], we introduce the term f(T)δr (where δr = 750 kg m-3) to reduce the density of the lava crust and to permit crustal pieces drifting with and not sinking into the lava. We consider the following conditions at the bound− ary Γ = Γ1 U Γ2 U Γ3 U Γ4 of the model domain. The velocity u1 and the temperature T1 are prescribed at the ~ k (T ) = 1.15 + 5.9 ⋅10−7 (T − T )2 , T < T 1.15 + 9.7 ⋅10−6 (T − T )2 , T > T ⎧ ⎨ ⎪ ⎩⎪ ~ ~ ~ ~ r=r m 1−f (T )( ) + rcf (T ) − (T )δrf ∇p = ∇ ⋅ η (T ) ∇u + ∇uT( )( ) − (T )ge2r ∇ ⋅ u = 0 ∇ ⋅ k (T ) ∇T( ) = u, c*∇ T( )r η (ϕ ) = 1+ ϕ δ( ) 1− F ( , , )⎡⎣ ⎤⎦ − B *~ ϕ ξ γ f ϕ = f / f * = 0.5 1− erf b 1 ( − 1 ) / b 2( )⎡⎣ ⎤⎦ / *f θ θ F (ϕ, ξ,γ ) = (1−ξ ) ⋅ erf π 2(1− ) (1+ ) ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ξ ϕ ϕγ η(T ) = * 12 (T ),η η η 12 (T ) = η 1 η 2 , 1 2 < 0 , η 0 , 1 2 > 0 , ⎧ ⎨ ⎪ ⎩⎪ η η η η η η η 2 (T ) = ( (T ) / * )η f f~ η 1 (T ) = exp n Tm − T( )( ) 5 LAVA CRUST DEVELOPMENT left boundary Γ1: , (12) No slip condition and zero heat flux (no heat lost to the base via conduction) or temperature T2 are pre− scribed at the lower boundary Γ2: , (13) where n is the outward unit normal vector at a point of the boundary. We assume at the right boundary Γ3: , (14) where σ = η(∇u+∇uT) is the deviatoric stress tensor. At the interface of lava flow with the atmosphere (boundary Γ4), zero normal flow and free−slip tangential conditions are prescribed. Lava cooling in the atmos− phere is the simultaneous action of radiation and ther− mal convection. Although the lava emissivity is such that radiation is the main mechanism of heat flow dur− ing the very first minutes of cooling, convection by the surrounding air becomes very important after just a few minutes [Neri, 1998]. A heat flow at the boundary Γ4 (the lava surface) can be then either linear: , (15) , or nonlinear (16) , where the first term in the right−hand side of Equa− tion (16) is related to the radiant heat flow, and the sec− ond term to the convective heat flow. Here Ta (= 300 K) is the air temperature; l ~ is the conductive heat transfer coefficient; ς (=5.668108x10-8 W m-2 K-4) is the Stefan–Boltzmann constant; ε is the effective emissivity of the lava surface; l is the dimensional (W m−2 K−M) constant; and M (= 1 or 4/3 in this study) is the pow− er−law exponent. In the dimensionless form the linear and nonlinear heat flow conditions can be represented as k〈∇T,n〉 = -Nu(T-1) and −k〈∇T,n〉 = a1 (T 4 -1) + a2(T -1) M respectively, where Nu = l ~ h / k* is the Nusselt number, a1 = hk* -1T3aςε, a2 = hk* -1T1/3λ, k* = rmc*ĸ* is the typ− ical heat conductivity, ĸ* (= 10 -6 m2 s-1) is the typical heat diffusivity of lava as it is a very poor thermal con− ductor, and h (= 2 m) is the typical lava thickness. The effective emissivity of the lava surface ε and the di− mensional constant λ vary from 0.6 to 0.95 and from 1 to 10, respectively [Neri, 1998; Harris et al., 1998]. The problem (1)−(16) is transformed to a dimensionless form assuming that length, temperature, viscosity and heat conductivity are normalised by h, Ta, η*, and k*. 3. SOLUTION METHOD To solve the problem (1)−(16), we use the numerical approach and code developed by Korotkii et al. [2016]. The OpenFOAM toolbox was used to develop a solver for this study. The finite volume method [e.g., Ismail−Zadeh and Tackley, 2010] is employed in the software to solve the numerical models on multiprocessor computers. The model domain Ω was discretized by about 104 polyhe− dral finite volumes. The SIMPLE method [Patankar and Spalding, 1972] was used to determine velocity and pressure at a given temperature (the relaxation parame− ters are chosen to be 0.7 and 0.3 for the velocity and pressure, respectively). The use of polyhedral finite vol− umes has some advantages compared to hexahedral fi− nite volumes employed earlier by Korotkii et al. [2016]; for example, it enhances the convergence rate of the SIMPLE method, in some cases by a factor of two. To implement the SIMPLE method, we employ the conjugate−gradient method [Ismail−Zadeh and Tackley, FIGURE 2. Dimensionless viscosity as a function of temperature. 1: (dashed−dotted line) temperature−dependent vis− cosity (η1, Equation 4); 2: (dashed line) crystal−vol− ume−fraction−dependent viscosity (η2, Equation 5); 3: (dotted line) combines viscosity (η1 · η2); and 4: (solid line) the viscosity used in the model (Equation 9). u,n = 0, σ n − σ n,n n = 0, − k ∇T ,n = ςε T 4 − Ta 4( ) + T − Ta( ) M l u = u 1 ( x 2 ), T = T 1 u = 0, ∂T ∂n = 0 σ n = 0, p = 0, ∂T ∂n = 0 u,n = 0, σ n − σ n,n n = 0 k ∇T ,n = − (T − Ta )l ~ 2010] to solve a set of linear algebraic equations (SLAE) with positive−definite and symmetric matrices, which are obtained after the discretization of the Stokes equa− tion. All computations were performed using one CPU Intel Core i3 2.13 GHz with 4GB memory, OS Kubuntu 14.4. In the case of the heat equation, SLAE were solved by the biconjugate gradient stabilized method [van der Vorst, 1992] with the pre−conditioner of incomplete LU−decomposition. The linear Gaussian scheme with a flow control was used to discretize the Laplace opera− tor. To approximate the convective operator, we em− ployed the total variation diminishing (TVD) method with the minmod limiter [Sweby, 1984; Ismail−Zadeh et al., 2007]. The nonlinear condition for temperature (16) at the upper model boundary is solved the simple iteration method, and at each iteration m+1 of the SIMPLE method the following scheme is used: . The descent step length 0 < γ < 1 is determined ex− perimentally to provide the convergence to the iterative process. This iterative process is stopped when the tem− perature, pressure and velocity residuals become less than 10-4. An average computational time for solution of the problem (1)−(16) was about 3 min for the chosen dimension of the discrete problem. 4. RESULTS: EVOLVING LAVA CRUST In this section, we present the results of several parametric case studies to analyse the influence of the heat flux at the interface between the lava flow and the atmosphere, and of the effusion rate on the lava flow temperature, viscosity, velocity, and the thickness of the lava’s crust. The crust is defined in the numerical mod− els as the area, where the lava viscosity is greater than η0/2. The parameters of the case study models are listed in Table 1, and η0=10 3. Numerical experiments have been performed for lin− ear (Equation 15) and nonlinear (Equation 16) heat flow conditions at the lava surface with the atmosphere for different injection rates |u1| (hereinafter we assume that the injection rate is the rate at which lava enters into the model domain, and the injection rate is proportional to the effusion rate per unit area). It is shown (see Sup- plementary Material S1, Figure S1) that in the case of the nonlinear heat flow, the lava cools down faster than in the case of linear heat flux, and hence the lava’s crust becomes thicker in the case of nonlinear heat flux. At higher injection rates, the lava is advected rapidly and hence cools slowly; whereas at lower rates, the crust becomes thicker (Figure 3a). The higher coefficient of the conductive heat transfer, the thicker the lava crust. When the injection rate is 0.023 m s-1 and the viscosity is 106 Pa s, a lava crust does not evolve at l ~ =1 (Figure 3b); at higher l ~ , a thin crust (about 6 cm) starts to develop at a certain TSEPELV ET AL. 6 FIGURE 3. Thickness of the lava crust with distance from the lava injection into the model domain depending on: (a) the injection rate |u1|, (b) conductive heat transfer l ~ , and (c) the typical lava viscosity η*. See Table 1 for case studies m1.1−m3.4. ∇T ( m+1) ,n = (1−γ ) ∇T ( m) ,n − − γ a 1 T ( m)( ) 4 − Ta 4⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + a2 T ( m) − Ta( ) 4/3⎡ ⎣ ⎢ ⎤ ⎦ ⎥ distance from the place of a lava injection into the model domain (boundary Γ1). With the increase of the lava viscosity (Figure 3c), the crust becomes thicker, and it reaches about 60−70 cm in the centre of the model domain. We note that the lava flow viscosity varies by several orders of magnitude depending on its composition. The typical viscosity is about 102-103 Pa s for basaltic lava flows, and increases up to 109- 1011 Pa s during emplacement of andesitic lava [Chevrel et al., 2016]. For all cases, the crust becomes thinner closer to the left model boundary because of a converging shape of the domain in the model. Tem− perature at the lava surface with the atmosphere de− creases with the increase in its effective emissivity and the power−law constant M, although the depen− dence on M is weak (Figure 4). At higher effective emissivity, the lava crust thickens (Figure 5). We have also performed numerical experiments for linear and nonlinear heat flow conditions with vary− ing slip conditions at the upper surface of the lava flow model, particularly, considering no−slip at the right part of the upper surface (a circular arc con− necting points C and D in Figure 1, which is about 1/5 of the model domain’s horizontal length). The no− slip condition can model the case, when the front of lava flow covered by a crust prevents lava to flow. The results show that the crust beneath the no−slip area becomes much thicker than beneath the free−slip area at both linear and nonlinear heat flow at the cooling surface of lava flow (see Supplementary Ma- terial S2, Figure S2). 7 LAVA CRUST DEVELOPMENT FIGURE 4. Dimensionless temperature at the interface between lava and the atmosphere with distance from the lava injection depending on the effective emissivity (case studies m4.1−m4.3) and the power−law exponent M (case studies m4.4 − m4.6). See Table 1 for the case studies. Model case |u1|, ms-1 η*, Pa s l ~ , W m-2 K-1 l, W m-2 K-M M ε m1.1 0.058 106 5.5 - - - m1.2 (reference) 0.023 106 5.5 - - - m1.3 0.012 106 5.5 - - - m2.1 0.023 106 1.0 - - - m2.3 0.023 106 10.0 - - - m3.1 0.023 105 5.5 - - - m3.3 0.023 107 5.5 - - 0 m3.4 0.023 108 5.5 - - - m4.1 0.023 106 - 5.5 1 0 m4.2 0.023 106 - 5.5 1 0.6 m4.3 0.023 106 - 5.5 1 0.95 m4.4 0.023 106 - 5.5 4/3 0 m4.5 0.023 106 - 5.5 4/3 0.6 m4.6 0.023 106 - 5.5 4/3 0.95 TABLE 1. Model parameters for case studies 5. THE DATA ASSIMILATION (INVERSE) MODEL PROBLEM Here we study the problem of determining thermal and flow characteristics of the lava and the thickness of its crust from measurements of temperature and heat flow at the interface with the atmosphere. We modify Equations (1)−(3) applying the Boussinesq approxima− tion [Boussinesq, 1903] to the equation by keeping the density constant everywhere except for the buoyancy term in the Stokes equation. In this approximation, the dimensionless Stokes, continuity, and heat equations take the form: , (17) , (18) , (19) where and the Rayleigh number is determined as Ra = α*grm ΔTh 3η* -1ĸ* -1, α* (= 10 −5 K−1) is the thermal expansivity; ΔT = Tm - Ts is the temperature contrast; ĸ(T) = k(T) / (rm c*) is the thermal diffusivity. The vis− cosity η(T) and the thermal conductivity k(T) are de− termined from Equations (9) and (10), respectively. We assume the following conditions for temperature and velocity at the model boundary: , (20) , (21) , (22) (23) . The principal problem is to find the solution to the problem (17)–(23), and hence to determine the velocity u = u(x), the pressure p = p(x), the temperature T = T(x), and hence the viscosity in the model domain Ω, when temperature T4 and heat flow φ are known at boundary (measured data). In this study measured data are obtained numerically by computing thermal heat flow and tem− perature at boundary Γ4 (see Supplementary Material S3). In addition to the principal problem, we state an auxiliary problem to find the solution to Equations (17)–(19) with the following boundary conditions: , (24) , (25) , (26) . (27) The auxiliary problem is a direct problem compared to the problem (17)–(23), which is an inverse problem, because the thermal condition is overdetermined at Γ4 and underdetermined at Γ2 .We note that the conditions at Γ1 and Γ3 are the same in the direct and inverse problems; temperature T2 is known at Γ2 and no heat flow is prescribed at Γ4 in the direct problem compared to the inverse problem. Consider “guess” temperature T2 = ξ at model boundary Γ2. Solving the auxiliary problem (17)–(19) and (24)–(27), we can determine the heat flow at model boundary Γ4 and compare it to the heat flow φ (the known synthetic data) at the same boundary. The fol− lowing cost functional for admissible functions ξ de− termined at T2 is to be then assessed: , (28) where k(Tξ) ∂T ξ / ∂n is the heat flow at Γ4 correspond− ing to temperature T2 = ξ at Γ2; and T ξ is the tempera− ture determined by solving the auxiliary problem. Therefore, we reduce the inverse problem to a minimiza− tion of the functional or to a variation of the function ξ at Γ2, so that, the model heat flow at Γ4 becomes closer TSEPELV ET AL. 8 FIGURE 5. Thickness of the lava crust with distance from the lava injection depending on the effective emissivity and the power−law exponent M (case studies m4.1, m4.3, and m4.6). See Table 1 for the case studies. Γ 4 : T = T 4 , − k ∇T ,n = ϕ, u,n = 0, σ n − n,n n = 0σ J (ξ ) = k (T ξ ) ∂T ξ ∂n − ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ 2 d Γ4 ∫ Γϕ Γ 1 : T = T 1 , u = u 1 Γ 2 : T = T 2 , u = 0 Γ 3 : T = T 3 , σ n = 0, p = 0 Γ 4 : T = T 4 , u,n = 0, σ n − n,n n = 0σ Γ 1 : T = T 1 , u = u 1 Γ 2 : u = 0 Γ 3 : T = T 3 , σ n = 0, p = 0 ∇p = ∇ ⋅ η (T ) ∇u + ∇uT( )( ) + Ra T e2 ∇ ⋅ u = 0 ∇ ⋅ κ (T ) ∇T( ) = u,∇T to the “observations” (known synthetic data) φ at Γ4. 6. SOLUTION METHOD Following Korotkii et al. [2016], we minimize the cost functional (28) using the Polak−Ribière conju− gate−gradient method [Polak and Ribière, 1969]. The gradient of the cost functional , (29) can be found as the solution (z, w, q) to the adjoint problem [see Appendix B in Korotkii et al., 2016] for the derivation of the adjoint problem): , (30) , (31) (32) , , (33) , (34) , (35) (36) , where σ ≡ ∇w+∇wT; the square brackets [A, B]=∑ aij bij denote the convolution of two m × m matrices A = (aij) and B = (bij); and sign ʹ means the derivation. The so− lution is a triplet (z, w, q) of quasi−temperature (z), quasi−velocity (w), and quasi−pressure q. The solution of the minimization problem is reduced to solutions of series of well−posed (direct and adjoint) problems. The algorithm for solving the data assimilation problem is based on solving iteratively the direct (aux− iliary) and adjoint problems and on the assessment J(ξ (i)), where i is the iteration number [Korotkii et al., 2016]. The numerical methods employed are described in sect. 3. An average computational time for 20 itera− tions was about 150 min on a single node; this included the time required for minimization of the cost func− tional by the conjugate gradient method, and the time to solve the auxiliary and adjoint problems. The performance of the algorithm is evaluated in terms of the number of iterations n required to achieve a prescribed relative reduction of the cost functional (Figure 6). The cost functional reaches a plateau after about 15 iterations. This behaviour is likely associated with errors in the model solution (see Section 8 for more detail). The quality of the gradient of the cost functional with respect to the control vari− able have been verified using the χ −test by Navon et al. [1992] (see Supplementary Material S4). 7. RESULTS OF THE ASSIMILATION In this section, two models of lava advancing down the slope on different topographic surfaces are analysed. The following values of the model param− eters are used h = 2 m, λ =5.5 W m-2 K-1, Ta = 300 K, ε = 0.95, ĸ* = 10 -6 m2 S-1, η* = 10 6 Pa s, and η0 = 1000. Initially we consider the model geometry presented in Figure 1. At Γ1 we prescribe the temperature T1(x1, x2) = 1320 K and the velocity u1(x2) = ϑ1U(x2)n1, x2 ∈ [x2 A, x2 E ], where n1=(√2/2, -√2/2), ϑ1 is a constant, and U(x2) is the parabola passing through the following three points: U(x2 E ) = 1, U(x2 A ) = 0, and U(0.5(x2 A + x2 E )) = 0.725. The cost functional is reduced to about 10−5 after about 15 iterations (Fig− ure 6). The number of iterations to get a given accu− racy in reduction of the cost functional depends on the initial “guess” temperature at Γ2. The closer the guess temperature is to the target temperature, the 9 LAVA CRUST DEVELOPMENT ~ m i,j=1 FIGURE 6. Relative reduction of the cost functional J with the number of iterations at three dimensionless rates of lava injection: ϑ = 10 (solid line), ϑ = 20 (dotted line), and ϑ = 40 (dashed line). Γ 4 : z = 2 k (T ξ ) ∂T ξ ∂n − ϕ ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ , w,n = 0, n − n,n n = 0σ~ σ~ ∇J (ξ ) = k (T ξ ) ∂ z ∂n ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ Γ2 ∇ ⋅ η (T ξ ) ∇w + ∇wT( )( ) = ∇q + z∇Tξ ∇ ⋅ w = 0 Γ 1 : z = 0, w = 0 Γ 2 : z = 0, w = 0 Γ 3 : z = 0, σ n = 0, q = 0~ ∇ ⋅ κ (T ξ )∇z( ) + uξ ,∇z + Ra e2 , w = = ʹκ (T ξ ) ∇T ξ ,∇z + + ʹη (T ξ ) (∇w + ∇wT ),∇u ξ ⎡ ⎣ ⎤ ⎦ fewer the number of iterations is needed. The target temperature at Γ2 is obtained from the solution of the direct problem (1)−(16). Figure 7 presents the target lava temperature, flow velocity, and viscosity as well as the residuals for vary− ing rates of the lava injection at the left boundary of the model domain. Here a residual is defined as the difference between the target physical parameter (e.g., temperature, velocity, or viscosity) and that obtained by the optimiza− tion (inverse) model. The results of this modelling show that at smaller rates of lava injection the residuals be− come small after 18 iterations within the almost entire model domain. At lower injection rates, the temperature, viscosity and velocity of a lava flow are well determined from the known thermal data on its surface with the at− mosphere (Figure 7a, b). Because of the topography over which the lava flows in the model case study, the crust thickens with the flow, reaches its maximal value in the middle of the channel, and diminishes toward the chan− nel’s end (the boundary Γ3). At higher rates, the recov− ery of the physical parameters worsens (Figure 7c, d); it can be explained by an advection−dominated lava flow with of a thin crust development (see Figure S1a), when the temperature of the lava flow approaches the temper− ature of the injected lava, and the thermal surface data are not properly assimilated into the lava flow. We consider a lava flow on an artificial volcanic to− pography in the second case study (Figure 8a). The lava flow from a vent on the topography has been modelled in a three−dimensional rectangular domain by Tsepelev et al. [2016]. As the optimization model is assumed to be in a steady state (that is, physical parameters of the model do not change with time), we select a portion of the lava flow on a two−dimensional profile along the lava flow; the profile is presented in the lower panel and its loca− tion in the upper panel of Figure 8a, where the selected portion of the lava (the model domain) is marked by red. In the model the lava thickness at the left and right TSEPELV ET AL. 10 FIGURE 7. Reconstruction of dimensionless viscosity, temperature, and velocity of lava flow in the model domain Ω. (a) and (c): Target functions at dimensionless injection rates ϑ = 10 and ϑ = 40, respectively; (b) and (d): the relevant residuals after the 1st iteration (the central panels) and several iterations (the right panels). 1 600 800200 400 1000 Viscosity 1 600 800200 400 1000 Viscosity 4.03.22.8 3.6 4.4 Temperature 4.03.22.8 3.6 4.4 Temperature 400 20 Velocity 400 20 Velocity 20 0 10 30 20 0 10 30 Velocity residual Velocity residual Viscosity residual 800 600 200 400 1000 0 800 600 200 400 1000 0 Viscosity residual Temperature residual 1 2 0 Temperature residual 1 2 0 a b c d 1st iteration 1st iteration 18th iteration 21st iteration sides are 14 m and 31 m, respectively, and the length of the model lava flow is 515 m. A lava is injected from the left boundary of the model domain. We solve the opti− mization problem in the selected domain to determine the lava’s viscosity, temperature, and flow rate based on the known thermal data on its surface with the atmo− sphere. Figure 8b presents the target viscosity, tempera− ture and flow rate, and Figure 8c their residuals. We see that after about 30 iterations the lava’s physical param− eters are recovered well enough from the surface thermal data. In this case study, a thin crust developing at the left end of the model domain becomes thicker toward its right end, and the flow velocity drops by a factor of about 3 with the lava advancement. 8. DISCUSSION Under relatively steady eruption conditions, a vis− cous lava flow cools at the interface with the atmo− sphere and develops a crust. A nonlinear heat flow at the interface, lava cooling, and crystallization of the uppermost layer of the moving melt results in the crust insulating the lava flow interior. The crust preserves the lava against rapid cooling and permits the lava flow to extend to substantial distances. If the lava supply ceases and the interior of the lava flow cools (so called “vol− ume−limited” flow regime) or a lava flow cools to such a degree that it can no longer moves even at a contin− uous slow supply of lava (“cooling−limited” flow regime), the lava will stop its further advance [Harris et al., 2007; Harris and Rowland, 2009; Rhéty et al., 2017]. In this paper, by solving direct problems we have studied how the lava crust forms at different conditions, namely, the changing effusion rates, and various pa− rameters of nonlinear heat flow and lava viscosity. Nu− merical experiments show that lava cools faster and its crust becomes thicker in the case of the nonlinear heat flow compared to the case of linear heat flow at the in− terface between lava and the atmosphere. At lower ef− 11 LAVA CRUST DEVELOPMENT FIGURE 8. Reconstruction of dimensionless viscosity, temperature, and velocity of lava flow on a synthetic topography. (a): A relief map (6000 m × 4000 m) of three−dimensional lava flow pattern, view from the top (upper panel) and the cross−section AB (lower panel) along the line indicated in the upper panel; the flow was computed by Tsepelev et al. [2016]; the red area marks the model domain in the inverse problem; (b): Target functions at dimensionless rates of lava injection ϑ = 10; (c): the relevant residuals after the 1st iteration (the left panel) and after 31st iteration (the right panel). 20 1 Velocity residual Viscosity residual Temperature residual 0.2 0.4 0.80 0.6 1st iteration 31st iteration c 6000 m 6000 m 4 0 0 0 m 3 0 0 m a 1 600400200 1000800 0 600400200 1000800 Viscosity 4.0 4.23.83.6 4.4 Temperature 4 60 2 Velocity b vent A A B B fusion rates of lava into the channel (the model do− main), cooling leads to development of a thick crust; the higher rates result in rapid lava advection and in a delay of a thick crustal formation . It has been also shown that a lava crust becomes thicker with either a higher coefficient of conductive heat transfer or higher lava viscosity. At a higher lava viscosity the flow be− comes slower, the lava cools and the crust develops more efficiently than in the case of a rapid lava flow due to its lower viscosity. Temperature at the lava sur− face with the atmosphere decreases and the lava crust thickens with the increasing its effective emissivity. If the front of lava flow covered by a crust prevents lava to flow easily, numerical experiments illustrate that the lava crust beneath the immobile region (the stuck area at the front of lava flow) becomes much thicker. If direct problems present models dealing mostly with understanding of phenomena, inverse models are related to finding either initial conditions or boundary conditions or model parameters based on geological, geophysical, geodetic and geochemical observations and measurements. Inverse problems in geodynamical models are mostly related to optimization of residuals between observed data and those obtained from rele− vant mathematical models. Earth orbiting radiometers can measure spectral ra− diance at a lava surface to be converted into thermal anomalies; lava temperature and heat flow can be then inferred from the detected anomalies. Results of nu− merical modelling illustrate that the lava’s physical pa− rameters, including the lava crust thickness, can be recovered from the surface thermal (known/measured) data well enough after a few dozens of iterations be− tween the adjoint and auxiliary problems. However, a spatial resolution of many satellites is coarse enough to not allow for high−resolution monitoring and precise measurements, and this gives a rise to uncertainties in thermal measurement as well as in the inferred param− eters [e.g. Zakšek et al., 2015]. Hence, if the measured temperature and heat flow data are biased, this infor− mation can be improperly assimilated into the lava flow models. And vice versa, if surface temperature and heat flow data are of high resolution and radiometric accu− racy, the temperature and velocity in the lava’s interior can be determined properly from measured data using the data assimilation approach. The data assimilation approach becomes be important in studies of natural lava flows, especially in the cases of thick lava flow. Synthetic Aperture Radar (SAR) satellite observations on lava thickness, volume, and flow extent [e.g., Kubanek et al., 2015], together with thermal measure− ment at the lava surface, could facilitate research and data−driven modelling of lava flow. The studied models describe steady−state flow, al− though lava flows are non−stationary. Meanwhile, as measurements on absolute temperature are discrete in time in most cases (e.g., depending on the location of Landsat satellites), a problem of non−steady state flow can be reduced to a series of steady−state flow problems with varying model domain and boundary conditions assimilating thermal data available at the discrete−in− time measurements. Also, airborne and space measure− ments of absolute temperature at the lava interface with the atmosphere, being almost instantaneous compared to the duration of the lava flow, allow searching for thermal conditions at the bottom of the lava flow using the cost functional (28). Once the boundary conditions at the lava bottom are determined, the steady−state problem can be replaced by a non−steady state problem, and the lava flow can be modelled forward in time to determine its extent, lava’s temperature and flow rate as well as backward in time using variational [Korotkii et al., 2016] or quasi−reversibility [Ismail−Zadeh et al., 2007] methods to search for the initial temperature of the lava flow and for the evolution of the effusion rate. 8.1 DATA ASSIMILATION LIMITATIONS The mathematical problem (17) − (23) does not possess the stability property of its solution [Korotkii and Kovtunov, 2011]. This means that if the problem is solved numerically, small errors in measured input data, e.g., the absolute temperature estimated from remote sensors measurements of thermal (electro− magnetic) radiation, and in the inferred heat flow at the interface of lava with the atmosphere, as well as computational errors, can result in significant errors in determination of temperature, viscosity, and flow velocity. Therefore, to solve the problem special numerical methods should be employed. Particularly, the prob− lem (17)−(23) is replaced by the optimization of the cost functional (28) which, in its turn, lead to the so− lution of the coupled auxiliary (17)−(19), (24)−(27) and adjoint (30)−(36) problems to find the gradient (29) of the cost functional. The Polak−Ribière stable iterative conjugate gradient method was used to min− imize the optimization problem. However, the itera− tive method provides the rapid and accurate solution TSEPELV ET AL. 12 to the adjoint and auxiliary problems only when the following condition is satisfied [Korotkii and Kov− tunov, 2011]: , (37) where C1 and C2 are positive constants. An iterative convergence slows at high Rayleigh numbers, and the iterations diverge at the Rayleigh numbers greater than 106. As the injection rate of lava into the model do− main increases, the minimization process slows down (Figure 6). A rapid injection results in advection of high temperature with flow and, hence, in a decrease of lava viscosity, and in slowing convergence of iterations and their divergence. 8.2 UNCERTAINTIES AND MEASUREMENT ERRORS The results of numerical modelling of the problem (17)−(23) show that the optimization works effec− tively: the temperature, viscosity and velocity resid− uals are very small already after a few dozens of iterations within the almost entire model domain. Normally any measurements (observations) are pol− luted by errors; for example, in the case of lava tem− perature measurements at the surface, the accuracy of the calibration curve of remote sensors and the noise of the sensors can influence measurements and con− tribute to the measurement errors [Short and Stuart, 1983]. Korotkii et al. [2016] performed numerical ex− periments introducing a noise on the “measured” heat flow data φδ(·) = φ(·) + δχ(·) and analysed the sensi− tivity of the model to the noise. Here, δ is the magni− tude of the noise, and χ(·) is the function generating numbers that are uniformly distributed over the in− terval [−1, 1]. They showed that the temperature and velocity residuals get larger with increase of the noise of the input data, but are still acceptable at the level of noise (or errors in measurements) of up to 10%. It should be noted that the error analysis of subpixel temperature retrieval from satellite infrared data showed that errors in measurements of the radiant heat flux are within about 5% to 10%, and can be re− duced [Lombardo et al., 2012]. In general, the cost functional related to the inverse (optimization) problem with noisy measurements at the surface of the lava flow with the atmosphere can be written in the following form: (38) . Substituting the solution ξ* to Equation (28) into Equation (38), we obtain Jδ(ξ*) ∼ δ 2, because the first and second terms of the right hand−side of Equation (38) turn to zero at ξ*. Therefore, at the minimization of the functional Jδ, the functional will be approaching the non−vanishing value equal to the square of the noise’s magnitude (δ2). In the case of synthetic thermal data prescribed at the upper model boundary (instead of real measurements), the ‘plateau’ in the curve illustrating the minimization of the functional (Figure 6) is likely to be associated with numerical errors. Forcing the solution to the functional to attain zero may lead to an unstable (or erroneous) solution. Moreover, some a priory infor− mation assists in solving the problem. For example, the temperature inside the model domain cannot be higher than that at the left boundary of the model domain (where the lava is injected into the model domain). This can serve as a control parameter for the computed tem− perature in the minimization problem. Rather accurate reconstructions of the model temper− ature, viscosity, and flow velocity in this study rely on the chosen method for minimization of the cost functional (28). In the general case, a Tikhonov regularization term should be introduced in the cost functional as: (39) , where α is a small positive regularization parameter, ᴧ > 0 is the operator accounting for a priori informa− tion on the problem’s solution (e.g., its monotony prop− erty, maximum and minimum values, and the total variation diminishing), and ξ0 is a priori known func− tion close to the solution of the problem. The introduc− tion of the regularization term in the cost functional makes the minimization problem more stable and less dependent on measurement errors. For a suitable regu− larization parameter α = α(δ), the minimum of the reg− ularized cost functional will tend to the minimum of the functional (28) at δ → 0 [Tikhonov and Arsenin, 1977]. The choice of the regularization parameter is a challenging issue as it depends on several factors, e.g., 13 LAVA CRUST DEVELOPMENT Ra η−1 C 1 ϕ 2 d Γ Γ4 ∫ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ 1/2 + C 2 ξ 2 d Γ Γ2 ∫ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ 1/2⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ < 1 J α (ξ ) = k (T ξ ) ∂T ξ ∂n − ϕ δ ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ 2 d Γ4 ∫ Γ+ + α Λ − 0( ) d Γ Γ2 ∫ ξ ξ 2δ k (T ξ ) ∂T ξ ∂n − ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ χd Γ4 ∫ Γ+ 2 2 δϕ χ+ J δ (ξ ) = k ( T ξ ) ∂T ξ ∂n −ϕ +δ χ ⎛ ⎝ ⎜⎜ ⎞ ⎠ ⎟⎟ 2 d Γ4 ∫ Γ = J (ξ ) + on errors of measured data [e.g., Kabanikhin, 2011]. Meanwhile, if there is a lack of the information on the solution of the problem, the better strategy is to mini− mize the functional (28) using a stable minimization method (e.g., the Polak−Ribière method). 9. CONCLUDING REMARKS Understanding the mechanisms influencing cooling and solidification of lava is essential in forecasting the lava flow advance. Numerical experiments of lava flows at different conditions, presented in this paper, provide an insight into lava cooling and its crust development and assist in bridging the gap to application for natu− ral lava flows. If the temperature and heat flow at the interface of lava with the atmosphere are of high reso− lution and radiometric accuracy, the temperature and velocity in the lava’s interior can be determined prop− erly from the measured data. Modelling of direct and inverse problems provides a knowledge of the thermal and dynamic characteristics of lava flow, and it be− comes important for lava flow hazard and disaster risk assessments [e.g., Wadge et al., 1994; Harris, 2015; Cut− ter et al., 2015]. Acknowledgements. We thank Augusto Neri, Donald Dingwell, and Frank Schilling for discussion related to nonlinear heat flow and lava rheology. We are grateful to Andrew Harris and Scott Rowland for very constructive comments, which improved the initial manuscript. This work is supported by the RSF grant 19- 17-00027. REFERENCES Boussinesq, J. (1903). Théorie Analytique de la Chaleur, vol. 2, Gauthier−Villars, Paris. Branca, S., De Beni, E. and Proietti, C. (2013). The large and destructive 1669 AD eruption at Etna vol− cano: reconstruction of the lava flow field evolu− tion and effusion rate trend, Bull. Volcanol., 75: 694, https://doi.org/10.1007/s00445−013−0694−5. Bryan, S.E., Peate, I.U., Peate, D.W., Self, S., Jerram, D.A., Mawby, M.R., Marsh, J.S. and Miller, J.A. (2010). The largest volcanic eruptions on Earth, Earth−Sci. Rev., 102 (3–4), 207−229. Calvari, S., Spampinato, L., Lodato, L., Harris, A. J. L., Patrick, M. R., Dehn, J., Burton, M. R. and An− dronico, D. (2005). Chronology and complex vol− canic processes during the 2002–2003 flank erup− tion at Stromboli volcano (Italy) reconstructed from direct observations and surveys with a handheld thermal camera, J. Geophys. Res., 110, B02201, doi:10.1029/2004JB003129. Castruccio, A. and Contreras, M.A. (2016). The influ− ence of effusion rate and rheology on lava flow dynamics and morphology: A case study from the 1971 and 1988–1990 eruptions at Villarrica and Lonquimay volcanoes, Southern Andes of Chile, J. Volcanol. Geother. Res., 327, 469−483. Chevrel, M.O., Platz, T., Hauber, E., Baratoux, D., Laval− lée, Y. and Dingwell, D.B. (2013). Lava flow rhe− ology: A comparison of morphological and petrological methods, Earth Planet. Sci. Lett., 384, 109−120. Chevrel, M.O., Guilbaud, M.N. and Siebe, C. (2016). The ∼AD 1250 effusive eruption of El Metate shield volcano (Michoacán, Mexico): magma source, crustal storage, eruptive dynamics, and lava rhe− ology, Bull. Volcanol., 78: 32, https://doi.org/10.1007/s00445−016−1020−9. Christiansen, R. (2001). The Quaternary and Pliocene Yellowstone Plateau volcanic field of Wyoming, Idaho, and Montana. US Geological Survey Pro− fessional Paper. G1−G145. Cordonnier, B., Lev, E. and Garel, F. (2015). Bench− marking lava−flow models. In: Harris, A.J.L., De Groeve, T., Garel, F. and Carn, S.A. (eds.) Detect− ing, Modelling and Responding to Effusive Erup− tions, Geological Society, London, Special Publications, 426, 425−445. Corsaro, R. A., Métrich, N., Allard, P., Andronico, D., Miraglia, L. and Fourmentraux, C. (2009). The 1974 flank eruption of Mount Etna: An archetype for deep dike−fed eruptions at basaltic volcanoes and a milestone in Etna’s recent history, J. Geo− phys. Res., 114, B07204, doi:10.1029/2008JB006013. Costa, A. and Macedonio, G. (2005). Computational modeling of lava flows: a review. In: Manga, M. and Ventura G (eds.) Kinematics and Dynamics of Lava Flows, Geol. Soc. Am. Spec. Publ. 396, Lon− don, pp. 209−218. Costa, A., Caricchi, L. and Bagdassarov, N. (2009). A model for the rheology of particle−bearing sus− pensions and partially molten tocks, Geochem. Geophys. Geosys., 10, Q03010, doi TSEPELV ET AL. 14 10.1029/2008Gc002138. Cutter, S., Ismail−Zadeh, A., Alcántara−Ayala, I. et al. (2015). Pool knowledge to stem losses from disas− ters, Nature, 522, 277−279. Dietterich, H.R., Downs, D.T., Stelten, M.E. and Zahran, H. (2018). Reconstructing lava flow emplacement histories with rheological and morphological analyses: the Harrat Rahat volcanic field, King− dom of Saudi Arabia, Bull. Volcanol., 80: 85. https://doi.org/10.1007/s00445−018−1259−4. Dragoni, M.A. (1989). A dynamical model of lava flows cooling by radiation, Bull. Volcanol., 51, 88−95. Fujita, E. and Nagai, M. (2015). LavaSIM: its physical basis and applicability. In: Harris, A.J.L., De Groeve, T., Garel, F. and Carn, S.A. (eds.) Detect− ing, Modelling and Responding to Effusive Erup− tions, Geological Society, London, Special Publications, 426, 375−386. Giordano, D. and Dingwell, D. B. (2003). Viscosity of hydrous Etna basalt: implications for Plinian−style basaltic eruptions, Bull. Volcanol., 65, 8−14. Griffiths, R.W. (2000). The dynamics of lava flows, Ann. Rev. Fluid. Mech., 32, 477−518. Harris, A.J.L. (2015). Basaltic lava flow hazard. In: Pa− pale, P. (ed.) Volcanic Hazards, Risks and Disas− ters, Elsevier, Amsterdam, pp. 17−47. Harris, A.J.L. and Allen, J.S. III (2008). One−, two− and three−phase viscosity treatments for basaltic lava flows, J. Geophys. Res., 113, B09212. Harris, A. J. L. and Rowland, S. K. (2001). FLOWGO: a kinematic thermo−rheological model for lava flowing in a channel, Bull. Volcanol., 63, 20−44. Harris, A. J. L., and Rowland, S. K. (2009). Effusion rate controls on lava flow length and the role of heat loss: A review. In: Thordarson, T., Self, S., Larsen, G., Rowland, S. K. and Höskuldsson, Á. (eds.) Studies in Volcanology: The Legacy of George Walker, Geol. Soc. London, Spec. Publ. IAVCEI, 2, pp. 33–51. Harris, A.J.L., Dehn, J. and Calvari, S. (2007). Lava ef− fusion rate definition and measurement: a review, Bull. Volcanol., 70, 1−22. Harris, A.J.L., Flynn, L.P., Keszthelyi, L., Mouginis− Mark, P.J., Rowland, S.K. and Resing, J.A. (1998). Calculation of lava effusion rates from Landsat TM data, Bull. Volcanol., 60, 52−71. Herault, A., Bilotta, G., Vicari, A., Rustico, E. and Del Negro, C. (2011). Numerical simulation of lava flow using a GPU SPH model, Ann. Geophys., 54(5), http://doi.org/10.4401/ag−5343. Hidaka, M., Goto, A., Umino, S. and Fujita, E. (2005). VTFS project: Development of the lava flow sim− ulation code LavaSIM with a model for three−di− mensional convection, spreading, and solidification, Geochem. Geophys. Geosyst., 6, Q07008, doi:10.1029/2004GC000869. Ismail−Zadeh, A. and Tackley, P. (2010). Computational Methods for Geodynamics. Cambridge University Press, Cambridge. Ismail−Zadeh, A., Korotkii, A. and Tsepelev, I. (2016). Data−Driven Numerical Modelling in Geodynam− ics: Methods and Applications. Springer, Heidel− berg. Ismail−Zadeh, A., Korotkii, A., Schubert, G. and Tse− pelev, I. (2007). Quasi−reversibility method for data assimilation in models of mantle dynamics, Geophys. J. Int., 170, 1381−1398. Jeffrey, D. and Acrivos, A. (1976). The rheological prop− erties of suspensions of rigid particles, AIChE J., 22, 417−432. Kabanikhin, S.I. (2011) Inverse and Ill−Posed Problems: Theory and Applications, De Gruyter, Berlin. Kelfoun, K. and Druitt, T. (2005). Numerical modeling of the emplacement of Socompa Rock Avalanche, Chile, J. Geophys. Res., 110, B12202, http://doi.org/10.1029/2005JB003758. Kelfoun, K. and Vargas, S. V. (2015). VolcFlow capabil− ities and potential development for the simulation of lava flows. In: Harris, A. J. L., De Groeve, T., Farel, F. and Carn, S. A. (eds.) Detecting, Modelling and Responding to Effusive Eruptions, Geol. Soc. London, Spec. Publ., 426, doi:10.1144/SP426.8. Kilburn, C.R.J. (2000). Lava flow and flow fields. In: Sigurdsson H. (ed.) Encyclopedia of Volcanoes, Academic Press, San Diego, pp. 291−305. Kirsch, A. (1996). An Introduction to the Mathematical Theory of Inverse Problems, Springer−Verlag, New York. Korotkii, A.I. and Kovtunov D.A. (2011). Optimal boundary control of a system describing thermal convection, Proc. Steklov Inst. Math., 272, suppl. 1, 74−100. Korotkii, A., Kovtunov, D., Ismail−Zadeh, A., Tsepelev, I. and Melnik, O. (2016). Quantitative reconstruc− tion of thermal and dynamic characteristics of lava from surface thermal measurements, Geo− phys. J. Int., 205, 1767−1779. Kubanek, J., Richardson, J.A., Charbonnier, S.J. and 15 LAVA CRUST DEVELOPMENT Connor, L.J. (2015). Lava flow mapping and vol− ume calculations for the 2012–2013 Tolbachik, Kamchatka, fissure eruption using bistatic Tan− DEM−X InSAR, Bull. Volcanol., 77(12), doi: 10.1007/s00445−015−0989−9. Lejeune, A. and Richet, P. (1995). Rheology of crystal− bearing silicate melts: An experimental study at high viscosity, J. Geophys. Res., 100, 4215−4229. Lombardo, V., Musacchio, M. and Buongiorno, M.F. (2012). Error analysis of subpixel lava tempera− ture measurements using infrared remotely sensed data, Geophys. J. Int., 191, 112−125. Macedonio, G., Costa, A. and Longo, A. (2005). A com− puter model for volcanic ash fallout and assess− ment of subsequent hazard, Computers & Geosciences, 31, 837–845. Marsh, B.D. (1981). On the crystallinity, probability of occurrence, and rheology of lava and magma, Contrib. Mineral. Petrol., 78, 85−98. Navon, I. M., Zou, X., Derber, J. and Sela, J. (1992). Variational data assimilation with an adiabatic version of the NMC spectral model, Monthly Weather Rev., 120(7), 1433−1446. Neri, A. (1998). A local heat transfer analysis of lava cooling in the atmosphere: application to thermal diffusion−dominated lava flows, J. Volcanol. Geother. Res., 81, 215−243 Patankar, S.V. and Spalding, D.B. (1972). A calculation procedure for heat and mass transfer in three−di− mensional parabolic flows, Int. J. Heat Mass Transfer, 15, 1787−1806. Patrick, M.R., Orr, T., Fisher, G., Trusdell, F. and Kauahikaua, J. (2017). Thermal mapping of a pa− hoehoe lava flow, Kilauea Volcano, J. Volcanol. Geother. Res., 332, 71−87. Pedersen, G.B.M., Belart, J.M.C., Magnússon, E., Vil− mundardóttir, O.K., Kizel, F., Sigurmundsson, F. S., et al. (2018). Hekla volcano, Iceland, in the 20th century: Lava volumes, production rates, and ef− fusion rates, Geophys. Res. Lett., 45, 1805−1813. Pinton, A., Giordano, G., Speranza, F., and Þórðarson, Þ. (2018). Paleomagnetism of Holocene lava flows from the Reykjanes Peninsula and the Tungnaá lava sequence (Iceland): implications for flow cor− relation and ages, Bull. Volcanol., 80: 10, https://doi.org/10.1007/s00445−017−1187−8. Polak, E. and Ribière, G. (1969). Note on the conver− gence of methods of conjugate directions. Revue Francaise d'Informatique et de Recherche Opera− tionnelle, 3(16), 35−43. Poland, M.P., Orr, T.R., Kauahikaua, J.P et al. (2016). The 2014−2015 Pāhoa lava flow crisis at Kīlauea Vol− cano, Hawai‘i: Disaster avoided and lessons learned, GSA Today, 26(2), 4−10. Rhéty, M., Harris, A., Villeneuve, N., Gurioli, L., Médard, E., Chevrel, E., and Bachélery P. (2017). A com− parison of cooling−limited and volume−limited flow systems: Examples from channels in the Piton de la Fournaise April 2007 lava−flow field, Geochem. Geophys. Geosyst., 18(9), 3270−3291. Rowland, S.K. and Walker, G.P.L. (1990). Pahoehoe and aa in Hawaii: volumetric flow rate controls the lava structure, Bull. Volcanol., 52, 615−628. Rumpf, M.E., Lev, E. and Wysocki, R. (2018). The influ− ence of topographic roughness on lava flow em− placement, Bull. Volcanol., 80: 63. https://doi.org/10.1007/s00445−018−1238−9. Self, S., Jay, A.E., Widdowson, M. and Keszthelyi, L.P. (2008). Correlation of the Deccan and Rajah− mundry Trap lavas: Are these the longest and largest lava flows on Earth? J. Volcanol. Geother. Res., 172, 3−19. Short, N.M. and Stuart, L.M. (1983). The Heat Capacity Mapping Mission (HCMM) Anthology. Scientific and Technical Information Branch, National Aero− nautics & Space Administration, Washington, D.C. Soldati, A., Harris, A.J.L., Gurioli, L., Villeneuve, N., Rhéty, M., Gomez, F. and Whittington, A. (2018). Textural, thermal, and topographic constraints on lava flow system structure: the December 2010 eruption of Piton de la Fournaise, Bull. Volcanol., 80: 74. https://doi.org/10.1007/s00445−018− 1246−9. Sweby, P.K. (1984). High resolution schemes using flux limiters for hyperbolic conservation laws, J. Numer. Anal., 21, 995−1011. Tikhonov, A.N. and Arsenin, V.Y. (1977). Solution of Ill− Posed Problems, Winston, Washington, D.C. Tsepelev, I., Ismail−Zadeh, A., Melnik, O. and Korotkii, A. (2016). Numerical modelling of fluid flow with rafts: An application to lava flows, J. Geodyn., 97, 31−41. Wadge, G., Young, P.A.V. and McKendrick, I.J. (1994). Mapping lava flow hazards using computer sim− ulation, J. Geophys. Res., 99, 489−504. Walker, G.P.L. (1973). Lengths of lava flows, Phil. Trans. R. Soc. London, 274, 107−118. Wright, R., Garbeil H. and Davies, A.G., (2010). Cooling TSEPELV ET AL. 16 rate of some active lavas determined using an or− bital imaging spectrometer, J. Geophys. Res., 115, B06205, doi:10.1029/2009JB006536. Wright, T.L. and Okamura, R.T. (1977). Cooling and crystallization of tholeiitic basalt, 1965 Makaop− uhi lava lake, Hawaii, US Geological Survey Pro− fessional Paper 1004, 78 p. Zakšek, K., Hort, M. and Lorenz, E. (2015). Satellite and ground based thermal observation of the 2014 ef− fusive eruption at Stromboli volcano, Remote Sens., 7, 17190−17211. *CORRESPONDING AUTHOR: Alik ISMAIL-ZADEH Karlsruhe Institute of Technology, Institute of Applied Geosciences, Adenauerring 20b, Karlsruhe 76131, Germany email: alik.ismail-zadeh@kit.edu © 2019 the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved 17 LAVA CRUST DEVELOPMENT << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.3 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.1000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize false /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile (None) /AlwaysEmbed [ true /AndaleMono /Apple-Chancery /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /CapitalsRegular /Charcoal /Chicago /ComicSansMS /ComicSansMS-Bold /Courier /Courier-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /GadgetRegular /Geneva /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Helvetica /Helvetica-Bold /HelveticaInserat-Roman /HoeflerText-Black /HoeflerText-BlackItalic /HoeflerText-Italic /HoeflerText-Ornaments /HoeflerText-Regular /Impact /Monaco /NewYork /Palatino-Bold /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman /SandRegular /Skia-Regular /Symbol /TechnoRegular /TextileRegular /Times-Bold /Times-BoldItalic /Times-Italic /Times-Roman /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Webdings ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 150 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.10000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 150 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.10000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.08250 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName (http://www.color.org) /PDFXTrapped /Unknown /CreateJDFFile false /SyntheticBoldness 1.000000 /Description << /ENU (Use these settings to create PDF documents with higher image resolution for high quality pre-press printing. The PDF documents can be opened with Acrobat and Reader 5.0 and later. These settings require font embedding.) /JPN /FRA /DEU /PTB /DAN /NLD /ESP /SUO /NOR /SVE /KOR /CHS /CHT /ITA >> >> setdistillerparams << /HWResolution [2400 2400] /PageSize [595.000 842.000] >> setpagedevice