A Method for Magma Viscosity Assessment by Lava Dome Morphology

Abstract

Lava domes form when a highly viscous magma erupts on the surface. Several types of lava dome morphology can be distinguished depending on the flow rate and the rheology of magma: obelisks, lava lobes, and endogenic structures. The viscosity of magma nonlinearly depends on the volume fraction of crystals and temperature. Here we present an approach to magma viscosity estimation based on a comparison of observed and simulated morphological forms of lava domes. We consider a two-dimensional axisymmetric model of magma extrusion on the surface and lava dome evolution, and assume that the lava viscosity depends only on the volume fraction of crystals. The crystallization is associated with a growth of the liquidus temperature due to the volatile loss from the magma, and it is determined by the characteristic time of crystal content growth (CCGT) and the discharge rate. Lava domes are modeled using a finite-volume method implemented in Ansys Fluent software for various CCGTs and volcanic vent sizes. For a selected eruption duration a set of morphological shapes of domes (shapes of the interface between lava dome and air) is obtained. Lava dome shapes modeled this way are compared with the observed shape of the lava dome (synthesized in the study by a random modification of one of the calculated shapes). To estimate magma viscosity, the deviation between the observed dome shape and the simulated dome shapes is assessed by three functionals: the symmetric difference, the peak signal-to-noise ratio, and the structural similarity index measure. These functionals are often used in the computer vision and in image processing. Although each functional allows to determine the best fit between the modeled and observed shapes of lava dome, the functional based on the structural similarity index measure performs it better. The viscosity of the observed dome can be then approximated by the viscosity of the modeled dome, which shape fits best the shape of the observed dome. This approach can be extended to three-dimensional case studies to restore the conditions of natural lava dome growth.

INTRODUCTION

Lava domes form as a result of the extrusion of highly viscous magma. It develops a solid surface layer (carapace) remaining mobile and undergoing deformations for days or even months. Several types of lava dome morphology are distinguished. In the endogenous regime, magma intrudes inside the dome without extrusion of fresh magma on the surface. In the exogenous regime, a fresh lava pours out over the surface forming various forms of domes, such as obelisks, lobes, pancake-shaped structures, and some others (Fig. 1). Lava dome collapse can cause explosive eruptions, pyroclastic flows, and lahars, and hence studies of the conditions of lava dome growth is important for hazard assessment and risk reduction.

Fig. 1.

Morphological structures of the lava domes at the Soufrière Hills Volcano, Montserrat (photos by Watts et al., 2002). Panels (a) lava obelisk (about 40 m high and 35 m wide); (b) lava whaleback structure; and (c) lobe structure (about 20 m thick) located on rock debris.

Detailed monitoring of lava domes has been conducted at several volcanoes, such as the Mount St. Helens in the United States (Swanson et al., 1987), Mount Pinatubo in the Philippines (Daag et al., 1996), Mount Unzen in Japan (Nakada et al., 1999), Santiaguito (Santa Maria) in Guatemala (Harris et al., 2003), Merapi and Sinabung in Indonesia (Voight et al., 2000; Nakada et al., 2019), Soufrière Hills on Montserrat (Watts et al., 2002; Wadge et al., 2014), and Volcán de Colima in Mexico (Zobin et al., 2015). Monitoring allows mapping the spatial and temporal development of lava domes, and determining the morphological changes during the dome growth as well as the changes in the lava volume over time (to assess the discharge rate).

The morphology of lava domes is influenced by the rheology of magma and lava discharge rate (DR). Magma viscosity depends on temperature and volume fraction of crystals, which in its turn is determined by crystallization kinetics (the characteristic time of crystal content growth, CCGT (Tsepelev et al., 2020)). At small CCGT values, i.e. fast lava crystallization, obelisk-type structures develop at lower DR and pancake-like structures at higher DR; at high CCGT values, the domes form either lava lobes or pancake-like structures. It was shown that cooling does not play a significant role in the development of the lava dome. If the crystals content is controlled only by the cooling, then the lava viscosity increases in the near-surface layer of the dome, and the thickness of the temperature boundary layer remains small compared to the dome height (Tsepelev et al., 2020). In the dome body, a significant increase in the viscosity occurs due to crystallization caused by a loss of volatiles. Thus, the evolution of the lava dome can be modeled using the rheology depending on CCGT and DR. Meanwhile, the following inverse problem is of an interest to volcanologists: to determine the lava dome viscosity (e.g., the rheological properties of the lava within the dome) by the observed shape of the lava dome for known discharge rate.

Here we propose an approach to solving this inverse problem based on minimizing the deviation between the observed and simulated lava dome shapes. We consider a two-dimensional axisymmetric model of lava dome evolution assuming that the lava viscosity depends only on the volume fraction of crystals, and this fraction, in its turn, depends on the CCGT. Lava domes are modeled numerically at different values of CCGT, DR and the conduit diameter. Using numerical experiments we develop a database of morphological shapes of modeled domes for specified extrusion durations. The results of the experiments (the elements of the database) and an observed dome (in the work, a synthetic dome is constructed to represent the observed dome) are analyzed in the form of two-dimensional images. To estimate the viscosity of the observed lava dome, we minimize the difference between the observed and simulated dome shapes using three different functionals used in computer vision and image processing theory. The viscosity of the observed lava dome is then assessed based on the rheological properties of the modeled lava dome, which shape fits best the shape of the observed dome.

PROBLEM STATEMENT AND THE METHOD OF NUMERICAL MODELING OF LAVA DOMES

We consider a two-dimensional axisymmetric model of two-phase immiscible incompressible fluid approximating the lava (one phase) and the air (another phase). The two phases are separated by a moving interface—the lava dome surface. The influence of the air phase on the lava dome growth is insignificant due to a large ratio between densities/viscosities of the air and the lava. In the model domain Ω (Fig. 2), the lava motion is described by the following set of equations supplemented by the initial and boundary conditions (Ismail-Zadeh and Tackley, 2010; Tsepelev et al., 2019, 2020).

Fig. 2.

Sketch of the model domain. The curved arrow points to the symmetry condition along the x₁=0 axis. ${{\Gamma }_{p}}$ (p = 1, 2, 3, 4, 5, 6) represents a part of the model boundary (see the text for the boundary conditions). The large black arrow shows the part of the boundary, through which magma enters the model domain. The interface between the lava and the air is indicated by a dotted line. The solid dashed line indicates the border of binary images stored in the database.

We use the Navier-Stokes equations with the initial condition u(t = 0, x) = 0 and the continuity equation to describe the lava dynamics

$$\frac{{\partial (\rho {\mathbf{u}})}}{{\partial t}} + \left\langle {{\mathbf{u}},\nabla } \right\rangle (\rho {\mathbf{u}})\, - \,\nabla \cdot \left( {\eta \left( {\nabla {\mathbf{u}}\, + \,\nabla {{{\mathbf{u}}}^{T}}} \right)} \right) = - \nabla p - \rho {\mathbf{g}},$$

(1)

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

(2)

where ${\mathbf{x}} = ({{x}_{1}},{{x}_{2}}) \in \Omega$ are the Cartesian coordinates; t$\in [0,\vartheta ]$ is the time; $\vartheta$ is the final time (the duration of the model experiments); ${\mathbf{u}} = ({{u}_{1}}(t,{\mathbf{x}}),{\text{ }}{{u}_{2}}(t,{\mathbf{x}}))$ is the velocity; $\rho$ is the density; $\eta$ is the viscosity; $p = p({\mathbf{x}})$ is the pressure; ${\mathbf{g}} = (0,g)$, g—is the acceleration due to gravity; $\nabla$, ^T, and $\left\langle {\,\bullet ,\bullet } \right\rangle$ denote the gradient vector, the transposed matrix, and the scalar product of vectors, respectively. We neglect the temperature dependence of the physical parameters of the model, and the surface tension forces. Model density and viscosity are represented as $\rho = {{\rho }_{L}}\alpha (t,{\mathbf{x}}) + {{\rho }_{A}}(1 - \alpha (t,{\mathbf{x}}))$ and $\eta = {{\eta }_{L}}\alpha (t,{\mathbf{x}}) + {{\eta }_{A}}(1 - \alpha (t,{\mathbf{x}}))$, respectively. Here ${{\rho }_{A}}$ is the air density; ${{\rho }_{L}}$ is the lava density; ${{\eta }_{A}}$ is the air viscosity; and ${{\eta }_{L}}$ is the lava viscosity. The function $\alpha (t,{\mathbf{x}})$ equals 1 for the lava and 0 for the air at each point x and at time t, and this function is transported with the velocity u according to the advection equation

$$\frac{{\partial \alpha }}{{\partial t}} + \nabla \cdot (\alpha {\mathbf{u}}) = 0$$

(3)

with the initial condition $\alpha (t = 0,{\mathbf{x}}) = 0$, which means that the entire model domain is filled with the air at the initial time.

We assume that the lava viscosity (measured in Pa s) depends on the volume fraction of crystals (Costa et al., 2009) as

$$\begin{gathered} {{\eta }_{L}}(\varphi ) = {{10}^{7}}\left( {1 + {{\varphi }^{\delta }}} \right) \\ \times \,\,{{\left[ {1 - (1 - \xi ){\text{erf}}\left( {\frac{{\sqrt \pi }}{{2(1 - \xi )}}\varphi (1 + {{\varphi }^{\gamma }})} \right)} \right]}^{{ - B{{\phi }_{*}}}}}, \\ \end{gathered} $$

(4)

where $\varphi = {\phi \mathord{\left/ {\vphantom {\phi {{{\phi }_{*}}}}} \right. \kern-0em} {{{\phi }_{*}}}},$ $\phi$ is the volume fraction of crystals; ${{\phi }_{*}}$ is the specific volume fraction of crystals; B is the Einstein coefficient’s theoretical value determined from the Einstein equation as $B = {{({{\eta }_{L}}(\phi ) - 1)} \mathord{\left/ {\vphantom {{({{\eta }_{L}}(\phi ) - 1)} \phi }} \right. \kern-0em} \phi }$ (Mardles, 1940) (it was experimentally determined that the Einstein coefficient varies from 1.5 to 5 (Jeffrey and Acrivos, 1976)); $\delta = 7.24$, $\gamma = 5.{\text{76}}$, and $\xi = 4.63 \times {{10}^{{ - 4}}}$ (Lejeune and Richet, 1995; Costa et al., 2009); ${\text{erf}}( \cdot )$ is the error function. The volume fraction of crystals is determined from the following evolutionary equation describing the simplified kinetics of crystal growth during crystallization due to magma degassing

$$\frac{{\partial \phi }}{{\partial t}} + \nabla \cdot (\alpha \phi {\mathbf{u}}) = - \alpha \frac{{\phi - {{\phi }_{{eq}}}}}{\tau },$$

(5)

with the initial condition $\phi (t = 0,{\mathbf{x}}) = 0$. Here ${{\phi }_{{eq}}}$ is the volume fraction of crystals at the equilibrium, which depends on the amount of water dissolved in the magma and on the temperature; $\tau$ is the CCGT. The smaller CCGT, the faster the crystallization process converges to its equilibrium state. CCGT is referred to as the relaxation time required to decrease the difference between actual ($\phi$) and equilibrium (${{\phi }_{{eq}}}$) values of the volume fraction of crystals by a factor of e (~2.72) relative to the initial difference (${{\phi }_{{in}}} - {{\phi }_{{eq}}}$), where ${{\phi }_{{in}}}$ is the initial volume fraction of crystals in the magma. For $\alpha = 1$ and u = 0, the CCGT can be found analytically as $\tau = - t{{\left( {\ln \left[ {{{({{\phi }_{{eq}}} - \phi )} \mathord{\left/ {\vphantom {{({{\phi }_{{eq}}} - \phi )} {({{\phi }_{{eq}}} - {{\phi }_{{in}}})}}} \right. \kern-0em} {({{\phi }_{{eq}}} - {{\phi }_{{in}}})}}} \right]} \right)}^{{ - 1}}}$. The relaxation time $\tau$ can vary from a few hours to several months depending on the temperature and water saturation of the magma, the number of pre-existing crystals in the magma and its composition (Tsepelev et al., 2020). Note that although the viscosity depends on the petrological (chemical) composition of the lava and the volatile content of the lava (its water saturation), these viscosity dependencies are not considered in this paper.

The following conditions are set on the boundary $\Gamma = {{\Gamma }_{1}} \cup {{\Gamma }_{2}} \cup {{\Gamma }_{3}} \cup {{\Gamma }_{4}} \cup {{\Gamma }_{5}} \cup {{\Gamma }_{6}}$ of the model domain (see Fig. 2). At the boundary ${{\Gamma }_{1}}$, the symmetry conditions are set, i.e., the impermeability condition $\left\langle {{\mathbf{u}},{\mathbf{n}}} \right\rangle = 0$ and the free slip condition $\left( {\nabla {\mathbf{u}} + \nabla {{{\mathbf{u}}}^{T}}} \right){\mathbf{n}} - \left\langle {\left( {\nabla {\mathbf{u}} + \nabla {{{\mathbf{u}}}^{T}}} \right)\,{\mathbf{n}},{\mathbf{n}}} \right\rangle {\mathbf{n}} = 0$. It is assumed that the lava of density ${{\rho }_{L}}$ and viscosity ${{\eta }_{L}}$ enters the model domain through the boundary ${{\Gamma }_{2}}$ at the given DR Q. At the boundaries ${{\Gamma }_{3}}$, ${{\Gamma }_{4}}$ and ${{\Gamma }_{5}}$, no-slip condition u = 0 is assumed. The air is removed from the model domain through the boundary ${{\Gamma }_{6}}$ according to the given DR. It is assumed that the volume fraction of crystals is equal $\phi = {{\phi }_{{in}}}$ at the boundary ${{\Gamma }_{2}}$ and $\phi = 0$ at ${{\Gamma }_{6}}$. The values of model parameters used in numerical simulations are presented in Table 1.

Table 1.   Model notations, parameters, and their values

The finite volume method implemented in Ansys Fluent software is employed to solve numerically the problem (1)-(5). To determine the position of the lava dome interface with the air, the volume of fluid (VOF) method is used (Hirt and Nichols, 1981). An implicit integration scheme is used to solve equations (1)–(5) with the appropriate boundary and initial conditions. The pressure is discretized by the second-order PRESTO! scheme (Peyret, 1996). We use a numerical scheme of the second-order accuracy to approximate the Laplace operator, and monotone schemes to discretize the convective terms in the equations (Ismail-Zadeh and Tackley, 2010). The SIMPLE method (Patankar and Spalding, 1972) is used to solve equations (1)–(2), where the relaxation parameters are chosen to be 0.01 and 0.3 for velocity and pressure, respectively. The relaxation parameter is chosen to be 0.5 for the function $\alpha$ and the volume fraction of crystals. Given the large jump between the lava and air viscosities, the choice of the relaxation parameters is critical. The time step is chosen in the range from 0.1 to 1 s to optimize the speed and to ensure the convergence of the solution to the system of linear algebraic equations (SLAE) obtained after the discretization of Eq. (1). The accuracy in solving the SLAE for the function $\alpha$ and the volume fraction of crystal is 10^-6.

DATABASE OF THE MORPHOLOGICAL SHAPES OF MODELED LAVA DOMES. CONSTRUCTION OF A SYNTHETIC DOME SHAPE

We develop a database of morphological shapes of lava domes using the following model of lava dome evolution. The model domain (see Fig. 2) is discretized by about 70000 hexagonal cells. We assign the following values for the CCGT $\tau = 5 \times {{10}^{4}}$ s, DR Q = 0.7 m³ s^–1, and the radius of the volcanic eruption vent r = 15 m. Using the parameters specified in Table 1, we numerically solve the problem (1)–(5) within the time intervals specified in Table 2 using the Ansys Fluent software.

Table 2.   Parameters of modeled lava domes

We consider a rectangular domain ${{\Omega }_{{ABCD}}} \subset \Omega$ with vertices A = (0,30), B = (100,30), C = (100,100), and D = (0,100), and denote $L = L(\tau ,Q,r,t) \subset {{\Omega }_{{ABCD}}}$ as the sub-domain containing the lava dome and $L{\kern 1pt} ' = {{\Omega }_{{ABCD}}}\backslash L$ as the sub-domain containing the air. We refer to the boundary $F = F(\tau ,Q,r,t)$ = $\left\{ {({{x}_{1}},{{x}_{2}}) \in \partial L\backslash \partial {{\Omega }_{{ABCD}}}} \right\}$ as the boundary (or the morphological shape) of the lava dome.

Using the results of the numerical experiments, we obtain the morphological shape ${{F}_{k}} = F(\tau ,Q,r,t = {{t}_{k}})$ for the time ${{t}_{k}}$ and place it in the database. To conduct the test study, we computed the morphological shapes for the following values of $\tau = \{ 1.8 \times {{10}^{4}},{\text{ }}6 \times {{10}^{4}},{\text{ }}7 \times {{10}^{4}},{\text{ }}8 \times {{10}^{4}},{\text{ }}5 \times {{10}^{5}}\}$ s and r = 5 m (see Table 2). The database can be replenished with modeled dome shapes for different parameters of CCGT, DR, and the vent’s radius. Figure 3 shows several morphological shapes of the domes obtained in numerical simulations for different parameters of the problem. Note that the area L (of each lava dome presented in Fig. 3) is approximately the same.

Fig. 3.

Morphological shapes of several modeled lava domes.

As the main goal of the study is to find a dome shape in the database that approximates a natural (observed) lava dome in the best way and to determine then its viscosity, we construct a synthetic dome as an example of a natural lava dome. To construct the synthetic dome, we choose a dome from the database, for instance, the dome with the morphological shape ${{F}_{{135}}} = F(\tau = 6 \times {{10}^{4}}\,\,{\text{s}},$ $Q = 0.7\,\,{{{\text{m}}}^{3}}\,{{{\text{s}}}^{{ - 1}}},$ $r = 15\,\,{\text{m}},$ $t = 5 \times {{10}^{4}}\,\,{\text{s}})$ (Fig. 4a) and introduce a random noise along the boundary of this dome to develop a synthetic dome shape F* (see Fig. 4b).

Fig. 4.

Morphological shapes of the modeled dome F = F₁₃₅ (a) and the synthetic dome (b).

The random noise is generated in the following way. Consider some point $a = ({{a}_{1}},{{a}_{2}}) \in F$ on the lava dome boundary. This point is randomly displaced as $a + \delta = ({{a}_{1}} + {{\delta }_{1}},{{a}_{2}} + {{\delta }_{2}})$, where ${{\delta }_{1}}$ and ${{\delta }_{2}}$ are normally distributed random variables with zero mathematical expectation and standard deviation equal to 1 (Wentzel, 1969). We consider the circle ${{C}_{{a + \delta }}}\left( {\left| \varepsilon \right|} \right)$ with the center at the point $a + \delta$ and radius |ε| for each displaced point, where $\varepsilon$ is a uniformly distributed random variable taking values on the interval $\left[ { - \omega ,\omega } \right]$, $\omega$ is a constant and equal to 1 in this case. Then we construct the following sub-domain:

$$L\text{*} = \left\{ {\begin{array}{*{20}{c}} {L \cup \left( {\bigcup\limits_{a \in F} {{{C}_{{a + {{\delta }}}}}\left( {\left| \varepsilon \right|} \right),\,\,\varepsilon \geqslant 0} } \right)} \\ {L\backslash \left( {\bigcup\limits_{a \in F} {{{C}_{{a + {{\delta }}}}}\left( {\left| \varepsilon \right|} \right),\,\,\varepsilon < 0} } \right)} \end{array}} \right..$$

Thus, the synthetic dome with boundary F* is obtained by combining the dome with boundary F with a set of random circles in the case of non-negative $\varepsilon$ and by truncating the synthetic dome with a set of random circles in the case of negative $\varepsilon$. The artificial distortion (noise) of the dome boundary allows for simulating real distortions caused by the growth of a natural lava dome and its partial collapse and/or by errors in measurements of the morphological shape of a lava dome.

PATTERN RECOGNITION METHODS

We compare the shape of the synthetic dome F* with the shapes of lava domes F_k from the database by the methods used in the theory of image processing and pattern recognition (e.g., Salomon, 2007). Here we represent lava domes by binary images and introduce a uniform rectangular partitioning of the area ${{\Omega }_{{ABCD}}}$ into cells $I \times J$ as ${{\Omega }_{{ij}}}\left( {{{\Omega }_{{ABCD}}} = \bigcup\nolimits_{i = 0,j = 0}^{I - 1,J - 1} {{{\Omega }_{{ij}}}} } \right)$. A rectangular matrix $P(F) = \{ {{p}_{{ij}}}\} _{{i = 0,j = 0}}^{{I - 1,J - 1}}$ of size $I \times J$ is assigned to each shape of the dome F, where the matrix element ${{p}_{{ij}}}$ equals to 0, if the corresponding cell Ω_ij contains more than 50 percent of the air, and equals to 1 in all other cases. We evaluate then the closeness of the synthetic dome shape F^* and the arbitrary dome shape F_k from the morphological shape database by means of the following quality functions.

1. The functional based on the symmetric difference:

$${{J}_{1}}(F\text{*},F) = {{k}_{1}}S((L{\text{*}} \cup L)\backslash (L{\text{*}} \cap L)),$$

where $S( \cdot )$ (m²) is the area of the region and k₁ (m^–2) is the scaling multiplier.

2. The functional based on peak signal-to-noise ratio (PSNR) measure (Salomon, 2007):

$$\begin{gathered} {{J}_{2}}(F\text{*},F) = \\ = {{k}_{2}}\left( {{{k}_{3}} + 10{{{\log }}_{{10}}}\left[ {{{\sum\limits_{i = 0,j = 0}^{I - 1,J - 1} {{{{\left( {{{p}_{{ij}}} - p_{{ij}}^{*}} \right)}}^{2}}} } \mathord{\left/ {\vphantom {{\sum\limits_{i = 0,j = 0}^{I - 1,J - 1} {{{{\left( {{{p}_{{ij}}} - p_{{ij}}^{*}} \right)}}^{2}}} } {(IJ)}}} \right. \kern-0em} {(IJ)}}} \right]} \right), \\ \end{gathered} $$

where k₂ is the scaling factor, and k₃ is a positive constant. In numerical implementations, if P(F) and P(F^*) match completely, the user receives a message containing the number of the modeled dome, where the condition of the complete match between the two matrices is reached.

3. The functional based on the SSIM (structure similarity index measure) (Wang et al., 2004):

$$\begin{gathered} {{J}_{3}}(F\text{*},F) = {{k}_{4}}\left( {2\mu (P)\mu (P\text{*}) + {{c}_{1}}} \right) \\ \times \,\,{{\left( {\sigma (P,P\text{*})\, + \,{{c}_{2}}} \right)} \mathord{\left/ {\vphantom {{\left( {\sigma (P,P\text{*})\, + \,{{c}_{2}}} \right)} {\left[ {\left( {{{\mu }^{2}}(P)\, + \,{{\mu }^{2}}(P\text{*})} \right)} \right.}}} \right. \kern-0em} {\left[ {\left( {{{\mu }^{2}}(P)\, + \,{{\mu }^{2}}(P\text{*})} \right)} \right.}} \\ \times \,\,\left. {\left( {{{\sigma }^{2}}(P)\, + \,{{\sigma }^{2}}(P\text{*})\, + \,{{c}_{2}}} \right)} \right], \\ \end{gathered} $$

where μ(P), μ(P*) is the mathematical expectation, σ(P,P*) is the covariance, σ²(P), σ²(P*) is the dispersion, k₄ is the scaling multiplier, c₁ = 0.01, and c₃ = 0.03. Here we consider a probabilistic model of image representation, namely, the image is considered as a field of random variables, and the value at each point of this field is a realization of a random variable.

The values of the functionals are calculated for each element of the database, and the obtained set of values is ordered in the descending order. In this case, the fact that one or another functional gives a more accurate estimate is based on the information about the synthetic dome. For example, Figure 5 showing the values of the three functionals on the elements presented in Fig. 3 illustrates that all three functionals reach a minimum at modeled dome #135. Note that functionals J₁, J₂, and J₃ estimate the quantitative deviation of the modeled and synthetic domes, while the functional J₃ also estimates the structural features of the morphological shapes of the domes, although it leads to complex and resource-consuming computations.

Fig. 5.

The values of functionals J₁, J₂, and J₃ related to the modeled domes (see Fig. 3) and presented in descending order.

RESULTS: DETERMINATION OF THE LAVA DOME VISCOSITY FROM THE KNOWN DOME SHAPE

Using the functionals described above, we compare the synthetic dome F* with all domes F_k from the database of morphological forms. The following scaling multipliers are prescribed to get the values of all functionals on the database’s elements within the interval [0, 1]: k₁ = 1/30, k₂ = 1/25, k₃ = 33, and k₄ = 1/14. Figures 6 and 7 present the values of the functionals J_m(F*, F), m = 1, 2, and 3 versus the k numbers of the domes from the database. The set of elements, on which the smallest values of the considered functionals are achieved, are almost identical: the modeled domes F₄₁, F₁₃₃, F₁₃₄, and F₁₃₅ are the closest to the synthetic dome F^* (see Fig. 7). Figure 8 shows the synthetic dome and four closest modeled domes from the database. The time of a lava dome formation, the discharge rate, and the vent’s size can help with a practical selection of the closest modeled dome.

Fig. 6.

The values of the functionals presented in descending order: ${{J}_{1}}(F\text{*},{{F}_{k}})$ (a), ${{J}_{2}}(F\text{*},{{F}_{k}})$ (b), and ${{J}_{3}}(F\text{*},{{F}_{k}})$ (c).

Fig. 7.

Fragment of Fig. 6 for small values of the functionals ${{J}_{1}}(F\text{*},{{F}_{k}})$ (a), ${{J}_{2}}(F\text{*},{{F}_{k}})$ (b), and ${{J}_{3}}(F\text{*},{{F}_{k}})$ (c).

Fig. 8.

Comparison of synthetic dome F* with modeled domes F₁₃₅, F₁₃₄, F₁₃₃ and F₄₁. Colors indicate the modeled lava viscosity distribution.

From the available information about the synthetic dome F^* (obtained by the introduction of a noise on the morphological shape of the modeled dome F₁₃₅) one can deduce that the SSIM-based functional provides a more accurate closeness estimate. Note that the functionals J₁ and J₂ give qualitatively similar results for the binary images, and hence only one of them can be considered. To consider qualitative and quantitative closeness of the synthetic and modeled domes simultaneously, a linear combination of the described functionals can be used, and in this case optimal qualitative and quantitative estimates can be achieved by a suitable choice of weighting coefficients. For example, considering the functional $J(F\text{*},F)$ = $0.5\left[ {{{J}_{1}}(F\text{*},F) + {{J}_{3}}(F\text{*},F)} \right]$, we see that the minimum of the functional is achieved for dome #135 (Fig. 9). Although the linear combination of the functionals can be used to choose the modeled dome that will optimally fit the natural dome, it should be noted that the choice of the weighting coefficients may be difficult in practice when the natural (and not synthetic) dome is to be assessed.

Fig. 9.

The values of functional $J(F{\text{*}},F)$ presented in descending order.

DISCUSSION AND CONCLUSION

Lava dome morphological diversity can be explained by changes in the magma viscosity caused by degassing and crystallization during the magma ascent through the volcanic conduit from the magma chamber (Melnik and Sparks, 1999, 2005). The rapid ascent of magma reduces the time available for magma crystallization, so magma behaves as a relatively low-viscosity fluid and, after extrusion to the surface, spreads out to form pancake-like morphological structures. With decreasing DR or at small CCGT, the crystalline magma undergoes a transition from a fluid to a quasi-solid state forming lava lobes and obelisks on the surface (Tsepelev et al., 2020).

The magma viscosity and the history of a lava dome growth can be determined from observations on the morphological shapes of the dome by solving inverse problems. This paper presents an approach to estimating the magma viscosity, based on the comparison of observed and modeled morphological shapes of lava domes, by three functionals used in the computer vision and image processing theory. Although each functional allows for determining the minimal difference between the modeled and observed morphological shapes of lava domes, the SSIM-based functional evaluates not only the quantitative deviation of the modeled and observed domes, but also their structural features, and thus performs a better assessment. The rheological parameters of the modeled dome, which morphological shape fits that of the observed lava dome, can be then estimated and adopted for the viscosity of the observed dome.

An observed lava dome was synthesized from a random modification of one of the modeled lava domes, and the functional estimates were based on the knowledge of the morphological shape of this synthetic dome. In the case of a natural lava dome, we would recommend considering the element of the database, at which the values of the functionals are minimal, along with a few other elements close to it. The obtained sample of the elements is to be subjected to an expert evaluation.

The inverse problem under consideration is a rather complicated object of research from both theoretical and computational points of view (Samarskii and Vabishchevich, 2004; Kabanikhin, 2009). As a rule, an inverse problem is ill-posed, that is, its solution is unstable and/or non-unique. The considered inverse problem is ill-posed, namely it is unstable (as small errors in the initial data or rounding errors can lead to significant errors in the solution of the problem), and also it has a non-unique solution (for example, different DR and different CCGT can result in a similar shape of the lava dome). However, the discharge rate can be practically determined by the dome volume, and therefore DR can be considered as a known characteristic of the process. Classical methods are not suitable for solving inverse problems, and various approaches are used to transform ill-posed problems into well-posed (when the solution is stable and unique) (Tikhonov and Arsenin, 1977). Many of these approaches may be satisfactory in the theoretical studies of the inverse problem, but unsatisfactory in its numerical solution. Other approaches involve various simplifications of inverse problems.

In this work, the inverse problem is related to a search for the lava dome viscosity distribution, which is presented as a certain function with a small number of model parameters. Actually, this search is reduced to determining CCGT, on which the viscosity depends, and this has been performed by minimizing the differences between the morphological shapes of the observed (synthetic, in this case) and modeled lava domes.

Prediction of a dome shape and stress distribution within the dome structure may help in assessments of its stability and possible collapse with the formation of pyroclastic flows or explosive eruptions. Ideally, a customized volcano-specific program should be used for short- and long-term forecasts by volcano observatory staff. A similar practice is used by oil companies to forecast oil production in the fields. A regular adjustment of the model by reproducing the history of oil production permits for reliable forecasts in this case and for determining optimal production strategies.

When applying this method to a real eruption, it is possible to consider deviations of the calculated dome shape from the measured one not only at a certain point in time, but also over the entire observation period, setting limitations on the crater shape, the conduit feeding the dome, and the magma petrology. Natural lava domes are three-dimensional objects. Although the approach presented in the paper is two-dimensional, it can be extended to the three-dimensional case and used to reconstruct the growth conditions of natural lava domes.