Ruz et al

Volcano tectonics of the Southern Andes

Most of the Chilean Andean margin (−18.4°S to −46.3°S) is currently characterized by dextral oblique convergence (ca. N78°E at a rate of ca. 6.6 cm/year) between the Nazca and South American plates, that has been consistently maintained during the last 20 Ma (Pardo-Casas and Molnar, 1987) (Fig. 1). Subduction related processes in the Chilean Andean margin have been ongoing since at least the Early Jurassic (~200 Ma) (Stern, 2004) when the volcanic arc was in the position of the present-day Coastal Cordillera. Combinations of several platemargin-scale processes, throughout this time, resulted in an overall eastward migration of the volcanic arc to its present location (e.g. Stern, 2004).

The Southern Volcanic Zone (SVZ) makes up the volcanic arc of the Andes between 33°S and 46°S and is divided into four segments (Stern, 2004). These segments form in close spatial relation to latitudinal differences in basement characteristics such as lithospheric elastic thickness (Tassara and Yañez, 2003;Stern, 2004), and petrographical, geochemical and isotopic variations of erupted products (Hildreth and Moorbath, 1988;López-Escobar et al., 1995;Dungan et al., 2001;Cembrano and Lara, 2009). The southernmost segments (37°S-46°S) are constructed over Cenozoic plutonic and metamorphic rocks (e.g. Charrier et al., 2002). Intra arc tectonics in the southernmost segments is dominantly governed by the margin parallel Liquiñe Ofqui Fault System (LOFS) (Fig. 1). The northernmost segments' (33°S-37°S) of the SVZ are characterized by a thick Meso -Cenozoic volcano sedimentary cover (e.g. Tapia et al., 2015;Cembrano and Lara, 2009 and references therein) and overlie margin-parallel Cenozoic folds and thrust faults that represent the westernmost expression of the Aconcagua and Malargüe fold and thrust belts (e.g. Giambiagi and Ramos, 2002).

Active tectonic features in the transitional segment of the SVZ (34.5°−37°S) are either poorly developed margin-parallel strikeslip faults or well-defined Andean Transverse Faults (ATF) (e.g. Piquer et al., 2018;Sielfeld et al., 2019b) (Fig. 1). Current NEtrending shortening at~36°S is suggested by the focal mechanism of the Mw 6.2 crustal earthquake on the~N10°E El Melado fault (Cardona et al., 2018) (Fig. 2). The geometry and nature of this crustal fault is constrained by >600 aftershocks and by a dextral strike-slip moment tensor (from NEIC and GCMT catalogues). Cardona et al. (2018) also registered crustal events along the WNW-and NE-striking Laguna Fea and Troncoso faults respectively, in the SW corner of the Laguna del Maule (Fig. 2). Active seismicity to the north of the study area (~35.5°S) also shows margin parallel and transverse-to-the-arc distribution (Pearce et al., 2020).

At this latitude (36°S), the Principal Cordillera consists mainly of Mesozoic and Cenozoic volcanoclastic and sedimentary units, as well as Late Mesozoic -Early Cenozoic granitoids (González and Vergara, 1962;Muñoz and Niemeyer, 1984). These units were deformed by regional scale thrusting and folding, with Oligocene -Miocene sedimentary sequences characterized by~NS-trending asymmetric folds (Sielfeld et al., 2019b), whereas Pliocene -Early Pleistocene volcanism is recognized as subhorizontal pyroclastic layers (Sielfeld et al., 2019b). Sedimentary sequences and granitoids are also intruded by pre-Quaternary dykes that are exposed in the Maule river-Los Cóndores canyon valley.

Active Quaternary volcanism in the study area is represented by the Descabezado Grande volcanic field, Tatara -San Pedro volcanic complex (TSPP), Laguna del Maule volcanic field (LMVF), and the Nevados de Longaví volcanic complex (e.g. Stern, 2004) (Fig. 2). This constitutes one of the widest segments of the SVZ, up to 60 km wide in the main arc and over 200 km wide in the back-arc region (Stern, 2004).

Tatara-San Pedro Pellado (TSPP) complex

The TSPP volcanic complex (36°S) is an ENE trending alignment of two stratovolcanoes (Tatara-San Pedro and Pellado-San Pablo summits from west to east), minor eruptive centres, lava flows, pyroclastic deposits, and dyke swarms that cover an area of over 387 km 2 (Singer et al., 1997). San Pedro volcano is the youngest and most westerly volcanic centre of the complex. Volcanic products erupted from at least three of the central vent regions include both lavas and pyroclastic deposits, and compositionally range from basalts to rhyolites, with mean compositions within the basaltic-andesitic range (ie. 52-56% SiO 2 ; Davidson et al., 1988;Singer et al., 1997, Dungan et al., 2001Costa and Singer, 2002). Volcanism in the TSPP commenced at~925 ka but preserved lavas (~54 km 3 ) represent only a fraction of this time interval (Singer et al., 1997). This implies that either volcanic quiescence, erosional mechanisms such as glacial abrasion, or a combination of both processes have removed significant parts of the eruptive record (Singer et al., 1997). Except for small neoglacial moraines, San Pedro volcano is unglaciated, and therefore is of Holocene age (Costa and Singer, 2002).

ENE-striking dyke swarms outcrop in the surroundings of the active San Pedro volcano and cut both pre-eruptive sequences and Pleistocene volcanic units of the TSPP complex Sielfeld et al., 2019aSielfeld et al., , 2019bEspinosa, 2019). The dykes vary in composition from andesite to rhyolite, and likely represent parts of the magmatic plumbing system that have fed the TSPP. Singer et al. (1997) mapped abundant ENE-striking dykes in the northern side of the complex. These dykes intrude Pleistocene dacitic lavas of the TSPP through normal faults (Singer et al., 1997). Moreover, Espinosa (2019) mapped ENE to EW-striking and NW-striking dykes to the south of the TSPP. Sielfeld et al. (2019b) suggest that Quaternary volcanism of the TSPP and the associated MGS have been constructed on top of a ENE-trending damage zone (namely the Tatara Damage Zone), a long-lived structural anisotropy of cross-cutting mesoscale ENE-WNW striking (Angermann et al., 1999). Rupture zones of the 1960 Valdivia Mw 9.5 and 2010 Maule Mw 8.8 megathrust earthquakes are outlined in lilac dashed lines. Moment tensors of the GCMT catalogue include subduction-type earthquakes (cyan-white beach-balls) and crustal events Mw > 5 and depth < 20 km (blue-white beach-balls). Crustal seismicity in the forearc region at 38°-39°S latitude extracted from Haberland et al. (2006). Size of these beach balls is scaled to MI. Red triangles indicate the location of active volcanoes (taken from the Smithsonian Institute Holocene Volcanic database) and yellow circles indicate the location of geothermal springs (taken from Aravena et al., 2016). The location of the TSPP is indicated with a green triangle. Crustal faults systems such as the Liquiñe-Ofqui Fault System (LOFS) and Andean Transverse Faults (ATF), and morphotectonic lineaments from Lara (2009), Pérez-Flores et al. (2016) and Sielfeld et al. (2019a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) transtensional oblique slip faults and dykes, characterized by a grabenlike geometry, that has developed before and during the construction of the TSPP. Strain field solutions from fault-slip data indicate a steeply to vertically plunging shortening axes and a horizontal NS-NNW trending stretching axes, thus indicating a transtensional to extensional regime (Fig. 2). This is in agreement with stress field solutions from Sielfeld et al. (2019b).

Physical volcanology/stratigraphy of the TSPP

The TSPP succession has been subdivided into 17 chronostratigraphic units (grouped in 8 sequences described here) based on geologic mapping, geochemistry and geochronology (K\ \Ar and 40 Ar\ \ 39 Ar) (Davidson et al., 1988;Singer et al., 1997;Dungan et al., 2001;Costa andSinger, 2002, andHildreth et al., 2010). From oldest to youngest, these sequences are: (1) Muñoz sequence (925-825 ka) composed of a silicic lava dome-flow complex, with phyric lavas and flow breccias, followed by basaltic andesitic and dacitic lavas, and rhyolites. The total preserved volume of the erupted material is 4.8 km 3 ;

(2) Quebrada Turbia (783-772 ka) consists of~1.8 km 3 basaltic andesitic to andesitic lava flows; (3) Estero Molina lavas (605-329 ka) composed of~4.2 km 3 preserved basaltic andesites and rare primitive basaltic flows; (4) Cordón Guadal lavas (512-359 ka) consist of basaltic to dacitic lavas; (5) Placeta San Pedro sequence (~230 ka) composed of 2 km 3 of a lower unit of glassy andesitic lava flows with mafic magmatic inclusions and an upper unit made of basalt flows with olivine and Singer et al. (1997) (for the detailed stratigraphy of the TSPP). Observed and inferred faults, folds and dykes are taken from Sielfeld et al. (2019b). Moment tensor of the 2012 Mw 6.2 crustal earthquake associated to the Melado fault is taken from the GCMT catalogue. Red line indicates the conductive anomaly of the Mariposa Geothermal System, taken from Hickson et al. (2011). Stereographic diagrams contain fault-slip data and strain field solutions for 6 selected sites around the TSPP (fault-slip data taken from Sielfeld et al. (2019b)). In each diagram, the maximum shortening and stretching axes (principal axes of average incremental strain) are shown with blue triangles and red stars respectively. Sites are indicated by white boxes and identified with numbers from 1 to 6. Location of the studied dyke swarm is indicated. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) clinopyroxene phenocrysts; (6) Pellado volcano sequence (188-83 ka) is made of lahars, lavas and pyroclastic rocks ranging in composition from basaltic andesite to andesite with a preserved volume of 12.5 km 3 , that erupted from vents in the vicinity of Cordón San Pablo and Cordón Pellado; (6) Tatara volcano sequence (90-19 ka) comprises 22 km 3 of basaltic andesitic to andesitic lava flows, and minor pyroclastic flows; (8) San Pedro volcano sequence (< 10 ka) is the youngest unit. Eruptive activity included early explosive eruptions with pumice, andesitic and dacitic lavas, followed by Strombolian eruptions producing basaltic andesitic scoria and lavas. Part of these lava flows fill a topographic depression generated by the collapse of the early Holocene cone.

Methods

The objectives of this research were addressed through a combination of structural field work, paleo-stress reconstruction and twodimensional finite element modelling. We took measurements from the well-exposed dyke swarm (0.4 km 2 in area) roughly 5 km to the east of Pellado volcano. Our data complements previous regional scale structural mapping from the surrounding area (Sielfeld et al., 2019b).

Field methods

Detailed mapping of the dyke segments at the TSPP plateau were conducted to diagnose variations in dyke geometry, mineralogy, attitude and texture. This mapping included systematic measurement of attitude (strike and dip), thickness and length on each dyke segment. Thickness and segment length were measured directly in the field using a measuring tape when possible. For longer segments, lengths were measured from Unmanned Aerial Vehicle (UAV) images and from a collection of marked GPS points at the start and end of each segment. Where it was not possible to collect field data due to access restrictions, data was recorded through remote sensing data (resolution of WorldView-1 image is 0.5 m/pixel). Sampling of representative dyke outcrops was carried out for further petrographical and fabric analyses at the thin section scale.

Stress reconstruction: paleostress analysis from dilatational fracture data

Here we combined the field measurements of dyke attitude with the GArcmB code (http://www.kueps.kyoto-u.ac.jp/~web-bs/tsg/software/ GArcmB/) to obtain the stress conditions (orientation of principal stresses and stress ratio) at the time of dyke formation (Faye et al., 2018;Yamaji and Sato, 2011;Jolly and Sanderson, 1997). The code clusters the distribution of poles of dykes (i.e. the vector normal to the dyke plane, assumed to be co-axial with σ 3 ) by fitting the dataset with a linear combination of Bingham distributions (mixed Bingham distribution), each with its own concentration axes (κ 1 , κ 2 ) and parameters, and interpreted as representing a stress state (Yamaji, 2016;Yamaji and Sato, 2011). Although the method makes a series of assumptions, that are likely a simplification of most natural cases, it provides a firstorder approximation of the stress state at the time of dyke formation (a detailed explanation of the method can be found in the supplementary material).

Finite element numerical method

We built a suite of plane-strain 2D finite element mechanical models to calculate the static stress field formed around both a volcanic complex and a dyke swarm under different regional loads in order to understand the observed pattern of dyke emplacement. We used the finite-element code ADELI (Hassani et al., 1997), which has previously been used to simulate a variety of geodynamic problems, including long-term deformation at subduction zones (Hassani et al., 1997;Cerpa et al., 2015), interseismic deformation (Contreras et al., 2016) and volcanic inflation (Novoa et al., 2019). This code uses a timeexplicit dynamic relaxation method (Cundall and Board, 1988) to solve the quasi-static equation of motion and has been thoroughly tested through a variety of elastic, visco-elastic and elasto-plastic rheology. In our models here we consider only elastic rheology to determine the conditions that influence the stress field which in turn drives dyke emplacement. More details on the ADELI code can be found in Chery et al. (2001), Cerpa et al. (2015), Gerbault et al. (2018) and Novoa et al. (2019).

We designed three types of model geometry that simulate 1) the volcanic complex as a sub-circular magmatic-hydrothermal reservoir and pressure source (Model-type 1) (Fig. 3a), 2) the dyke swarm as an elongated pressure source on the margins of the complex (Modeltype 2) (Fig. 3b), and 3) an amalgamation of both geometries with the emplacement of a younger dyke swarm at a greater distance from the complex (Model-type 3) (Fig. 3c). Since the modelled domain covers a plane-view section of the tectono-magmatic system at depth, the initial lithostatic equilibrium can be approximated by an initially null stress tensor, hence the effect of gravity is equivalent at constant depth and does not affect the elastic behaviour. All models assume an elastic domain where the crustal section has a density (ρ) of 2750 kg/m 3 , a Poisson's ratio (ν) of 0.27 and a Young's modulus (E) of 5 GPa, which are all reasonable values for crustal segments hosting volcanoes (Gudmundsson, 2011;Heap et al., 2020).

The specific model geometries are as follows:

• Model type 1: The reservoir is modelled as a circular cavity of radius 1 km. The walls of the cavity are assigned a magmatic overpressure of 10 MPa. The dimensions of the cavity and distance to the elliptical cavities (ie. dyke swarm) are roughly estimated from the location of a conductive anomaly associated with a magmatic reservoir beneath the TSPP (Reyes-Wagner et al., 2017). • Model type 2: The active dyke swarm is modelled as an elliptical cavity with a length of 1 km and a thickness of 150 m (diameter of long and short axes), and a strike of N80°E. A magmatic overpressure of 10 MPa was applied to all of the walls of the cavity. • Model type 3: The chilled dyke swarm is modelled as a meshed ellipse with a stiffness (Young's moduli) of 50 GPa, i.e. 5 times higher than the rest of the modelled domain. The new active dyke swarm is modelled as an elliptical cavity with the same dimensions as those of Model-type 2, and a strike of N60°E. This second dyke injection is modelled from the likely propagation path derived from Model-type 2 and field observations, thus considers the remnant stress generated from Model-type 2 over a time interval that allows the cooling of the first dyke injection (the results of a model considering two active dyke swarms can be found in Section 3 of the supplementary material).

The crustal properties are kept constant throughout all of the model simulations. The size of the magmatic-hydrothermal reservoir is constrained by the shape of conductive bodies obtained from a magnetotelluric profile (Reyes-Wagner et al., 2017), while the orientation of dykes was measured directly in the field. The lengths of the dykes in the models are deliberately one order of magnitude larger than measured in the field for purposes of mesh resolution. As we assume plane-strain deformation the model is symmetrical in the third, vertical dimension, hence the cavities and ellipses are cylindrical features, similar to previous assumptions such as by Gudmundsson and Andrew (2007), Andrew and Gudmundsson (2008), Gerbault et al. (2012).

The model geometries are meshed with the Gmsh software (http:// gmsh.info/) using tetrahedral elements (the third dimension is meshed to a thickness of 0.75 km but plane-strain conditions leads to effectively 2D models). The mesh resolution is highest in areas near the cavity boundaries. The total number of elements is in the order of tens of thousands and slightly varies from one model to the other.

All of the three model geometries are loaded using two different boundary conditions. They are as follows and shown schematically in Fig. 3 (d, e):

1. Loading condition 1: Internal magmatic excess pressure of 10 MPa applied to the boundaries of both the circular reservoir and elongated dyke, in order to simulate fluid overpressure. As we assume that each cavity is filled with melt, the tangential stress at the cavity boundaries is zero. The other boundaries of the model are set free slip, i.e. they are only fixed in their normal direction (i.e. avoiding rigid body rotation) (Fig. 3d).

2. Loading condition 2: A combination of internal magmatic excess pressure of 10 MPa applied on the boundaries of the reservoir and dykes, and a horizontal compression (Sx) of 1 MPa is applied at the horizontal boundaries of the model domain (X direction, mimicking E-W compression). This horizontal compression is applied as a velocity normal to the vertical borders, so that the total displacement (shortening) at the end of the model run reaches Dx~4.3 m. With free-slip horizontal borders, this loading condition builds uniaxial strain, and yields a background stress state of 1 MPa within the homogeneous domain.

Results

We now present results from all of the three types of analyses that were carried out: field observations (Section 4.1), paleo-stress reconstructions (Section 4.2), and the numerical modelling (Section 4.3). The results are further synthetized together in Section 5.

Field observations and analyses

The dykes are exposed over a narrow area (approximately 400 m in width and 1 km in length) (Figs. 4; 5a) and cross cut 103 ka old basaltic lavas of the Pellado volcano stratigraphic sequence (Singer et al., 1997). To the south of the dyke swarm exposure extends a plateau at~2500 m. a.s.l. where eroded lava flows and scoria of various hues crop out. Some dyke segments of the swarm outcrop to the north yet are inaccessible due to steep-sided slopes. Dyke segments were easily identified because they rise sharply over the eroded lavas which they intrude.

In the central part of the dyke swarm we observe the remnants of a steep-sided scoria cone of approximately 358 m length and 175 m width (Fig. 5b). The cone is made of distinct units of multidirectional jointed lavas, topped by near-horizontal units of scoria (Fig. 5c). On the upper most surface of the cone we observe glacial striations trending ENE and NNW (Fig. 5b). Within a section of wall of the eastern most side of the cone we observed three dyke segments. One of the dykes is~5 m-thick and transitions into a lava feeder near the surface of the cone (Fig. 5d). The other two dyke segments are~1 m in thickness and are arrested within the cone, making it possible to identify their tips. These segments did not feed an eruption (Fig. 5e).

In total, we measured 35 dyke segments which range in thickness from 0.56 m to 7.15 m (Table 1) (Fig. 6). The total length of the dyke swarm is~1 km, but individual dyke segments range in length between 4.3 m to 60.6 m (Fig. 6). In the field it was not possible to identify dyke tips along their strikes, and therefore we do not rule out the possibility that these dyke lengths are underestimated, covered by soils, unconsolidated sediments and modern debris. Based on the orientation, thickness and mineralogy we interpret that the 35 segments belong to only nine dykes, each comprised of between 2 and 5 segments. This results in total outcropping dyke lengths between 40 m and 170 m (Fig. 4b).

Each dyke is distinguished with a number in Table 1 and Fig. 4. The segment of each dyke is also recorded.

The dyke segments exposed in the area all dip steeply betweeñ 70-90°, and strike predominantly between N50°E and N90°E, with some segments striking WNW (n = 8) (Fig. 6). We can distinguish two preferential orientations at N50°-70°E (n = 10) and N70°-90°E (n = 12) indicating an apparent rotation in strike orientation of about 20°. The distance over which segmentation and changes in strike orientation occur, are within the same order of magnitude as the dyke thickness, i.e. cm to m. Such changes also relate to dyke segments with similar textures and mineralogy (Fig. 7b, d). Poor cross-cutting relations between these predominant orientations suggests that the N50°-60°E striking dykes are younger than the E-W-striking dykes (Fig. 7c). There does not appear to be a clear en-échelon segmentation of the dyke segments but the~20°rotation in strike orientation of the dyke segments could reflect an incipient en-échelon pattern.

The composition of all the dykes is porphyric basalts with a plagioclase, clinopyroxene and orthopyroxene phenocrystic phase and small Table 1). Segmented black boxes indicate the location of the dyke segments shown in Figs. 5 and 7. (b) Interpretative map of the studied dyke swarm. Colours indicate to which "total dyke length" each segment is associated to. Numbers indicate the dyke segment ID. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) amounts of olivine in a crypto crystalline groundmass of plagioclase, pyroxenes and glass. Textural differences, such as grain size and vesicularity, are observed among dyke segments with different orientations. We did not observe evidence of more evolved magmatic compositions in the swarm. In some samples, micro plagioclase crystals in the groundmass show a preferred orientation, parallel to the dyke's edge. Most samples possess vesicles but only one dyke exhibited vesicles elongated parallel to the margin (a more detailed petrographical description of some dyke segments can be found in the supplementary material). Moreover, some dyke segments (5.1-5.5) have a~40 cm chilled margin and brecciated clastic edges, with clast size of 2-6 cm. They also display quench textures in the border, indicative of rapid cooling. There is no clear breccia texture between the~40 cm margin and the centre of the dyke, but a smooth textural transition.

Paleo stress reconstruction

To reconstruct the paleo stress conditions during dyke emplacement we performed an inversion from the orientation (strike and dip) of 35 dyke segments. The data set shows an optimal K value of 1, meaning that one single stress condition can be used to explain the orientation of all the dyke segments. This stress state corresponds to a stress regime with a horizontal NW-trending minimum compressive principal stress (σ 3 ) (336°/0.4°), a horizontal NE-trending intermediate compressive principal stress (σ 2 ) (066°/0.8°), and a vertical maximum compressive principal stress (σ 1 ) (219°/89.1°). The resulting stress shape ratio is Φ = 0.1, indicating that σ 2 ≈ σ 3 (Fig. 8).

Model results

A total of 6 numerical models were run by varying the geometrical setups (Model types 1, 2, and 3) and loading conditions (1 and 2). Here we discuss the results after applying 1) 10 MPa of excess pressure only along the cavity boundaries, and 2) applying both 10 MPa excess pressure along the cavity boundaries and a 1 MPa E-W compression applied from the model edges. To examine possible propagation paths and interactions between dykes, we focus on the magnitude of σ 3 , or tensile stress, and the orientation of σ 1 . We also display the volumetric strain, which is the first invariant of the strain tensor and expresses its isotropic component. Through analysis of principal stresses and volumetric strain, zones of both contraction and dilation can be discerned.

Model-type 1: circular cavity

In Fig. 9 we show the results of the numerical simulations using the geometry of Model type 1 and the two end-member loading conditions. In loading condition 1, with only the application of an excess pressure of 10 MPa from the circular cavity walls, the tensile stress (σ 3 ) decreases radially from a maximum value of 10 MPa at the walls of the cavity to almost zero in areas >2 km away. The orientation of σ 1 trajectories is radial (Fig. 9a). These results reproduce analytical solutions in which for plane-strain, the stress field decays as 1/R 2 (e.g. Timoshenko and Goodier, 1970, see also benchmarks in Gerbault et al., 2012). Fig. 9 (c, d) shows the result of loading the same geometry with an overpressure in the circular cavity and an E-W compression of 1 MPa (loading condition 2). The effect is an overall increase in the magnitude of σ 3 in the bedrock, decreasing from a maximum value of~25 MPa at the cavity walls to~1 MPa in areas away from the cavities. The Table 1 Dyke segment parameters as measured in the field, from UAV images and from a collection of marked GPS points. Dykes are identified by two numbers: the first indicates dyke number and the second indicates the segment number. Last column indicates the aspect ratio of 33 dyke segments. (*) Dykes are exposed in a vertical section, thus length is not an observable feature.

Segment

Strike ( trajectories of σ 1 still rotate around the circular cavity as in the model with loading condition type 1, but now align more closely along the orientation of the applied horizontal compression in areas away from the cavity (Fig. 9c). Zones of contraction and dilation develop to the north-south and east-west of the cavity respectively (Fig. 9d), with volumetric strain decaying away from the cavity walls.

These results indicate that regions to the east-west of the cavity would favour magma propagation through dykes striking ENE and WNW at locations close to the cavity (~2 km) and striking E-W at locations further away from the cavity (>2 km). Regions to the north and south of the cavity undergo uniform compression, and also exhibit unfavourable orientations of σ 1 , thus dyke propagation would be inhibited to the north and south of the magmatic cavity. These results are in agreement with classic mechanical studies by Griffith (1921) (among others) and recent experimental results obtained by Li et al. (2017) for marble specimens with circular holes, where the formation of tensile cracks occurs in the loading direction due to the increase of micro crack density. Fig. 10 shows the results of the model geometry that considers a plane-view section of both the volcanic complex and a nearby radial dyke swarm striking N80°E (Model-type 2). This model-type represents the first stage of a dyke swarm injection as discerned from the field measurements. When an excess pressure of 10 MPa is applied in both the circular and elliptical cavities, the stress field develops and decays away from these cavities. Tensile stress (σ 3 ) concentrates at the tips of the dyke swarm (Fig. 10a), as well as dilatational strain (Fig. 10b). In contrast, along the elongated sides of the dyke swarm we note contractional strain. The trajectories of σ 1 are distributed radially around the circular and elliptical cavity; particularly, at the tip of the dyke they are oriented parallel to the direction of the major semi-axis of the dyke, while along its long edges, the trajectories are oriented perpendicular. This is in agreement with the volumetric strain pattern. Fig. 10 (c, d) shows the result from loading the same geometrical configuration with loading condition (2) (e.g. applied cavity overpressure plus lateral compression). The region located in between the inflating cavities experiences a larger amount of tensile stress and hence dilatational strain (Fig. 10 c, d). The σ 3 stress contour of 2-3 MPa links both cavities. The trajectories of σ 1 rotate to a near EW orientation at the east of the dyke tip. Ahead of the dyke tip, the trajectories of σ 1 strike within ±~13°to the orientation of the dyke strike, rotating within 4°to 20°(clockwise and anticlockwise rotations) with respect to the trajectories of σ 1 obtained with loading condition (1). Within 1.5-2 km (~dyke length) to the north and south of the dyke, trajectories of σ 1 rotate by 7°to 30°(clockwise and anti-clockwise), while areas further away than 2 km, they rotate by up to ±75°compared to the trajectories of σ 1 obtained with loading condition 1, as they align more closely to the loading direction. When considering an order of magnitude smaller dyke swarm under the same loading condition, the results are essentially the same (these results can be found in the supplementary material). Fig. 11 shows the results of a model geometry that considers together the volcanic complex, a cooled or stiff dyke swarm, and a new radial swarm striking N60°E (Model-type 3). Now, the dyke that was used in Model-type 2, is modelled as a stiff inclusion, i.e. a unit with a Young's modulus five times that of the host-rock (50 GPa). A second elliptical cavity striking N60°E is added to the east of the previous dyke, simulating a later dyke injection.

Model-type 2: elliptical cavity

Model-type 3: elliptical cavity and stiff inclusion

When applying an over pressure of 10 MPa (loading condition (1), Fig. 11 a, c) to both the circular and second elliptical cavity boundaries, the stiff dyke has the effect of increasing the magnitude of the tensile stress at the closest dyke tip, i.e. the westernmost area of the pressurized dyke. The bedrock in this area also undergoes dilatation. On the other hand, the magnitude of σ 3 in the bedrock to the east of the easternmost dyke tip is less, most likely due to its broader distance from the circular cavity compared to the westernmost dyke tip. This broader distance also inhibits linkage between the two over-pressured cavities. σ 1 trajectories in the region between the inflating cavities strike ENE, and in the region east of the easternmost inflating dyke tip, the trajectories strike NE.

The influence of the stiff dyke becomes more pronounced when applying loading condition (2). In this case the stiff dyke acts partly as a stress barrier to the σ 3 stress generated by the new pressurized dyke (Fig. 11c). The region located in between the inflating cavities and the tips of the inflating dyke experiences dilatational strain, still without interacting much (Fig. 11 c, d). The regions to the NW/SE of the inflating dyke undergo contractional strain. The trajectories of σ 1 in the region between the inflating cavities are ENE to EW striking, rotating slightly (up to 15°) with respect to σ 1 trajectories obtained from loading condition (1). In the region east of the easternmost inflating dyke tip they rotate by 7°to 24°with respect to σ 1 trajectories obtained from loading condition (1), from predominantly NE striking to ENE striking trajectories, as again, they align more closely to the loading direction.

Discussion

In this section, we integrate field observations, paleo-stress analysis and numerical models to consider the emplacement of a dyke swarm of TSPP volcanic complex. Our field observations indicate that the observed section of the dyke swarm was likely emplaced at very shallow depths and contributed to feeding a series of fissure eruptions. Our numerical models indicate that the orientation of principal stresses can rotate in the presence of previous dyke injections, which is due in part to the geometry of the dykes interacting with both local and far-field loading sources. These results are discussed and interpreted following four main topics to facilitate their understanding: (1) the dyke segment geometries; (2) the temporal evolution of the dyke swarm and the role of previous dykes in influencing propagation paths; (3) estimation of reservoir conditions from the dyke geometry; and (4) a conceptual model and regional implications associated with dyke emplacement.

Dyke segment geometry and evidences of a dyke-fed eruption

The 35 dyke segments studied here, form a~1 km long swarm that fed a sub-areal or sub-glacial eruption. This is indicated by a feeder dyke that transitions into a lava feeder near the surface of ã 358 × 175 m eruptive cone. Striations in lavas at the top of the eruptive cone could be related to glacial erosion (Singer et al., 1997).

The measured thickness of the feeder dyke (> 4 m) is similar to that of other typical feeder dyke thicknesses measured, for example by Geshi et al. (2010Geshi et al. ( , 2020 in the caldera walls of Miyakejima volcano. The thickest segments measured in the swarm have maximum lengths <26 m, while the longest segments (L > 50 m) are no more than 3.5 m thick. The thickness: length aspect ratio of individual segments ranges between 0.03 and 0.45 (mean and median of 0.20 and 0.16 respectively), values that are abnormally high for dykes, which generally have thicknesses: length aspect ratios in the range of 10 −2 and 10 −4 (Rubin, 1995). However, when we consider the entire dyke lengths rather than individual segments, the aspect ratio falls much closer to this range, with ratios between 0.016 and 0.090 (mean and median of 0.04). This highlights the importance of considering, where possible, full dyke lengths rather than dyke segments which may lead to erroneous interpretations (e.g. Daniels et al., 2012).

However, segmentation of dykes as observed in the studied swarm is a common and inherent feature that occurs due to lateral variations in crustal stresses, crustal discontinuities, variable in-situ tensile strength and/or rock layering, among others (e.g., Delaney and Pollard, 1981;Baer and Beyth, 1990;Gudmundsson, 2020). Therefore, it should be assumed that the dyke segments are connected to a deeper parent dyke and are hence continuous at some greater depth below the surface. We hypothesize that this dyke swarm reached near-surface crustal levels and generated an eruption, therefore segmentations could have occurred through changes in the crustal properties at shallow depths. The changes in strike observed in the exposed dyke segments result from a rotation of the stress field around a vertical axis (usually assumed to be parallel to the propagation direction), yet the expected en-échelon pattern is not clearly observed.

Based on Hoek's classification of dyke-fracture geometry (1991), the exposed dykes can be described as a combination of both irregular and braided patterns. The dykes present little or no overlap between segments and occur in offsets, both characteristic of irregular segmented fracture systems. The dykes in the swarm are also characterized by the side-by-side occurrence of segments that defines braided dykes (Hoek, 1991).

Overpressures and depth of magma chamber of the feeder dyke swarm

Aspect ratios (thickness:length) of the dykes were used as input to analytically estimate the magmatic overpressure p o and depth to the dyke source (h) of this feeder dyke system. We thus undertake the methodology adopted by e.g. Gudmundsson (2011); Becerril et al. (2013), as a first order insight on the depth source of our dyke systems. The aspect ratio of a dyke (w/L) is related to its internal fluid overpressure (p 0 ) and elastic parameters (E, ν) by Eq. (1). Moreover, p 0 is also related to the magma's buoyancy (ρ r − ρ m ), (where ρ r and ρ m are the host rock and magma densities), the fluid excess pressure in the source (p e ) (in excess of the lithostatic stress), the differential stress state at the surface (σ d ) and the depth to source (h) by Eq. (2).

From the first equation, we can calculate p o , and by combining Eqs. (1) and (2) we can estimate the depth to the source as follows: Table 2 shows the resulting overpressure and depth of origin of magmas for 9 dykes, obtained by assuming constant standard values for p e and σ d (2.5 and 1 MPa respectively). We are aware that there could exist uncertainty in these values that could thus affect our estimations. Nevertheless, using the aspect ratios of "total dykes" (0.016-0.09), we estimate magmatic overpressures to range between 8.7 and 28.7 MPa. The calculated overpressures are similar to those obtained from dyke aspect ratios in other regions and numerical models from Laguna del Maule (e.g. Geshi et al., 2010;Gudmundsson, 2011;Becerril et al., 2013, Novoa et al., 2019Masoud, 2020;Geshi et al., 2020), and could be associated to mid crustal reservoirs. The estimated depths to the source reservoir of the dykes range between 2.7 and 22.8 km.

The range of source depths calculated here could be related to (i) shallow intrusion(s) that likely act as a heat source to MGS and/or (ii) a mid-crustal magma reservoir beneath the volcanic complex (as obtained by Reyes-Wagner et al. (2017) through 2D magnetotelluric inversion). Within this range of depths, small magnitude (< 3 M L ) volcano tectonic events near the volcanic edifice have been reported (OVDAS, Observatorio Volcanológico de los Andes del Sur). The values of source depths calculated, represent the long-term location of the plumbing system of the volcanic complex, and should be complemented with other geophysical (seismic tomography, active seismicity, gravity) and geochemical methods (mineral chemistry, geobarometry) to reduce uncertainty.

Several studies have already attempted to estimate the overpressure and source depths of magmas using the aforementioned technique in other volcano tectonic settings (e.g. Ray et al., 2006;Becerril et al., 2013;Browning et al., 2015;Elshaafi and Gudmundsson, 2016;Karaoğlu et al., 2018;Drymoni et al., 2020). Nonetheless, analytical calculation of magma overpressure and magma source depth is also highly sensitive to dyke aspect ratios and the appropriate choice of crustal stiffness (Young's modulus). Further studies are needed to confirm general relations between magmatic driving pressures and aspect ratios of dykes and other fluid filled fractures in Andean settings.

Dynamic paleo stress reconstruction and the role of dyke injections in local stress field rotations

The paleo stress reconstruction from dyke orientation data indicates an extensional stress regime at the time of dyke formation, with a horizontal NW-trending σ 3 and steeply dipping σ 1 (Fig. 8a). This paleostress state partly explains the overall ENE orientation of the dyke swarm and we suggest that this result describes a "bulk" tectonic stress field, since it does not necessarily explain all the dyke orientations observed in the dyke swarm. Dyke strike ranges almost 90°, with a high concentration of dykes striking between N30°E and N90°E (n = 27), with two preferential orientations observed that appear to cross-cut each other: (1) N50°-70°E (n = 10) and (2) N70°-90°E (n = 12) (Figs. 6, 7c), which we expected to represent two stress states. The stress inversion method used here is highly sensitive to outlier orientations (e.g. dyke 1.1), and therefore the resulting stress field could be subject to debate. The stress ratio obtained from the inversion and the 95% confidence ellipse of the least and intermediate principal stress axes (Fig. 8a) (as seen in the stereoplot) indicate that the least and intermediate principal stresses could have rotated and swapped.

Dykes are dynamic features that influence their surrounding stress field by locally compressing the crust in orientations perpendicular to their strike to facilitate opening (e.g. Roman and Cashman, 2006;Tibaldi and Bonali, 2017). Hence, the emplacement of any one dyke influences the stress field and the orientation of later dyke emplacement events. In Figs. 10 and 11 (a, c) of our numerical results, we note that tensile stress is generated at the tips of the pressurized dykes of a sufficient magnitude to promote hydrofracture, dilation and the continuation of the dyke propagation path. At the tips of the dykes, the trajectories of σ 1 , have a~60°range (±30°) in the absence of a regional tectonic load, but that range is reduced to nearly~40°(±20°) when we add a horizontal compression. These findings suggest that dykes do not necessarily emplace following a "bulk" stress field, but instead form over a wider range of orientations that takes account of both the bulk stress field and local stress anisotropy, as for example, that generated by a previous dyke injection. Nonetheless, the effect of a far field stress regime remains the most important control on regional dyke swarms, as we observe from our numerical models that σ 1 rotates in the direction of the regional stress field and therefore dyke propagation paths seek to maintain their orientation parallel to a plane of principal stress.

From our field observations and the numerical models, it is possible to consider the evolution of the studied dyke swarm as follows:

Injection of a basaltic magma through ENE-striking sub-vertical dykes from either shallow or mid-crustal magmatic reservoirs beneath the TSPP volcanic complex. Injection of these dykes alters the stress field, rotating the σ 1 trajectories to a~60°range of orientations between N50°E and N110°E. Preferential injection of NE-striking vertical dykes again rotates the trajectories of σ 1 and influences the propagation path of later dykes. Injection of these dykes also, in turn, alters the stress field, and σ 1 trajectories rotate to a~60°range between N25°E and N95°E (Fig. 12). The effect of regional tectonic compression is to reduce the range of angles over which any new dyke segment can become emplaced (Fig. 12b). This result is supported by observations from laboratory true-triaxial deformation experiments which indicate that the range of angles over which new mode I fractures grow is reduced (a) σ 3 and σ 1 respectively colour contour and trajectories. Both magnitude and trajectories decay radially around the circular cavity, while around the dyke, stress concentrates at the tips. (b) Volumetric strain. The bedrock to the WSW and ENE of the dyke undergoes dilatational strain. (c) σ 3 colour contour and σ 1 trajectories. The magnitude of σ 3 increases in the bedrock compared to loading condition (1); σ 1 trajectories align to the loading direction (d) Volumetric strain. A region undergoing dilatational strain links both cavities, as well as to the ENE of the dyke. (2). The stiff dyke has Young's modulus equal to 50 GPa. (a) σ 3 and σ 1 respectively colour contour and trajectories. Tensile stress concentrates at the pressurized dyke tips and is intensified at the westernmost dyke tip. (b) Volumetric strain. The bedrock to the WSW and ENE of the pressurized dyke undergoes dilatational strain. (c) Contour plot of the σ 3 and σ 1 trajectories. The magnitude of σ 3 increases in the bedrock compared to case (a). The stiff chilled dyke acts as a stress raiser (concentrator), where σ 3 reaches~4 MPa. (d) Volumetric strain. The chilled dyke undergoes constriction while the surroundings undergo greater dilatational strain than in case d) (by~0.00002).

Table 2

Magmatic overpressure and depth of origin of magmas for 9 total length dykes. Constant values used for these calculations are the following: Young's modulus (E) and Poisson's ratio (ν) of the host rock is 1 GPa and 0.25 respectively; the average density of the host rock (ρ r ) and magma (ρ m ) is 2650 kg/m 3 and 2450 kg/m 3 respectively; internal excess magmatic pressure in the reservoir (p e ) is 2,5 MPa; differential stress (σ d ) is 1 MPa; and acceleration due to gravity (g) is 9,81 m/s 2 . (*) Thickness is taken as the arithmetic mean of dyke segment thickness; (**) Total length is taken as the sum of dyke segment lengths.

Dyke

Thickness ( when the level of applied horizontal stress is greater (Browning et al., 2017), in that case the horizontal stress is the intermediate principal stress.

Emplacement model and regional geological implications

The ENE oriented Tatara-San Pedro-Pellado volcanic chain is characterized by abundant ENE-striking dykes and minor eruptive vents, such as the feeder dyke swarm described here. In the Southern Andes, this orientation of volcanism can also be observed at the NE oriented Callaqui stratovolcano and dyke swarms, which is formed in an overall dextral transpressional setting (Sielfeld et al., 2017) and at the Llaima-Sierra Nevada volcanic chain (Naranjo and Moreno, 2005). These examples suggest that Pleistocene-Holocene magma migration has occurred in ENE oriented pathways, driven by a combination of both Quaternary intra arc ENE trending shortening (Lavenu and Cembrano, 1999;Cembrano and Lara, 2009) and optimally oriented fault-fracture systems.

The paleo stress reconstruction from dyke orientation agrees with the fault-slip and dyke orientation paleo-stress inversion in the same area (Sielfeld et al., 2019b) and aids in the definition of the Tatara-Damage Zone, a ca. 4 km wide and 9 km long graben-like transtensional block over which the TSPP is constructed. Faulting activity within this structure is long-lived (late Miocene up to late Pleistocene) and could have helped develop the construction of the volcanic complex by providing active pathways for interseismic eruptive activity, as well as the heat source and permeability structures for the emplacement of a geothermal reservoir (MGS) Sielfeld et al., 2019b).

Resistivity models from magnetotelluric data inversion show a wide conductive region from the volcanic front at this latitude to the east (Reyes-Wagner et al., 2017). Three conductive bodies, related to three distinct magmatic-hydrothermal bodies, are associated with TSPP, MGS and LDM. Reyes-Wagner et al. (2017) correlated the deepest and largest anomaly (~8-10 km deep) to a magma reservoir at the volcanic front (TSPP), and associated the shallow body (~3 km deep) with the heat source for the high enthalpy MGS. The source depths obtained from the analytical method used in this study agree with Reyes-Wagner et al. (2017) models: shallow depths (< 5 km) obtained from our calculations could be correlated to a magmatic heat source for MGS, while depths >7 km could be associated with deeper crustal reservoirs. The latter could also be supported with the dykes more primitive composition (Fig. 13).

It is reasonably assumed that most dykes possess a strike dimension significantly smaller than their dip dimension, indicating vertical propagation. The geometry of the dyke swarm, made up of several offset segments, is also indication that the dykes likely propagated predominantly vertically (Gudmundsson, 1990;Ray et al., 2006), although we do not rule out the possibility of some lateral propagation. However, whether the dykes propagated vertically or laterally bears no impact on the described mechanical controls on the strike orientations. The exposed feeder dyke(s) of this swarm are the near-surface representation of likely only a mere 10-20% of the dykes that may have initially propagated from a deeper magmatic source (e.g. Menand et al., 2010). The remaining majority of 80-90% of the rest of the dykes would have be arrested or deflected into sills at depth (Gudmundsson, 2020). We hypothesize that the dykes in the TSPP propagated vertically from a deep magmatic reservoir and that the arrest and deflection of dykes into sills at~2-3 km has served to mature the heat source and generation of the active (MGS) hydrothermal system (Fig. 13).

Conclusions

1. We studied a near surface section of a feeder dyke system,~5 km to the east of the Pellado volcano (Southern Volcanic Zone of Chile). This area lies directly above a conductive anomaly which has previously been interpreted as an active geothermal system (Mariposa Geothermal System; MGS). In total, 35 dyke segments were measured, some of which likely fed a fissure eruption, as an intrusion into a minor eruptive centre. 2. Analytical estimates of the depth of magma origin from aspect ratios of 33 dyke segments indicate a range of depths between 2.7 and 22.8 km for the origin of the magma. This range of depths represents the long-term location of the plumbing system of the volcanic complex, and could be related to shallow crustal reservoirs or to deeper reservoirs associated with a~8-10 km deep conductive anomaly previously observed with 2D magnetotelluric inversion (Reyes-Wagner et al., 2017) in the volcanic arc region. From these observations, we hypothesize that magma propagated vertically from a magmatic reservoir at around 8-10 km depth directly to the surface to feed eruptions, as was the case with the dyke swarm studied. However, it is possible that magma from the deep reservoir has been transported to shallower crustal levels where it is arrested and acts as a heat source for the generation of the previously identified MGS. 3. Our two-dimensional finite element models show that a combination of both internal magma pressure and regional tectonic compressive loading locally rotate crustal principal stresses. Rotation is exacerbated when the crust hosts heterogeneities, such as previous dyke injections. The effect is to locally rotate the principal stress orientations and locally increase the magnitude of the bedrock's tensile stress. These effects combine to influence the trajectories over which new dyke injections will propagate and hence the arrangement of shallow magmatic plumbing systems, as well as the distribution of heat sources for hydrothermal systems. 4. We find that individual dyke injections themselves can locally rotate the orientation of the principal stresses if the regional tectonic loading direction is somewhat offset from the dyke strike. The effect again, is that previous dyke injection events modify the propagation paths of future dyke injections. This finding should be considered during periods of unrest, and for volcanic hazard assessments, when considering the location and orientation of secondary fissure eruptions in active volcanic areas such as the Southern Volcanic Zone of Chile.

CRediT authorship contribution statement