# Quantitative reconstruction of thermal and dynamic characteristics of lava flow from surface thermal measurements

Alexander Korotkii, Dmitry Kovtunov, Alik Ismail-Zadeh, Igor Tsepelev, Oleg Melnik

► To cite this version:

Alexander Korotkii, Dmitry Kovtunov, Alik Ismail-Zadeh, Igor Tsepelev, Oleg Melnik. Quantitative reconstruction of thermal and dynamic characteristics of lava flow from surface thermal measurements. Geophysical Journal International, 2016, 205, pp.1767-1779. 10.1093/gji/ggw117 . insu-03581657

HAL Id: insu-03581657

https://insu.hal.science/insu-03581657

Submitted on 20 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantitative reconstruction of thermal and dynamic characteristics of lava flow from surface thermal measurements

Alexander Korotkii,<sup>1,2,3</sup> Dmitry Kovtunov,<sup>1,2</sup> Alik Ismail-Zadeh,<sup>1,4,5</sup> Igor Tsepelev<sup>1,2</sup> and Oleg Melnik<sup>1,6</sup>

Accepted 2016 March 29. Received 2016 January 31; in original form 2015 August 26

## SUMMARY

We study a model of lava flow to determine its thermal and dynamic characteristics from thermal measurements of the lava at its surface. Mathematically this problem is reduced to solving an inverse boundary problem. Namely, using known conditions at one part of the model boundary we determine the missing condition at the remaining part of the boundary. We develop a numerical approach to the mathematical problem in the case of steady-state flow. Assuming that the temperature and the heat flow are prescribed at the upper surface of the model domain, we determine the flow characteristics in the entire model domain using a variational (adjoint) method. We have performed computations of model examples and showed that in the case of smooth input data the lava temperature and the flow velocity can be reconstructed with a high accuracy. As expected, a noise imposed on the smooth input data results in a less accurate solution, but still acceptable below some noise level. Also we analyse the influence of optimization methods on the solution convergence rate. The proposed method for reconstruction of physical parameters of lava flows can also be applied to other problems in geophysical fluid flows.

**Key words:** Numerical solutions; Inverse theory; Numerical approximations and analysis; Effusive volcanism.

## 1 INTRODUCTION

Modern remote sensing technologies (e.g. air-borne or space-borne infrared sensors) allow for detecting the absolute temperature at the Earth's surface (e.g. Flynn et al. 2001; Fig. 1). The absolute temperature is a measure of the average kinetic energy of charged particles within matter in motions. Collisions of the particles result in changes of energy emitted as thermal (electromagnetic) radiation, which can be detected by remote sensors. The Stefan–Boltzmann law relates the total energy radiated per unit surface area of a body across all wavelengths per unit time to the fourth power of the absolute temperature of the body. Hence the absolute temperature can be determined from the measurements by remote sensors (e.g. Harris et al. 2004). The heat flow could be then inferred from the Stefan–Boltzmann law using the temperature.

Is it possible to use the surface thermal data so obtained to constrain the thermal and dynamic conditions beneath the surface? In this paper we propose a quantitative approach to reconstruct temperature and velocity in the steady-state lava flow. The knowledge

of the thermal and dynamic characteristics of lava is important, particularly, for lava flow hazard and risk assessments (Wadge et al. 1994) and disaster risk reduction (Cutter et al. 2015).

Many thermal/dynamic problems can be described by mathematical models, that is, by a set of partial differential equations, boundary and initial conditions defined in a specific domain. If a direct mathematical problem concerns an analysis of the effects of surface dynamic processes, an inverse problem intends to determine causes of the processes from their effects.

The problem of reconstruction of lava thermal and flow characteristics is considered here in the case when the temperature and the heat flow are known on the lava surface, but the lava temperature and velocity are unknown. The problem is reduced to determination of temperature and velocity as the solution to the model of steady-state flow of viscous fluid with prescribed conditions for velocity and temperature at the boundary $\Gamma = \partial\Omega$ of the model domain $\Omega$. At a part of the model boundary the conditions are abundant (e.g. both temperature and heat flow are known), and at another part of the boundary there is a lack of information on the temperature

<footnote n="1"><sup>1</sup>Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Moscow, Russia</footnote>

<footnote n="2"><sup>2</sup>Institute of Mathematics and Mechanics, Russian Academy of Sciences, Yekaterinburg, Russia</footnote>

<footnote n="3"><sup>3</sup>Institute of Mathematics and Computer Sciences, Ural Federal University, Yekaterinburg, Russia</footnote>

<footnote n="4"><sup>4</sup>Institute of Applied Geosciences, Karlsruhe Institute of Technology, Karlsruhe, Germany. E-mail: Alik.Ismail-Zadeh@kit.edu</footnote>

<footnote n="5"><sup>5</sup>Institut de Physique du Globe de Paris, Paris, France</footnote>

<footnote n="6"><sup>6</sup>Institute of Mechanics, Lomonosov Moscow State University, Moscow, Russia</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Surface temperature of lava can be extracted from the satellite measurements, for example, from LANDSAT 7 ETM+ thermal and infrared bands (the satellite image is courtesy of USGS).</figcaption></figure>

(because of no direct measurements at this part of the boundary). This mathematical problem is reduced in its turn to solving the inverse problem for determination of the temperature at the bottom of the lava and for subsequent search for the temperature and velocity of the lava.

For clarity of subsequent discussion, we introduce a few mathematical definitions used in the paper. A mathematical model for a geophysical problem has to be *well-posed* in the sense that it has to have the properties of existence, uniqueness and stability of a solution to the problem (e.g. Hadamard 1923; Tikhonov & Arsenin 1977; Kirsch 1996). Problems for which at least one of these properties does not hold are called *ill-posed*. If a problem lacks the property of stability then its solution is almost impossible to compute because numerical computations are polluted by unavoidable errors. If the solution of a problem does not depend continuously on the initial data, then, in general, the computed solution may have nothing to do with the true solution.

The problem of temperature and flow reconstruction in a lava from the temperature and heat flow data at the lava surface is an ill-posed because of a lack of the solution stability as a small perturbation of the conditions at the lava surface may lead to significant errors in the solution to the problem. To solve ill-posed problems, special methods, sometimes called data assimilation techniques, are required (e.g. Ismail-Zadeh & Tackley 2010). The techniques can be used to constrain the condition at the lower boundary of the lava from observations at the lava surface. Assimilation of data in this case can be defined as the incorporation of observations at relevant boundaries of a model domain in an explicit model to provide coupling among the physical fields (e.g. velocity, temperature). The basic principle of this assimilation is then to consider the condition at the lava bottom as a control variable and to optimize this condition in order to minimize the misfit between the observations at the lava surface and the model solution at the same surface.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Model geometry.</figcaption></figure>

In this study, we assimilate measured temperature and heat flux from the upper boundary into the interior of the model domain to constrain the thermal condition at the lower boundary. This type of data assimilation belongs to optimal boundary control problems (e.g. Zou *et al.* 1995). Thacker & Long (1988) suggested to investigate open-boundary control in ocean models from a data assimilation perspective, pointing out that open-boundary conditions are analogous to initial conditions and should be determined as part of fitting dynamics to data. Le Dimet (1988) formulated open-boundary control in a general mathematical framework. With a limited-area shallow-water equations model, Zou *et al.* (1993) examined the performance of variational data assimilation via both open-boundary and initial condition controls.

The goal of this research is to develop a numerical approach and algorithm for stable numerical solution of the optimal boundary control problem in lava flow models. In Section 2, we present a mathematical statement of the 2-D problem for reconstruction of temperature and flow pattern in a model of lava. In Section 3, a variational (adjoint) method to solve the problem is discussed. The method is based on the reduction of the problem to minimization of a specially constructed cost functional describing the difference between heat flow inferred from measurements and that from the model (similar techniques were developed, e.g. by Ismail-Zadeh *et al.* 2004, 2006; Korotkii & Kovtunov 2006; and Korotkii & Starodubtseva 2015a). The algorithm for a search of the minimum of the cost functional is described in the section, and the numerical approach to solve the problem is presented in Section 4. In Section 5, we show how this approach can be applied to a lava model and how the temperature and flow velocities can be reconstructed from the surface measurements as well as how the errors in observations can influence the results of the reconstruction. We discuss the challenges of the reconstructions and present conclusions in Section 6.

## 2 STATEMENT OF THE MATHEMATICAL PROBLEM

We study the problem of lava flow in model domain $\Omega \subset \mathbb{R}^2$ (Fig. 2) and assume that the lava behaves as a Newtonian incompressible fluid with a temperature-dependent viscosity and thermal conductivity. This flow is governed by the physical laws of the conservation of momentum, mass and energy, and is described by motion, continuity and heat equations (e.g. Hidaka *et al.* 2005). Although a lava flow is non-stationary depending on an effusive rate, for simplicity of the mathematical problem's description we assume a steady-state lava flow in the modelling. Considering a constant effusion rate, the lava flow in the vicinity of the volcanic vent can be approximated
as a steady-state (see Appendix A). In this assumption and in the Boussinesq approximation, the dimensionless Stokes, continuity and heat equations take the following form:

$$ \nabla \cdot (\mu(T)(\nabla \mathbf{u} + \nabla \mathbf{u}^T)) = \nabla p - Ra T \mathbf{e}_2, \quad (1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2) $$

$$ \nabla \cdot (\kappa(T)\nabla T) = \langle \mathbf{u}, \nabla T \rangle, \quad (3) $$

where $\mathbf{x} = (x_1, x_2) \in \Omega$ are the Cartesian coordinates; $\mathbf{u} = (u_1(\mathbf{x}), u_2(\mathbf{x}))$ is the vector velocity; $p = p(\mathbf{x})$ is the pressure; $T = T(\mathbf{x})$ is the temperature; $\mu = \mu(T)$ is the viscosity; $\kappa = k(T)/(\rho_{\text{ref}}c_p)$ is the thermal diffusivity; $k = k(T)$ is the heat conductivity; $\rho_{\text{ref}}$ is the typical density; and $c_p$ is the specific heat capacity. The Rayleigh number is defined as $Ra = \alpha g \rho_{\text{ref}} \Delta T h^3 \mu_{\text{ref}}^{-1} \kappa_{\text{ref}}^{-1}$, where $\alpha$ is the thermal expansivity; $g$ is the acceleration due to gravity; $\mu_{\text{ref}}$ and $\kappa_{\text{ref}}$ are the typical viscosity and thermal diffusivity, respectively; $\Delta T$ is the temperature contrast; $h$ is the typical length; $\mathbf{e}_2 = (0, -1)$ is the unit vector; $\nabla, \mathbf{T}$ and $\langle \cdot, \cdot \rangle$ denote the gradient vector, the transposed matrix and the scalar product of vectors, respectively. Length and temperature are normalized by $h$ and $\Delta T$, respectively.

We assumed the following conditions for temperature and velocity at the model boundary $\Gamma = \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$. The temperature $T_1$ and the velocity $\mathbf{u}_1$ are prescribed at the left boundary $\Gamma_1$:

$$ T = T_1, \quad \mathbf{u} = \mathbf{u}_1. \quad (4) $$

No slip condition is prescribed at the lower boundary $\Gamma_2$ (unknown temperature is to be found):

$$ \mathbf{u} = 0. \quad (5) $$

At the right boundary $\Gamma_3$, the temperature $T_3$ is prescribed (a strong assumption, which might be omitted, but will complicate the problem solution), the velocity gradient and the pressure are vanishing:

$$ T = T_3, \quad \sigma \mathbf{n} = 0, \quad p = 0. \quad (6) $$

At the upper surface $\Gamma_4$, the temperature $T_4$ and heat flow $\varphi$ are given, and no normal flow and free-slip tangential conditions are used:

$$ T = T_4, \quad k \langle \nabla T, \mathbf{n} \rangle = \varphi, \quad \langle \mathbf{u}, \mathbf{n} \rangle = 0, \quad \sigma \mathbf{n} - \langle \sigma \mathbf{n}, \mathbf{n} \rangle \mathbf{n} = 0, \quad (7) $$

where $\sigma = \mu(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ is the deviatoric stress tensor, and $\mathbf{n}$ is the outward unit normal vector at a point on the model boundary.

The principal problem is to find the solution to eqs (1)–(3) with the boundary conditions (4)–(7), and hence to determine the velocity $\mathbf{u} = \mathbf{u}(\mathbf{x})$, the pressure $p = p(\mathbf{x})$ and the temperature $T = T(\mathbf{x})$ in the model domain $\Omega$ when temperature $T_4$ and heat flow $\varphi = k \partial T / \partial \mathbf{n}$ are known at boundary $\Gamma_4$.

In addition to the principal problem, we define an auxiliary problem, which is formulated as, to find solution to eqs (1)–(3) (i.e. to determine $\mathbf{u}, p$ and $T$ in $\Omega$) with the following boundary conditions,

$$ \Gamma_1: \quad T = T_1, \quad \mathbf{u} = \mathbf{u}_1, \quad (8) $$

$$ \Gamma_2: \quad T = T_2, \quad \mathbf{u} = 0, \quad (9) $$

$$ \Gamma_3: \quad T = T_3, \quad \sigma \mathbf{n} = 0, \quad p = 0, \quad (10) $$

$$ \Gamma_4: \quad T = T_4, \quad \langle \mathbf{u}, \mathbf{n} \rangle = 0, \quad \sigma \mathbf{n} - \langle \sigma \mathbf{n}, \mathbf{n} \rangle \mathbf{n} = 0. \quad (11) $$

The auxiliary problem (1)–(3) and (8)–(11) is a direct problem compared to the problem (1)–(7), which is an inverse problem. We note that the conditions at $\Gamma_1$ and $\Gamma_3$ are the same in the direct and inverse problems, but the temperature $T_2$ is known at $\Gamma_2$ and no heat flow is prescribed at $\Gamma_4$ in the auxiliary problem compared to the inverse problem (1)–(7). The well- and ill-posedness of the similar problems have been studied by Ladyzhenskaya (1969), Lions (1971), Temam (1977), Korotkii & Kovtunov (2006) and Korotkii & Starodubtseva (2015a).

Let us assume now the (measured) heat flow $\varphi = k(T) \partial T / \partial \mathbf{n}$ at model boundary $\Gamma_4$ be related to some (unknown as yet) temperature $T = T_2 = \xi^*$ at model boundary $\Gamma_2$. Let the temperature $T^*$ be a component of the solution $(T^*, \mathbf{u}^*, p^*)$ to the auxiliary problem, when the temperature $T = T_2$ at $\Gamma_2$ equals to $\xi^*$ (see eq. 9), and hence $\varphi = k(T^*) \partial T^* / \partial \mathbf{n}$ at $\Gamma_4$.

Now consider the following cost functional for admissible functions $\xi$ determined at $\Gamma_2$:

$$ J(\xi) = \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right)^2 d\Gamma, \quad (12) $$

where $T_{\xi}$ is the component of the solution $(T_{\xi}, \mathbf{u}_{\xi}, p_{\xi})$ of the auxiliary problem with the condition $T = \xi$ at $\Gamma_2$ in eq. (9). The functional has its global minimum at value $\xi = \xi^*$ and $J(\xi^*) = 0$, that is, temperature $\xi = \xi^*$ attains a minimal value of the functional

$$ J(\xi) \rightarrow \min: \quad \xi \in \Xi, \quad (13) $$

where $\Xi$ denotes a set of admissible temperatures at boundary $\Gamma_2$. Therefore, we reduce the inverse problem to a minimization of the functional or to a variation of the function $\xi$ at $\Gamma_2$, so that heat flow $k \partial T / \partial \mathbf{n}$ at $\Gamma_4$ becomes closer to the prescribed value $\varphi$ at $\Gamma_4$.

### 3 SOLUTION TO THE MINIMIZATION PROBLEM

To minimize the cost functional (12), we use the Polak–Ribière conjugate-gradient method (Polak & Ribière 1969; Polak 1997):

$$ \xi^{(n+1)} = \xi^{(n)} + \gamma^{(n)} d^{(n)}, \quad n = 1, 2, 3, \dots, \quad (14) $$

$$ d^{(n)} = \begin{cases} -\nabla J(\xi^{(n)}), & n = 1 \\ -\nabla J(\xi^{(n)}) + \beta^{(n)} d^{(n-1)}, & n = 2, 3, \dots \end{cases} \quad (15) $$

$$ \beta^{(n)} = \int_{\Gamma_2} \nabla J(\xi^{(n)}) (\nabla J(\xi^{(n)}) - \nabla J(\xi^{(n-1)})) d\Gamma / \int_{\Gamma_2} (\nabla J(\xi^{(n-1)}))^2 d\Gamma, \quad n = 2, 3, \dots, \quad (16) $$

and the descent step length $\gamma^{(n)}$ can be found from the Wolfe conditions (Wolfe 1968, 1969, 1971):

$$ \left\{ \begin{array}{l} J(\xi^{(n)} + \gamma^{(n)} d^{(n)}) \le J(\xi^{(n)}) + c_1 \gamma^{(n)} \int_{\Gamma_2} \nabla J(\xi^{(n)}) d^{(n)} d\Gamma, \\ \int_{\Gamma_2} \nabla J(\xi^{(n)} + \gamma^{(n)} d^{(n)}) d^{(n)} d\Gamma \ge c_2 \int_{\Gamma_2} \nabla J(\xi^{(n)}) d^{(n)} d\Gamma, \end{array} \right. \quad (17) $$

where the search for the descent step length is based on a number of iterations (e.g. Fletcher 2000). Here $\nabla J$ is the gradient of the cost functional; $\xi^{(n)}$ is the $n$-iteration of the admissible function $\xi$; and $0 < c_1 < c_2 < 1$. We assume $c_1 = 0.001$ and $c_2 = 0.01$ in the case of the conjugate-gradient method (and $c_1 = 0.01$ and $c_2 = 0.9$ in the case of the limited-memory Broyden–Fletcher–Goldfarb–Shanno, L-BFGS method; see Section 6).
The gradient of the cost functional

$$ \nabla J(\xi) = \left( k(T_{\xi}) \frac{\partial z}{\partial \mathbf{n}} \right) \Big|_{\Gamma_2}, \qquad (18) $$

can be found as the solution $(z, \mathbf{w}, q)$ to the adjoint problem (see Appendix B for derivation of the adjoint problem)

$$ \nabla \cdot (\mu(T_{\xi}) (\nabla \mathbf{w} + \nabla \mathbf{w}^T)) = \nabla q + z \nabla T_{\xi}, \qquad (19) $$

$$ \nabla \cdot \mathbf{w} = 0, \qquad (20) $$

$$ \begin{aligned} & \nabla \cdot (\kappa(T_{\xi})\nabla z) + \langle \mathbf{u}_{\xi}, \nabla z \rangle + Ra \langle \mathbf{e}_2, \mathbf{w} \rangle \\ & = \kappa'(T_{\xi}) \langle \nabla T_{\xi}, \nabla z \rangle + \mu'(T_{\xi}) [\langle \nabla \mathbf{w} + \nabla \mathbf{w}^T \rangle, \nabla \mathbf{u}_{\xi}], \end{aligned} \qquad (21) $$

with the following boundary conditions

$$ \Gamma_1: z = 0, \quad \mathbf{w} = 0, \qquad (22) $$

$$ \Gamma_2: z = 0, \quad \mathbf{w} = 0, \qquad (23) $$

$$ \Gamma_3: z = 0, \quad \tilde{\sigma} \mathbf{n} = 0, \quad q = 0, \qquad (24) $$

$$ \Gamma_4: z = 2 \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right), \quad \langle \mathbf{w}, \mathbf{n} \rangle = 0, \\ \tilde{\sigma} \mathbf{n} - \langle \tilde{\sigma} \mathbf{n}, \mathbf{n} \rangle \mathbf{n} = 0, \qquad (25) $$

where $\tilde{\sigma} \equiv \nabla \mathbf{w} + \nabla \mathbf{w}^T$; the square brackets $[\mathbf{A}, \mathbf{B}] = \sum_{i,j=1}^{m} a_{ij} b_{ij}$ denote the convolution of two $m \times m$ matrices $\mathbf{A} = (a_{ij})$ and $\mathbf{B} = (b_{ij})$; and $\mathbf{sign}'$ means the derivation. The solution is a triplet $(z, \mathbf{w}, q)$ of quasi-temperature ($z$), quasi-velocity ($\mathbf{w}$) and quasi-pressure $q$.

The algorithm for solving the principal problem can be presented using the following steps (at the initial iteration we assume some guess function $\xi^{(1)} = \xi^{(1)}(\mathbf{x}) \in \Xi$ determined at $\Gamma_2$):

Step 1. Consider $\xi^{(i)} = \xi^{(i)}(\mathbf{x})$, $\mathbf{x} \in \Gamma_2$ ($i = 1, 2, \dots$) as the boundary condition (9) of the auxiliary problem (eqs 1–3 and 8–11) and determine the solution $(T_{\xi^{(i)}}, \mathbf{u}_{\xi^{(i)}}, p_{\xi^{(i)}})$ of this problem in $\Omega$.

Step 2. Insert the components $T_{\xi^{(i)}}$ and $\mathbf{u}_{\xi^{(i)}}$ of the solution into the adjoint problem (eqs 19–25) and determine the solution $(z = z_{\xi^{(i)}}, \mathbf{w} = \mathbf{w}_{\xi^{(i)}}, q = q_{\xi^{(i)}})$ of this adjoint problem in $\Omega$.

Step 3. Determine the gradient of the cost functional $\nabla J(\xi^{(i)})$ from eq. (18) and then $d^{(i)}, \beta^{(i)}$ and $\gamma^{(i)}$ from the conditions (15)–(17), respectively.

Step 4. Determine the value $\xi^{(i+1)}$ from eq. (14).

Step 5. If $J(\xi^{(i+1)}) + \|\nabla J(\xi^{(i)})\|^2 < \varepsilon$, where $\varepsilon > 0$ is a given small number, terminate the minimization problem. Otherwise, the procedure is repeated until the inequality is satisfied.

The performance of the algorithm is evaluated in terms of the number of iterations $n$ required to achieve a prescribed relative reduction of $\xi^{(n)}$. Fig. 3 presents the evolution of the cost functional $J(\xi^{(n)})$ and the norm of the gradient of the objective functional $\|\nabla J(\xi^{(n)})\| = (\int_{\Gamma_2} (\nabla J(\xi^{(n)}))^2 d\Gamma)^{1/2}$ versus the number of iterations.

Implementation of the minimization algorithm requires the evaluation of both the cost functional (12) and its gradient (18). Each evaluation of the objective functional requires an integration of the model eqs (1)–(3) with the appropriate boundary conditions (8)–(11), whereas the gradient is obtained through the integration of the adjoint problem (eqs 19–25). We note that information on the properties of the Hessian $(\nabla^2 J)$ is important in many aspects of minimization problems (Le Dimet et al. 2002). To obtain sufficient

<figure><img src="image_4.png" /><figcaption>Figure 3. Relative reductions of the objective functional (dashed line 1) and the norm of the gradient of the objective functional (solid line 2) as functions of the number of iterations.</figcaption></figure>

conditions for the existence of the minimum of the problem, the Hessian operator must be positive definite.

A viscous dissipation has been neglected in the modelling as the adjoint problem becomes more complicated otherwise. However, the dissipation number, $Di = \mu_{\text{ref}}\kappa_{\text{ref}}/(c_p\rho_{\text{ref}}\Delta T h^2)$, is small enough (about $10^{-7}$ for the lava), so that the viscous dissipation term can be neglected. Thus, the solution of the minimization problem is reduced to solutions of series of well-posed (direct and adjoint) problems.

## 4 NUMERICAL APPROACH

To implement the algorithm for solving the minimization problem, a numerical code was developed using OpenFOAM (http://www.openfoam.org). The finite volume method (e.g. Ismail-Zadeh & Tackley 2010; chapter 3) is used in this software to solve a wide spectrum of fluid dynamic models using multiprocessor computers. Particularly, it includes the codes for numerical solution of the Stokes and advection-diffusion problems with various boundary conditions and model geometries.

The model domain $\Omega$ was discretized by 1500 hexahedral finite volumes (60 and 25 volumes in the horizontal and vertical directions). The SIMPLE method (Patankar & Spalding 1972; Ismail-Zadeh & Tackley 2010, chapter 6.5.2) was used to determine velocity and pressure at a given temperature (the relaxation parameters are chosen to be 0.7 and 0.3 for the velocity and pressure, respectively). To implement the SIMPLE method, we employ the conjugate-gradient method (Ismail-Zadeh & Tackley 2010, chapter 6.3.3) to solve a set of linear algebraic equations (SLAE) with positive-definite and symmetric matrices, which are obtained after the discretization of the Stokes equation. 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 relative accuracy of the numerical solutions to the derived SLAE (i.e. the ratio between the norm of the residual to the norm of the right-hand-side of the SLAE) is $10^{-3}$ in the case of the model domain discretization by 1500 finite volumes. The
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Reconstruction of the temperature at the boundary $\Gamma_2$ (a). The red curve corresponds to the target temperature, the green curve to the guess temperature, the brown curve to the temperature after 5 iterations and the blue curve to the temperature after the 10 iterations. The reconstructed temperature after 10 iterations (b) in the case of no noise in the heat flow at the upper boundary of the lava (solid line; the blue curve in panel a) and in the case of the noise magnitude $\delta = 0.1$ in the heat flow (dashed line).</figcaption></figure>

linear Gaussian scheme with a flow control was used to discretize the Laplace operator. To approximate the convective operator, we employed the total variation diminishing method with the minmod limiter (Sweby 1984; Wang & Hutter 2001; Ismail-Zadeh et al. 2007).

All computations were performed using one CPU Intel Core i5 2.6 GHz with 16GB memory, OS X 10.10. An average computational time for 80 iterations in the inverse problem was 75 min: this included the time required for minimization of the cost functional by the conjugate gradient method, and the time to solve the direct and adjoint problems (normally 4–5 times) to determine the descent step length.

## 5 MODEL RESULTS

We consider a model of lava advancing down the slope (Fig. 2) and assume that temperature and heat flow are available from remote thermal measurements. The boundary of the model domain consists of the following parts: $\Gamma_1$ is a line segment connecting points $\mathbf{x}^A = (x_1^A, x_2^A) = (0, 2.5)$ and $\mathbf{x}^B = (x_1^B, x_2^B) = (0, 1.5)$; $\Gamma_2$ is a circular arc connecting points $\mathbf{x}^B, \mathbf{x}^C = (x_1^C, x_2^C) = (1.5, 0.5)$ and $\mathbf{x}^D = (x_1^D, x_2^D) = (3.0, 0.0)$; $\Gamma_3$ is a line segment connecting points $\mathbf{x}^D$ and $\mathbf{x}^E = (x_1^E, x_2^E) = (3.0, 0.5)$; and $\Gamma_4$ is a circular arc connecting points $\mathbf{x}^E, \mathbf{x}^F = (x_1^F, x_2^F) = (1.5, 1.2)$ and $\mathbf{x}^A$.

The following dimensional parameters are used in the modelling: $\alpha = 10^{-5} \text{ K}^{-1}$, $g = 9.8 \text{ m s}^{-2}$, $h = 1 \text{ m}$, $\rho_{\text{ref}} = 3000 \text{ kg m}^{-3}$, $\mu_{\text{ref}} = 3.5 \times 10^6 \text{ Pa s}$, $T_{\text{ref}} = 300 \text{ K}$, $T_* = 1473 \text{ K}$, $\Delta T = T_* - T_{\text{ref}}$, $\kappa_{\text{ref}} = 10^{-6} \text{ m}^2 \text{ s}^{-1}$, $c_p = 1200 \text{ J kg}^{-1} \text{ K}^{-1}$, $k_{\text{ref}} = \rho_{\text{ref}} c_p \kappa_{\text{ref}} = 3.6 \text{ W m}^{-1} \text{ K}^{-1}$, and therefore, the Rayleigh number is $Ra = 100$. We consider the dimensionless temperature-dependent viscosity (Griffiths 2000)

$$ \mu(T) = \exp(0.039(4.91 - T)) \quad (26) $$

and the dimensionless temperature-dependent thermal conductivity (Hidaka et al. 2005)

$$ k(T) = \begin{cases} 0.32 + 4.92 \cdot 10^{-5}(T - 4.91)^2, & T < 4.91, \\ 0.32 + 8.08 \cdot 10^{-4}(T - 4.91)^2, & T > 4.91. \end{cases} \quad (27) $$

At $\Gamma_1$ we prescribe the temperature $T_1(x_1, x_2) = 5.0 - 0.5(x_2 - x_1^B)$, $x_2 \in [x_2^B, x_2^A]$, and the velocity $\mathbf{u}_1(x_2) = U(x_2)\mathbf{n}_1$, where $\mathbf{n}_1 = (\sqrt{2}/2, -\sqrt{2}/2)$ and $U(x_2)$ is the parabola passing through the following three points: $U(x_2^A) = 10$, $U(x_2^B) = 0$ and $U(0.5(x_2^A + x_2^B)) = 7.25$.

The temperature is $T_3(x_1, x_2) = 3.5 - 2(x_2 - x_1^D)$, $x_2 \in [x_2^D, x_2^E]$ at $\Gamma_3$ and $T_4(x_1, x_2) = 4.5 - 2(x_1 - x_2^A)/3$, $x_1 \in [x_1^A, x_1^E]$ at $\Gamma_4$. Considering guess temperature $\xi^{(1)} = \xi^{(1)}(\mathbf{x})$ at $\Gamma_2$, we use the algorithm described in Section 3 (steps 1–5) to find the temperature at $\Gamma_2$. Doing so, we assimilate the thermal data from the boundary $\Gamma_4$ to $\Gamma_2$ through the model domain by solving the inverse problem (1)–(7).

The cost functional is reduced to about $10^{-5}$ after 30 iterations (Fig. 3). The reconstruction of the temperature at the boundary $\Gamma_2$ versus the number of iterations is presented in Fig. 4(a). We note that the number of iterations to get a given accuracy in reduction of the cost functional depends on the initial ‘guess’ temperature at $\Gamma_2$. The closer is the guess temperature to the target temperature, the less number of iterations is needed.

Fig. 5 shows the reconstruction process of the lava temperature and flow velocity from the initial iteration to the 80th iteration. The temperature and velocity residuals, that is, the difference between the temperature and velocities predicted by the forward model (with the target temperature at $\Gamma_2$) and those reconstructed, are also presented in Fig. 5. The results of this modelling show that the restoration works quite well: the temperature residuals are very low already after 80 iterations within the almost entire model domain.

The comparison between ‘measured’ (modelled) and reconstructed lava temperature is quite natural from the computational point of view, but not from the geophysical point of view, because the measurements (observations) are polluted by errors. The accuracy of temperature measurements and inferred heat flux density can be attributed to the accuracy of the calibration curve of remote sensors and the noise of the sensors. Considering these sources of errors of measured temperatures, the errors would range from 0.1 K to 1 K (Short & Stuart 1983). The heat flow errors inferred from the Stefan–Bolzmann law can be then estimated between 0.6 and 6 W m$^{-2}$ at the reference temperature $T_{\text{ref}} = 300 \text{ K}$, which are related to dimensionless error values from 0.0013 to 0.013 (normalized with respect to heat flow at the reference temperature).

Hence, we perform numerical experiments introducing a noise on the ‘measured’ data and study how well the problem can be resolved. Particularly, we introduce a disturbance on the heat flow $\varphi(\cdot)$ at the boundary $\Gamma_4$ as $\varphi_\delta(\cdot) = \varphi(\cdot) + \delta\gamma(\cdot)$, where $\delta$ is the magnitude of the disturbance; $\gamma(\cdot)$ is the function generating numbers that are uniformly distributed over the interval $[-1, 1]$; and $\varphi(\cdot)$ is obtained from the solution of eqs (1)–(3) with the conditions at the boundaries
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Reconstruction of the lava temperature (a) and the flow velocity (c) after 20 and 80 iterations. The relevant residuals of the temperature (b) and the velocity (d) indicate the quality of the reconstruction.</figcaption></figure>

(8)–(11) for $T_2 = \xi^{(1)}$ at $\Gamma_2$. We choose three values for the noise magnitude $\delta$ (0.001, 0.01 and 0.1) to approximate the possible noise level of the remote thermal measurements.

We analyse the influence of the noise on the reconstruction of the temperature at the boundary $\Gamma_2$ (Fig. 4b), on the temperature and velocity residuals (Fig. 6) and on the reduction of the cost functional and the norm of its gradient (Fig. 7). The computations show that the errors (temperature and velocity residuals, Fig. 6) get larger with increase of the noise of the input data. Meanwhile for some range of the noise ($\delta \le 0.01$) the reconstructions are still reasonable as the temperature and velocity residuals are not high (Fig. 6). Namely, if we consider $M_T = \max_{x \in \Omega} |T_{30}(x) - T^0(x)|$ and $M_u = \max_{x \in \Omega} \|u_{30}(x) - u^0(x)\|_{\mathbb{R}^2}$, where $T^0(x)$ and $u^0(x)$ are the solution of the direct problem (1)–(3) and (8)–(11), then $M_T = 0.095$, 0.096, 0.099 and 0.265, and $M_u = 0.0073, 0.0074, 0.0075$ and 0.01526 for $\delta = 0, 0.001, 0.01$ and $0.1$, respectively.

# 6 DISCUSSION AND CONCLUSION

We have proposed a numerical approach to solving the optimal boundary control problem arising in studies of lava dynamics. This approach is stable to numerical errors in the input data and permits reconstructing the thermal state of the lava flow based on the measured temperature and inferred heat flow information at the lava's surface. The measured data can be assimilated to the lava's lower part using direct and adjoint lava flow models. The efficiency of two optimization techniques—the conjugate gradient and the L-BFGS algorithms—has been compared, and the convergence to the solution has been analysed.

A rapid development of ground-based thermal cameras, drones and satellite data allows getting repeated thermal images of the sur-

face of the lava flow (Calvari *et al.* 2005). Available instrumentation allows getting a large amount of data during a single lava flow eruption. For example, lava flow emplacement between 2010 December 4 and 6 at Tungurahua volcano (Ecuador) was recorded with a thermal Forward Looking InfraRed camera; about 90 recorded thermal images of the lava flow (Kelfoun & Vargas 2015). These data require development of appropriate quantitative approaches to link subsurface dynamics with observations. Our approach permits to link the observations with the dynamics of lava flow and is a step forward to determine a heat budget of the active lava flows.

Earth orbiting radiometers can measure spectral radiance at a lava surface to be converted then into thermal anomalies. Lava temperature and heat flow can be inferred from the detected anomalies. However, a spatial resolution of many satellites is coarse enough to allow for high-resolution monitoring and precise measurements. This gives a rise to uncertainties in thermal measurement as well as in the inferred parameters. Hence, if the measured temperature and heat flow data are biased, this information can be improperly assimilated into the lava flow models. According to Zakšek *et al.* (2015) to reduce major sources of uncertainties in thermal anomaly monitoring, satellite instruments should allow for (i) fine spatial resolution (e.g. less than 50 m), (ii) short revisit time (less than 15 min), (iii) multiple spectral bands, (iv) high radiometric accuracy (<0.1 K) and (v) observations at different instrument gain settings to determine high and low temperature anomalies. If surface temperature and heat flow data are of high resolution and radiometric accuracy, the temperature and velocity in the lava's interior can be determined properly from measured data using the proposed numerical approach.

The proposed approach can become important in studies of natural lava flows, especially in the cases of thick lava flow. Synthetic Aperture Radar satellite observations on lava thickness, volume
<figure><img src="image_7.png" /><figcaption>Figure 6. Residuals of temperature (the left panels) and velocity (the right panels) at different noise magnitudes δ.</figcaption></figure>

and flow extent (e.g. Kubanek et al. 2015), together with thermal measurement at the lava surface, facilitate research and data-driven modelling of lava flow. Moreover, the proposed approach to assimilate measured data into the model brings a unique opportunity to estimate thermal budget of the lava flow. Widely used assumptions that basal heat flux is much smaller than the flux to the atmosphere can be validated if the whole temperature field inside the lava flow is known.

To model numerically a lava flow, boundary and initial conditions should be known. Meanwhile the temperature or heat flux at the lava's bottom is unknown as it is almost impossible to measure it. Airborne and space measurements of temperature (heat flux) at the surface of lava flows, being almost instantaneous compared to the duration of lava flows, allow to search for thermal conditions at the lava's bottom. Once the boundary conditions at the lava bottom are determined, the steady-state problem can be replaced by a non-stationary 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 or quasi-reversibility methods (e.g. Ismail-Zadeh et al. 2009) to search for the initial temperature of the lava flow and for the evolution of the effusion rate.

In this study, two thermal conditions at the upper surface of the model domain have been required to assimilate these data to the lower boundary of the model domain. We have used the known temperature at the upper surface to solve the direct problem and the known heat flow at the same surface to solve the adjoint problem. Meanwhile the thermal conditions at the upper surface can be permuted, namely: heat flow can be prescribed (instead of temperature) and temperature (instead of heat flow) at the upper surface to solving the direct and the adjoint problems, respectively (Korotkii & Starodubtseva 2015b).

The effusion rate, at which lava is erupted, controls the way in which a lava body grows and extends influencing its dimensional properties. Estimations of effusion rates based on radiated heat flux (Harris et al. 2004) is at the moment based on a very crude model that does not account for basal heat flux. To test a sensitivity of the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Relative reductions of the objective functional (dashed lines) and the norm of the gradient of the objective functional (solid lines) as functions of the number of iterations for several values of noise magnitude $\delta$: black lines stand for $\delta = 0$, red lines for $\delta = 0.001$, blue lines for $\delta = 0.01$ and green lines for $\delta = 0.1$.</figcaption></figure>

approach to changes in flow patterns, we varied the magnitude of the velocity $|U|$ at the left-side boundary $\Gamma_1$ of the model domain between 1 and 25, and the Rayleigh number, $Ra$, between 1 and 10000. The approach is rather robust to changes in the velocity magnitude and in $Ra$. (Note that in the case of a lava flow the Rayleigh number is small enough and in many applications considered to be close to zero. Higher $Ra$ were used here to test the stability of the numerical algorithm.)

A performance of the numerical approach depends on optimization methods employed. In the presented approach we used the conjugate-gradient method. To compare its performance with other optimization methods, we have employed the L-BFGS method of Liu & Nocedal (1989), based on the L-BFGS algorithm described by Nocedal (1980). To minimize the cost functional (12) using the L-BFGS method, components $d^{(n)}$ in eq. (14) are determined as $d^{(1)} = -\nabla J(\xi^{(1)})$ and $d^{(n)} = -\mathbf{B}^{(n)} \nabla J(\xi^{(n)})(n = 2, 3, ...)$, where $\mathbf{B}^{(n)}$ is the approximated inverse Hessian operator.

When the L-BFGS method is used, the average computational time to perform 80 iterations for minimization of the cost functional is reduced to 15 min (by the factor of 5) compared to the case of the conjugate-gradient method used. The computational time reduction is achieved because the descent step length in the iteration scheme is determined much faster. The reduction of the objective functional and the norm of the gradient of the objective functional with the number of iterations is faster than in the case of the conjugate-gradient method (Fig. 8), although the dependence of the solution on the noise magnitude is similar in the both cases (compare Figs 7 and 9).

Limited-memory quasi-Newton (LMQN) methods represent a class of algorithms, which use a low amount of storage to accelerate the convergence rate, which is important in cases of large-scale problems. Several comparative studies were performed to clarify the best-performance methods/algorithms among LMQN and conjugate-gradient methods, for example, Navon & Legler (1987), Gilbert & Lemarichal (1989), Liu & Nocedal (1989) and Zou *et al.* (1993) indicated that the L-BFGS method, belonging to the LMQN

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Relative reductions of the objective functional (dashed lines) and the norm of the gradient of the objective functional (solid lines) in the case of the conjugate gradient method (black lines) and in the case of the L-BFGS method (red lines).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The influence of noise on model solutions in the case of the L-BFGS method used: relative reductions of the objective functional (dashed lines) and the norm of the gradient of the objective functional (solid lines) versus the number of iterations. Red lines stand for $\delta = 0$, blue lines for $\delta = 0.001$, brown lines for $\delta = 0.01$ and green lines for $\delta = 0.1$.</figcaption></figure>

family (e.g. Nocedal & Wright 1999), was one of the best algorithms for several problems examined. Meanwhile the conjugate gradient method with guaranteed descent developed later (Hager & Zhang 2006) shows in some cases higher performance compared to the L-BFGS and other gradient methods (Alekseev *et al.* 2009).

There are several simplifications in the presented model of lava flow that can be overcome in future, but require further development of the algorithm and increase in computational resources. For example, the proposed numerical approach allows also for reconstructing the temperature at the right boundary of the model domain (if heat flux is negligible at its lower boundary) or at lower
and right boundaries simultaneously. In the present model, the shape of the lava is prescribed (the position of the right boundary $\Gamma_3$ is constant). In reality, the lava flows and its right boundary should also be reconstructed during the solution with an extra boundary condition related to the absence of tangential stress at the free surface.

The problem can be extended to the non-steady-state flow, but this will complicate the mathematical and computational approaches. Meanwhile, as the measurements on absolute temperature are discrete in time (e.g. depending on the location of Landsat satellites), a problem of non-stationary flow can be reduced to a number of steady-state flow problems with varying boundary conditions at the upper model surface (where the discrete-in-time measurement are available).

A more complicated lava rheology with formation and disintegration of solid crust (e.g. Tsepelev et al. 2016) should be considered. The influence of the shape of the crust and the degree of its disintegration on the radiated heat flux (Neri 1998) can be significant. Proposed algorithm and its numerical implementation have a wide range of applications in other problems of reconstruction of the flows of fluids with strongly temperature dependent viscosity, for example, in chemical technology or oil industry.

## ACKNOWLEDGEMENTS

We thank S. Fedotov and A. Harris for discussions on lava flows and Landsat measurements, and two reviewers and the Editor for their constructive comments, which improved the initial manuscript. This research was supported by the grant of the Russian Science Foundation (RSF project 14-17-00520).

## REFERENCES

Alekseev, A.K., Navon, I.M. & Steward, J., 2009. Comparison of advanced large-scale minimization algorithms for the solution of inverse ill-posed problems, *Optim. Methods Softw.*, **24**(1), 63–87.

Calvari, S., Spampinato, L., Lodato, L., Harris, A.J.L., Patrick, M.R., Dehn, J., Burton, M.R. & Andronico, D., 2005. Chronology and complex volcanic processes during the 2002–2003 flank eruption at Stromboli volcano (Italy) reconstructed from direct observations and surveys with a handheld thermal camera, *J. geophys. Res.*, **110**, B02201, doi:10.1029/2004JB003129.

Cutter, S et al., 2015. Global risks: pool knowledge to stem losses from disasters, *Nature*, **522** (7556), 277–279.

Fletcher, R., 2000. *Practical Methods of Optimization*, 2nd edn, Wiley & Sons.

Flynn, L.P., Harris, A.J.L. & Wright, R., 2001. Improved identification of volcanic features using Landsat 7 ETM+, *Remote Sens. Environ.*, **78**, 180–193.

Gilbert, J.C. & Lemarichal, C., 1989. Some numerical experiments with variable-storage quasi-Newton algorithms, *Math. Program.*, **45**, 407–435.

Griffiths, R.W., 2000. The dynamics of lava flows, *Annu. Rev. Fluid Mech.*, **32**, 477–518.

Hadamard, J., 1923. *Lectures on the Cauchy Problem in Linear Partial Differential Equations*, Yale University Press.

Hager, W.W. & Zhang, H., 2006. Algorithm 851: CG DESCENT, a conjugate gradient method with guaranteed descent, *ACM Trans. Math. Softw.*, **32**(1), 113–137.

Harris, A.J.L., Flynn, L.P., Matias, O., Rose, W.I. & Cornejo, J., 2004. The evolution of an active silicic lava flow field: an ETM+ perspective, *J. Volcanol. Geotherm. Res.*, **135**, 147–168.

Hidaka, M., Goto, A., Umino, S. & Fujita, E., 2005. VTFS project: Development of the lava flow simulation code LavaSIM with a model for three-dimensional convection, spreading, and solidification, *Geochem. Geophys. Geosyst.*, **6**, Q07008, doi:10.1029/2004GC000869.

Huppert, H.E., 1982. The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface, *J Fluid Mech.*, **21**, 43–58.

Ismail-Zadeh, A. & Tackley, P., 2010. *Computational Methods for Geodynamics*, Cambridge Univ. Press.

Ismail-Zadeh, A., Schubert, G., Tsepelev, I. & Korotkii, A., 2004. Inverse problem of thermal convection: Numerical approach and application to mantle plume restoration, *Phys. Earth planet. Inter.*, **145**, 99–114.

Ismail-Zadeh, A., Schubert, G., Tsepelev, I. & Korotkii, A., 2006. Three-dimensional forward and backward numerical modeling of mantle plume evolution: Effects of thermal diffusion, *J. geophys. Res.*, **111**, B06401, doi:10.1029/2005JB003782.

Ismail-Zadeh, A., Korotkii, A., Schubert, G. & Tsepelev, I., 2007. Quasireversibility method for data assimilation in models of mantle dynamics, *Geophys. J. Int.*, **170**, 1381–1398.

Ismail-Zadeh, A., Korotkii, A., Schubert, G. & Tsepelev, I., 2009. Numerical techniques for solving the inverse retrospective problem of thermal evolution of the Earth interior, *Comput. Struct.*, **87**, 802–811.

Kelfoun, K. & Vargas, S.V., 2015. VolcFlow capabilities and potential development for the simulation of lava flows, in *Detecting, Modelling and Responding to Effusive Eruptions*, vol. 426, eds Harris, A.J.L., De Groeve, T., Farel, F. & Carn, S.A., Geological Society, London, Special Publications, doi:10.1144/SP426.8.

Kirsch, A., 1996. *An Introduction to the Mathematical Theory of Inverse Problems*, Springer-Verlag.

Korotkii, A.I. & Kovtunov, D.A., 2006. Reconstruction of boundary regimes in the inverse problem of thermal convection of a high-viscosity fluid, *Proceedings of the Steklov Institute of Mathematics*, **295**(Suppl. 2), S81–S92. (Translated and published by the Pleiades Publishing, Inc. from the original Russian text published in *Proceedings of the Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences*, **12**(2), 88–97, 2006).

Korotkii, A.I. & Starodubtseva, Yu.V., 2015a. Direct and inverse boundary value problems for models of stationary reaction–convection–diffusion, *Proceedings of the Steklov Institute of Mathematics*, **291**(Suppl. 1), S96–S112. (Translated and published by the Pleiades Publishing, Inc. from the original Russian text published in *Proceedings of the Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences*, **20**(3), 98–113, 2014).

Korotkii, A.I. & Starodubtseva, Yu.V., 2015b. Modelling of Direct and Inverse Boundary Problems for Steady-State Models of Heat and Mass Transfer, The Ural State University Press (in Russian).

Kubanek, J., Richardson, J.A., Charbonnier, S.J. & Connor, L.J., 2015. Lava flow mapping and volume calculations for the 2012–2013 Tolbachik, Kamchatka, fissure eruption using bistatic TanDEM-X InSAR, *Bull. Volcanol.*, **77**(12), doi:10.1007/s00445-015-0989-9.

Ladyzhenskaya, O.A., 1969. *The Mathematical Theory of Viscous Incompressible Flow*, Gordon and Breach.

Le Dimet, F.X., 1988. Determination of the adjoint of a numerical weather prediction model. Tech. Rep. FSU-SCRI-88-79, The Florida State University, Tallahassee, FL, 22 pp.

Le Dimet, F.X., Navon, I.M. & Daescu, D.N., 2002. Second order information in data assimilation, *Mon. Weather Rev.*, **130**(3), 629–648.

Lions, J.L., 1971. *Optimal Control of Systems Governed by Partial Differential Equations*, Springer-Verlag.

Liu, D.C. & Nocedal, J., 1989. On the limited memory BFGS method for large scale optimization, *Math. Program.*, **45**, 503–528.

Navon, I.M. & Legler, D.M., 1987. Conjugate-gradient methods for large-scale minimization in meteorology, *Mon. Weather Rev.*, **115**(4), 1479–1502.

Navon, I.M., Zou, X., Derber, J. & Sela, J., 1992. Variational data assimilation with an adiabatic version of the NMC spectral model, *Mon. 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. Geotherm. Res.*, **81**, 215–243.
Nocedal, J., 1980. Updating quasi-Newton matrices with limited storage, *Math. Comput.*, **35**, 773–782.

Nocedal, J. & Wright, S.J., 1999. *Numerical Optimization*, Springer.

Patankar, S.V. & Spalding, D.B., 1972. A calculation procedure for heat and mass transfer in three-dimensional parabolic flows, *Int. J. Heat Mass Transfer*, **15**, 1787–1806.

Polak, E., 1997. *Optimization: Algorithms and Consistent Approximations*, Springer-Verlag.

Polak, E. & Ribière, G., 1969. Note on the convergence of methods of conjugate directions, *Rev. Fr. Inform. Rech. Oper.*, **3**(16), 35–43.

Short, N.M. & Stuart, L.M., 1983. *The Heat Capacity Mapping Mission (HCMM) Anthology*, Scientific and Technical Information Branch, National Aeronautics & Space Administration.

Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws, *SIAM J. Numer. Anal.*, **21**, 995–1011.

Temam, R., 1977. *Navier-Stokes Equations: Theory and Numerical Analysis*, North-Holland.

Thacker, W.C. & Long, R.B., 1988. Fitting dynamics to data, *J. geophys. Res.*, **93**(C2), 1227–1240.

Tikhonov, A.N. & Arsenin, V.Y., 1977. *Solution of Ill-Posed Problems*, Winston.

Tsepelev, I., Ismail-Zadeh, A., Melnik, O. & Korotkii, A., 2016. Numerical modelling of fluid flow with rafts: an application to lava flows, *J. Geodyn.*, doi:10.1016/j.jog.2016.02.010.

van der Vorst, H.A., 1992. BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems, *SIAM J. Sci. Stat. Comput.*, **13**(2), 631–644.

Wadge, G., Young, P.A.V. & McKendrick, I.J., 1994. Mapping lava flow hazards using computer simulation, *J. geophys. Res.*, **99**, 489–504.

Wang, Y. & Hutter, K., 2001. Comparison of numerical methods with respect to convectively dominated problems, *Int. J. Numer. Methods Fluids*, **37**, 721–745.

Wolfe, P., 1968. The secant method for simultaneous nonlinear equations, *Commun. ACM*, **2**, 12–13.

Wolfe, P., 1969. Convergence conditions for ascent methods, *SIAM Rev.*, **11**, 226–235.

Wolfe, P., 1971. Convergence conditions for ascent methods. II: Some corrections, *SIAM Rev.*, **13**, 185–188.

Zakšek, K., Hort, M. & Lorenz, E., 2015. Satellite and ground based thermal observation of the 2014 effusive eruption at Stromboli volcano, *Remote Sens.*, **7**, 17190–17211.

Zou, J., Hsieh, W.H. & Navon, I.M., 1995. Sequential open-boundary control by data assimilation in a limited area model, *Mon. Weather Rev.*, **123**(9), 2899–2909.

Zou, X., Navon, I.M., Berger, M., Phua, M.K., Schlick, T. & Le Dimet, F.X., 1993. Numerical experience with limited-memory, quasi-Newton methods for large-scale unconstrained nonlinear minimization, *SIAM J. Optim.*, **3**, 582–608.

## APPENDIX A: FLUID FLOW HEIGHT ASSESSMENT

In this study we use a steady-state condition for the temperature field. In order to validate this assumption, we employ the self-similar solution for a fluid flow by Huppert (1982). A temperature field can be considered as steady-state if the effusion rate is constant, and the lava flow height and its velocity do not change significantly with time. Consider an intruding fluid (lava) of density $\rho$ and viscosity $\mu$ (Fig. A1). In the case of the constant discharge rate of lava $q$, the flow height $z(x, t)$ in point $x$ and time $t$ can be found from the self-similar solution as:

$$z(x, t) = 1.247 \left( \frac{q^2 \mu}{\rho g} \right)^{1/5} t^{1/5} \phi(\xi), \quad \xi = \frac{x}{x_N}, \quad x_N = 0.804 \left( \frac{\rho g q^3}{\mu} \right)^{1/5} t^{4/5}, \qquad (A1)$$

$$\phi(\xi) = 1.3387(1 - \xi)^{1/3} (1.0104 - 0.0104\xi), \qquad (A2)$$

where $g$ is the acceleration due to gravity. All other flow characteristics are functions of the height and its special derivatives, thus, small changes in the flow height in a particular part of the flow with time on the timescale of temperature equilibration guarantees steady-state conditions. We normalize equations in (A1) by the following change of variables $x = \tilde{x}h$, $z = \tilde{z}h$, $\rho = \tilde{\rho}\rho_{\text{ref}}$, $\mu = \tilde{\mu}\rho_{\text{ref}}gh t_{\text{ref}}$, $t = \theta t_{\text{ref}}$ and $q = \tilde{q}h^2/t_{\text{ref}}$. For validity of the lubrication-theory approximation we assume $\tilde{\rho} = 1$, $\tilde{\mu} = 10^{-4}$, and $\tilde{q} = 50$. Omitting tilde from the dimensionless variables, we obtain the dimensionless equations: $z(\xi, \theta) = 0.9478\theta^{1/5}\phi(\xi)$, $x_N = 53.05\theta^{4/5}$.

Considering that the thermal equilibrium is attained at time $t_{\text{ref}} = h^2/\kappa_{\text{ref}}$, we analyse the lava height at dimensionless times $\theta \ge 3$. The lava tip $x_N = 128$ at $\theta_1 = 3$ and 280 at $\theta_2 = 8$. Fig. A2 shows the evolution of the lava height at some distances from the vent ($x = 0$). The relative change of the height expressed by the formula $(z(\theta_2) - z(\theta_1))/(z(\theta_1)(\theta_2 - \theta_1))$ is about 4.6 per cent at $x = 6.84$, 5.9 per cent at $x = 34.22$ and 16.1 per cent at $x = 102.66$. With a desired accuracy, the steady-state approach for the temperature field can be then valid in the area close to the vent.

<figure><img src="image_11.png" /><figcaption>Figure A1. A sketch of the flow field and coordinate system (modified after Huppert 1982).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A2.** Evolution of the lava height with time at three distances $x$ from the vent.</figcaption></figure>

## APPENDIX B: COST FUNCTIONAL J AND ITS GRADIENT

Here we derive the adjoint problem. Let the triplet $(T_{\xi+\chi}, \mathbf{u}_{\xi+\chi}, p_{\xi+\chi})$ be the solution of the auxiliary problem (1)–(3), (8)–(11) for the prescribed condition $T = T_2 = \xi + \chi$ at the boundary $\Gamma_2$ (see eq. 9) and the triplet $(T_{\xi}, \mathbf{u}_{\xi}, p_{\xi})$ be the solution of the same problem for the prescribed condition $T = T_2 = \xi$ at the same boundary, where $\chi$ is an admissible increment of the boundary element $\xi$. The difference of the two solutions $T = T_{\xi+\chi} - T_{\xi}$, $\mathbf{u} = \mathbf{u}_{\xi+\chi} - \mathbf{u}_{\xi}$ and $p = p_{\xi+\chi} - p_{\xi}$ should satisfy the following boundary value problem for $\mathbf{x} \in \Omega$:

$$ \nabla \cdot (\delta\mu(T_{\xi})(\nabla\mathbf{u} + \nabla\mathbf{u}^T)) + \nabla \cdot (\mu(T_{\xi})(\nabla\mathbf{u} + \nabla\mathbf{u}^T)) + \nabla \cdot (\delta\mu(T_{\xi})(\nabla\mathbf{u}_{\xi} + \nabla\mathbf{u}_{\xi}^T)) = \nabla p - Ra T \mathbf{e}_2, \quad (B1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (B2) $$

$$ \nabla \cdot (\delta\kappa(T_{\xi})\nabla T) + \nabla \cdot (\kappa(T_{\xi})\nabla T) + \nabla \cdot (\delta\kappa(T_{\xi})\nabla T_{\xi}) = \langle \mathbf{u}, \nabla T \rangle + \langle \mathbf{u}_{\xi}, \nabla T \rangle + \langle \mathbf{u}, \nabla T_{\xi} \rangle, \quad (B3) $$

with the following boundary conditions

$$ \Gamma_1: \quad T = 0, \quad \mathbf{u} = 0, \quad (B4) $$

$$ \Gamma_2: \quad T = \chi, \quad \mathbf{u} = 0, \quad (B5) $$

$$ \Gamma_3: \quad T = 0, \quad \sigma \mathbf{n} = 0, \quad p = 0, \quad (B6) $$

$$ \Gamma_4: \quad T = 0, \quad \langle \mathbf{u}, \mathbf{n} \rangle = 0, \quad \sigma \mathbf{n} - \langle \sigma \mathbf{n}, \mathbf{n} \rangle \mathbf{n} = 0, \quad (B7) $$

where $\delta\mu(T_{\xi}) = \mu(T_{\xi+\chi}) - \mu(T_{\xi})$ and $\delta\kappa(T_{\xi}) = \kappa(T_{\xi+\chi}) - \kappa(T_{\xi})$. We note that

$$ \begin{align*} J(\xi + \chi) - J(\xi) &= \int_{\Gamma_4} \left( k(T_{\xi+\chi}) \frac{\partial T_{\xi+\chi}}{\partial \mathbf{n}} - \varphi \right)^2 d\Gamma - \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right)^2 d\Gamma \\ &= 2 \int_{\Gamma_4} \left( k(T_{\xi+\chi}) \frac{\partial T_{\xi+\chi}}{\partial \mathbf{n}} - k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + \int_{\Gamma_4} \left( k(T_{\xi+\chi}) \frac{\partial T_{\xi+\chi}}{\partial \mathbf{n}} - k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right)^2 d\Gamma \\ &= 2 \int_{\Gamma_4} \left( k(T_{\xi+\chi}) \frac{\partial T_{\xi+\chi}}{\partial \mathbf{n}} - k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|), \end{align*} $$
and accounting for $k(T_{\xi+\chi}) = k(T_{\xi}) + k'(T_{\xi})T + o(\|T\|) = k(T_{\xi}) + k'(T_{\xi})T + o(\|\chi\|)$, we obtain

$$ 
\begin{aligned}
& 2 \int_{\Gamma_4} \left( k(T_{\xi+\chi}) \frac{\partial T_{\xi+\chi}}{\partial \mathbf{n}} - k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|) \\
&= 2 \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T}{\partial \mathbf{n}} + o(\|\chi\|) \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|) \\
&= 2 \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} + o(\|\chi\|) \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|) \\
&= 2 \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|),
\end{aligned}
$$

and hence

$$ J(\xi + \chi) - J(\xi) = \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) 2 \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right) d\Gamma + o(\|\chi\|). \quad (B8) $$

We assume that a test function $\mathbf{w} = \mathbf{w}(\mathbf{x})$, $\mathbf{x} \in \Omega$ satisfies the incompressibility condition

$$ \nabla \cdot \mathbf{w} = 0 \quad (B9) $$

and the following boundary conditions

$$ \Gamma_1, \Gamma_2: \mathbf{w} = 0, \quad (B10) $$

$$ \Gamma_3: \tilde{\sigma} \mathbf{n} = 0, \quad (B11) $$

$$ \Gamma_4: \langle \mathbf{w}, \mathbf{n} \rangle = 0, \quad \tilde{\sigma} \mathbf{n} - \langle \tilde{\sigma} \mathbf{n}, \mathbf{n} \rangle = 0. \quad (B12) $$

Now we multiply eq. (B1) by a test function $\mathbf{w} = \mathbf{w}(\mathbf{x})$ and integrate the resultant equation over $\Omega$. Considering eqs (B9)–(B12) and after integrating by parts, we obtain

$$ \int_{\Omega} \langle \mathbf{u}, \nabla \cdot (\mu(T_{\xi}) (\nabla \mathbf{w} + \nabla \mathbf{w}^T)) \rangle dx - \int_{\Omega} \mu'(T_{\xi}) T [\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi}] dx + \int_{\Omega} Ra T \langle \mathbf{w}, \mathbf{e}_2 \rangle dx = o(\|\chi\|), \quad (B13) $$

where the relation $[\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi}]$ can be represented in a symmetric form as $[\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi} + \nabla \mathbf{u}_{\xi}^T]/2$. We multiply eq. (B2) by a test scalar function $q = q(\mathbf{x})$, $\mathbf{x} \in \Omega$, and integrate by parts the resultant equation over $\Omega$. Assuming that the function $q = 0$ at $\Gamma_3$ and considering boundary conditions (B4)–(B7) for the vector function $\mathbf{u}$, we obtain

$$ \int_{\Omega} \langle \mathbf{u}, \nabla q \rangle dx = 0. \quad (B14) $$

We multiply eq. (B3) by a test scalar function $z = z(\mathbf{x})$, $\mathbf{x} \in \Omega$, and integrate by parts the resultant equation over $\Omega$. Considering boundary conditions (B4)–(B7) for the function $T$ and assuming that the function $z$ satisfies the following boundary conditions: $z = 0$ at $\Gamma_1$, $\Gamma_2$, $\Gamma_3$, and $z = 2(k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi)$ at $\Gamma_4$, the modified equation can be presented as

$$ \begin{aligned} & \int_{\Omega} T \{ \nabla \cdot (\kappa(T_{\xi}) \nabla z) - \kappa'(T_{\xi}) \langle \nabla T_{\xi}, \nabla z \rangle + \langle \mathbf{u}_{\xi}, \nabla z \rangle \} dx - \int_{\Omega} \langle \mathbf{u}, \nabla T_{\xi} \rangle z dx \\ & \quad + \int_{\Gamma_4} \left( \kappa(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + \kappa'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) z d\Gamma - \int_{\Gamma_2} \kappa(T_{\xi}) \frac{\partial z}{\partial \mathbf{n}} \chi d\Gamma = o(\|\chi\|). \end{aligned} \quad (B15) $$

Adding eq. (B15) to eq. (B13) and deducting eq. (B14), we obtain

$$ \begin{aligned} & \int_{\Omega} \langle \mathbf{u}, \{ \nabla \cdot (\mu(T_{\xi}) (\nabla \mathbf{w} + \nabla \mathbf{w}^T)) - z \nabla T_{\xi} - \nabla q \} \rangle dx + \int_{\Omega} T \{ \nabla \cdot (\kappa(T_{\xi}) \nabla z) - \kappa'(T_{\xi}) \langle \nabla T_{\xi}, \nabla z \rangle + \langle \mathbf{u}_{\xi}, \nabla z \rangle \} \\ & \quad - \mu'(T_{\xi}) [\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi}] + Ra \langle \mathbf{w}, \mathbf{e}_2 \rangle \} dx + \int_{\Gamma_4} \left( \kappa(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + \kappa'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) z d\Gamma - \int_{\Gamma_2} \kappa(T_{\xi}) \frac{\partial z}{\partial \mathbf{n}} \chi d\Gamma = o(\|\chi\|). \end{aligned} \quad (B16) $$
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Verification of the calculation of the gradient of the cost function $J$.</figcaption></figure>

Now assuming that the expression in braces in eq. (B16) equals to zero, we obtain two equations for $\mathbf{w}$ and $z$, $\mathbf{x} \in \Omega$, as well as the equality for two boundary integrals (to be used for determination of the increment and gradient of the functional):

$$ \nabla \cdot (\mu(T_{\xi}) (\nabla \mathbf{w} + \nabla \mathbf{w}^T)) - z \nabla T_{\xi} - \nabla q = 0, $$

$$ \nabla \cdot (\kappa(T_{\xi})\nabla z) - \kappa'(T_{\xi})\langle\nabla T_{\xi}, \nabla z\rangle + \langle\mathbf{u}_{\xi}, \nabla z\rangle - \mu'(T_{\xi})[\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi}] + Ra \langle\mathbf{w}, \mathbf{e}_2\rangle = 0, $$

and

$$ \int_{\Gamma_4} \left( k(T_{\xi}) \frac{\partial T}{\partial \mathbf{n}} + k'(T_{\xi}) T \frac{\partial T_{\xi}}{\partial \mathbf{n}} \right) z \, d\Gamma = \int_{\Gamma_2} k(T_{\xi}) \frac{\partial z}{\partial \mathbf{n}} \chi \, d\Gamma + o(\|\chi\|). \quad (B17) $$

Finally, we obtain the adjoint problem ($\mathbf{x} \in \Omega$)

$$ \nabla \cdot (\mu(T_{\xi}) (\nabla \mathbf{w} + \nabla \mathbf{w}^T)) = \nabla q + z \nabla T_{\xi}, $$

$$ \nabla \cdot \mathbf{w} = 0, $$

$$ \nabla \cdot (\kappa(T_{\xi})\nabla z) + \langle\mathbf{u}_{\xi}, \nabla z\rangle + Ra \langle\mathbf{w}, \mathbf{e}_2\rangle = \mu'(T_{\xi}) \frac{1}{2} [\nabla \mathbf{w} + \nabla \mathbf{w}^T, \nabla \mathbf{u}_{\xi} + \nabla \mathbf{u}_{\xi}^T] + \kappa'(T_{\xi}) \langle\nabla T_{\xi}, \nabla z\rangle, $$

with the boundary conditions

$$ \Gamma_1 \text{ and } \Gamma_2: \quad z = 0, \quad \mathbf{w} = 0, $$

$$ \Gamma_3: \quad z = 0, \quad \tilde{\mathbf{n}} = 0, \quad q = 0, $$

$$ \Gamma_4: \quad z = 2 \left( k(T_{\xi}) \frac{\partial T_{\xi}}{\partial \mathbf{n}} - \varphi \right), \quad \langle \mathbf{w}, \mathbf{n} \rangle = 0, \quad \tilde{\mathbf{n}} - \langle \tilde{\mathbf{n}}, \mathbf{n} \rangle \mathbf{n} = 0. $$

The derived adjoint problem provides the formulae for increment and gradient of the functional (see eqs A8 and A17):

$$ J(\xi + \chi) - J(\xi) = \int_{\Gamma_2} \chi \nabla J(\xi) d\Gamma + O(\|\chi\|^2), $$

where $\nabla J(\xi) = k(T_{\xi}) \frac{\partial z}{\partial \mathbf{n}}|_{\Gamma_2}$. We have performed the $\chi$-test by Navon et al. (1992) to verify the quality of the gradient of the cost functional with respect to the control variable. For this aim we choose the following increment $\chi = \varepsilon \nabla J(\xi) / \|\nabla J(\xi)\|$, where $\varepsilon$ is small. We rewrite then the last equation introducing a function of $\varepsilon$ as

$$ \varphi(\varepsilon) = \frac{J(\xi + \varepsilon \nabla J(\xi) / \|\nabla J(\xi)\|) - J(\xi)}{\varepsilon \|\nabla J(\xi)\|} = 1 + O(\varepsilon). $$

For values of $\varepsilon$ that are small but not too close to the machine zero, one should expect to obtain a value for $\varphi(\varepsilon)$ that is close to 1. For $\xi = \xi^{(1)} \in \Xi$ (see the main text) the values of $\varphi(\varepsilon)$ are shown in Fig. B1. It is clear that for a value of $\varepsilon$ between $10^{-3}$ and $10^{-10}$, a near unit value of $\varphi(\varepsilon)$ is obtained. This validates the quality of the adjoint model for use in obtaining the gradient of the cost function with respect to the control variable.