“CRUST DEVELOPMENT INFERRED FROM NUMERICAL MODELS OF LAVA FLOW AND ITS SURFACE THERMAL MEASUREMENTS„

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.

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

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 Ω ⊂ R 2 (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 (x 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: , . (3) Here x = (x 1 , x 2 ) ∈ Ω are the Cartesian coordinates; u = (u 1 (x),u 2 (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 temperature-dependent heat conductivity; g is the acceleration due to gravity; c * (= 1000 J kg -1 K -1 ) is the specific heat capacity; e 2 = (0,1) is the unit vector; ∇, T , and 〈•,•〉 denote the gradient vector, the transposed matrix, and the scalar product of vectors, respectively.
Finally, we use in the modelling the following lava viscosity restricting its exponential growth with tem− perature by the prescribed value η 0 (Figure 2): where η * (=10 6 Pa s) is the typical lava viscosity, and η 0 takes values from 10 3 to 10 6 in the modelling.
In addition to these equations, we consider the fol− lowing relationships for the thermal conductivity and the density: where, T = 1473 K; r m (= 2750 kg m -3 ) and r c (= 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 u 1 and the temperature T 1 are prescribed at the No slip condition and zero heat flux (no heat lost to the base via conduction) or temperature T 2 are pre− scribed at the lower boundary Γ 2 : where n is the outward unit normal vector at a point of the boundary.We assume at the right boundary Γ 3 : where σ = η(∇u+∇u T ) 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: ) 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, T a , η * , and k * .

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 10 4 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.
u,n = 0, σ n − σ n,n n = 0, 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.

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 |u 1 | (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 Supplementary 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 10 6 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.
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 10 2 -10 3 Pa s for basaltic lava flows, and increases up to 10 9 -10 11 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 Material S2, Figure S2).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.

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) where and the Rayleigh number is determined as 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: , ( 21) 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 T 4 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: . ( 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 T 2 is known at Γ 2 and no heat flow is prescribed at Γ 4 in the direct problem compared to the inverse problem.
Consider "guess" temperature T 2 = ξ 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 T 2 is to be then assessed: where k(T ξ ) ∂T ξ / ∂n is the heat flow at Γ 4 correspond− ing to temperature T 2 = ξ 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.
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.

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) , (36) , where σ ≡ ∇w+∇w T ; the square brackets [A, B]=∑ a ij b ij denote the convolution of two m × m matrices A = (a ij ) and B = (b ij ); 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).
Initially we consider the model geometry presented in Figure 1.At Γ 1 we prescribe the temperature T 1 (x 1 , x 2 ) = 1320 K and the velocity u where n 1 =(√2/2, -√2/2), ϑ 1 is a constant, and U(x 2 ) is the parabola passing through the following three points: and Γ 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.

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

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 to the adjoint and auxiliary problems only when the following condition is satisfied [Korotkii and Kov− tunov, 2011]: , (37) where C 1 and C 2 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.

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

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].
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 T a (= 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〉 = a 1 (T 4 -1) + a 2 (T -1) M respectively, where Nu = l ~h / k * is the Nusselt number,

FIGURE 3 .
FIGURE 3. Thickness of the lava crust with distance from the lava injection into the model domain depending on: (a) the injection rate |u 1 |, (b) conductive heat transfer l ~, and (c) the typical lava viscosity η*.See Table 1 for case studies m1.1−m3.4.

FIGURE 6 .
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).

FIGURE 7 .
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 1 st iteration (the central panels) and several iterations (the right panels).
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.Figure8bpresents the target viscosity, tempera− ture and flow rate, and Figure8ctheir 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.

FIGURE 8 .
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 1 st iteration (the left panel) and after 31 st iteration (the right panel).

TABLE 1 .
Model parameters for case studies