Lava Dome Morphology and Viscosity Inferred From Data-Driven Numerical Modeling of Dome Growth at Volcán de Colima, Mexico During 2007-2009

Magma extrusion, lava dome growth, collapse of domes, and associated pyroclastic flow hazards are among important volcanological studies. In this paper, we analyze the influence of the magma viscosity and discharge rates on the lava dome morphology at Volcán de Colima in Mexico during a long dome-building episode lasting from early 2007 to fall 2009 without explosive dome destruction. Camera images of the lava dome growth together with recorded volumes of the erupted lava have been used to constrain numerical modeling and hence to match the history of the dome growth by nudging model forecasts to observations. Our viscosity model incorporates crystal growth kinetics and depends on the characteristic time of crystal content growth (or CCGT) and the crystal-free magma viscosity. Initially, we analyze how this viscosity, CCGT, and the rate of lava extrusion influence the morphology of the growing dome. Several model scenarios of lava dome growth are then considered depending on the crater geometry, the conduit location, the effective viscosity of dome carapace, and the extrusion rates. These rates are determined either empirically by optimizing the fit between the morphological shape of modeled domes and that of the observed dome or from the recorded lava dome volumes. The maximum height of the modeled lava dome and its horizontal extent are in a good agreement with observations in the case of the empirically-derived extrusion rates. It is shown that the topography of the crater at Volcán de Colima is likely to be inclined toward the west. The viscosity of the modeled lava dome (∼1012 Pa s) is in a good agreement with the effective viscosity estimated experimentally from lavas of Volcán de Colima. Due to the interplay between the lava extrusion and the gravity forces, the dome reaches a height threshold, and after that a horizontal gravity spreading starts to play an essential role in the lava dome evolution. The model forecasts that the dome carapace of higher viscosity (∼1014 Pa s) influences the dome growth and its morphology during long dome-building episodes by retarding horizontal advancement and developing steep-sided eastern edge of the dome at the volcano. The developed model can be used in assessments of future effusive eruptions and lava dome growth at Volcán de Colima or elsewhere. History matching modeling of lava dome growth sheds a light on dynamic processes inside the dome and may assist in assessing stress state in the dome carapace and in forecasting the dome failures.


INTRODUCTION
Lava domes grow by the extrusion of viscous magma from a volcanic conduit (e.g., Calder et al., 2015). The lavas of domebuilding volcanoes have low average eruption rates and high viscosities that are commonly associated with high groundmass crystallinity and, consequently, substantial yield strength (Lavallée et al., 2008;Lavallée et al., 2008;Sheldrake et al., 2016). Nonlinear dynamics of lava dome growth is caused by crystallization and outgassing of the highly viscous lava in the volcanic conduit (Melnik and Sparks, 1999). Through the intermittent buildup of gas pressure, growing domes can often experience episodes of explosive activity (e.g., Voight and Elsworth, 2000;Heap et al., 2019).
The morphology of lava domes is controlled by the magma rheology, the crater topography, the ascent dynamics, and the mechanism of dome growth (e.g., Fink and Griffiths, 1998;Tsepelev et al., 2020). Domes may reach heights of several hundred meters, and they may grow rapidly for some days or slowly and steadily for months to years. Dome growth in the crater may alternate between endogenous to exogenous styles. Endogenous growth refers to the enlargement of a lava dome due to expansion caused by intrusion of new magma. Exogenous growth refers to dome enlargement as a result of magma forcing its way through a preexisting lava dome carapace to the surface or flowing directly from the vent and forming discrete lobes of lava that pile on top of, or adjacent to, each other (Calder et al., 2002;Watts et al., 2002;Harris et al., 2003;Simmons et al., 2005;Calder et al., 2015;Rhodes et al., 2018;Harnett et al., 2019).
Numerical models of lava dome growth clarify how the rheological properties of magma within a lava dome influence morphological features of the dome and what are the effects of crystal content, temperature, pressure, and the discharge rate on the morphology (Hale and Wadge, 2003;Husain et al., 2014;Harnett et al., 2018;Tsepelev et al., 2020). Numerical modeling of lava dome growth allows reconstructing the process of filling of the crater with magma and estimating the lava viscosity, the characteristic time of crystal content growth, and the extrusion rate (Tsepelev et al., 2020;Starodubtseva et al., 2021). This modeling requires information about the dynamics of lava dome growth collected during monitoring of the eruption. Here, we present a quantitative modeling study of the 2007-2009 endogenous stage of lava dome growth at Volcán de Colima, México, based on data collected by the University of Colima during the monitoring of the dome growth process. The goal of this study is to understand how the viscosity of magma/lava, the lava discharge rate, the conduit location, and the topography of the crater influence the dome morphology at Volcán de Colima. Compared to other numerical models of dome growth at this volcano, which are characterized by a rapid dome building (e.g., Walter et al., 2019), this modeling study deals with a slow domebuilding process without explosive dome destruction.
VOLCÁN DE COLIMA AND ITS LAVA DOME GROWTH DURING 2007-2011 The andesitic 3860-m-high stratovolcano Volcán de Colima is the most active volcano in Mexico (Luhr and Carmichael, 1990). It is located in the western part of the Trans-Mexican Volcanic Belt ( Figure 1A), and together with the Pleistocene volcano Nevado de Colima, they form the Colima Volcanic Complex. The recent dome-building eruption at Volcán de Colima began in November 1998 and was preceded by five swarms of volcano-tectonic earthquakes recorded from November 1997. The earthquake swarms indicated the pathway used by juvenile magma to ascend into the shallow magma storage zone, which is located at depths between 4 and 12 km beneath the crater (Zobin et al., 2002). This batch of magma fed seven lava dome-building episodes occurring during the 1998-2017 eruptions. The duration of the growth of lava domes varied from 14 days (rapid growth) to 33 months (slow growth) (Zobin and Tellez, 2019). In this paper we present the results of the numerical study of the long lava dome-building episode during 2007-2009, which is characterized by an uninterrupted dome extrusion. This removes complications associated with lava dome destruction by explosions or collapses, which are difficult to incorporate in a single mathematical/numerical model, and permits to analyze a dome extrusion evolving over several years.
The photos in Figures 1B-D illustrate the lava dome growth within the crater during [2007][2008][2009]. The appearance of a new lava dome was observed at the beginning of February 2007 ( Figure 1B; Zobin et al., 2017). As shown in Figure 1D, the crater was filled by a lava dome by the end of November 2009. In February 2010, the initially endogenous process of the dome growth was accompanied by exogenous growth of a small lobe of lava within the western sector of the crater of the volcano. The simultaneous endogenous and exogenous growth of the dome continued until June 2011 (Zobin et al., 2017). Video monitoring of lava dome growth during 2007-2011, as well as photos taken with a camera located at a distance of about 5,700 m from the center of the lava dome (Bretón-González et al., 2013), allowed for outlining the morphology of the growing dome in detail ( Figure 1E).
Dome growth was monitored by photogrammetry tools (aerial photos) and by the related seismic recording (Zobin et al., 2017). Visual inspection of aerial photos was used to calculate the volume of the lava dome (Figure 2A), and the average discharge rates were then estimated ( Figure 2B). The average discharge rate until the beginning of 2009 was about 0.02 m 3 s −1 .
By March 2009, lava dome volume reached about 1,400,000 m 3 , filling roughly 80% of the total crater. In fall 2009, the dome reached the western rim of the crater as shown in Figure 1D ( Zobin et al., 2017). Figure 2C shows the temporal variations in the number of seismic signals generated by small gas-and-ash explosions (Zobin et al., 2015). Appearance of explosive events may serve as an indicator of magma degassing during the lava dome growth. Three periods of the degassing process can be distinguished. During February-September 2007, the number of explosions sharply increased demonstrating the activation of degassing at the initial stage of the lava dome growth. Then, during the ongoing stage of endogenous growth of the dome between October 2007 and January 2009, the level of degassing stayed at rather high level. Later the level of degassing continuously decreased up to the end of the lava dome growth in June 2011 (Zobin et al., 2017). The maximum likelihood regression lines shown in Figure 2C illustrate the tendency in the variation of the number of explosions during three periods of lava dome growth.

STATEMENT OF THE PROBLEM AND MODEL SETUP Geometry
A two-dimensional model geometry presented in Figure 3A approximates the observed geometry of the crater rim of Volcán de Colima across the E-W profile ( Figure 1E) with an assumption that the lava dome base is located 20 m below the central part of the visible crater's rim (Bretón-González et al., 2013). Although the crater geometry varies in numerical experiments, the following parameters of the model domain remain the same: the width of the crater is 257 m, the width and the depth of the conduit are 14.5 and 30 m, respectively (Zobin et al., 2015). In model domain Ω ( Figure 3A), we consider a two-dimensional, viscous, incompressible, two-phase, immiscible fluid flow approximating an extrusion of lava (one fluid phase) into the air (another fluid phase) on the surface of the crater of Volcán de Colima.

Governing Equations With the Initial and Boundary Conditions
A lava dome growth is described mathematically by the Navier-Stokes equations (Eq. 1) with the initial condition u(t 0, x) u 0 , the continuity equation (Eq. 2), and the advection equation (Eq. 3) for the interface between the air and the lava: where x (x 1 , x 2 ) ∈ Ω are the Cartesian coordinates; t is the time; u (u 1 (t, x), u 2 (t, x)) is the velocity; ρ is the density; η is the viscosity; p p(x) is the pressure; g (0, g), and g ( 9.81 m s −2 ) is the acceleration due to gravity; α(t, x) takes the value of one for the lava and zero for the air at each point x and time t; ∇, ∇·, T , and 〈·, ·〉 denote the gradient operator, the divergence operator, the transposed matrix, and the scalar product of vectors, respectively. We consider that the model density and viscosity are presented as ρ ρ L α(t, x)), respectively, where ρ A 1.225 kg m −3 and η A 10 −3 Pa s are the density and viscosity of the air; ρ L 2,500 kg m −3 and η L are the lava density and viscosity, respectively. At the initial time, the sub-domain Ω 1 (the model conduit) is filled by an older magma with high volume fraction of crystals (ϕ 0.8), and hence α(t 0, x) 0 for x ∈ Ω 2 and 1 for x ∈ Ω 2 . The initial velocity field u 0 is chosen so as to generate a laminar flow in the entire model domain, including the part filled by the air (see Appendix A1). In the numerical modeling, we neglect the dependence of the density on temperature and/or phase transformation due to crystallization or volatile exsolution. The following conditions are imposed on the model boundary Γ Γ 1 ∪ Γ 2 ∪ Γ 3 (Figure 3). A new lava (with the volume fraction of crystals ϕ ϕ in ) enters into the model domain at a specified extrusion rate u ext (0, u ext ) through the part Γ 1 of the model boundary. The lava extrusion rate (m s −1 ) in the modeling is determined either 1) empirically from camera images of the morphological shapes taken during lava dome growth from 2007 to 2009 ( Figure 1E) or 2) from the observed volumes of extruded lava and calculated discharge rates (m 3 s −1 ) (see Appendix A2). In case (i), we vary the model extrusion rate and choose one that provides the best fit between the morphological shapes of modeled and observed domes at specified times. A no-slip condition u 0 is prescribed at Γ 2 . The outflow conditions are determined at Γ 3 (the blue curve) by removing the air from the model domain proportional to the lava discharge rate and to guarantee the condition of incompressibility: u out u out n and u out −|Γ 3 | −1 Γ 1 〈u ext , n〉dΓ, where n is the outward unit normal vector at a point on the model boundary.

Lava Viscosity
Degassing-induced (i.e., temperature independent) crystallization is the main driving force in the bulk viscosity increase in conduits and lava domes compared to the coolinginduced viscosity increases of the melt and crystallization (e.g., Chevrel et al., 2013;Tsepelev et al., 2020). Meanwhile, a coolinginduced increase in the viscosity of the dome carapace becomes important during long episodes of lava dome growth. Hence, although the lava dome rheology is more complex, we assume here that the lava viscosity η L depends on the volume fraction of crystals only (Costa et al., 2009): where η p is the crystal-free magma viscosity; φ ϕ/ϕ p , ϕ is the volume fraction of crystals; ϕ p is the specific volume fraction of crystals, representing the critical solid fraction at the onset of the exponential increase of the lava viscosity (here we assume ϕ p 0.384); the empirical parameters δ 7.24, ξ 4.63 × 10 −4 , and c 5.76 are taken from Lejeune and Richet (1995) and Costa et al. (2009); and erf (·) is the error function. The coefficient B is determined from the Einstein equation as B (η L (ϕ) − 1)/ϕ (Mardles, 1940). It was experimentally determined that this coefficient varies from 1.5 to 5 (Jeffrey and Acrivos, 1976); and we assume B 2.5 (Costa et al., 2009). The volume fraction of crystals ϕ is determined from the evolutionary equation describing the simplified crystal content growth kinetics of degassing-induced crystallization (e.g., Tsepelev et al., 2020): with the initial condition for the volume fraction of crystals ϕ(t 0, x) 0 for x ∈ Ω 2 and ϕ(t 0, x) 0.8 for x ∈ Ω 1 . Here, ϕ eq is the volume fraction of crystals at the equilibrium state, which depends on the fraction of water dissolved in magma and temperature; τ is the characteristic time of crystal content growth (CCGT) needed by crystals to reach ϕ eq . We assume that ϕ eq 0.83, which is in the range of experimental constraints (Riker et al., 2015;Cashman, 2020). The smaller CCGT, the faster the crystallization process converges to its equilibrium state. CCGT is referred to as the relaxation time (Tsepelev et al., 2020), which is required to reduce the difference between the actual (ϕ) and equilibrium (ϕ eq ) values of the volume fractions of crystals by a factor of e with respect to the difference (ϕ en − ϕ eq ), where ϕ en is the volume fraction of crystals in the magma entering the model conduit at Γ 1 (in the modeling we assume that ϕ en 0.4).

METHOD
We consider the problem of lava dome viscosity determination from observations of the dome morphology (e.g., the observed height and horizontal extent of the dome). The relevant mathematical problem belongs to the class of inverse problems. There are several methods to solve the problem, that is, to determine the distribution of lava dome viscosity Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 735914 from relevant observations. One of the methods is based on optimization using data assimilation techniques, and helps solving the inverse problem (and determining the viscosity) by organizing an iteration process to minimize the difference between the model solution and observations (e.g., Ismail-Zadeh et al., 2016;Korotkii et al., 2016). This method involves an analytical derivation of the optimization (adjoint) problem, and hence it is restricted to simple cases. Another method is based on the use of image processing and machine learning to determine the viscosity of a lava dome by analyzing images of the observed and modeled morphological shapes of the dome at each time of lava dome growth (Starodubtseva et al., 2021). An inverse problem can be replaced by the direct problem conjugated with the inverse problem, and the direct problem can be solved by varying model parameters and fitting observations. Although this method is simpler than the previous two approaches, it requires a significant number of numerical experiments to fit the observations. In this work we use the latter method, namely, we solve numerically Eqs 1-5 with the relevant initial and boundary conditions and vary the model parameters to get the best fit to the height and width of the lava dome during its growth.
To solve this problem numerically, we employ the Ansys Fluent software, where the finite volume method is used to solve numerical models on multiprocessor computers. The software is based on the volume of fluid (VOF) method (Hirt and Nichols, 1981) allowing for computationally inexpensive treatment of a moving interface between two fluid phases, e.g., the lava and the air. The cells containing the interface between the lava and air have α values between zero and one depending on the lava proportion in the cells. Because of the large discontinuity between the lava and air viscosity, the interface of the two fluids does not represent a sharp boundary, and some smearing can be observed during computations. More detail about the numerical approach can be found in Appendix A1.

Sensitivity Analysis
We develop initially a series of numerical experiments of lava dome growth varying the model parameters to understand how each of them influences the morphological shape of the lava dome (e.g., its height and width). Particularly, we have performed a sensitivity analysis with respect to the crystal-free magma viscosity η p , CCGT τ, and extrusion rate u ext . Experiments 1: We study the influence of the crystal-free magma viscosity η p on the dome morphology. In the experiments, this factor varies from 1.6 × 10 5 to 8 × 10 6 Pa s. At higher viscosity, the erupted magma tends to remain close to the vent, forming obelisks, rather than to advance horizontally forming lava flows. Figure 4A shows that the dome height increases with the crystal-free magma viscosity while the lateral extent of the dome decreases. Therefore, the crystal-free magma viscosity controls the height and extent of lava domes.
Experiments 2: These experiments explore the influence of CCGT on the morphology of lava domes. The crystal-free magma viscosity is prescribed to be η p 1.6 × 10 5 Pa s. All other parameters are the same as in experiments 1. Figure 4B shows the dome morphology for the four CCGT. The greater CCGT the more the lava dome advances horizontally and less vertically. To better fit the shape of the observed lava dome, we choose 5 days for CCGT in our modeling.
Experiments 3: These numerical experiments explore the influence of the lava extrusion rate on the dome morphology. Figure 4C illustrates the morphological shapes of lava domes at various effusion rates. The higher the rate the larger the dome as it accommodates the larger volume of the extruded lava.

Lava Dome Morphology
Based on the sensitivity analysis, we have constrained the model parameters to get a better fit between the observed and modeled heights of the lava dome at Volcán de Colima. We present here the results of three numerical experiments, which are consistent with the observations on the dome growth from February 1, 2007 to October 20, 2009. The model parameters are listed in Table 1; the model experiments differ from each other by their geometry.
In these model experiments, the lava extrusion rate u ext is obtained empirically from the analysis of the lava dome images at different times and from the results of the sensitivity analysis with respect to varying rates of extrusion (experiments 3) This empirical approach is supported theoretically (see Eq. A4 in Appendix A2; we note that the calculated values of the rate are in close vicinity to the rates obtained empirically). To match the observed and modeled heights of the lava dome, the following extrusion rates u ext have been assigned (see Figure 5): 5 × 10 −6 m s −1 for the first 300 days, 7.5 × 10 −6 m s −1 from day 301 to day 455, and 2 × 10 −5 m s −1 from day 456 until day 1,020, the end time of the numerical experiments.
The crystal-free magma viscosity is assumed to be 1.6 × 10 5 Pa s for 150 days, and it has been increased by a factor of 5 (8 × 10 5 Pa s) after that and kept unchanged until the end of the experiments. This increase of the model viscosity can be attributed to the increased degassing at the initial stage of the lava dome growth evident from seismic records ( Figure 2B). We note that the effects of volatiles and bubbles on the magma rheology are not considered explicitly in the model. This would complicate the model by additional equations, describing the bubble dynamics. Also, volatile exsolution is not equilibrium in high viscosity magmas, and faster extrusion rates lead to higher amounts of dissolved volatiles. This increase of the model viscosity favors upward growth of the dome compared to horizontal spreading of the extruded lava.
Experiment 4: The morphological shapes of the modeled lava dome are presented in Figures 6A,D. The results show that the extruded lava develops a dome structure, which fills the volcano crater ( Figure 6A). Due to the interplay between the extrusion rate and the lava viscosity, the dome grows up and advances horizontally. While dome dynamics is dominated by vertical growth for at least 480 days, a horizontal advancement  becomes dominant later due to gravitational spreading of the erupted lava. The maximum heights of the modeled and observed lava domes show a good agreement (Figures 6D, 7). Meanwhile, the morphological shapes of both domes do not fit each other well on the left side of the model domain (or in the eastern part of the crater). Although the photo of the dome initiation in February 2007 ( Figure 1B) shows that it started to build nearly in the center of the crater, with time the extruded lava seems to have been moving toward the west ( Figure 1D). This westward slope in the dome building could be due to the location of the lava vent, which is likely to be opened closer to the western side of the crater, and/ or due to a westward inclination of the crater topography promoting the lava advancement in the direction. Moreover, it was proposed that the upper part of the conduit at Volcán de Colima may be bending westward close to the surface (Lavallée et al., 2012;Salzer et al., 2014). To test the hypotheses, we have modified the geometry of the crater and performed two other experiments.
Experiment 5: This experiment differs from exp. 4 by the geometry of the computational domain ( Figure 3B). The topography of the crater is assumed to be flat on its right side and inclined on its left side. The location of the conduit is shifted to the right of the model domain (to the west) by the size of the conduit's diameter (14.5 m). The model parameters are the same as in exp. 4. Figures 6B,E illustrates the model dome evolution in exp. 5. Since initiation of the lava dome growth, the dome moves easily toward the west (to the right side of the model domain) because of the slope of the crater's topography on the left side of the model domain. Once the lava reaches the western border of the crater, the dome builds up in the west and extends horizontally to the east ( Figure 6B). At the later stage of the dome growth, the morphological shapes of the domes in experiments 4 and 5 show similarity. Figure 6E shows that the morphological shapes of the observed and modeled lava domes fit well enough at least until May 1, 2008 (day 480), and the maximum heights of both domes are close to each other (see also Figure 7). Although the shapes of both domes on the western part of the crater converge better in exp. 5 compared to exp. 4, there are still discrepancies between the shapes on the eastern side of the crater, namely, the eastern side of the observed dome is steeper than that  of the modeled dome. More complex geometry of the crater on its eastern side as well as an increasing viscosity, as the lava fully solidifies at the edges, could explain the observed morphology of the dome. Experiment 6: In this experiment, the crater base is flat ( Figure 3B), and the conduit is shifted (by 14.5 m) toward the west compared to the geometry of the crater in the case of exp. 4. The model parameters are the same as in exp. 4. Figures 6C,F illustrates the modeled dome evolution in this experiment. The flat topography of the model crater promotes the lava dome to advance in both direction (to the west and to the east), although the lava reaches the western side of the crater faster than its eastern side due to the location of the conduit in the model ( Figure 6C). The morphological shapes of the modeled and observed domes fit each other rather good at least until May 2008 ( Figure 6F). However, later the maximum heights of the modeled lava dome are smaller than those of the observed lava dome. Comparing the maximum heights of the modeled domes, one can see that the dome heights in exps. 4 and 5 agree well with the observed data especially for the time interval from 480 to 704 days. Meanwhile, the heights of the modeled dome in exp. 6 do not match well those of the observed dome except the time closer to the end of 2009. As all parameters of the model in exps. 4-6 are the same except the model geometry, we consider that the best fit of the modeled dome height and width to the observations is achieved in exp. 5, where the crater floor is inclined westward on the left of the vent and the floor is horizontal on its right.

Lava Dome Viscosity
We illustrate here the distribution of the modeled crystal content, viscosity, and velocity for times of November 1, 2007, May 1, 2008, and April 1, 2009 in the case of experiments four to six. The crystal content within the model conduit and the lava dome is presented in Figure 8. A new magma enters the model conduit with the initial volume of the crystal content equal to 0.4, and it reaches almost its equilibrium value 0.83 at the surface near the vent, because the magma crystalizes rapidly at small extrusion rates and small CCGT. We note the magma viscosity remains finite at ϕ eq 1 in terms of viscosity formulation (see Eq. 4) by Costa et al. (2009). Moreover, 100% crystal content is never observed in eruptive products. According to Eq. 4, more rapid growth in viscosity occurs at lower crystal contents of 0.4-0.7 depending on the crystal shape and size distribution. The volume fraction of crystal and melt content in the magma uplift controls its effective (bulk) viscosity (Hess and Dingwell, 1996;Melnik and Sparks, 2005). This volume fraction changes with pressure and flow velocity, and the resulting morphology of a lava dome can be affected (Watts et al., 2002;Melnik and Sparks, 2005).
The viscosity distribution in the lava dome mimics the distribution of the crystal content ( Figure 9). The lava dome viscosity reaches 10 12 Pa s, while the viscosity of magma in the conduit varies from about 10 9.8 to about 10 11.9 Pa s. For lavas from Volcán de Colima, Lavallée et al. (2008) estimated the melt viscosity and the effective viscosity to be 10 8.58 and 10 11.28 Pa s, respectively, where the effective viscosity was calculated using the Einstein-Roscoe equation (Einstein, 1906;Roscoe, 1952). The viscosity obtained in the model is in a good agreement with the estimated effective viscosity of lavas from Volcán de Colima, although direct comparison between calculated and measured viscosities has some difficulties including strain localization, shear heating and other effects. For example, the apparent viscosity of lavas from Volcán de Colima estimated experimentally by Lavallee et al. (2008) is between about 10 9 and 10 9.4 Pa s at stress of 8 MPa. At small extrusion rates and the prolonged duration of the dome growth, the higher viscosity provides the lava dome to grow upward (see exp. 1), while the lower viscosity leads to a horizonal advancement of the lava dome and an insignificant vertical growth. The higher viscosity could be associated with magma degassing during the lava dome building and rapid crystallization as recorded by the increased seismic activity ( Figure 2B) during 2007-2009. The velocity field shows higher flow rates in the conduit and near the vent oriented almost vertically, and the rates decay with the increase of the dome height.

Extrusion Rates
In the following experiments, we employ the extrusion rates calculated from the observed lava dome volumes and calculated discharge (effusion) rates using Eq. 5 in Appendix A2 (see also Figure 2B). In the experiments, the initial extrusion rate is u ext 5.6 × 10 −6 m s −1 ; it increases to 2.25 × 10 −5 m s −1 from day 368 to day 434, decreases to 1.07 × 10 −5 m s −1 from day 435 to day 704 (see Figure 5). Starting from day 705 and until day 1,020, the extrusion rate decreases to 2.71 × 10 −6 m s −1 . The numerical experiments here differ from each other by the geometry of the model domain: in exp. 7 the model geometry is the same as in exp. 4 ( Figure 3A), in exp. 8 as in exp. 5, and in exp. 9 as in exp. 6 ( Figure 3B).
Initially the morphological shape and the maximum height of the lava domes in experiments 7-9 ( Figure 10) follow closely those in experiments 4-6, respectively, as the initial extrusion rate is the same in both sets of the experiments (Table 1; Figure 5). However, the subsequent evolution of the lava dome in experiments 7-9 differs from that in experiments 4-6 starting from day 368, and it is associated with the increase of the extrusion rate by a factor of about four in experiments 7-9. The decrease of the extrusion rate after day 434 results in the gravitational spreading of lava dome instead of its growth upward. The dome height reaches its maximum value at day 704 as after that the extrusion rate decreases by a factor of about 4, and the dome starts to advance horizontally. Unfortunately, the morphology of the modeled domes does not fit that of the observed dome, and the dome heights are low compared to the observations (Figure 10).

Lava Dome Carapace
Lava domes usually have a solid surface layer (carapace) but remain mobile and undergo deformation for days to months (Anderson and Fink, 1990). Iverson (1990) highlighted the importance of the thermal carapace in determining the morphology of domes. The role of the surface cooling and the carapace formation were analysed in laboratory and numerical experiments (e.g., Fink and Griffiths, 1998;Griffiths, 2000;Hale and Wadge, 2003;Hale et al., 2007;Hale, 2008;Husain et al., 2014Husain et al., , 2018Husain et al., , 2019Harnett et al., 2018). A thickness of the thermal boundary layer ( κt √ ), which is associated with a rigid part of the lava dome carapace, ranges from ∼1 to ∼10 m for dome growth time t from 10 days to 3 years, respectively (assuming the coefficient of the thermal diffusivity to be κ 10 −6 m 2 s −1 ). Tsepelev et al. (2020) showed that at high discharge rates and for a short time duration of lava extrusion, this lava generates a thin carapace due to cooling, which influence only slightly the morphology of a dome promoting a steeper slope similar to that observed during the dome growth at Volcán de Colima.
To simulate a thermal carapace development during lava dome growth, our mathematical/numerical model should be supplemented by the heat equation, which would complicate the model and computations. Therefore, to keep the same numerical model, we modify the lava dome viscosity by introducing a higher viscous layer at the interface between the dome and its surrounding in order to simulate a dome carapace (see Appendix A1).
In experiment 10, we assume that the crater geometry and the model parameters are the same as in exp. 5, except for the extrusion rate. This numerical experiment starts at day 481 of the model dome growth as in exp. 5. A higher viscosity carapace of ∼6 m thick is introduced in the model at day 481 simultaneously with a decrease of the extrusion rate by 15% compared to the rate in exp. 5 (see Figure 5). This decrease in the extrusion rate has been introduced to reduce the volume of erupted lava and to better fit the observed morphological shapes. Figure 11 presents several stages of the dome growth in exp. 10. The dome morphology differs significantly compared to that in exp. 5 ( Figure 11A). Particularly, 1) the dome develops a steep flank on its left side, which fits observations at Volcán de Colima rather well, and 2) the advancement of lava dome at its left side is retarded by about 50 m for 334 days, compared to dome advancement in exp. 5. Meanwhile, the dome heights in both experiments ( Figure 11B) are in a good agreement with the observations. The modeled viscosity and velocity in the lava dome are presented in Figure 12 together with the images of the dome at Volcán de Colima for three stages of its growth. It is seen that exp. 10 provides much better fit between the model and observations.

Model Outcomes
We have presented a plausible numerical model for lava dome growth at Volcán de Colima during the long dome-building episode lasting from early 2007 to fall 2009 without explosive Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 735914 dome destruction. This numerical study has allowed for estimating the influence of the model geometry (the crater topography and the conduit's location), extrusion rates, and a high-viscous carapace on the morphology of the lava dome. Three geometries of the crater have been explored with the different locations of the vent, and it has been shown that the crater topography is likely to be inclined toward the west as shown in experiments 5 and 10. The numerical results show that the extrusion rates calculated empirically from 2-D images of the morphological shapes of the lava dome at different times lead to a better fit between the observed and modeled shapes. The best fit provides the reduced extrusion rates as shown in exp. 10; the extrusion rates calculated using the erupted lava dome volumes yield to underestimation of the maximum height of the observed lava dome. As a lava advancement within the dome depends on the surface topography of the crater as well as the conduit geometry, two-dimensional (2-D) models cannot capture complexity of the three-dimensional (3-D) lava dome growth.
The developed model has shown that high viscosity of the dome is required to maintain a vertical growth. As the viscosity depends on the volume fraction of crystals in our model, and the latter depends on CCGT, the crystal content grows rapidly in the uppermost conduit at low discharge rates ( Figure 2B), and lava becomes highly viscous. The modeled lava viscosity within the dome core varies from 10 9.8 to 10 12 Pa s and is consistent with the viscosity measurements of lavas from Volcán de Colima.
Using recorded dataset from Volcán de Colima we have shown that solidification and rheological stiffening within a lava dome are controlled in part by cooling and degassing-induced crystallization. While the crystallization due to degassing is a dominant process in dome-building eruptions (e.g., Melnik and Sparks, 1999), especially during short dome-building episodes (e.g., Tsepelev et al., 2020), the model forecasts indicate that a high-viscous carapace significantly influences the dome growth process during long lava dome-building episodes. A carapace of higher viscosity (about 10 14 Pa s) prevents a rapid dome advancement to the east and promotes the development of a steep slope on the dome's eastern side.
Cooling influences the lava viscosity at the interface between the dome and its surrounding making the dome carapace more viscous and promoting a development of lobe-shaped lava dome (Watts et al., 2002;Tsepelev et al., 2020). In the cases of short episodes of lava dome growth, the thermal carapace of the lava dome is not thick enough to present significant resistance to lava dome horizontal advancement (e.g., Bourgouin et al., 2007;Tsepelev et al., 2020). The dome carapace becomes thicker in the case of radiative and convective heat transfer at the lava/air interface (Tsepelev et al., 2019) or in the cases of long domebuilding episodes, as shown in this work. Further studies related to the influence of temperature on the viscosity of the carapace during long lava dome-building episodes can refine the current results. Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 735914

Model Limitations
The numerical model is unlikely to represent exact dynamics of lava dome growth due to the complexity of magmatic processes at depths and unknowns such as degassing, initial crystal volume content, and viscosity. Any model comes with its simplifications (e.g., less complex viscosity, and 2-D or axisymmetric instead of 3-D geometry) and assumptions (e.g., temperature-independent density and viscosity, and no density change due to crystallization). For example, difficulties related to 3-D datadriven numerical modeling (e.g., data assimilation, inversions, and computational cost) are also associated with data issues. Several techniques for volcano monitoring are available to collect and use data on lava dome growth, for example, video recording (Bretón-Gonzalez et al., 2013), infrared camera imaging (Stevenson and Varley, 2008), terrestrial photogrammetry (James and Varley, 2012), and satellite radar observations (e.g., Salzer et al., 2014;Walter et al., 2019). Meanwhile, it is still challenging to determine precise 3-D morphological changes during dome growth, especially when the data come from the records using a single video (photo) camera as 2-D images have a limited information related to a 3-D object. Other challenging issues are related to the determination of the detailed crater geometry prior a new dome starts to evolve, and the location and geometry of the uppermost conduit (e.g., Pinel et al., 2011). The problem of a search for the best fit between the morphological shapes of observed and modeled lava domes tuning the parameters of the lava viscosity belongs to the class of inverse problems. The solution to the inverse problem can be non-unique; for example, different discharge rates and different CCGT can produce similar lava dome shapes (Starodubtseva et al., 2021). As for the discharge rate, it can be determined from observations, and, therefore, this rate can be considered as a known characteristic of the lava dome dynamics constraining the model. In our numerical modeling the optimization of the misfit between the modeled and observed morphologies has been based on tuning CCGT, extrusion rates, and the crystal-free magma viscosity at the times of available observations.
Matching the timing of lava dome growth at Volcán de Colima to its morphology is challenging, because the estimation of the volume of the erupted lava is restricted to analysis of camera images and hence uncertain; the existing 3-D temporal coverage of the extrusion process is not complete; the location of the vent is not precisely documented; and the determination of the extrusion rate in 2-D models is complicated. Therefore, with more observations 3-D numerical models of lava domes should provide additional information on lava dome growth and their eventual collapse.

Model Applications
The developed model in this work and similar models can be used to analyze future effusive eruptions and lava dome growth at Volcán de Colima or elsewhere after proper calibration based on history matching of dome growth by nudging model forecasts to observations (i.e., minimizing misfits between the modeled and observed morphological shapes of domes). The model can be used to assess a stress localization in the dome carapace and its potential failure, which may lead to pyroclastic flow hazards.
In the future, every potentially hazardous volcanic eruption should be accompanied by its virtual numerical model that is constantly tuned by a history matching procedure and gives short-and long-term forecasts of the eruption dynamics and associated hazards. This will require increasing the accuracy of monitoring techniques and considerable investment in geophysical studies of volcanic systems.

DATA AVAILABILITY STATEMENT
The data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AI-Z, OM, and VZ contributed to the conceptualization of the paper. VZ provided data on the dome growth and also an insight into the development of the dome at Volcán de Colima during 2007-2009. NZ processed the data, developed a numerical model, and performed all model experiments. AIZ supervised the study, and assisted NZ in numerical modeling and data processing. IT contributed to the numerical statement of the problem, data processing, and the development of numerical codes. All authors contributed to the interpretation of numerical results. AIZ and NZ performed writing-original draft preparation and figures. All co-authors performed writing-review and editing.

ACKNOWLEDGMENTS
This is a contribution to the research project funded by the Russian Science Foundation's grant . NZ thanks also Deutsche Forschungsgemeinschaft (project DFG IS203/14-1) for support. We are grateful to Frank Schilling for providing helpful insight into physical properties of lavas. Also, we are thankful to Valerio Acocella, Shanaka de Silva, and two reviewers for careful reviews of the initial manuscript and constructive comments.

APPENDIX A1 Numerical Approach
The numerical solver uses an implicit scheme of integration for the governing equations (Eqs. 1-5) with the relevant boundary and initial conditions. The model domain is discretized by about 54,000 hexahedral cells. The pressure is discretized by the secondorder staggered scheme PRESTO! (e.g., Peyret, 1996). To approximate the Laplacian, we use the numerical scheme of the second-order accuracy, and the monotonic schemes are used to discretize convective terms in the equations (e.g., Ismail-Zadeh and Tackley, 2010). The pressure-velocity coupling is handled by the SIMPLE method (Patankar and Spalding, 1972), where the relaxation parameters are chosen to be 0.015 and 0.3 for the velocity and pressure, respectively. We assign 0.5 to the relaxation parameter for the function α, and the volume fraction of crystals. Considering a discontinuity between the lava and air viscosities, the choice of the relaxation parameters is critical and sometimes it was lowered to ensure the solution's stability. A time step is chosen in the range of 0.1-40 s depending on the stability and optimization of the velocity to assure a convergence of a set of linear algebraic equations (SLAE), which is obtained after the discretization of the Navier-Stokes equations. The implicit scheme results in stable computations with a relatively large timestep. We note that an employment of explicit schemes in the model led to unstable numerical results even for the small Courant number. In the modelling, we employ the conjugate-gradient method to solve SLAE. The numerical accuracy attains 10 −3 for the velocity and pressure and 10 −6 for the function α and the volume fraction of crystals. Numerical experiments were carried out on the multi-processor highperformance computer (bwHPC) of the Karlsruhe Institute of Technology. In the modelling, 1 day of the lava dome evolution takes about 4 min of the computational time depending on the time step size.
The problem of magma extrusion and lava dome growth involves the interface between the lava and the air. Numerical schemes usually produce errors, originating from discontinuities of physical parameters at the interface (e.g., the density and viscosity), and related smearing of the parameters along the interface. To minimize such errors and smearing, special numerical methods have been developed: the methods are rather simple for the case of discontinuous density (e.g., Christensen, 1992;Naimark and Ismail-Zadeh, 1995), but more difficult for the case of discontinuous viscosity (e.g., Naimark et al., 1998). In this modeling the viscosity jump at the lava-air interface is significant (many orders of the magnitude), and this generates indeed smearing at the interface. To reduce the smearing, a finer mesh can be introduced around the interface.
At the initial time t 0 the conduit (sub-domain Ω 1 ) is filled by the old magma. To determine the initial velocity in the entire model domain, we consider u(t 0, x) 0 and small time step (10 −1 -10 −3 s), and perform a numerical experiment for 100 to 1,000 time steps. The velocity obtained at the end of this experiment we take as the initial velocity u 0 to run further computations with larger time steps. We note that the Reynolds number for the lava dome dynamics is small for the slow flow and the high lava viscosity, and hence, the influence of the inertial terms in the Navier-Stokes equations is insignificant. Therefore, the initial condition for velocity influence insignificantly the dynamics of lava flow and dome growth evolution. On the contrary, the Reynolds number is large for the air dynamics, and the inertial terms play a significant role. The described approach yields to a laminar flow in the air subdomain ensuring the stability of further calculations. The influence of the air phase on lava dome growth is negligible because of a large density and viscosity contrast between the magma and the air.
To simulate numerically the development of a highly viscous dome carapace due to cooling, a 'thermal boundary layer' is introduced at the interface between the lava dome and its surrounding. This allows to avoid introducing and solving the heat equation, and hence do not complicate the existing model. In experiment 10, the viscosity of the boundary layer is assumed to be by a factor of 200 higher than the viscosity of the dome core. The thickness of the layer is about 5-6 computational cells adjacent perpendicularly to the interface, and each cell is of about 1 m thick. Figure A1 presents the model carapace.

A2 Determination of the Model Extrusion Rates Using the Extruded Lava Volume Data
Since January 2007 a slow endogenous growth of the lava dome continued until the end of 2009. Figure 2 presents the values of the lava dome volume and the discharge rates U dis (m 3 s −1 ) for the time intervals, where the data on the dome lava volumes have been available (Figure 2A). When solving numerically 2-D models of the lava extrusion, difficulties arise with the FIGURE A1 | A layer (red) of the higher viscosity introduced in numerical experiment 10 to simulate a dome carapace.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 735914 calculation of the model extrusion rate, which will correspond to the observed effusion rate. Below we derive the formulas for conversion of the observed discharge rate U dis into the model extrusion rate u ext (measured in m s −1 ), which have been used in our 2-D model. This approach is applicable in the case of endogenous lava dome growth. We approximate a lava dome shape by a half-spheroid with the equal a-and b-semi-axis (a b) and height c ( Figure A2). In Figure A2, r is the radius of the conduit, z* is the height of the lava dome, r + x* is the radius of the lava dome, d x and d z are the increments of the radius and height of the lava dome for the time t 2t 1 , respectively. We calculate the difference between the volume V 2 of the outer half spheroid (the volume of lava at time t 2 ) and the volume V 1 of the inner half-spheroid (the volume of lava at time t 1 ). This difference will characterize the volume of the lava extruded for time interval t 2t 1 .
The volume of the extruded lava for time interval t 2 -t 1 can be found using the discharge rate as U dis (t 2 − t 1 ). Hence Equation (A1) can be re-written as: Considering the 2-D case study, we calculate the difference between the area of the cross-section through the center of the outer half spheroid S 2 at time t 2 and the area of the cross-section through the center of the inner half-spheroid S 1 at time t 1 . This difference will characterize the area of the lava, which was extruded during time interval t 2t 1 : The lava area extruded for time interval t 2 -t 1 equals to 2u ext r(t 2 − t 1 ), and Equation (A3) can be re-written as: Inserting (t 2t 1 ) from Eq. (A2) into Eq. (A4) we can get the conversion formula from the observed discharge rate to the model extrusion rate: The parameters x*, z*, d x , and d z are accessed from the images of the 2-D morphological shapes of the growing lava dome (see Figure 1E; Bretón-González et al., 2013). The extrusion rate u ext is then calculated using Equation (A5), where the discharge rate U dis is presented in Figure 2B. This extrusion rate is used in numerical experiments 7-9. Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 735914