Introduction
Geothermal energy exploits the Earth's natural heat for power generation. It is a low-carbon, base-load alternative energy that, despite its low cost relative to other energy sources, is not yet widely used (Goldstein et al., 2011): financial risk associated with locating and developing the resource - which means understanding the underground shape of permeable features such as faults, for their exploitation - have represented so far, together with environmental impacts, main barriers (Subir and Morrow, 2011; IFC, 2013). Geophysical techniques are often useful for discovering unknown subsurface conditions and can be utilized as preliminary screening before performing direct investigations. Non-invasive potential methods - such as gravity and magnetics - are often appropriate for locating structures (e.g. Glen et al., 2008) and in the general characterization of an underground (geothermal) system, aiming to its exploitation (Pipan et al., 2010; Kipsang, 2015). During the last five decades, a variety of methods based on the use of vertical and horizontal gradients of potential-field anomalies have been developed and implemented for the determination of geologic boundaries, such as contacts and faults - hence also valuable for the exploration of geothermal resources (Nabighian, 1972, 1974; Keating and Pilkington, 1990; Ferreira et al., 2013).
In the trans-tensional sector of the Great Basin, Basin and Range (B&R), western USA, the importance of characterizing known geothermal resources is heightened by the fact that most of the still-undiscovered geothermal systems (~75%) may show no surface hydrothermal springs, thus defined as blind or hidden (Coolbaugh et al., 2006; Faulds and Hinz, 2015). Faulds et al. (2010; 2012) categorized almost 250 geothermal fields in the Great Basin and found that the most favorable tectonic setting for the existence of a system is associated with the presence of a (major) normal fault: geothermal fluids flowing preferentially at fault tips and in fault interaction zones. These zones are characterized by a high density of fractures and dynamically-maintained permeability (Curewitz and Karson, 1997; Rowland and Sibson, 2004; Cashman et al., 2012).
Through the study of aerial data collected along the Surprise Valley, CA (North-West Basin & Range), Glen et al. (2013) identified an intra-basin magnetic high, running a significant length of the valley, which they interpreted as a buried, faulted mafic dyke. All the geothermal manifestations in the valley seem to be related to this body as they lie where the high is cut by perpendicular faults (Glen et al., 2013).
In this paper we examine an area within the central Trans-Mexican Volcanic Belt (TMVB, Figure 1), in the northeastern part of the Michoacán-Guanajuato Volcanic Field, a sector interested by, approximately, N-S, E-W, NW-SE and NE-SW structures, repeatedly involved in the formation of small, monogenetic and larger, polygenetic volcanoes (Hasenaka and Carmichael, 1985; Suter et al., 1999; Cebriá et al., 2010). We use magnetic data recorded through detailed surface surveying for characterizing a localized positive magnetic anomaly corresponding to a stretch of the central Cuitzeo Lake (Michoacán, Mexico). Southernly, on the lake shore, hydrothermal manifestations reveal the presence of a geothermal system of medium enthalpy at depth (Arredondo-Fragoso, 1983; Campos-Enríquez et al., 1988).
In the following, we first give a broad review of the regional geology, describing the different tectonics interesting the area. We then apply filters based on combinations of directional derivatives to magnetic data surveyed in the study area, for structural interpretation and with the aim of identifying the heat source and the local structures that allow the rise of the geothermal fluid. We eventually try to give a volcanological explanation for the presence of such heat source at depth.
Geologic setting
The study area is located within the central segment of the TMVB, between 99o and 102o W longitude, Figures 1 and 2. Regionally, the area is interested by the northward extension of the arc, carried by ~E-W striking normal faults with a minor but consistent left-lateral slip component (Suter et al., 1995a; Ego and Ansan, 2002; Garduño-Monroy et al., 2009). From Oligocene times, this territory has been affected by two main trends of subduction-arc related volcanism and the related sin- and post-volcanic extensions, nearly perpendicular to each other. These are the NNW-trending Sierra Madre Occidental volcanic province and the roughly East-West-trending Trans-Mexican Volcanic Belt (TMVB). The two volcanic belts overlap between the Pacific coast and the longitude of Mexico City (Figure 1a) and have at least two characteristic features in common: the broad orientation of the arc and the dominant composition of the rocks (Ferrari et al., 1999).
1 Sierra Madre Occidental volcanic and Basin and Range extensional provinces
The Sierra Madre Occidental (SMO) of western Mexico is a huge Silicic plateau which runs a NNW-trending direction for over 2000 km, from south of the TMVB to the MX/USA boarder (Ferrari et al., 1999; Ferrari et al., 2000). It is the result of Cretaceous-Cenozoic magmatic and tectonic episodes related to the subduction of the Farallon plate beneath the North American Plate (NAP). The rapid increase of the subduction angle during the removal of the plate generated a flux of uprising hotter astenospheric magma, underplating the NAP and intruding within the lower crust (Ferrari et al., 2002), allowing for the extrusion of >300,000 km3 of silicic material (Huppert and Sparks, 1988; Ferrari et al., 2007), apparently erupted from fissure vents corresponding to Basin & Range faults traces (Henry and Aranda-Gomez, 1992; Aguirre-Díaz and Labarthe-Hernández, 2003).
The extensional phases that followed the main volcanic peaks and aged about 30, 23, 10 and 5 Ma (Ferrari et al., 1999) have been sufficiently intense (up to 100%) in the northern part of the SMO (northern Mexico, Figure 1a) to exhume part of the Proterozoic crystalline basement (Dickinson and Lawton, 2001; Ferrari et al., 2007), whereas in the reminder of the arc, and particularly within the southern B&R province, extension has not exceeded 20% (Henry and Aranda-Gomez, 1992; Aranda-Gomez and McDowell, 1998), achieving in the Mesa central, north of TMVB (Fig. 1a), a maximum of about 8% (Nieto-Samaniego et al., 1999).
While an ignimbrite from the upper part of a 300 m thick rhyolitic succession 15 km south of the city of Morelia has been dated by K-Ar at 21±1 Ma (Pasquaré et al., 1991), B&R faults allowing for the extrusion of SMO lavas and aged between 38 and 25 Ma are widespread south of the TMVB (Alba-Aldarve et al., 1996; Morán-Zenteno et al., 1999). Garduño-Monroy and Gutiérrez-Negrín (1992) noted how the regional NNW-SSE B&R Taxco-Queretaro fault, passing some tens of kilometers east of Morelia, represents the eastern limit of SMO volcanism and of the Guerrero Terrain (Campa-Uranga and Coney, 1983).
Seismicity, high heat flow and recent volcanism indicate that the Basin and Range province is actively extending, within its northern and southern parts (Parson, 1995).
2 Trans-Mexican Volcanic Belt and subduction of Fracture Zones
The TMVB is an east-west, Late-Miocene to Quaternary mostly calc-alkaline continental volcanic arc of more than 1000 km longitudinal length and up to 150 km N-S width, originated by the subduction of the Rivera and Cocos plates under the North American plate (Nixon, 1982; Garduño-Monroy et al., 1993; Ferrari et al., 2012), Figure 1b. The mean elevation of the arc is >1500 m with respect to the southern forearc region; summit elevations are between 3000 and 5700 m (Suter, 1991). TMVB presents high mean heat flow (>80 Wm-2; Prol, 1991) and the combination of this with features such as a negative regional Bouguer anomaly (Urrutia-Fucugauchi and Flores-Ruiz, 1996), volcanic activity, sin- and post-volcanic shallow extensional faulting and seismicity, suggests that its central E-W Chapala-Tula fault system may correspond to a young active continental rift (Luhr et al., 1985; Marquez et al., 1999). To the point that Verma et al. (2009) coined the terminology 'Mexican Volcanic Rift'.
Along the arc, volcanoes are distributed with a ~15o oblique trend relative to the Middle American Trench (MAT, Figure 1): this suggests that their location is controlled by the slab geometry (Pardo and Suarez, 1995; Ego and Ansan, 2002; Ferrari et al., 2012). Particularly, the area of Mexico inland of the MAT can be split into several sections with different subduction angles (from W: Jalisco block, JB, 50o; Michoacán block, MB, 30o; Guerrero block, GB, 0o until the latitude of Mexico City and Oaxaca block, 30o. Pardo and Suarez, 1995; Stubailo et al., 2012. Fig. 1b).
The Rivera Fracture Zone (RFZ) is generally recognized as the Rivera-Cocos Plates Boundary (RCPB, DeMets et al., 1990; Bandy et al., 1995), it can be dated at 10 Ma as the Rivera plate (Bandy, 1992; Stubailo et al., 2012) and is interested by higher heat fluxes than the surrounding oceanic crust. Gravity, heat-flow and sea-morphology analyses suggest that the subducted part of the RCPB lies directly beneath and is oriented parallel to the Southern Colima graben (Bandy, 1992). Bandy et al. (1995) showed that the RCPB lies east of the Central and Northern Colima graben and is related to the Colima Gap (CG, Figure 1b), coinciding in plan with the area where the Wadati-Benioff zone undergoes a sharp decrease in dip, eastward. These authors suggested the presence of a NE-SW oriented tear-fault between the subducting plates, explaining the low density zone beneath the granitic highlands (CG) by thermal convective movements in the upper mantle, providing a conceptual model consistent with roughly E-W extension within the subducting lithosphere (Nixon, 1982; Bandy, 1992; Ferrari et al., 1994).
Easterly of the RFZ, Stubailo et al. (2012) confirmed previous results (e.g., Pardo and Suarez, 1995; Blatter et al., 2007) that the angle of subduction varies substantially at the longitude of the Tzitzio Gap (TG, figure 1b), consistent with the location of the Orozco Fracture Zone (OFZ, Fig. 1), which separates older (17.6 Ma), cooler and denser oceanic crust at the NW side from younger crust (12.3 Ma) at the SE (Manea et al., 2005). Slab rollback (retreat) below the MB occurs at about 5 cm yr-1 (Singh and Pardo, 1993; Bandy et al., 2000; Suter et al., 2001): this process should displace mantle asthenosphere, and, in the presence of a slab tear, the mantle material would flow through it. Tear faults occur where the plate is weakened by the presence of the fracture zone.
During subduction, fracture zones are buoyant enough to usually lower the angle of descent of the subducting plate: the arc volcanism is consequently displaced landward, creating indentations in the front of the volcanic arc (McCann and Habermann, 1989). Through vulcanological studies, Blatter and Hammersley (2010) concluded that the OFZ is being subducted under the TG and that the buoyancy of the fracture zone would cause this section of the Cocos slab to subduct at a lower angle than the Cocos slab on either side of the OFZ. This causes the subducting material to reach the melting depth of about 100 km at higher distance from the trench, in correspondence with the OFZ (and TG), spatially 'delaying' volcanism at the surface towards the north. Dougherty et al. (2012) proposed that the Cocos slab is currently fragmenting into a N Cocos and a S Cocos plates along the projection of the OFZ, in agreement with other authors (Bandy, 1992; Bandy et al., 2000). They identified this tearing event as a younger (~0.9 Ma) analogy of the 10 Ma old Rivera-Cocos plate boundary.
3 Structural geology, regional seismicity and state of stress of central tmvb
The central part of the TMVB (99o-102oW) is characterized by major E-W intra-arc basins: from W to E, the Chapala, Cuitzeo, Acambay and Mezquital grabens, 8-10 Ma old (Garduño-Monroy et al., 1993; Rosas-Elguera et al., 1997), whose E-W to ENE-WSW bordering normal faults have mean Quaternary slip rate of <0.1 mm yr-1, characterizing the area with a ~N-S oriented extension of <5% (Suter et al., 1995a; Ferrari and Rosas-Elguera, 2000). Recent studies have proposed, for some of these structures, an age of 18 Ma (Mendoza-Ponce et al., 2018). Figure 2 shows the Acambay and Cuitzeo extensional areas. The Acambay graben is a ~30 km long E-W structure with delimiting faults dipping 50o-70o, with a mean slip rate of 0.17 mm yr-1 during Holocene times (Langridge et al., 2000). Cuitzeo graben (and half-graben) can be traced for over 45 km and its master faults dip between 45o and 75o (Ferrari et al., 1991; Garduño-Monroy et al., 2009. See Figure 3), having slip rates between 0.09 and 0.18 mm yr-1, which is approximately three times higher than in the Morelia area, 20 km south (Suter et al., 2001). E-W fault planes' striations indicate an extensional (N-S) dip-slip with a left-lateral strike-slip component (Pasquaré et al., 1988; Garcia-Palomo et al., 2000), which was explained by Ego and Ansen (2002) with slip-partitioning at the trench.
Central TMVB is seismically active (Suter et al., 1996; Garduño-Monroy et al., 2009) and during the last century, eight shallow events with magnitude 4.1<M<6.9 occurred along E-W faults; among these, the Acambay, 1912 (M=6.7) and the Maravatio, 1979 (M=5.3) earthquakes (Suter et al., 2001).
Southern B&R NNW-SSE structures present a normal-right-lateral sense of motion (Garduño-Monroy and Gutierrez-Negrin, 1992; Garcia-Palomo et al., 2000). Historical seismic events of Mw>3 comprehend the Pinal de Amoles (1887) and the Sanfandila (1998) events (Suter et al., 1996; Zuñiga et al., 2003) and Peñamiller sequence (2010-2011; Clemente-Chavez et al., 2013). Although we do not have data on mean slip rates, Alaniz-Alvarez et al. (1998) determined how, on average, north-south B&R faults only have about one-third of the horizontal displacement rate of the east-west TMVB faults. The contemporaneous activity of TMVB and B&R faults in the area was explained with a recurrent permutation of σ1 and σ2 (see figure 2, where σ1=σV) which would allow for the activation of NNW-striking B&R structure as right-lateral during the deposition of the Cuitzeo (and Acambay) basin sedimentation (Suter et al., 1995b; Ego and Ansan, 2002).
3.1 The Michoacán-Guanajuafo Volcanic Field
A monogenetic volcanic field, covering an area of some 40,000 km2 is located at the west-central sector of TMVB, extending over northern Michoacán and southern Guanajuato and interesting the study area (Figures 1 and 4). A unique part of TMVB for its lack of young large composite volcanoes, dominant in the rest of the arc, the Michoacán-Guanajuato Volcanic Field (MGVF) has been active from Pliocene to present (Hasenaka and Carmichael, 1985; Ferrari et al., 2007). It contains more than 1000 small volcanic centers, as shown in Figure 4, and its main eruptive products are calc-alkaline olivine basalt and basaltic andesites (Hasenaka and Carmichael, 1987): rocks crystallized from these magmas can hold relatively high abundances of ferromagnetic minerals. The majority of the small-sized volcanoes are cinder cones (~900) but the field also comprehends lava domes, maars, tuff rings and thick lava flows not associated with cones (Williams, 1950; Hasenaka, 1994), coexisting in time and space with over 300 medium-sized shield volcanoes distributed throughout the volcanic field and whose origin may not be univocal (monogenetic or polygenetic. Hasenaka et al., 1994; Verma and Hasenaka, 2004). While small volcanic centres in MGVF are spread all over the extension of the field, authors have identified alignments of cones following geological structures (Cebriá et al., 2010; Gomez-Vasconcelos et al., 2015), implying that, although volcanism may be a consequence of contemporaneous extensional regimes, the distribution of monogenetic vents is actually controlled by reactivation of older fractures, producing space for magma ascent at near surface levels (Hasenaka and Carmichael, 1985).
Through geochemical and isotopic studies of rock samples from the whole extension of the MGVF, Verma and Hasenaka (2004) concluded that no slab input is present in the basic magmas of the field and that different degree of partial melting of a heterogeneous mantle source could explain the origin of most magmas in the field, allowing for some crustal assimilation for the most evolved - supported by the presence of frequent granitic xenoliths (Verma and Hasenaka 2004, and references therein).
4 Regional geophysics and the study area
Figure 5 presents the relative Bouguer
anomaly with respect to the WGS84 across the area within and around the Cuitzeo
Basin. This map is based on data of Petroleum of Mexico (Pemex) from 1980 (Arredondo-Fragoso, 1983). The lowest values
of the gravity field within the lake are linked with relatively light lacustrine
sediments (density, ϱ
Figure 6 shows the reduction to the pole of the regional magnetic field around the study area, from aeromagnetic data of the Mexican Geological Survey (Servicio Geológico Mexicano, 1988). The magnitude of an anomaly in the magnetic field recorded at the surface, with respect to the International Geomagnetic Reference Field (IGRF; Thebault et al., 2015) mean value, is dependent on the concentration of ferromagnetic minerals in the rocks of the upper crust at the site (Jackson and Bowies, 2014). Strong positive anomalies may then be related with large volumes of rock of basic composition (andesitic/basaltic solidified magmatic bodies), below Curie temperature for magnetite (around 575oC) and crystallized after the Matuyama-Brunhes reversal, about 780 ka (or during any other period of normal polarity), outcropping at the surface or trapped within the shallow crust. Cuitzeo basin is bordered by different magnetic highs linked with volcanic edifices at the surface (Figure 6), particularly the SSA, NE of the lake.
A relatively strong positive magnetic anomaly lies in the middle of the actual lake, some hundreds of meters NW of the gravity low, in plan. The origin of this anomaly may be found in a volcanic structure older than the lake and consequently covered by the lake sediments, or it may be correlated to an ascending magma body, blocked by the plastic lacustrine sediments.
In Figure 7a principal lithologies and main structures are presented: massive andesitic and ignimbritic highly fractured banks outcrops in the area (Fig. 7 b and c). Our structural survey showed a statistical prevalence of ~E-W fault planes on outcrop, with dips <75o. Nearly vertical ~N-S striking fault planes (showing horizontal striae and with fractures completely sealed by hydrothermal calcite) outcrop in an area south of San Juan Tararameo (SJT, Figure 7d), denoting the existence of a paleo-geothermal system. In the image, inferred continuations of some of these structures below the lake are shown with dotted lines and an andesitic island can be seen in the center of the image. Shown are also the areas interested by geothermal manifestations, with more than 100 springs between San Agustin del Mais (SAM) and SJT, laying near faults intersection (Sibson, 1981; Rawland and Sibson, 2004) and marking the areas with long-lived siliceous sinter terraced-deposits, associated with high temperature hydrothermal reservoirs (Fournier and Rowe, 1966).
Water temperatures of the springs reach 93oC (our work), water-boiling T at the location (1820 m a.s.l.). A geochemical study of different geothermal springs around lake Cuitzeo (Segovia et al., 2005) found at SAM and SJT steam heated waters in partial equilibrium, observing a mixing trend among the samples - indicating reservoir temperatures between 150 and 220oC. Their Radon analysis results indicated a highly efficient fluid flow transport in the zones where higher temperatures were estimated. Lastly, the 8 NE-SW traces over which magnetic measurements were taken, as described in the following, are also depicted in Figure 7.
Magnetic survey: data acquisition, processing and interpretation
Survey lines, spaced 1 km, were 15 km long but the northern one (reduced to 8 km, Fig. 7a) and were covered on foot and by rowboat (made of plastic and wood for not affecting measurements), from May 2015 to July 2016. Measurements of magnetic field where taken at stations spaced 100 m on each line, using a Geometrix G-857 proton-precision magnetometer, connected with a Garmin Oregon 450 GPS that guaranteed each measure to be taken within 2 meters from the station location, on each line. Different stations could not be covered due to the presence of anthropic barriers, corrupting measures in their proximities (e.g., buildings and other constructions, power lines), leaving a total of 975 point-values for analysis, over an area of some 110 km2. Instrumentation also comprehended an independent base-station for calculating daily variations of the magnetic field due to solar wind interferences during working hours, recording values at a fixed point. Due to theft of the base-station, we used the values of the Magnetic Observatory of Teoloyúcan to make diurnal corrections to our data, during the last part of the survey.
We worked on the magnetic raw database (after correcting for daily variations) using the software Encom Discover PA 2010 (®RockWare). We applied to the mapped dataset computational filters based on operations with directional derivatives of the data, as described in the following. Figure 8a maps the anomaly values of the Total Magnetic Field, relative to the IGRF, displaying a positive anomaly in the middle of the surveyed area.
Surface data suffer from high-frequency noise due to shallow sources and probably also due to higher-than-optimum stations separation, for some stations. One of the tasks of data processing is to simplify the complex information provided in the original data: we have applied to our data an Upward Continuation (UC) filter, a low(wave-frequency)-pass filter which allows to avoid considering the contribution to the total magnetic field recorded at the surface given by rock layers above a certain depth (e.g. Gianiyu et al., 2013, Ferreira et al., 2013). Specifically, the filter was applied to the data-grid for the depth of 500 m in order to discard the contribution of high-frequency shallower potential sources to the recorded magnetic field values, as shown in Figure 8b. The filter Reduction-to-Pole (RP) is utilized to image measured values as if the measurements were taken at the magnetic pole: it reduces bipolar anomalies to monopolar ones, considering magnetic Inclination and Declination at the site (respectively 47° and 7° for the study area), basically placing anomalous values of the magnetic field on the vertical of the geological bodies causing them (Baranov, 1957). The application of the RP filter to UC data is shown in Figure 8c: these data have been used as the basis for the application of the mathematical filters described below.
Local magnetic field of the study area, upward continued and reduced-to-pole, presents a principal positive anomaly (reaching values higher than 400 nT with respect to IGRF) sited near the center of the area (Figure 8c) which is potentially due to the underground geologic body we aim to characterize. Also present is a smaller peak, south of the principal one, which might be related with the volcanic stratigraphy of the shore.
1 Data filtering: operations with directional gradients
In the last decades, a variety of methods based on vertical and horizontal derivatives of surveyed potential field data have been developed as efficient tool for the determination of geometric parameters of causative bodies (such as location of boundaries and depth-to-top) and structures (i.e. major faults. Nabighian, 1972; Roest et al, 1992; Miller and Singh, 1994; Ravat et al., 2002). The first vertical derivative (VDr) of the magnetic field is the rate of change of its intensity (M in Table 1) in the vertical direction. Its computation removes long wavelength features of the magnetic field and, while it amplifies signals from shallower sources, it significantly improves the resolution of closely spaced anomalies (Hood, 1965). VDr (Equation 1 in Table 1), on a map, has its zero values over the vertical edges of thick source bodies and positive and negative values over positive and negative anomalies (Cooper and Cowan, 2004). In Figure 9a the first order VDr of the UC and RP data is shown, together with the representation of a likely anomaly' source edge, drawn approximately following the zero contour. In the image, the magnetic response from the southern anomaly seems slightly heightened.
Equation # | Formula | Measure |
---|---|---|
1 |
|
|
2 |
|
|
3 |
|
|
4 |
|
|
5 |
|
|
6 |
|
|
7 |
|
|
The second vertical derivative (2VDr, Equation 2 in Table 1) is used for improving resolution of anomalies and to delineate geological discontinuities in the subsurface. Lineaments in a 2VDr map would lie at value 2VDr = 0, following abrupt change in magnetization due to geologic structures (Rebolledo-Vieyra et al., 2010; Lopez-Loera et al., 2010). In Figure 9b this filter is applied to our data, showing some imaged potential fault-lines following zero contours on the map, where horizontal gradients are highest. Being a second order filter, 2VDr enhances near surface effects at the expenses of deeper anomalies, it amplifies noise and may produce artificial second derivative anomalies (Wahyudi et al, 2017).
Total Horizontal Derivative (THDr) is considered as the simplest approach to estimate contact locations and has been used as edge detector (Cordell & Grauch, 1985, Cooper & Cowan, 2008). It displays maxima above nearly vertical borders of source bodies (or faults with decent offset) and relative minima at the center and outside of sources. Figure 9c shows the application of this filter to the VDr grid (Equation 3), as in Fedi and Florio (2001), where they used the filter for defining boundaries of a calderic collapse in southern Italy. Superimposed to the map are some interpreted structures.
1.1 Analytic Signal
Mathematics, an Analytic Signal (AS) is a complex-valued function that has no negative frequency components. Developed by Nabighian (1972, 1974) and also known as total gradient, the AS amplitude is defined as the square root of the sum of the squared directional derivatives, as described by eq. 4. AS presents maxima above near vertical magnetic contrasts, identifying faults or geologic bodies of different susceptibility or induced magnetization.
Potential field data correspond to the superposition of all causative sources and nearby source interference yields mislocations (Grauch and Cordell, 1987). Authors have used the filter AS on grids of vertical derivatives of n order (for n = 1, 2) for sharpening nearby boarders (Nabighian, 1972, Roest et al., 1992). Figure 9d displays the AS filter applied to the VDr grid (eq. 5). The signal due to the southern positive anomaly appears strongly enhanced and, while this filter distinguishes between the two anomalies, it only marginally helps in structures identification. From the image we could recognize a N-S and a E-W striking structures crossing near the centre of the picture and a NW-SE structure which we could not recognize through the application of other filters but that may also belong to B&R deformation. Application of AS to higher order vertical derivatives of our data amplifies high-frequency noise to intolerable levels, further attenuating deep source signal.
Since the strength of a magnetic field decays proportionally to source-receptor distance cubed, a handicap of these simple derivative-based filters is that they tend to magnify shallow sources at the expenses of deeper ones, whose gradients measured at surface are weaker. Particularly for our data, the southern anomaly recognized on the RP map (fig. 8c) seems enlarged by vertical derivations (fig. 9a, b and d).
2 Tilt Angle
The Tilt Angle (TA), introduced by Miller and Singh (1994) for profiled data and improved by Verduzco et al. (2004) for 3D cases, overcomes the problem of amplifying shallow sources, which characterizes simple derivative-based filters. This is accomplished by dealing with the ratio of the vertical derivative to the horizontal one, as in Eq. 6. Both VDr and THDr will be smaller for deep sources than for shallower ones, so that their ratio will still be large over the source, pass through zero over or near the edges (where VDr is zero and THDr presents maxima), and will be negative outside of the magnetic body (where VDr < 0), in 2D. Due to the nature of the arctan trigonometric function, all amplitudes are restricted to values between π/2 and -π/2 (+90° and -90°) regardless of the amplitudes of VDr and THDr (Miller and Singh, 1994). ffacts make this filter function like an automatic gain control (AGC) sieve, tending to equalize the amplitude of the signal for shallow and deeper sources (Verduzco et al, 2004; Salem et al, 2007). The application of the TA to our data is shown in Figure 10. This filter gives the two sources (northern and southern, in the image) balanced leverages: the dimensions of the southern anomaly appear downsized and can be linked with a shallow source (shore' volcanic stratigraphy?).
Salem et al. (2007) developed the 'Tilt-depth method' for estimating the depth of a magnetic source from 2D TA maps. They demonstrated how, under certain assumptions such as when contacts (between differently magnetized lithologies) are nearly vertical and the magnetic field is vertical (or RP), contours of TA on a map can help identifying both the edges of magnetic structures (TA = 0) and their depth, as half the physical distance between ±π/2 contours. They applied the method on aeromagnetic data over the Karoo sedimentary rift structures of southeast Tanzania (Salem et al., 2007). Applied to our data (Figure 10), the methodology estimates a depth for the roof of the magnetic body (which we argue to equal the depth of the bed of Cuitzeo' lacustrine sedimentation, at the site) between about 700 and 900 m (depth = AB/2, on the TA map).
2.1 Tilt angle of the horizontal gradient
Ferreira et al. (2013) presented an edge detection method that is based on the enhancement of the THDr of magnetic anomalies using the TA. Referred to as Tilt angle of the (total) horizontal gradient, (TH)TA was applied by the authors to 3D synthetic models, displaying balanced maxima over the edges of magnetic prisms located at different depths with outstanding precision, particularly when compared to other edge detection methods. Also, the filter was tested for the detection of edges of superimposed sources (Ferreira et al., 2013), obtaining promising results. The authors applied the filter to aeromagnetic data of an area in the central portion of the Ribeira belt (state of Rio de Janeiro, southeastern Brazil), a Neoproterozoic range consisting of four tectonostratigraphic terranes. Through its maxima, the filter helped in defining major regional faults and shear zones and in characterizing dykes and other intrusions (Ferreira et al., 2011).
We applied (TH)TA to our data (eq., 7. Figure 11) and it revealed very useful in locating potential faults, giving support to our structural interpretation. Main TMVB (E-W) and B&R (N-S) faults in our study area, identified by this method through positive sharp peaks, affect the volcanic bedrock (pre-TMVB) and the >8 Ma old lacustrine sedimentation, thicker than 700 m above the magnetic source in study. There are recurrent E-W structures affecting the study area, one of which interests the Cuitzeo peninsula (Figures 7, 9a, b, c, 11) that, being located on top of the magnetic anomaly in maps, could be directly related to the once hotter magma body, now causing the magnetic anomaly. Particularly, the southern half of the peninsula is made up of lava flows of basic composition (andesite) and is separated from the northern half, constructed with ignimbritic products, by an E-W structure lying approximately in the middle (figure 7a): the two lithologies might have been erupted from the fissure itself, potentially during the same event. On the other hand, we believe that the volcanic island in the middle of the lake may be linked with a N-S structure, as explained in the conceptual model of the system, in the following.
Discussion and Conclusions
The application of filters based on combinations of directional derivatives to magnetic data surveyed over the middle section of Lake Cuitzeo helped in the identification of structures potentially affecting the underground, below and within the cap-rock layer consisting of the basin sedimentation. While the application of the vertical derivatives and tilt angle filters to the data may indicate the presence of a cylindrically shaped magnetic body at a depth shallower than 1 km (Figures 9a and 10), the interpretation of results is based on restricting assumptions (e.g. verticality of lateral borders) which would indeed characterize the magnetic body based on its inherent geometry, related to its nature and origin. The Authors of this work have been criticized on a previous publication (Mazzoldi et al., 2016) in that the magnetic anomaly in the study area may be due to a volcanic structure older than the Cuitzeo lake (8-10 Ma) and submerged by the lake sediments afterward, thus not related with the geothermal system. Transect in Figure 3 shows the stratigraphy below the sediments of lake Cuitzeo, displaying massive amounts of andesitic material constituting the TMVB basement. Authors agree that a TMVB 'basement volcanism', with the same general geometry of the present arc and developed during a time-span comprehending the last activities of SMO and the beginning of proper TMVB volcanism, should be considered as an independent volcanic province (Gutiérrez-Negrín, 1988; Ferrari et al., 1999). With an approximate age between 15 and 7 Ma, its volcanic products would still belong to a calc-alkaline series, but more acidic than TMVB magmas (Venegas et al., 1985; Garduño-Monroy and Gutiérrez-Negrín, 1992), somehow representing a middle member between SMO and TMVB - which contrast with the thesis of the magnetic high being due to a volcanic edifice of ultrabasic composition (needed to justify the anomaly under the lake).
The highly basic and ultrabasic magmas of the MGVF, solidified underground below the lake sediments and forming the studied magnetic anomaly, allowed to shed some light on the structural condition of the geothermal reservoir at Cuitzeo, through the application of directional filters to our raw magnetic data. Our interpretation of results images E-W and N-S to NW-SE strikes (e.g. Figure 11) as principal directions of the faults present in the area. These would be transtensive structures belonging respectively to the TMVB tectonics, extending to the north, and to the older B&R, ENE-directed extensive tectonics. We have also defined dimensions and depth of the magnetic body that we believe could be the heat source of the geothermal system in existence, described below.
1 Geothermal system under Lake Cuitzéo
Thinned crust is often characterized by listric normal faults which tend to flatten at a depth proximal to the brittle-ductile transition (Jackson and McKenzie, 1983; Brogi et al., 2003). The top of the reflective zone is generally located at about 15±5 km (McCarthy and Thompson, 1988; Mayer et al. 1997) but it tends to be shallower beneath regions with higher heat flow (Ranalli, 1995). An example can be found in southern Tuscany, central Italy, in the geothermal areas of Lardarello and Amiata (Batini et al., 1978; 1983; Borgia et al., 2014), where the regional average heat-flow is very high (120 mW/m2; Brogi et al., 2003) and the brittle-ductile transition reach depths as shallow as 3-6 km (Batini et al., 1985; Liotta and Ranalli, 1999), together with its shear-decollement zone properties (Borgia et al., 2014; Mazzoldi et al., 2015). At this depth main listric normal faults flatten.
Along the length of the E-W Chapala-Tula fault system - a proto-continental rift for some authors (see above) - the general northward extension, the frequency of geothermal areas, the high average heat flow (>80 mW m-2) and the average low dip of the domino fault-systems (lower than 70o between Morelia and lake Cuitzeo), make us suppose that the main E-W structures, coeval to the TMVB (8-10 Ma, and potentially already active during TMVB basement volcanism, >15 Ma, Mendoza-Ponce et al., 2018), find a decollement shear-zone at a depth of no more than 8-10 km (probably shallower) - at which level they flatten. On the other hand, being active since the construction of the TMVB, NW- to N-trending Basin&Range faults, older than 20 Ma and with a higher angle of dip, could reach deeper than 15 km: the imposition of the regional average high heat-flow in the area, some 10 Ma (or 15-18 Ma, depending on whether the TMVB basement volcanism is considered as a main player for heat-flow and E-W faults activation), only relatively affected their geometry.
From these observations we derive our essential conceptual model of the geothermal system. In Figure 12a an idealized section cutting along the longitude of the TMVB is presented. The image summarizes the most recent models for the subduction of the Rivera and Cocos plates, previously described, with the two main gaps in Quaternary volcanism at the surface (Colima Gap and Tzitzio Gap) coinciding with the presence of tear faults at depth. This geometry involves a relatively strong vertical heat-flow for those sectors interested by the tears. Figure 12b depicts our conceptual model of the geothermal system. A soaring basaltic/andesitic magmatic body, exploiting a B&R structure of the T-VdS fault system for its rising, got trapped at near-surface level under a thick layer of plastic clay, identifiable with Cuitzeo lacustrine sedimentation, which interrupted its vertical motion, sometime during the last 500 ka. Hydrothermal calcite interesting B&R fault planes on the southern shore of the lake (figure 7d), and thick sinter deposits near hydrothermal springs, suggest that the system might have been bigger and more vigorous in the past. Activity of the N-S faults would have made it possible for a portion of the magma to reach the surface, through the clayey caprock, forming the andesitic island in the middle of the lake. Also, as shown by geological and geophysical analyses (e.g., figures 7a, 9b, 11), the whole peninsula of Cuitzeo may be related to this same intrusion through a E-W structure, more active than B&R, interesting the main part of its extension. Along the depth of the crystallized, still hot intrusion, water would heat up and consequently ascend to the surface, through listric E-W fault planes, more active and highly fractured (permeable) along their intersections with almost vertical B&R planes. These permeable intersections control the location of geothermal springs at the surface (e.g. Sibson, 1990; Figure 12b).
At San Agustin del Mais, there is already some experience in geothermal development, thanks to different thermal resorts in the area, and geothermal fluid at about 130oC is exploited at some 300 m depth (personal communication from spa owner). The estimated reservoir temperatures (150-220oC, par. 2.4) and the mixing trend of the geothermal fluid (Segovia et al., 2005) lead to the hypothesis that fluid at T>150oC can be economically harvested somewhere around 500-600 m depth, provided the fractured zone of the main E-W fault bordering the southern bank of the Cuitzeo lake is targeted, particularly near its intersection with B&R structures, to avoid permeability issues during production. Geothermal exploitation at the site may be profitable to the industry and probably also to the local economy, although, on top of the presumed environmental stress delivered to an already polluted area, less heat may be available to local spas after beginning of production.
2 Alignments of cones in the MGVF and origin of the plume
Within MGVF authors have detected small-cones alignments following geologic lineaments, confirming that distribution of volcanoes can be the expression of stress conditions in the crust during activity (Hasenaka, 1984; Hasenaka and Carmichael, 1985; Cebriá et al., 2010). In the southern half of the field, NE-striking alignments of cinder cones are evident (Ban et al., 1992; Figure 13), coinciding with the orientation of the Tenochtitlán Fault System whose mean faults strike parallels the relative motion vector of the Cocos and North American plates (40o; DeMets et al., 1990), defining an approximately NW-SE σ3. These structures controlled the spatial distribution of many monogenetic and polygenetic volcanoes within the TMVB, hold evidences of Holocene reactivation, and disciplined the eruption of >70 cinder cones in MGVF, younger than 40,000yr (Ban et al., 1992; Cebriá et al., 2010; Figure 13). NE-striking lineaments contain Paricutin and Jorullo volcanoes (recently active) at their SW edges, and, if prolonged towards the NE, would interest our study area. Although results of our geophysical survey do not support the existence of NE-SW faults below Cuitzeo lake, these structures exist some km south of the SW shore of the lake (see Figures 2, 4, 7 and 13).
In the northeastern and northwestern parts of MGVF, authors described alignments of cinder cones with E-W strikes (Hasaka and Carmichael, 1985; Ban et al., 1992), clearly related to the Chapala-Tula fault system. These faults, potentially delineating a young rift zone over the axis of TMVB, may also be invoked as cause for the ascent of the plume, below the lake, although their relatively low angle of dip would make this option less preferable.
The last lineament of interest is located at the northeastern sector of the MGVF, where some 13 maars (and many cinder-cones, not represented in figure 13), aged between 1.8 and 0.075 Ma (Aranda-Gomez and Carrasco-Nuñez, 2014) are distributed within an elongated 7 km by 50 km stripe, oriented NNW-SSE. This has been taken as evidence of a pre-existing buried fault system, a zone of deep crustal discontinuity enabling magmas to reach the surface (Murphy, 1982; Uribe-Cifuentes and Urrutia-Fucugauchi, 1999): ascending magma encountered the regional aquifer, favoring the generation of phreatomagmatic eruptions (Aranda-Gomez and Carrasco-Nuñez, 2014). Some of these maars hold lower crustal xenoliths, which is proof of their deep origins (Uribe-Cifuentes and Urrutia-Fucugauchi, 1999; Ortega-Gutiérrez et al., 2014).
Following the most recent models (e.g. Blatter and Hammersley, 2010; Stubailo et al., 2012), this zone of crustal weakness occupy part of the projection on the surface of the tear fault caused at depth by the subduction of the Orozco Fracture Zone (Figure 12a), which has a high heat flow, favoring the ascent of astenospheric magma through the lower crust. Figure 13 highlights how, on a map, the alignment of maars strikes parallel to the continuation of the subducting OFZ and, on the ground, to the Tzitzio-Valle de Santiago fault system, already active during Basin & Range tectonics and still today.
We can further our hypothesis based on the volume of magnetic rock needed to create such an anomaly on the surface. If we assign to the magnetic source a diameter of say 3 km in 2D (Figure 9a and 10) and estimate for the plume a similar thickness (3 km, low appraisal), we would have a volume somehow higher than 20 km3, higher than characteristic values of erupted material for monogenetic cinder cones (Max registered 5 km3, Hasenaka and Carmichael, 1985), let alone maars volcanoes. Alaniz-Alvarez et al. (1998) observed how, over the TMVB, monogenetic small cones are preferentially oriented parallel to east-west normal faults with high slip rates; on the other hand, polygenetic volcanoes align along faults with low displacement rates (north-south B&R faults). The crystallized magmatic body under Cuitzeo lake, affected by different extensional tectonics, may be more relatable to a shield volcano of the MGVF and, in turn, might have enjoyed different recharge of magma, piling below the clay sediments, during its history.
2.1 New results
In 2017 researchers from Universidad Michoacana de San Nicolas de Hidalgo (UMSNH) dated some of the E-W normal faults near Cuitzeo lake at about 18 Ma and the island in the middle of the lake (considered in this study) at 17 Ma, through 40Ar/39Ar analysis. As we detailed through the text, we believe this age for E-W faults of the Morelia-Acambay fault system reflects the imposition of an early TMVB-basement volcanism (and related heat-flow) which generated the right stress-field for the creation of these structures, probably with a pure extensional motion, during early Miocene times. This new insight might have scientists think again on the age of lake Cuitzeo, currently believed to be around 10 Ma, but probably older. The age of the island is another issue, having been related to the still-hot, crystallized, basic plume under the lake, in this work. An age of 17 Ma would not give any chance to a magmatic body to still be at an interesting temperature for feeding a hydrothermal system at present. A volume of magma as the one imaged in this work would cool down in no longer than say 500 ka after underground emplacement and we may fit this new finding (the age of the island) in our model in basically two ways. First, as we explained in the previous paragraph, the specific structural background of the study area (main N-S and E-W faults intersection) can influence the magma supply and provide recharge over time. Hence, while the island may be an early expression of the phenomenon, repeated magma top-ups from below might have influenced the strength of the geothermal system along time. In this case, a direct affinity with MGVF volcanism is less unambiguous. A second loophole in accounting for the age of the volcanic island is that it might not be related at all with the geothermal arrangement at the site. We believe new geophysical and geological studies should focus on the characterization of perspective geothermal sites laying along the TMVB, with the aim of easing future explorations and broadening general knowledge.