Analysing stress field conditions of the Colima Volcanic Complex (Mexico) by integrating finite-element modelling (FEM) simulations and geological data

In recent decades, finite-element modelling (FEM) has become a very popular tool in volcanological studies and has even been used to describe complex system geometries by accounting for multiple reservoirs, topography, and heterogeneous distribution of host rock mechanical properties. In spite of this, the influence of geological information on numerical simulations is still poorly considered. In this work, 2D FEM of the Colima Volcanic Complex (Mexico) is provided by using the Linear Static Analysis (LISA) software in order to investigate the stress field conditions with increasingly detailed geological data. By integrating the published geophysical, volcanological, and petrological data, we modelled the stress field considering either one or two magma chambers connected to the surface via dykes or isolated (not connected) in the elastic host rocks (considered homogeneous and non-homogeneous). We also introduced tectonic disturbance, considering the effects of direct faults bordering the Colima Rift and imposing an extensional far-field stress of 5 MPa. We ran the model using the gravity in calculations. Our results suggest that an appropriate set of geological data is of pivotal importance for obtaining reliable numerical outputs, which can be considered a proxy for natural systems. Beside and beyond the importance of geological data in FEM simulations, the model runs using the complex feeding system geometry and tectonics show how the present-day Colima volcanic system can be considered in equilibrium from a stress state point of view, in agreement with the long-lasting open conduit dynamics that have lasted since 1913.

Abstract. In recent decades, finite-element modelling (FEM) has become a very popular tool in volcanological studies and has even been used to describe complex system geometries by accounting for multiple reservoirs, topography, and heterogeneous distribution of host rock mechanical properties. In spite of this, the influence of geological information on numerical simulations is still poorly considered. In this work, 2D FEM of the Colima Volcanic Complex (Mexico) is provided by using the Linear Static Analysis (LISA) software in order to investigate the stress field conditions with increasingly detailed geological data. By integrating the published geophysical, volcanological, and petrological data, we modelled the stress field considering either one or two magma chambers connected to the surface via dykes or isolated (not connected) in the elastic host rocks (considered homogeneous and non-homogeneous). We also introduced tectonic disturbance, considering the effects of direct faults bordering the Colima Rift and imposing an extensional far-field stress of 5 MPa. We ran the model using the gravity in calculations. Our results suggest that an appropriate set of geological data is of pivotal importance for obtaining reliable numerical outputs, which can be considered a proxy for natural systems. Beside and beyond the importance of geological data in FEM simulations, the model runs using the complex feeding sys-

Introduction
Magmatism and tectonism are strongly related to regional and local stress fields, affecting both the orientation of faults and the location of volcanic vents (Geyer et al., 2016). The stress field around a magmatic source originates from three main contributions: (1) the background stress, composed of a vertical gravitational load, a lateral horizontal load (lithostatic confinement), and tectonic regime; (2) the stress field caused by the loading of the volcano edifice; and (3) the stress field generated by the magmatic pressure (e.g. Marti and Geyer, 2009;Currenti and Williams, 2014). In recent years, a large number of semi-analytical and numerical methods have been proposed for the solution of the stress field of natural systems (e.g. Cayol and Cornet, 1998;Simms and Graven, 2004;Manconi et al., 2007;Long and Grosfils, 2009;Currenti et al., 2010;Currenti and Williams, 2014;Zehner et al., 2015), taking into account the static-elastic deforma-Published by Copernicus Publications on behalf of the European Geosciences Union. 2516 S. Massaro et al.: Analysing the stress field at the CVC (Mexico) tion in a multilayered half-space (e.g. Dieterich and Decker, 1975;Bonafede et al., 2002;Wang et al., 2003;Gudmundsson and Brenner, 2004;Zhao et al., 2004;Pritchard and Simons, 2004;Gottsmann et al., 2006;Geyer and Gottsmann, 2010;Zhong et al., 2019). Following successful application in mechanical engineering, the use of finite-element modelling (FEM) has been extensively introduced in volcanology in order to investigate the effects of topography, lithologic heterogeneities, tectonic stresses, and the gravity field on the stress of volcanic systems (e.g. Cailleau et al., 2003Cailleau et al., , 2005Buchmann and Conolly, 2007;Manconi et al., 2009;Masterlark et al., 2012) including volcanoes (e.g. Fujita et al., 2013;Charco and Galán del Sastre, 2014;Coco et al., 2016;Ronchin et al., 2015;Hickey et al., 2015;Cabaniss et al., 2018;Rivalta et al., 2019). There are several examples of the use of FEM for volcanic areas, which span from the influence of layered materials on the surface deformation process during volcanic inflation (e.g. Darwin volcano, Galapagos Islands; Manconi et al., 2007;Albino et al., 2010) to processes affecting chamber rupture (e.g. Grosfils, 2007;Long and Grosfils, 2009). The local stress around a volcanic feeding system depends on the geometry of the magma plumbing system, including the chamber(s) and dykes forming it, and on the mechanical properties of the host rock around it (e.g. Marti and Geyer, 2009), especially on changes in the Young modulus (e.g. Gudmundsson, 2011;Jeanne et al., 2017;Heap et al., 2020). For instance, limestone, lava flows, welded pyroclastic units, and intrusive rocks can be very stiff (high Young modulus from ca. 1.7 to 27 GPa for limestones, Touloukian et al., 1989;ca. 5.4 GPa for volcanic rocks, Heap et al., 2020), but young and non-welded pyroclastic units may be very soft (low Young modulus ca. 1.7-3.1 GPa, Margottini, 2013). Therefore, the local stress may abruptly change from one layer to another (e.g. Gudmundsson, 2006). Irrespective of the scope of the numerical investigation, the importance of applying accurate physical constraints to FEM has already been discussed in many studies (e.g. Folch et al., 2000;Newman et al., 2001;Fernández et al., 2001;Currenti et al., 2010;Geshi et al., 2012). However, in the last decade, few investigations have been carried out to assess the influence of the amount and quality of geological data on FEM computations (Kinvig et al., 2009;Norini et al., 2010Norini et al., , 2019Cianetti et al., 2012;Ronchin et al., 2013;Chaput et al., 2014). To bridge this gap, in this work we use the Linear Static Analysis (LISA) software (version 8.0; https://www.lisafea.com/, last access: 10 December 2020) to study the subsurface field state in the Colima Volcanic Complex (CVC, Mexico) at increasing geological detail. The CVC area is a good candidate for testing the response of FEM software against different geological conditions, being constituted by a large volcanic complex (Lungarini et al., 2005) within a tectonic graben filled with volcano-clastic material ( Fig. 1a; Norini et al., 2010Norini et al., , 2019. The FEM was run starting from simple homogeneous vs. stratified lithology of the subsurface in successively more detail through the ad-dition of single and double magma chambers, feeder dykes, faults, and extensional far-field tectonic stress (Fig. 1b).

Geological framework
The Pleistocene-Holocene CVC is one of the most prominent volcanic edifices within the Trans-Mexican Volcanic Belt (TMVB) (Macías et al., 2006;Capra et al., 2016;Norini et al., 2019, Fig. 1a). In this area, the Rivera microplate and the Cocos plate subduct beneath the North America Plate along the Middle American Trench, forming a triple junction that delimits the tectonic units known as the Jalisco Block (JB) and the Michoacán Block (MB) (Luhr and Carmichael, 1985;Allan, 1986;Rosas-Elguera et al., 1996Ferrari et al., 1999;Rosas-Elguera et al., 2003;Frey et al., 2007). The three rifts of this system are the Tepic-Zacoalco (TZR), the Chapala-Tula (CTR), and the Colima Rift (CR). The still active N-S-trending Colima Rift (CR) was formed during an extensional phase that occurred after the Late Cretaceous-Paleogene compressive and transpressive phase (Allan, 1986;Bandy et al., 1995;Cortés et al., 2010). While opening, the CR was gradually filled with Pliocene-Quaternary lacustrine sediments, alluvium, and colluvium (e.g. Allan, 1986Allan, , 1991Norini et al., 2010). The geometry, kinematics, and dynamics of the CR have been studied on the basis of field, seismic, and geodetic data, mainly collected in its northern and central sectors (see Fig. 1 in Norini et al., 2010). The magnitude of vertical displacement of the northern and central sectors is ca. 2.5 km by adding the topographic relief of the bounding fault scarps (1.5-1.6 km) to the calculated sediment depth (Allan, 1985;Serpa et al., 1992). Field data and focal mechanism solutions are consistent with a direction of opening of the northern and central sectors oriented from E-W to NW-SE, with mainly normal and minor rightlateral displacements of the bounding faults (Barrier et al., 1990;Suárez et al., 1994;Rosas-Elguera et al., 1996;Norini et al., 2010Norini et al., , 2019Garduño-Monroy et al., 1998). In contrast to field and seismic evidence of long-term slightly dextral oblique extension, recent GPS geodetic measurements suggest a possible sinistral oblique extension of the CR (Selvans et al., 2011). In both cases, the stress regime is mainly extensional, with an approximately E-W orientation of the minimum horizontal stress in the basement of the CVC (Barrier et al., 1990;Suárez et al., 1994;Rosas-Elguera et al., 1996;Selvans et al., 2011;Norini et al., 2010Norini et al., , 2019. The CVC stands within the central sector of the CR on top of the Cretaceous limestones, late Miocene-Pleistocene volcanic rocks, and Pliocene-Holocene lacustrine sediments, alluvium, and colluvium (Allan, 1985(Allan, , 1986(Allan, , 1991Cortés et al., 2010;Norini et al., 2010;Escudero and Bandy, 2017) considering zero displacement along the bottom and left and right sides. Note that for case (vi) in (b) the zero displacement is removed from the lateral sides; (d) sketch of the Fuego de Colima feeding system composed of a 15 km deep magma chamber connected to the surface via a 6 km deep magma chamber and dykes. P chs and P chd are the magmatic overpressures in the shallow and deep chambers, respectively (modified from Massaro et al., 2019). the youngest and active Volcán de Colima (3763 m a.s.l.) (Norini et al., 2019).

Eruptive activity
The eruptive history of the CVC started in the northeast area with the formation of the Cantaro volcano at ca. 1-1.5 Ma, followed by Nevado de Colima starting at ca. 0.53 Ma, which is composed of voluminous andesitic lava domes and deposits associated with caldera-forming eruptions and partial sector collapses (Robin et al., 1987;Roverato et al., 2011;Roverato and Capra, 2013;Cortés et al., 2019). The youngest -Volcán de Colima -comprises the Paleofuego edifice, which suffered several sector collapses that formed a horseshoe-shaped depression where the new active cone (also known as Volcán de Fuego) grew up. Its activity was 2518 S. Massaro et al.: Analysing the stress field at the CVC (Mexico) characterized by dome growths and collapses, extrusion of lava flows, and Vulcanian and occasionally sub-Plinian explosive eruptions Massaro et al., 2018Massaro et al., , 2019.

The CVC plumbing system
Seismic tomography (Spica et al., 2017) highlights a 15 km deep low-velocity body (LVB) which was interpreted as a deep magmatic reservoir. It is confined within the CR, suggesting a structural control of the normal fault system on it (Spica et al., 2014). The LVB has an extent of ca. 55 km × 30 km in the N-S and E-W directions, respectively, showing an averaged thickness < 8 km. Escudero and Bandy (2017) obtained a higher-resolution tomographic image of the CVC subsurface area, showing that the most active magma generation zone is now under the Fuego de Colima edifice. The ambient seismic noise tomographic study of Spica et al. (2014) indicates a shallow magma chamber above ca. 7 km of depth, in agreement with petrological studies (Medina-Martínez et al., 1996;Luhr, 2002;Zobin et al., 2002;López-Loera et al., 2011;Reubi et al., 2013Reubi et al., , 2019Macías et al., 2017). Cabrera-Gutiérrez and Espíndola (2010) suggested that the shallow active magma storage has a volume of ca. 30 km 3 . It is connected to the surface by conduits, whose path is facilitated by the presence of the CR fault zone, which provides a natural pathway for fluids (e.g. Allan, 1986;Norini et al., 2010Norini et al., , 2019. The arrangement of dykes and the alignment of volcanic centres of the CVC suggest that the dykes swarm draining the magma chambers developed along the NNE-SSW-trending, steep, eastwarddipping normal fault exposed on the northern CVC flank (Norini et al., 2010(Norini et al., , 2019. Massaro et al. (2018) provided a first-order geometrical reconstruction of the Fuego de Colima feeding system during the 1913 sub-Plinian eruption by using volcanological data Bonasia et al., 2011) as inputs and constraints for numerical simulations. Results showed good matches for a hybrid configuration of the shallow conduit feeding system composed of a ca. 5500 m long, 200-2000 m wide, and 40 m thick dyke passing into a shallower (500 m long, 40 m diameter) cylindrical conduit. The shallow magma chamber top was set at 6 km of depth with a dyke-cylinder transition at 500 m below the summit as inferred from geophysical data (Salzer et al., 2014;Arámbula-Mendoza et al., 2018).

Methods
In this study, we used the commercial 8.0 version of LISA (https://www.lisafea.com/, last access: 10 December 2020), a general-purpose finite-element analysis (FEA) software developed in the 1990s based on the formulations proposed by Rao (1989) and successively integrated from other sources (Michaeli, 1991;Schwarz, 1991;Babuška et al., 1995;Bathe , 1999). Despite LISA originally being used for structural analysis (Rao, 1989(Rao, , 2013, it successfully predicts the stress-strain behaviour of rock masses in elastic models, in particular the deformation mechanisms, even in layered rock masses (Gabrieli et al., 2015).

Modelling approach
The stress field of the CVC plumbing system is simulated considering an E-W cross section, parallel to the extension associated with the active CR (Norini et al., 2010), as shown in Fig. 1a-b (a-a ). Since the extent of the CVC magma chambers in the NNE-SSW direction is typically much longer than the dimensions of the E-W cross section (Spica et al., 2017), 2D solutions of either numerical or analytical models describing E-W elongated magma chambers in the crust can be reasonably adopted (Jaeger et al., 2009;Costa et al., 2011). A topographic profile and 2D plane along the chosen E-W cross section of the CVC area were obtained in ESRI ArcGIS from a digital elevation model (DEM; resolution 50 m; Instituto Nacional de Estadística y Geografía -INEGI, https://en.www.inegi.org.mx/, last access: 10 December 2020) and imported into Autodesk Auto-Cad R13 using a third-degree spline approximation. The IGES file was then imported into LISA for the mesh discretization. The investigated domain extends 60 × 30 km in an x-z Cartesian coordinate system with the three-and four-node finiteelement discretization (Table 1). Zero normal displacements are assigned at the bottom and the lateral boundaries, while the upper boundary represents the free-stress ground surface (Fig. 1c). The FEM is carried out by using a plane strain approximation, implying that the deformation in the third direction is assumed to be negligible. As reported by Zehner et al. (2015), FEM of geological structures requires accurate discretization of the computational domain. It follows that the unstructured tetrahedral meshes have to fulfil the following requirements: (i) sufficient mesh quality, with the tetrahedrons not too acute-angled, since numerical instabilities can occur; (ii) incorporation of geometry for defining boundary conditions and constraints; and (iii) local adaption, which is a refinement of the mesh in the vicinity of physical sources in order to avoid numerical errors during the simulation. In this work, we adopted a mesh composed of 4660 plane continuum elements, which have been refined in the regions of higher gradients (i.e. near the contours of the magmatic feeding system). In our simulations, the extent of the rock layers ( Table 2) refers to that used by Norini et al. (2010Norini et al. ( , 2019. The configuration of the CVC feeding system (i.e. depth, shape, and dimensions of the magma chambers and feeder dykes) derives from the literature (Spica et al., 2014(Spica et al., , 2017Massaro et al., 2018Massaro et al., , 2019, and it is simplified in Fig. 1d. In particular, magma chambers and dykes are considered to be pressurized finite-size bodies in an elastic crustal segment acting as fluid-filled holes. The boundary condition (pressurization) is provided by applying internal forces that act on the walls. This approach has been used extensively in several analytical and numerical models that treat magma reservoirs as internally pressurized ellipsoidal cavities within an elastic half-space in order to gain insight into the behaviour of magma plumbing systems (Pinel and Jaupart, 2004;Gudmundsson, 2006;Grosfils, 2007;Andrew and Gudmundsson, 2008;Hautmann et al., 2009;Currenti and Williams, 2014). Previously published studies indicate that differences between, and problems with, elastic models derive principally from the key role played by gravity (e.g. Lister and Kerr, 1991;Watanabe et al., 2002;Gerbault et al., 2012;Albino et al., 2018). Some authors have argued regarding whether it is appropriate or not to account for the gravity body force in models of volcanic systems (e.g. Currenti and Williams, 2014;Grosfils et al., 2015). When the gravitational loading is not included in the model, the volcanic deformation results from a change with respect to a stage previously at equilibrium (e.g. Gerbault et al., 2018). In this work, we carried out simulations considering the effect of the gravitational loading in the host rock, implemented via body forces. The model initial condition has a pre-assigned lithostatic stress, whose computation, in the presence of topography and material heterogeneities, is not trivial because it requires applying the gravity load preserving the original non-deformed geometry of the mesh (Cianetti et al., 2012). Due to the presence of a lithostatic stress field, the load applied at the reservoir boundaries represents a superposition of the magmatic overpressure and lithostatic component. We define the magmatic pressure as either excess pressure ( P e , magmatic minus lithostatic pressure but below the tensile strength of wall rocks) or overpressure (or driving pressure, P o ), which is the magmatic pressure exceeding the tensile strength of wall rocks (Gudmundsson, 2012). The first pertains to FEM using isolated magma chambers (single or double), while the second is used for models with connected magma chambers (with a conduit feeding system).
We also took into account the effect of the existing faults of the CR system even though LISA cannot include a frictional law to represent the fault movement (Chaput et al., 2014). As reported in Jeanne et al. (2017), the damage induced by faults increases from the host rocks to the fault core, implying a reduction in the effective elastic moduli. In this light, we represented the faults bordering the CR as two damage zones (ca. 70 • of inclination, ca. 1 km thick, and down to 10 km of depth) showing reduced elastic properties with respect to the surrounding rocks. To take into account the effect of the far-field extensional regime, we applied a uniform stress of 5 MPa to the lateral boundaries of the domain, as reported by Marti and Geyer (2009). Considering the E-W cross section (a-a ; Fig. 1a), we provide six domain configurations: (i) the homogeneous lithology model in which the volcanic domain is only composed of andesite rocks; (ii) the non-homogeneous lithology model wherein different geological units are considered; (iii) the single magma chamber model composed of a non-homogeneous lithology and a 15 km deep magma chamber; (iv) the dual magma chamber model composed of a non-homogeneous and 6 and 15 km deep magma chambers; (v) a conduit feeding system model composed of non-homogeneous lithology, 6 and 15 km deep magma chambers connected through a deep dyke, and a shallow conduit connecting to the surface; (vi) an extensional model with a 5 MPa horizontal extensional stress (far field); and (vii) a faulted model, in which two damaged zones mimicking the CR faults (local stress) are also added. The number of nodes is set at 4426 for the only substratum and single magma chamber models, at 4161 for the dual magma chamber model, and at 3737 for the conduit feeding system and faulted models. It is important to note that simulation outputs are shown using different colour scales. Although such a choice may make the visual comparison of the different runs difficult, it preserves the necessary details of stress distribution, which would have been lost using a common colour scale. Finally, in the following we refer to σ 1 as the greatest compressive stress and σ 3 as the least compressive stress.
4 Geological data

Stratigraphy and rock mechanics
Four units forming the CVC system were defined from the available geological data (  Cortés et al., 2010;Norini et al., 2010Norini et al., , 2019. We assumed constant mechanical characteristics within each unit using the typical rock mass properties of density (ρ), the Young modulus (E), and the Poisson ratio (ν) ( Table 2). The rock masses are considered dry for (eventual) pore pressure to be neglected. Only for Unit GF was a higher value for the Poisson ratio used close to the surface in order to mimic high water content in the graben sediments. The maximum thickness of the graben fill (about 1 km) is assumed from the literature (Allan, 1985;Serpa et al., 1992; Hoek and Brown (1997) and Marinos and Hoek (2000), while for volcanic materials (units FC and VD; Table 2) they are estimated according to the approach proposed by Del Potro and Hurlimann (2008). In order to describe the effects of the CG faults on stress field distribution, the mechanical properties are locally degraded in proximity to the faults themselves.

Geometry of the plumbing system
In our 2D model, we assumed the CVC is composed of two magma chambers connected by dykes and to the surface by a conduit (Fig. 1d). The shape of the magma chambers and dykes is represented by elliptical cross sections with the major (2a) and minor (2b) axes. Generally, magma chambers have a sill-like shape that is often imaged in seismic studies of volcanoes and rift zones (Macdonald, 1982;Sinton and Detrick, 1992;Mutter et al., 1995;MacLeod and Yaouancq, 2000;Canales et al., 2009). Most of them are not totally molten but rather a mixture of melt and crystal mush (i.e. Parfitt and Wilson, 2008 Gudmundsson, 2011). During ascent to the surface, the dykes align themselves with the most energyefficient orientation, which is roughly perpendicular to the least compressive stress (e.g. Gonnermann and Taisne, 2015;Rivalta et al., 2019), providing the magma-driving pressure remains small compared to the deviatoric stress (Pinel et al., 2017;Maccaferri et al., 2014a). This behaviour, however, can be modulated in the presence of significant variations in fracture toughness of the surrounding rock due to stratification (Maccaferri et al., 2010) or to old and inactive fracture systems (Norini et al., 2019). Although for oblate magma chambers the propagation of dykes is most probable from the tip areas, in our simulations the orientation of dykes is assumed to be vertical because of the preferential pathways represented by the CR fault planes (Spica et al., 2014). We set the dimensions of feeder dykes in agreement with Massaro et al. (2018): deep dyke 2ad = 2 km; shallow dyke 2a varies from 1 km at the bottom to 500 m in the upper part of the volcano; width of both deep and shallow dyke 2bd = 2b = 100 m (Fig. 1d). It is worth noting that it is outside the scope of this work to provide the conditions for rupture of the magma chamber, since LISA accounts only for the elastic regime. For these reasons, we fixed P e and P o (for isolated and connected magma chambers, respectively) in the range of 10-20 MPa for the 15 km deep chamber and 5 MPa for the 6 km deep one. For the dykes and conduit, P o is set to 10 MPa in the deeper dyke and 5 MPa in the shallower dyke, while in the upper 500 m of conduit it is 0.4 MPa.

Results
In this section we report the sensitivity analysis carried out to quantify the approximation of the Young modulus variation in FEM outputs and a description of the model outputs when adding complexity to the input geological-geophysical data.  Table 2.

Sensitivity analysis of Young modulus
Using the single magma chamber model as a reference case, we quantified the influence of Young modulus variation in each geological unit. Taking into account the mechanical properties of rocks (Table 2) as reference values, we compared the stress state of the computational domain on changing the Young modulus by (±) 1 order of magnitude. This sensitivity analysis, although incomplete, may lead to increased awareness of the selection of input data when running an FEM. The sensitivity analysis was carried out on a reduced simulation domain (the x axis was set to 35 km) in order to diminish the influence of binding effects along the domain borders. We applied the Euclidean norm (L 2 ) method to illustrate the results. The L 2 norm applied on a vector space x (having components i = 1, . . .n) is strongly related to the Euclidean distance from its origin and is equal to In our case, the vector space x is composed of all nodes of the computational domain (Table 1). We defined x ref as the vector containing the results for the maximum and mini-mum principal stress when using the selected values of material properties (Table 1), and x (−) and x (+) are the vectors upon varying the Young modulus by 1 order of magnitude in each unit. In Fig. 2 the global relative variations in L 2 of σ 1 and σ 3 caused by the variation of the Young modulus in each unit are reported for each model configuration (i.e. nonhomogeneous lithology, single magma chamber, dual magma chamber, and dual magma chamber model with conduits) as follows: All the models show variability of less than 15 %, with a few exceptions within Unit B that have variability of over 30 % (Fig. 2). In this light, the spatial distribution of the major variations seems not to significantly affect the final stress distributions because (i) they are located near the mesh borders ( Fig. 3a and b), and (ii) when not at the mesh borders, the variations are limited to a few percent (Fig. 3c and d). It means that a 1 order of magnitude variation in the Young modulus produces variation in FEM outputs distributed over a large domain, and the change affecting the single nodes is limited to a few percent.

Homogeneous and non-homogeneous lithology
In Fig. 4 we report gravity-loaded models with homogeneous lithology composed of andesitic lavas (Fig. 4a) and nonhomogeneous lithology composed of carbonates (Unit B) and alluvial, volcaniclastic, and pyroclastic deposits (Units GF and VD; Fig. 4b). It is important to stress that the x-z zero displacement assigned at the bottom and at the lateral boundaries of the domain created substantial artefacts in the results (i.e. curved patterns of stress), especially considering σ 3 (Fig. 4, panels i and ii) whereby the boundary effect on the x axis is amplified by the presence of the upper free surface. It follows that the only unperturbed area extends ca. 30 km horizontally and ca. 15 km vertically (within the blue contour in Fig. 4). It is worth noting that the homogeneous and non-homogeneous models show quite similar stress pattern results (Fig. 4).

Gravitational modelling using inferred feeding system geometry
In Figs. 5 and 6 we show three cross-sectional profiles describing the feeding system starting from a single magma chamber to two chambers, then adding the conduits, and finally considering the full complexity by adding the effects of far-field stress and CR faults. Figure 5a describes the σ 3 (panel i) and σ 1 (panel ii) stress distribution for the single magma chamber model and P e = 10 MPa. No significant differences in the magnitude and pattern of stresses are visible using P e = 20 MPa (Appendix Fig. A1). The addition of the shallow magma chamber significantly changes the values and pattern of both σ 3 and σ 1 (Fig. 5b). In particular, σ 3 and σ 1 stresses describe a typical inflation pattern produced by excess pressure in the magma chamber(s) (Anderson, 1936;Gudmundsson, 2006Gudmundsson, , 2012, producing well-defined stress arches of σ 3 (red dotted lines in Fig. 5bi) and divergent strong gradients of σ 1 around the deep magma chamber (Fig. 5bii). Very slight differences in the magnitude and pattern of stresses appear when using P o = 10 MPa (Fig. 5b) or 20 MPa (Appendix Fig. B1). Looking at Fig. 6, it is evident how the insertion of the conduits in the CVC feeding system dramatically changes the stress distribution, with the disappearance of the stress arch and a nearly constant stress in the computational domain except around the deep magma chamber tips.

Application of an extensional stress field
In order to explore the influence of the extensional far-field stress on stress patterns (Fig. 1a), we ran simulations applying 5 MPa (a typical low value for rift zones; Turcotte and Schubert, 2002;Moeck et al., 2009;Maccaferri et al., 2014b; along the lateral boundaries of the computational domain (Fig. 7). In the case of a single magma chamber ( P e = 10 MPa; Fig. 7, panels i-ii), the addition of the far-field stress reduces the confinement effect due to the no-displacement condition imposed along the xz directions (plane strain approximation). When considering the double magma chamber configuration ( P o = 10 MPa in the deep chamber and P o = 5 MPa in the shallower one), the presence of the far-field stress produces slight changes in the stress magnitude and pattern for both σ 3 and σ 1 (Fig. 7, panels iii and iv) with respect to Fig. 5b. Very similar effects appears in the complete feeding system configuration model (Fig. 7, panels v-vi). Also, in this case using P o = 20 MPa in the deep magma chamber does not significantly affect the model outputs (Appendix B).

Faults bordering the Colima Rift
The effect of faults bordering the CR on the final feeding system configuration is simulated through two damage zones by degrading their elastic properties. Adding these elements does not significantly alter the stress distribution observed in Fig. 7v and vi but only provides a slight reduction in both σ 1 and σ 3 intensities around their edges ( Fig. 7vii and viii). The different distance of the two damage zones from the feed- ing system produces a small asymmetry in both σ 1 and σ 3 patterns with respect to simulations without damage zones, especially near the deep magma chamber (Fig. 7v-viii).

FEM analysis at increasing geological detail
The study highlights some important features of crustal stress distribution in terms of changing the geological and geophysical constraints used as input conditions (Spica et al., 2014(Spica et al., , 2017Massaro et al., 2018). Although the results have to be considered a first-order approximation, the changes in stress distribution are appreciable and useful for a better understanding of the advantages of FEM. Under the assumptions of plane strain and gravitational loading, the use of homogeneous or non-homogeneous lithology provides negligible effects on stress intensity and pattern (Fig. 4). This is likely due to the limited thickness of the shallow units (Units FC, VD, and GF; Table 2) in the simulated domain, the results of which are dominated by Unit B (Table 2). However, this does not mean that the influence of the upper units may still be negligible using smaller scales of the simulated domain. Analysing the single magma chamber model outputs, it emerges that P e limits the effects of gravitational loading. In contrast, the dual magma chamber geometry better describes the inflation induced by P e within magma chambers, with the formation of the stress arch in the σ 3 plot. It is worth noting that for both the single and dual magma chamber models, changing P e from 10 to 20 MPa slightly affects the magnitude of the stress but not its general pattern (Appendices A and B). The presence of dykes in the magma feeding system dramatically changes the σ 3 and σ 1 patterns (Fig. 6), which become quite homogeneous throughout the computational domain, with the only exception of sidewall effects induced by the zero displacement conditions. The addition of extensional field stress of 5 MPa reduces the sidewall effects and produces an almost homogeneous stress dis-tribution in the upper part of the computational domain above the top of the deep magma chamber. This, along with the additional inclusion of the damage zones introduced to mimic the effects of CR faults, describes a volcanic system close to equilibrium, in which pressure within the volcano feeding system almost equilibrates the lithostatic stress .

Some implications of the stress state of the CVC inferred from FEM
The results from the most complete FEM runs highlight an almost homogeneous stress distribution in the CVC area. This means the dual magma chamber model and the application of far-field stress provide a stable geometry, which limits the stress changes to a few megapascals (MPa). The majority of stress variations are located at the tips of the magma chambers, as expected for pressurized or under-pressurized cavities in the lithosphere (Marti and Geyer, 2009), implying that the whole feeding system is in a quasi-equilibrium state. Even if we consider the scenario of complete emptying of the upper conduit and part of the shallow magma chamber, as occurred occasionally during the past sub-Plinian and Plinian eruptions (Luhr, 2002;Saucedo et al., 2010;Massaro et al., 2018), this would result in the restoration of the stress arch, which is still a stable stress configuration. Even complete emptying of the shallow magma chamber would probably be ineffective for triggering a large collapse (caldera-forming) of the feeding system. Beside and beyond the limitations due to the first-order approximation of the FEM analysis, other sources of uncertainty in the discussion about the present and future stress state of the CVC come from not considering gravity-driven processes, such as volcano spreading due to plastic deformation of Unit GF (Norini et al., 2010(Norini et al., , 2019) and detailed regional tectonics (Norini et al., 2010(Norini et al., , 2019. The effect of the two fault systems bordering the CR is simulated here by degrading the mechanic properties of rocks in an area of about 1 km width up to a depth of 10 km. Although the effects are negligible at the scale of the computational domain, the possibility cannot be excluded that there are some local significant effects that cannot be resolved using the described approach.

Summary and conclusion
The presented study highlights the importance of using complete and detailed geological and geophysical data when dealing with FEM of volcanic areas. The different geological detail used in the model runs showed how the stress pattern critically depends on the geometry of the volcano feeding system, with huge differences in having a single or double magma chamber system and, in particular, whether or not the magma chamber(s) are connected to the surface by feeder dykes and conduits. The geometry of the feeding system is Figure 7. E-W gravitational modelling of the CVC domain with a non-homogeneous stratigraphy considering the extensional field stress. The magnitude and pattern of the principal stresses are shown for the single magma chamber model (i-ii), the dual magma chamber model (iii-iv), and the dual magma chamber model with conduits (v-viii). Note that in (vii-viii) the faults bordering the Colima graben are shown. For all configurations an extensive far-field stress of 5 MPa is applied at the lateral boundaries of the domain. In panels (vii-viii) the additional effect of the local extensive field is simulated using a reduced value of material properties (Table 2). P o is set to 10 and 5 MPa for the 15 km deep and 6 km deep magma chambers, respectively. Black dotted lines highlight the passage from different stress values. The red arrows indicate the direction of the applied far-field stress. Note that the scales of stress values are different for each panel in order to maximize the simulation details. prevalent in model outputs with respect to varying rock properties (i.e. Young modulus) of 1 order of magnitude. In the case of CVC, the use of subsurface homogeneous or stratified lithology does not influence the FEM outputs much, with the subsurface geology of the computational domain being dominated by carbonates (Unit B). Beside and beyond the results obtained by analysing the influence of detailed geological and geophysical data, the presented modelling confirms the close-to-equilibrium state of the volcano, which is the expected stress distribution induced by a feeding system directly connected to the surface. The complete emptying of the upper conduit and part of the shallow magma chamber has in the past led to sub-Plinian and Plinian eruptions. This process permits the restoration of the stress arch and maintains a stable subsurface stress configuration. It follows that large-magnitude, caldera-forming eruptions are possible only if the bigger deep magma chamber is also involved and significantly emptied during an eruption. Code availability. The LISA code is available at https://lisafea. com/ (last access: 10 December 2020; Sonnenhof Holdings (Canada), 2013).
Author contributions. SM, RS, AC, GN, and GG conceived the study. SM and RS wrote the bulk of the paper with the input of all the co-authors. SM and GL compiled the numerical simulations and formulated the adopted methodology. MP and SM carried out the sensitivity analysis. RS, AC, SM, GN, GG, LC, GL, MP, and AG worked on the interpretation of the results.
Competing interests. The authors declare that they have no conflict of interest.