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.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 2Rock mass and mechanical properties of the geological units used in FEM from Norini et al. (2010, 2019).

Download Print Version | Download XLSX

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

Download

Figure 3Spatial 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.

Download

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

Figure 4E–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.

Download

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:

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

Download

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 6E–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.

Download

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

Download

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.