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.

1 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 deformation 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., 2003, 2005; Buchmann 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., 2010, 2019; Cianetti 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., 2010, 2019). The FEM was run starting from simple homogeneous vs. stratified lithology of the subsurface in successively more detail through the addition of single and double magma chambers, feeder dykes, faults, and extensional far-field tectonic stress (Fig. 1b).

Figure 1 (a) Morphotectonic map of the Colima Volcanic Complex (NC: Nevado de Colima volcano; FC: Fuego de Colima volcano) and Colima Rift with the main tectonic and volcano-tectonic structures (NCG: northern Colima graben; CCG: central Colima graben; from Norini et al., 2019). In the inset, the location of the Colima Volcanic Complex (CVC) within the Trans-Mexican Volcanic Belt (TMVB) is shown in the frame of the subduction-type geodynamic setting of Central America (from Davila et al., 2019); (b) general sketch of the geometrical configurations used in LISA; (c) example of mesh of the investigated area for the dual magma chamber model with conduits (case v in b), 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).

2 The Colima Volcanic Complex, Mexico

2.1 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., 1996, 1997; Ferrari 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, 1986, 1991; Norini 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 right-lateral displacements of the bounding faults (Barrier et al., 1990; Suárez et al., 1994; Rosas-Elguera et al., 1996; Norini et al., 2010, 2019; Garduñ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., 2010, 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, 1986, 1991; Cortés et al., 2010; Norini et al., 2010; Escudero and Bandy, 2017). It is formed by three andesitic stratovolcanoes: Cantaro (2900 m a.s.l.), Nevado de Colima (4255 m a.s.l.), and, in the southern part, the youngest and active Volcán de Colima (3763 m a.s.l.) (Norini et al., 2019).

2.2 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 characterized by dome growths and collapses, extrusion of lava flows, and Vulcanian and occasionally sub-Plinian explosive eruptions (Saucedo et al., 2010; Massaro et al., 2018, 2019).

2.3 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., 2013, 2019; Mací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³. 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., 2010, 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, eastward-dipping normal fault exposed on the northern CVC flank (Norini et al., 2010, 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 (Saucedo et al., 2010, 2011; 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).

3 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 et al., 1999). Despite LISA originally being used for structural analysis (Rao, 1989, 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).

Table 1 Element types used in LISA considering the final conduit feeding system configuration and the E–W cross section (a–a^′ – Fig. 1d, panel vi). The total number of elements is 4209.

3.1 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 finite-element 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. (2010, 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, 2017; Massaro et al., 2018, 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 σ₁ as the greatest compressive stress and σ₃ as the least compressive stress.

Table 2 Rock mass and mechanical properties of the geological units used in FEM from Norini et al. (2010, 2019).

4 Geological data

4.1 Stratigraphy and rock mechanics

Four units forming the CVC system were defined from the available geological data (Table 2): (i) a basement (Unit B) of Cretaceous limestones and intrusive rocks forming the bedrock underlying the CVC; (ii) graben fill deposits (Unit GF): Quaternary alluvial, colluvial, and lacustrine deposits filling the graben; (iii) Fuego de Colima deposits (Unit FC) of andesitic lavas and pyroclastic deposits forming the Paleofuego–Fuego de Colima edifices; and (iv) volcaniclastic deposits (Unit VD) covering the southern flank of the CVC (e.g. Cortés et al., 2010; Norini et al., 2010, 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; Norini et al., 2010, 2019). For Units B and GF rock mass proprieties are derived from 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.

Figure 2 Results of the sensitivity analysis carried out on the Young modulus variations within each rock layer of the domain considering different configurations (stratified substratum model – nodes: 4426; single magma chamber model – nodes: 4426; dual magma chamber model – nodes: 4161; dual magma chamber model with conduits – nodes: 3737). For each geological unit (B, FC, GF, VD), the relative global variation in L₂ (%) is provided for σ₁ and σ₃. The x₍₋₎ and x₍₊₎ vectors indicate the Young modulus variation by an order of magnitude with respect to the x_ref vector, containing the stress values calculated by using the values of material properties indicated in Table 2.

Figure 3 Spatial variation (%) of the L₂ norm's components at varying Young modulus values for the selected cases of Units B and VD: (a) Unit B in the stratified substratum model (nodes: 4426); (b) Unit B in the single magma chamber model (nodes: 4426); (c) Unit B in the dual magma chamber model (nodes: 4161); (d) Unit VD in the dual magma chamber model with conduits (nodes: 3737). Symbols x₍₋₎ and x₍₊₎ have the same meaning as in Fig. 2.

4.2 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). Various estimates have been made to infer the actual amount of melt in a magmatic body, showing that it is only ca. 10 % of the total chamber volume (Gudmundsson, 2012, and references therein). After Spica et al. (2017), the 15 km deep low-velocity body (LVB) is ca. 7000 km³. Therefore, if we assume the melt as 10 %, the deep magma chamber volume would be ca. 700 km³. Simplifying this volume in an elliptical sill-like geometry, the magma chamber dimensions (i.e. $\mathrm{2}a,\mathrm{2}b,\mathrm{2}c$ axes) have to be scaled according to those of the LVB (55 × 30 × 8 km; Spica et al., 2017) using 2a=14 km, 2b=3.6 km, and 2c=26 km, with 2c being elongated in the NW–SE direction. For the shallow part of the feeding system, we have no detailed geophysical constraints. However, Massaro et al. (2019) reproduced the nonlinear cyclic eruptive activity at Fuego de Colima in the last 20 years through numerical modelling using a shallow magma chamber volume in the range of 20–50 km³, also according to the estimation of Cabrera-Gutiérrez and Espíndola (2010). Here we assume a volume of 30 km³ using 2a=3.5 km, 2b=2 km, and 2c=8 km as the dimensions of the shallow magma chamber. Numerous theoretical and field studies have established that host rock stresses dictate the magma pathways (e.g. Maccaferri et al., 2011; Gudmundsson, 2011). During ascent to the surface, the dykes align themselves with the most energy-efficient 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 $\mathrm{2}bd=\mathrm{2}b=\mathrm{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.

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

5.1 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₂) method to illustrate the results. The L₂ norm applied on a vector space x (having components $i=\mathrm{1},\mathrm{\dots }n$) is strongly related to the Euclidean distance from its origin and is equal to

$$\begin{array}{}\text{(1)}& \left|\right|x|{|}_{\mathrm{2}}=\sqrt{\sum _{i=\mathrm{0}}^n{x_i}^{\mathrm{2}}}.\end{array}$$

Figure 4 E–W gravitational modelling of the CVC domain. The scale of the mesh is expressed in units of design (1 UD = 1 km). The domain extends 60 km along the x axis and 30 km along the z axis. The number of nodes used in the mesh is set to 4426. The magnitude and pattern of the principal stresses (dotted black lines) are reported for (a) the homogeneous stratigraphy (Unit FC: andesitic lavas and pyroclastic deposits) and for (b) the non-homogeneous stratigraphy (Unit FC; Unit B: Cretaceous limestones and intrusive rocks forming the bedrock underlying the CVC; Unit GF: Quaternary alluvial, colluvial, and lacustrine deposits filling the graben; Unit VD: volcaniclastic deposits covering the southern flank of the CVC). The blue line contours the unperturbed part of the domain, which extends ca. 30 km horizontally and ca. 25 km vertically. Note that the scale of stress values is the same for the all simulations.

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 minimum 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₂ of σ₁ and σ₃ caused by the variation of the Young modulus in each unit are reported for each model configuration (i.e. non-homogeneous lithology, single magma chamber, dual magma chamber, and dual magma chamber model with conduits) as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }\left|\right|x|{|}_{\mathrm{2}(-)}={\displaystyle \frac{x_{\text{ref}}-x_{(-)}|{|}_{\mathrm{2}}}{\left|\right|x_{\text{ref}}|{|}_{\mathrm{2}}}}\text{(3)}& {\displaystyle }\left|\right|x|{|}_{\mathrm{2}(+)}={\displaystyle \frac{x_{\text{ref}}-x_{(+)}|{|}_{\mathrm{2}}}{\left|\right|x_{\text{ref}}|{|}_{\mathrm{2}}}}\end{array}$$

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.

5.2 Homogeneous and non-homogeneous lithology

In Fig. 4 we report gravity-loaded models with homogeneous lithology composed of andesitic lavas (Fig. 4a) and non-homogeneous 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 σ₃ (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).

Figure 5 E–W gravitational modelling of the CVC domain with a non-homogeneous stratigraphy. The magnitude and pattern of the principal stresses are reported for (a) the single magma chamber model represented by a magma chamber (2a=14 and 2b=3.6 km) at 15 km of depth and (b) the dual magma chamber model composed of a 15 km deep magma chamber (2a=14 and 2b=3.6 km) and a shallow 6 km deep one (2a=3.5 and 2b=2 km). The magma chambers are not connected. ΔP_o is set to 10 and 5 MPa for the 15 km deep and 6 km deep magma chambers, respectively. The number of nodes is set to 4426 and 4161 for the single and dual magma chamber models, respectively. Black dotted lines highlight the passage from different stress values. The red dotted line in (b) (panel i) indicates the formation of the stress arch. Note that the scales of stress values are different for each panel in order to maximize the simulation details.

5.3  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 σ₃ (panel i) and σ₁ (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 σ₃ and σ₁ (Fig. 5b). In particular, σ₃ and σ₁ stresses describe a typical inflation pattern produced by excess pressure in the magma chamber(s) (Anderson, 1936; Gudmundsson, 2006, 2012), producing well-defined stress arches of σ₃ (red dotted lines in Fig. 5bi) and divergent strong gradients of σ₁ 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.

Figure 6 E–W gravitational modelling of the CVC domain with a non-homogeneous stratigraphy accounting for a dual magma chamber system connected by dykes via the surface (deep magma chamber, 2a=14 and 2b=3.6 km at 15 km of depth; shallow magma chamber, 2a=3.5 and 2b=2 km at 6 km of depth). The magnitude and pattern of the principal stresses are shown. The number of nodes used is set to 3737. ΔP_o is set to 10 and 5 MPa for the 15 km deep and 6 km deep magma chambers, respectively. The black dotted lines in (ii) the passage from different stress values. Note that the scales of stress values are different for each panel in order to maximize the simulation details.

5.4 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; Sulpizio and Massaro, 2017) 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 x–z 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 σ₃ and σ₁ (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).

5.5 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 σ₁ and σ₃ intensities around their edges (Fig. 7vii and viii). The different distance of the two damage zones from the feeding system produces a small asymmetry in both σ₁ and σ₃ patterns with respect to simulations without damage zones, especially near the deep magma chamber (Fig. 7v–viii).

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.

6 Discussion

6.1 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, 2017; Massaro 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 σ₃ 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 σ₃ and σ₁ 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 distribution 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 (Sulpizio et al., 2017).

6.2 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, 2019) and detailed regional tectonics (Norini et al., 2010, 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.

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

Appendix A

Figure A1 E–W gravitational modelling of the CVC domain (stratified lithology) for all configurations investigated. The magnitude and pattern of the principal stress account for the (i–ii) single magma chamber model (number of nodes: 4426), (iii–iv) dual magma chamber model (number of nodes: 4161), and (v–vi) dual magma chamber model with conduits (number of nodes: 3737). The dimensions of the deep magma chamber are 2a=14 and 2b=3.6 km at 15 km of depth; shallow magma chamber: 2a=3.5 and 2b=2 km at 6 km. ΔP_e and ΔP_o are equal to 20 MPa for the deep chamber and 5 MPa for the shallower. Black dotted lines highlight the passage from different stress values. Note that the scales of stress values are different for each panel in order to maximize the simulation details.

Appendix B

Figure B1 E–W gravitational modelling of the CVC domain (stratified lithology) considering a far extensional stress field of 5 MPa for all configurations investigated. The magnitude and pattern of the principal stress account for the (i–ii) single magma chamber model (number of nodes: 4426), (iii–iv) dual magma chamber model (number of nodes: 4161), (v–vi) dual magma chamber model with conduits (number of elements: 3737). The dimensions of the deep magma chamber are 2a=14 km and 2b=3.6 km at 15 km of depth; shallow magma chamber: 2a=3.5 km and 2b=2 km at 6 km. ΔP_e and ΔP_o are equal to 20 MPa for the deep chamber and 5 MPa for the shallower. 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.