The impact of hydrothermal alteration on the physiochemical characteristics of reservoir rocks: the case of the Los Humeros geothermal field (Mexico)

Abstract

Hydrothermal alteration is a common process in active geothermal systems and can significantly change the physiochemical properties of rocks. To improve reservoir assessment and modeling of high-temperature geothermal resources linked to active volcanic settings, a detailed understanding of the reservoir is needed. The Los Humeros Volcanic Complex, hosting the third largest exploited geothermal field in Mexico, represents a natural laboratory to investigate the impact of hydrothermal processes on the rock properties through andesitic reservoir cores and outcropping analogs. Complementary petrographic and chemical analyses were used to characterize the intensities and facies of hydrothermal alteration. The alteration varies from argillic and propylitic facies characterized by no significant changes of the REE budget indicating an inert behavior to silicic facies and skarn instead showing highly variable REE contents. Unaltered outcrop samples predominantly feature low matrix permeabilities (< 10^–17 m²) as well as low to intermediate matrix porosities (< 5–15%), thermal conductivities (0.89–1.49 W m⁻¹ K⁻¹), thermal diffusivities (~ 0.83 10^–6 m² s⁻¹), and sonic wave velocities (V_P: ~ 2800–4100 m s⁻¹, V_S: ~ 1600–2400 m s⁻¹). Average magnetic susceptibility and specific heat capacity range between 2.4–7.0 10^–3 SI and 752–772 J kg⁻¹ K⁻¹, respectively. In contrast, the hydrothermally altered reservoir samples show enhanced porosities (~ 7–23%), permeabilities (10^–17–10^–14 m²), and thermal properties (> 1.67 W m⁻¹ K⁻¹; > 0.91 10^–6 m² s⁻¹), but a significant loss of magnetic susceptibility (10^–3–10^–6 SI). In particular, this latter characteristic appears to be a suitable indicator during geophysical survey for the identification of hydrothermalized domains and possible pathways for fluids. The lack of clear trends between alteration facies, alteration intensity, and chemical indices in the studied samples is interpreted as the response to multiple and/or repeated hydrothermal events. Finally, the proposed integrated field-based approach shows the capability to unravel the complexity of geothermal reservoir rocks in active volcanic settings.

Introduction

Unconventional geothermal systems, such as the high-temperature or so-called super-hot geothermal systems (SHGS, > 350 °C), are important targets for the generation of electric power (Reinsch et al. 2017). The majority of the currently exploited high-temperature geothermal resources are related to magmatic settings, predominantly to active volcanic systems (Heřmanská et al. 2019). The petrophysical characterization of volcanic systems is often challenging due to their complex structural and geological evolution, resulting in various rock types that exhibit complex characteristics regarding mineralogy, petrophysical and rock mechanical behavior (Pola et al. 2012). Additionally, hydrothermal alteration is a widely observed process in volcanic settings, which significantly changes the intrinsic physiochemical properties of pristine rocks (Pola et al. 2016; Mordensky et al. 2019a, b; Durán et al. 2019; Pandarinath et al., 2020; Heap et al. 2019, 2022). The circulation of hydrothermal fluids leads to mineral dissolution, replacement, or precipitation (Browne 1978) causing partial to complete modification in the mineralogical composition and subsequent physical behavior of the affected rocks (Wyering et al. 2014). Thereby, the impact of hydrothermal alteration on the rock properties strongly depends on the original rock types and their mineralogical and petrophysical characteristics, pressure, and temperature conditions within the reservoir, chemical composition of the reservoir fluid, as well as the duration of the fluid–rock interaction (e.g., Sillitoe 2010; Frolova et al. 2014; Rabiee et al. 2019; Sillitoe and Brogi 2021).

The influence of hydrothermal alteration on rock properties of volcanic rocks has been investigated in the past to assess not only the longevity of geothermal reservoirs (Villeneuve et al. 2019) but also regarding slope instability (Sánchez-Núñez et al. 2021), heat flux, volcanic activity, and the possible impact on phreatic eruptions (Mayer et al. 2016; Heap et al. 2019, 2022). Previous studies aimed to identify general trends in altered rocks, e.g., increased porosity and permeability along with reduced rock strength (Pola et al. 2014; Wyering et al. 2014), reduced permeability due to mineral precipitation (Mordensky et al. 2018), reduced porosity and permeability due to silicification (Dobson et al. 2003), or reduced porosity and increased thermal conductivity associated with an increased degree of hydrothermal alteration (Mielke et al. 2015), respectively. However, such changes are highly variable, even within a single reservoir, and are not fully understood yet. Comprehensive datasets are rarely available and most reservoir assessment studies or models are based on assumed or generalized data, and thus, effects of fluid-rock interactions on rock properties and their spatial extent are commonly neglected for simplification reasons. Since hydrothermal alteration can both increase or decrease key reservoir properties, e.g., matrix porosity, permeability, rock strength, it also impacts the economic potential of the geothermal reservoir, and the controlling factors need to be understood and considered during reservoir assessment.

To investigate the influence of hydrothermal alteration on the physiochemical properties of rocks, we investigated borehole core samples retrieved from an active hydrothermal geothermal reservoir and compared the results with outcrop samples that are stratigraphically equivalent to the reservoir units. For this purpose, the Los Humeros Volcanic Complex (LHVC) hosting one of the largest active silicic calderas located in the north-eastern part of the Trans-Mexican Volcanic Belt (TMVB) was selected as case study. The LHVC hosts a steam-dominated hydrothermal system in predominantly andesitic sequences with temperatures above 380 °C below 2 km in the northern part of the caldera (Pinti et al. 2017; Deb et al. 2021a, b). The geothermal reservoir has been operated by the Comisión Federal Electricidad (CFE) since 1990 (Romo-Jones et al. 2020) and 65 boreholes have been drilled so far. However, a sustainable utilization of these so-called “super-hot” regions for power production has not yet been possible due to aggressive reservoir fluids causing corrosion and scaling problems (Flores-Armenta et al. 2010). Various geological (Ferriz and Mahood 1984; Carrasco-Núñez et al. 2017a, b), geochemical (Prol-Ledesma and Browne 1989; Martı́nez-Serrano 1993, 2002; Izquíerdo et al. 2011, 2015), geophysical (Lermo et al. 2008; Arzate et al. 2018), and hydrological studies (Tello 2005) have been performed in the past and most recently within the framework of the GEMex project (EU-H2020, GA Nr. 727,550; Jolie et al. 2018; Weydt et al. 2018, 2021b; Lucci et al. 2020; Urbani et al. 2020, 2021) to improve reservoir understanding and to create conceptual (Cedillo 2000; Arellano et al. 2003) and 3D geological models (Calcagno et al. 2020; Deb et al. 2019). However, even after more than 40 years of exploration, information on rock properties of the different geological units in the study area were scarce or not available. To overcome the lack of suitable data for economic assessments, numerical reservoir simulations, and the interpretation of geophysical surveys, an extensive outcrop analog study has been performed to characterize all relevant key units regarding their geochemistry, mineralogy, petrophysical, and mechanical rock properties (Weydt et al. 2021a, 2022).

This study focuses on the geochemical and petrophysical characterization of mainly andesitic lavas collected from outcrops and borehole core samples (from here on defined as reservoir samples) that are considered to represent the reservoir units hosting the hydrothermal system in the Los Humeros geothermal field. Therefore, 66 reservoir samples were drilled from 37 core sections covering 16 boreholes (H7 until H40). Additionally, 24 samples of Miocene to Pleistocene lavas of the Teziutlán and Cuyoaco andesite units representing the stratigraphically equivalent analogs to the reservoir formations in the deep subsurface were collected from 14 outcrops located in the surrounding area of the LHVC. The samples were analyzed for bulk and particle density, porosity, permeability, thermal conductivity, thermal diffusivity, specific heat capacity, as well as ultra-sonic wave velocities and magnetic susceptibility. Polarized light and scanning electron microscopes (PLM and SEM, respectively) investigations were performed to evaluate the intensity of hydrothermal alteration and to define alteration facies. Furthermore, whole-rock compositions were determined and chemical weathering indices (CWI) were applied to quantify elemental concentration changes, and thus, the impact of fluid-rock interactions (Pola et al. 2012; Pandarinath et al. 2020) and the post-alteration physiochemical properties. Correlation of reservoir samples to outcropping analogs were done following the approach proposed by Carrasco-Núñez et al. (2017b) and using selected elements considered to be unaffected by hydrothermalism and capable to act as geochemical indices (Winchester and Floyd 1977; MacLean and Barrett 1993).

The dataset presented in this study aims to provide new insights on the impact of hydrothermal alteration on reservoir rocks in high-temperature hydrothermal systems. It will help future reservoir assessments define more accurate estimations of petrophysical properties at reservoir depths—both at the Los Humeros geothermal field and also geologically similar super-hot geothermal systems elsewhere in the world.

Geological setting

The LHVC is located in the eastern sector of the Trans-Mexican Volcanic Belt (TMVB; Fig. 1, inset box). The TMVB is an E–W trending Miocene-Holocene continental volcanic arc some 1000 km long, which is linked to the subduction of the Rivera and Cocos plates beneath the North American plate along the Middle-American Trench (López-Hernández et al. 2009; Ferrari et al. 2012). The LHVC occupies an area approximately 21 × 15 km (Ferriz and Mahood 1984) and is the largest active caldera of the TMVB, predominantly comprising Pleistocene to Holocene basaltic to rhyolitic volcanic rocks (Carrasco-Núñez et al. 2017a, 2017b, 2018; Cavazos-Álvarez et al. 2020; Lucci et al. 2020). The geology in the study area (Fig. 1) is characterized by (1) pre-volcanic basement units, (2) pre-caldera volcanism units, (3) caldera stage units, and (4) post-caldera stage units (Carrasco-Núñez et al. 2017a, 2018).

Fig. 1

Geological map of the Los Humeros Volcanic Complex (LHVC) modified from Carrasco-Núñez et al. (2017a). The red points mark the sampling locations of the outcrop samples collected from the Cuyoaco andesites (Tm, with purple color) and the Teziutlán andesites (Tpa, pink color). Inset map showing the location of the LHVC and extension of the TMVB

The pre-volcanic basement is represented by metamorphic (greenschists) and intrusive (granites and granodiorites) rocks dated at 246–131 Ma mainly exposed in the Teziutlán Massif and partially covered by up to 3 km thick Mesozoic sedimentary rocks belonging to the Sierra Madre Oriental (Yáñez and García 1982). The Mesozoic sedimentary successions comprise Jurassic and Cretaceous limestones, dolomites, shales, marls, and sandstones. NW–SE striking thrusts and folds formed during the Late Cretaceous–Eocene compressive Laramide Orogeny and subordinate NE-striking normal faults developed during the Eocene–Pliocene extensional tectonic deformation phase (Fitz-Díaz et al. 2017; López-Hernández et al. 1995). Andesitic to basaltic dykes as well as granitic to syenitic plutons (Oligocene–Miocene) intruded into the sedimentary sequences and led to local thermal metamorphism producing marble, hornfels, and skarn (Ferriz and Mahood 1984).

The pre-caldera volcanism in the study area started in the Late Miocene at ~ 10.5 ± 0.7 Ma (Yáñez and García 1982) with the emplacement of the Cuyoaco and Alseseca lavas comprising mainly fractured andesitic and dacitic lava flows with a cumulative thickness of 800–900 m. These lava flows can be related to the activity of the Cerro Grande volcanic complex dated between 8.9 and 11 Ma (Carrasco-Núñez et al. 1997; Gómez-Tuena and Carrasco-Núñez 2000). The Pliocene to Pleistocene pre-caldera volcanism is instead represented by the Teziutlán andesites dated between 1.44 ± 0.31 and 2.65 ± 0.43 Ma (Carrasco-Núñez et al. 2017a), mainly comprising fractured dark gray massive, porphyritic andesitic lavas composed of plagioclase and two-pyroxene phenocrysts as well as olivine-bearing basaltic lavas (Carrasco-Núñez et al. 2017a). These late Miocene to Pleistocene fractured pre-caldera andesites show a thickness of up to 1500 m as reported in the lithostratigraphic profiles of the geothermal wells (Fig. 2; Carrasco-Núñez et al. 2017b; López-Hernández et al. 1995) and constitute the currently exploited reservoir of the Los Humeros geothermal field. In previous studies, the pre-caldera reservoir samples were described as augite andesites, tuff, hornblende andesites, and basalts (from top to base; Cedillo 2000).

Fig. 2

a Location of the geothermal wells drilled within the Los Humeros caldera and b the lithostratigraphic profile of well H43 (modified from Carrasco-Núñez et al., 2017b). c Lithostratigraphic profiles and in-depth correlation of the main stratigraphic groups modified from Urbani et al. (2020, 2021) and Carrasco-Núñez et al. (2017b) with new data from Cavazos-Álvarez et al. (2020)

Despite the variability in chemical composition, the Cuyoaco and Teziutlán lavas are commonly referred to as pre-caldera andesites or Teziutlán and Cuyoaco andesites for simplicity (e.g., Arellano et al. 2003; Calcagno et al. 2020). In this article, we will also use these terms for simplification reasons.

The onset of magmatic activity of the LHVC caldera stage is represented by partially buried rhyolitic lavas and abundant rhyolitic domes mainly located at the western side of the LHVC, which were dated between 270 ± 17 and 693 ± 1.9 ka (Carrasco-Núñez et al. 2018).

Two main caldera-forming eruptive events separated by a sequence of several large plinian and sub-plinian eruptive phases took place during the caldera stage (Carrasco-Núñez et al. 2021). The Los Humeros caldera collapse is associated with the emplacement of the Xáltipan ignimbrite (164 ± 4.2 ka; Carrasco-Núñez et al. 2018), which is well exposed in the surrounding area of the LHVC, but has a highly variable thickness (90–780 m) within the reservoir (Cavazos-Álvarez et al. 2020). Afterward, a sequence of explosive episodes occurred at 70 ± 23 ka (Carrasco-Núñez et al. 2018) depositing thick rhyodacitic Plinian deposits named as the Faby Tuff (Ferriz and Mahood 1984; Willcox 2011). The second caldera forming event occurred shortly after the Faby Tuff emplacement and is related to the rhyodacitic to andesitic Zaragoza ignimbrite at 69 ± 16 ka (Carrasco-Núñez et al. 2018) forming the smaller Los Potreros nested caldera (8 to 10 km in diameter).

The LHVC post-caldera stage represents the most recent volcanic activity in the study area and can be divided into two eruptive phases, which are a Late Pleistocene resurgence phase followed by a Holocene reactivation phase (Carrasco-Núñez et al. 2021). The Late Pleistocene phase is characterized by rhyolitic and dacitic domes within the center of the Los Humeros caldera center (44.8 ± 1.7 ka, Carrasco-Núñez et al. 2018) and north of the Los Humeros caldera (55.7 ± 4.4 ka, Carrasco-Núñez et al. 2018) followed by a sequence of explosive eruptions producing dacitic pumice fall units (Xoxoctic Tuff; Ferriz and Mahood 1984), volcaniclastic breccias, and pyroclastic flow deposits (Llano Tuff, Ferriz and Mahood 1984; Willcox 2011). The Holocene includes several basaltic to trachyandesite lava flows as the result of alternating episodes of effusive and explosive eruptions (8.9 ± 0.03 ka, Carrasco-Núñez et al., 2017a), followed by a rhythmic alternation of contemporaneous basaltic and trachyandesite fall out deposits constituting the Cuicuiltic Member (7.3 ± 0.1 ka, Dávila-Harris and Carrasco-Núñez 2014), a ring-fracture episode erupting trachyandesite and olivine basalt lava flows (3.9 ± 0.13 ka, Carrasco-Núñez et al. 2017a) and trachytic lava flow at 2.8 ± 0.03 ka (Carrasco-Núñez et al. 2017a).

In particular, the Holocene LHVC post-caldera volcanic activity is associated with a complex magma plumbing system vertically distributed within the whole crust (Lucci et al. 2020) and shallow intrusions in the Los Potreros caldera within the Los Humeros geothermal field (< 1 km depth; Urbani et al. 2020, 2021; Deb et al. 2019), which causes several very localized temperature anomalies mainly along the NNW–SSE trending Maxtaloya and Los Humeros lineament (Jentsch et al. 2020; Deb et al. 2019). Except for the water-producing well H1, the produced fluids are predominantly steam with an enthalpy of more than 2000 kJ kg⁻¹ (Romo-Jones et al. 2020). The geothermal field contains low-saline NaCl to H₂CO₃-SO₄ fluids, which are oversaturated in quartz and calcite and locally contain high boron, ammonia, and arsenic concentrations (Izquíerdo et al. 2009). Most of the geothermal fluids show a near neutral pH between 7 and 8 (Tello 2005), while a few wells in the central to northern part of the geothermal field feature very low pH brines (< 5, Flores-Armenta et al. 2010). The Los Humeros geothermal field is characterized by a trap-door collapse structure (Bonini et al. 2021) and the topographic top of the sedimentary basement shows a high displacement between the northern-central part of the geothermal field (~ 1720 m bgl [below ground level] in H43; north of Los Humeros in Fig. 2) and the southern sector close to the Xalapazco crater (deeper than 3300 m bgl in H24).

Up to now, hydrothermal alteration was predominantly investigated on cuttings retrieved from boreholes (Prol-Ledesma and Browne 1989; Prol-Ledesma 1990; Martı́nez-Serrano 2002; González-Partida et al. 2022). Prol-Ledesma and Browne (1989) defined four laterally distributed areas through the identification of secondary mineral assemblages and defining temperatures ranging from 250 to > 300 °C in the northern-central part of the geothermal field and lower temperatures of ~ 120 °C moving to the caldera rim. Martı́nez-Serrano (2002) defined three main alteration zones from top to base based on cutting analyses of five boreholes in the central collapse zone: (1) a shallow argillic zone which is present in the upper levels of the volcanic sequences and characterized by kaolinite–montmorillonite, chlorite, and zeolites as well as calcite; (2) a propylitic zone between 500 and 1800 m depth mainly comprising epidote, quartz, amphibole, calcite, chlorite, montmorillonite–illite, sulfides, and iron oxides; and (3) a skarn zone at depths of greater than 1800 m composed of garnet, wollastonite, clinopyroxene, and biotite. González-Partida et al. (2022) redefined these units and proposed a high-temperature high-silica advanced argillic alteration zone in the deeper sections of the pre-caldera andesites associated with subcritical brines.

Material and methods

During the field campaigns 24 samples with a dimension of ~ 30 × 30 × 20 cm were collected from 14 outcrops (Table 5) representing the Teziutlán and Cuyoaco andesite units (Fig. 1). Cylindrical sub-samples with diameters ranging from 25 to 64 mm were drilled from the outcrop samples. In addition, 66 sub-samples with a diameter of 40 mm and an approximately length of 75 mm were drilled from 37 wellbore cores of the Los Humeros geothermal field (Fig. 2, Table 6). Following the international standard ASTM D4543 (2019), 207 and 99 plugs were prepared from the outcrop and reservoir core samples (length: ~ 25–128 mm), respectively. Given the sample size used in this study, the petrophysical measurements presented here provide matrix properties of the rocks, eventually including small-scale or single fractures.

The methodology of the petrophysical and geochemical measurements is described in detail in Weydt et al. (2021a) and thus is only mentioned briefly in the following sections. All petrophysical measurements were performed under ambient laboratory conditions (0.1 MPa and ~ 20 °C). To ensure reproducibility of the results, the plugs were oven-dried and cooled down to room temperature in a desiccator (20 °C). A vacuum desiccator (approx. − 1 bar) filled with de-ionized water was used to saturate the samples.

Particle and bulk density measurements were accomplished by using a helium pycnometer (AccuPyc 1330) and a powder pycnometer (GeoPyc 1360); thereby, the particle and bulk volume were measured five times, respectively. Matrix porosity was then calculated from the resulting differences in volume and represents the gas-effective or so-called connected porosity (measurement reproducibility = 1.1%; Micromeritics 1997, 1998).

Matrix permeability was determined with a column permeameter constructed according to ASTM D4525 (2013) standard using dried compressed air at five air pressures (1–3 bar) to obtain the apparent permeability. The apparent permeability was then corrected by applying the Klinkenberg method (Klinkenberg 1941) to determine the intrinsic matrix permeability. According to Bär (2012), the measurement error varies from 5% for high permeable rocks (K > 10^–14 m²) to 400% for impermeable rocks (K < 10^–16 m²).

A thermal conductivity scanner was used to measure thermal conductivity and diffusivity after Popov et al. (2016). Both parameters were simultaneously measured four to six times on each plug under saturated and dry conditions, respectively. According to the manufacturer (Lippman and Rauen 2009) the measurement error is 3% for thermal conductivity and 5% for thermal diffusivity.

A heat-flux differential scanning calorimeter was used to determine specific heat capacity (measurement error = 1% according to Setaram Instrumentation 2009). Over a period of 24 h, crushed sample material was heated from 20 up to 200 °C and subsequently cooled down to room temperature. Thereby, specific heat capacity was calculated from the resulting temperature curves through heat flow differences. Volumetric heat capacity was calculated by multiplying the specific heat capacity with the associated bulk density of each sample.

An ultrasound generator (UKS-D) from Geotron-Elektronik (2011) was used to determine compressional and shear wave velocities. Continuous measurements were performed with a contact pressure of 0.1 MPa and a frequency of 80 kHz to 250 kHz. P- and S-wave velocities were measured simultaneously four to six times on each plug under saturated and dry conditions, respectively, whereby the arrival times were picked manually.

A magnetic susceptibility meter SM30 from ZH Instruments (2008) was used to measure magnetic susceptibility directly on the plane surface of each plug. To account for mineralogical heterogeneities each plane surface of a plug was measured five times.

The classification of Bär (2012) was used to evaluate the rock properties regarding their geothermal potential (ranging from very low to very high).

To analyze the bulk chemical composition of the samples whole-rock X-ray fluorescence spectroscopy (XRF) analyses were performed using Panalytical Axios spectrometers (Max WD-XRF at TU Delft and Advanced at GFZ) in combination with the software SuperQ for data evaluation as well as a Bruker S8Tiger 4 WD-XRF spectrometer using the Quant Express method (at TU Darmstadt). Measurement error is < 5% for the major elements and < 10% for the trace elements. The limit of detection varies between 400 ppm (Na) and < 10 ppm (e.g., Rb, Sr, Nb). Rare-earth elements (REE) of the reservoir samples and selected outcrop samples were analyzed at TU Darmstadt using an Analytik Jena Plasma Quant MS Elite^® ICP-MS based on an application note by Analytik Jena and Yu et al. (2001). Method validation and recovery experiments were carried out using two certified reference standards (Basalt, Columbia River [BCR-2] and Andesite, Guano Valley [AGV-2]-Fluxana Bedburg-Hau, Germany) as well as sample material (cuttings) from Carrasco-Núñez et al. (2017b). Geochemical data were analyzed using a nonparametric statistical approach (Kruskal–Wallis test with p = 0.05, Vargha and Delaney 1998; XLSTAT Premium V2020.5.1). A detailed analysis of clay minerals or aluminosilicates and their distribution using X-ray Diffraction methods has already been performed in previous studies (e.g., Prol-Ledesma 1990; Martı́nez-Serrano 2002; González-Partida et al. 2022) and thus was not repeated in this study.

Furthermore, the collected samples were prepared as polished thin sections and studied by polarized light microscope (PLM). Selected reservoir samples were then analyzed with a scanning electron microscope (SEM; further details are included in the supplementary material) to gain information on the primary assemblages, intensity of alteration, mineralogy and alteration facies, texture, porosity, (micro-)fractures, and grain size and shape. Mineral abbreviations used in this study follow Whitney and Evans (2010). Samples were classified into three groups based on the intensity of alteration and macroscopic features (e.g., change of structure, color, intensity of mineral replacement, or neoformed minerals; Pola et al. 2012; Navelot et al. 2018), including “none-to-weak alteration” (pristine mineralogy to first signs of alteration), “moderate alteration” (hydrothermal alteration impacts the groundmass or phenocrysts to a variable extent), and “strong alteration” (alteration is pervasive and strongly affects the groundmass and phenocrysts).

Several chemical weathering indices (CWI) were tested, such as chemical index of weathering (CIW; Harnois 1988), chemical index of alteration (CIA; Nesbitt and Young 1982), plagioclase index of alteration (PIA; Fedo et al. 1995), Ishikawa alteration index (AI; Mathieu, 2018), and chlorite-carbonate-pyrite index (CCPI; Mathieu 2018) to investigate the relationship between chemical changes and alteration intensity as described in Pola et al. (2012) or Lee et al. (2021). These indices have been developed to better characterize the alteration intensity using weight percentage of major elements and assume that the (re-)distribution of these elements is controlled by the alteration intensity.

Furthermore, variations of the major, trace, and rare-earth elemental concentrations of the reservoir samples were calculated according to Pandarinath et al. 2020 using the outcrop analogs as preserved protoliths.

Results

Petrographic analyses

Outcrop samples

The Teziutlán andesite unit comprise different massive to porous, predominantly porphyritic lavas with a dark gray to medium gray as well as brownish-gray color (Fig. 3a–h; Table 1). The primary assemblage is dominated by plagioclase with pyroxene (clinopyroxene and orthopyroxene), olivine and subordinate (Ti-)magnetite and ilmenite. The groundmass is made up of a comparable microcrystal assemblage. The samples LH13-LH17, RLM1, RLM2, and RLM9 show no macroscopically visible pores. The samples LH8-LH12 contain a vesicular matrix with abundant circular to elongated pores (1–5 mm in diameter), while LH29, LH30, and RLM3 contain only a few pores with an irregular distribution. Occasionally fractures occur. RLM3 shows alteration rims along fine cracks and fissures and contains some secondary alteration minerals (most likely clay minerals). The lavas related to the late Miocene Cuyoaco andesite unit show a variable gray to beige and reddish color with a porphyritic to glomerophyric texture comprising plagioclase, clinopyroxene (augite), olivine, and ilmenite (Fig. 3i–n, Table 1). The nonporous lavas often show small-scaled fractures. LH40, LH41, LH72, as well as OLH1 samples show a weak alteration characterized by overprints at outer rims of pyroxenes and appearance of secondary oxides such as hematite (Fig. 3i–j). LH50 and LH51 are significantly different with a beige, nonporous microcrystalline matrix containing plagioclase, quartz, pyroxene, hornblende, and biotite (Fig. 3o–p). These samples were interpreted as subvolcanic/hypabyssal rocks.

Fig. 3

Scans of polished plugs and corresponding thin section images (XPL) from selected outcrop samples. a–h represent the Teziutlán andesite unit and i–p represent the Miocene Cuyoaco andesite unit. Cpx, Clinopyroxene; Hbl, hornblende; p, pores; Pl, plagioclase; Opx, Orthopyroxene

Table 1 Summary of petrographic features for the outcrop samples

Reservoir samples

The investigated reservoir core samples include a highly variable suite of magmatic and metamorphic rocks from aphanitic to porphyritic basaltic andesitic to rhyolitic lavas and ignimbrites as well as skarns and marbles collected from 37 core sections drilled at 353 m to 2900 m bgl (Table 2, Figs. 4, 5, 6, 16, and 17). The core sections cover 16 boreholes, thus providing good spatial coverage of the inner part of the caldera. Most of the samples were affected by hydrothermal alteration of different intensities, brecciation, and fracturing, thus resulting in high matrix heterogeneity (Fig. 16). Hydrothermal alteration varies from weak to strong and predominantly occurs along fractures and cracks with a limited lateral extension, often of a few centimeters only. For example, H39-2-C2 shows intense hydrothermal alteration and bleaching along a fracture (Fig. 6a), while the second plug (H39-2-C1), which was drilled only a few centimeters away, shows only a weak to moderate hydrothermal overprint. In some cases, the rock matrix was completely replaced by secondary minerals, thus preventing a clear identification of the original protolith.

Table 2 Summary of petrographic features for the reservoir core samples

Fig. 4

Scans and photographs of selected reservoir samples and corresponding thin sections with PPL = (b) and (d) and XPL = (f), (h), (j), and (l). a to d represent the upper Xáltipan ignimbrite and andesite/ignimbrite breccias of the caldera group, e to l represent mafic to rhyolitic lavas of the upper part of the pre-caldera group predominantly affected by argillic alteration. Cal, calcite; Chl, chlorite; Cpx, Clinopyroxene; Hem, hematite; Mag, magnetite; p, pores; Pl plagioclase; Qz, quartz

Fig. 5

Low-temperature (a–d), moderate-temperature (e–h) to high-temperature propylitic alteration (i to l) as well as advanced silicification (m–p) with relicts of original phenocrysts and secondary porosity (black patches) observed in thin sections with PPL = (b) and (f) and XPL = (c), (d), (g), (k), and (n) as well as SEM images of selected altered andesitic lavas displayed as (h), (l), (o), and (p). Cal calcite; Chl chlorite; Ep, epidote; Czo, clinozoisite; Fsp, feldspar; Kfs, K-feldspar; p, pores; Pl, plagioclase; Py, pyrite; Qz, quartz

Fig. 6

Scans and photographs of selected samples and corresponding thin sections (PPL = b, e, g, and XPL = i) and SEM images (c, f, k) representing a to c andesitic lava affected by propylitic alteration and silicification, d–g high-temperature propylitic alteration in a fracture of an andesitic lava, h, i skarn, and j–k marble. Ap, apatite; Cal, calcite; Chl, chlorite; Di, diopside; Ep, epidote; FeOX, iron oxides; Grt, garnet; Grs, grossular garnet; Kfs, K-feldspar; Mag, magnetite; p, pores; Pl, plagioclase; Py, pyrite; Qz, quartz; Ttn, titanite; TiOx, titanium oxide; Wol, wollastonite

The shallowest sample H26-1 retrieved from 353 m bgl comprises a beige to light gray vitroclastic tuff that most likely represents the upper part of the Xáltipan ignimbrite unit (Fig. 4a, b). The sample shows only a weak hydrothermal overprint.

The samples retrieved from 633 to 1400 m bgl were identified as mafic, andesitic, trachyandesite, and rhyolite lavas as well as andesite–ignimbrite breccias (at 633–673 m bgl, Fig. 4c, d). The samples have a highly variable appearance and texture again characterized by strongly variable alteration intensities, from weak to strong. While some samples represent relatively fresh lavas showing unaltered plagioclase, olivine, and clinopyroxene (H27-1 at 1110 m bgl, Fig. 4e, f), other lavas retrieved from a similar depth level are strongly altered (H18-1 at 796 m bgl). This section of the reservoir is dominated by argillic alteration as well as the beginning of propylitic alteration with the main alteration assemblages made up of calcite, quartz, chlorite, hematite, and clay minerals, such as kaolinite and mixed layers smectite–illite. Calcite is very abundant and occurs in pores, fractures, and pervasively in the matrix replacing also both plagioclase and K-feldspar phenocrysts. Allotriomorphic to idiomorphic quartz crystals predominantly grow at the outer rim of pores and fractures. Chlorite appears as light green, very fine-grained patches or very fine radial aggregates in pores and fractures, especially in the upper depth levels together with other clay minerals. Epidote occurs locally as small subhedral to euhedral crystals overgrowth on plagioclase. Pores are usually smaller than 0.5 mm in diameter. However, one sample of strongly altered andesitic lava (H20-4, 1400 m bgl, Fig. 5a–d) contains largely elongated, partially filled pores of up to 10 × 30 mm in size, leading to a tuffaceous appearance.

Samples retrieved from ≥ 1500 m depths (bgl) predominantly comprise basaltic to andesitic lavas (Table 2) and are dominated by propylitic alteration with an alteration assemblage composed of epidote, chlorite, albite, K-feldspar, quartz, calcite, pyrite, hematite, titanite, Ti-Fe-oxides, as well as occasionally clinozoisite. The abundance of chlorite and epidote led to the typical greenish-yellow color and patchy appearance of the rock matrix (Figs. 5, 16). The primary plagioclase is completely transformed to albite, calcite, and epidote with also the appearance of secondary K-feldspar, chlorite and subordinate apatite, and Ti-bearing minerals, such as titanite and Ti oxides. Likewise, clino- and orthopyroxene phenocrysts are completely altered. Amygdales are filled with chlorite, epidote, and calcite rimmed with quartz (H19-3, Figs. 5f, 17). However, in several samples some pores remained unfilled (Fig. 5). Fractures, cracks, and brecciated sections are abundant and are also typically filled with an assemblage made up of chlorite, epidote, calcite, quartz, subordinate pyrite, and Fe-oxides (Figs. 5, 16). Magnetite occurs in some samples in the most preserved portions (H39-2, Fig. 6b), but is commonly substituted by secondary titanite, hematite, and other Fe-oxides. Clinozoisite locally occurs associated with epidote s.s. and chlorite. In samples with reduced Fe₂O₃ concentrations, clinozoisite occurs as nests of euhedral radiating crystals growing in vugs or on plagioclase (Fig. 5i–l). The H38-4 sample, which originated from 1950 m bgl represents an exceptional case comprising ≤ 5 mm euhedral to subhedral crystals of andradite garnet intergrown with elongate and blocky calcite crystals in fractures and in the matrix along with notable secondary porosity (Lacinska et al. 2020). The garnets are rich in mineral inclusions, predominantly clinopyroxene (diopside–augite) with a lesser amount of calcite, pyrite, and wollastonite. The latter was also observed as acicular crystals on the garnet–calcite intergrowths. Furthermore, rare amphibole phenocrysts altered to chlorite–epidote assemblage were observed (Fig. 6d–g).

Two skarn samples were retrieved from the deepest section of the reservoir, close to the Xalapazco crater (2844–2900 m bgl). The samples are composed of quartz, calcite, apatite, and possibly acicular wollastonite with very small garnet crystals growing on the crystal surfaces as well as actinolite and clays. The significant leaching of elements and replacement of the matrix with quartz and other secondary minerals led to a white to grayish color of the samples (also called bleaching, Fig. 6h). However, pervasive silicification was also observed at shallower depth levels (~ 2000 m bgl, H26-4 (Fig. 16) and H28-2). SEM investigations and X-ray mapping revealed that the samples underwent intensive chemical changes, including the dissolution of the original phenocrysts (Figs. 5o–p, 17), the development of secondary porosity, and extensive silicification of the groundmass with locally mica, chlorite, as well as some radial nests of pore-filling biotite.

Intensively fractured marble samples were retrieved from cores in the central part of the caldera (1970–2414 m bgl), which represent the upper section of the pre-volcanic basement (Figs. 6j–k, 17). They predominantly consist of a mosaic of interlocking coarsely crystalline calcite cross-cut by numerous calcite-filled fractures and veins that occasionally include wollastonite, garnet and subordinate diopside, apatite, pyrite, epidote, and sulfides. Due to the dense array of calcite fractures and veins, the marbles are very friable.

Whole-rock chemistry

The major, trace, and rare-earth elemental chemistry of studied samples is presented in Figs. 7, 8, and 9 as well as in Additional file 1: Table S1. The outcrop samples show basaltic andesitic to dacitic composition (Fig. 7a) with a general calc-alkaline character (Rickwood 1989). The SiO₂ content of the Teziutlán and Cuyoaco andesite units ranges from 53 to 62 wt% and from 61 to 65 wt%, respectively with Na₂O + K₂O ranging from 4.8 to 6.6 wt% and from 5.4 to 6.1 wt%, respectively. Likewise, the reservoir samples consist of varying SiO₂ concentrations ranging between ~ 52 wt% and ~ 77 wt% (basaltic to rhyolitic composition). Harker variation diagrams of selected major elements using SiO₂ wt% as differentiation index are shown in Fig. 7b–i. Increasing SiO₂ concentrations are typically correlated with decreasing concentrations of CaO, MgO, Fe₂O_3, Al₂O_3, P₂O₅, and TiO₂ as well as increasing K₂O concentrations. In contrast, Na₂O vs. SiO₂ shows no distinct trend. Samples affected by advanced silicification and skarns show increased SiO₂ concentrations and plot in the same range as the rhyolitic lavas and Xáltipan ignimbrite (samples H19-1 and H26-1, respectively). Furthermore, they are characterized by lower Fe₂O_3, P₂O₅, and TiO₂ contents as well as highly variable CaO, MgO, and Na₂O values. Similar compositional changes were observed for trace elements and REE (Figs. 8 and 9). Trace elements, such as Nb, Sr, and Zr, and related ratios (e.g., Zr/Al₂O_3, Ti/1000, Sr/Nb, and Nb/La), which are usually believed to be immobile and constant, respectively, during hydrothermal alteration (Winchester and Floyd 1977; Floyd and Winchester 1978; MacLean and Barrett 1993; Carrasco-Núñez et al. 2017b; Pandarinath et al. 2020) show instead great variability and scattering.

Fig. 7

a Total alkali versus silica (TAS) diagram (Le Maitre et al. 2002) of the pre-caldera andesite outcrop samples. Selected major-element Harker diagrams of the outcrop and reservoir samples in (b–i). The pre-caldera reservoir samples were classified regarding their degree of hydrothermal alteration observed during petrographic analyses (weak, moderate, and strong)

Fig. 8

Selected trace and REE element concentrations of the reservoir and outcrop samples in (a–h). i Nb/La vs. Sr/Nb plot after Carrasco-Núñez et al. (2017b). Chemical indices such as (j) CIW, (k) CIA and (l) PIA values were applied to investigate the relationship between the intensity of alteration and the physiochemical changes of the rocks

Fig. 9

REE/Chondrite-normalized values of (a) the outcrop and (b, c) reservoir samples after Sun and Mc Donough (1989)

Chondrite-normalized (Sun and McDonough 1989) spider diagrams for the REE are presented in Fig. 9. The studied andesitic lavas are characterized by fractionated patterns with LREE (light rare-earth elements) enrichment and gently sloping profiles moving to HREE (heavy rare-earth elements) with no evidence of Eu anomaly (Fig. 9a). Most of the analyzed Teziutlán andesites are characterized by higher REE concentrations compared to the Cuyoaco andesites. To note, a relevant number of samples from both Teziutlán and Cuyoaco andesites show a positive Sm anomaly. The reservoir samples are characterized by REE patterns overlapping with those of Teziutlán lavas. The rhyolite sample H19-1 and the shallow-depth cuttings from H20 and H43 boreholes show fractionated REE patterns with a well-developed Eu negative anomaly (Fig. 9b). A slight Ce positive anomaly is instead reported in samples H13-1, H18-1, and H39-1 characterized by a strong alteration overprint dominated by calcite forming and collected between ~ 800 and 1200 m bgl, from boreholes located in the southern sector of the geothermal field, close to the Los Humeros-Maxtaloya lineament. Skarns and silicified samples are characterized by fractionated REE patterns (Fig. 9c) showing LREE enriched gently dipping profiles generally overlapping with both reservoir samples and outcropping Teziutlán and Cuyoaco andesites.

The results of the chemical indices are presented in Table 3 and Additional file 1: Table S2. In Table 3, the results of the reservoir core samples are classified with respect to their alteration facies and alteration intensity, respectively. Samples affected by strong silicification show increasing CIW, CIA, PIA, and AI values, whereas skarns are again characterized by strongly variable CIW, CIA, and PIA indices values. Both silicified andesites and skarns are characterized by decreasing CCPI values with respect to non-altered and analog outcropping samples (Fig. 8).

Table 3 Chemical indices calculated for the outcrop and reservoir samples with respect to alteration facies and intensity

Petrophysical properties

The results of the petrophysical measurements are presented in Figs. 10, 11, 12, and 13 and in Table 4. Particle densities of the outcrop and reservoir samples are within the same range and vary between 2.59 and 2.88 g cm⁻³ (Fig. 10e and f). Exceptions form the highly porous ignimbrites (2.43 g cm⁻³) and sample H38-4 (2.99 g cm⁻³), which comprises a high amount of calcite and large garnet crystals. Bulk densities of the outcrop samples range between 2.00 and 2.78 g cm⁻³ and bulk densities of the reservoir samples range between 1.63 and 2.78 g cm⁻³, in ignimbrites and skarn, respectively. Thereby, bulk density is negatively correlated with porosity, and thus, is positively correlated with (i) thermal conductivity, (ii) thermal diffusivity, and (iii) to a lesser extent with P-wave and S-wave velocities. While particle and bulk densities of the outcrop samples show location specific clusters (Fig. 10a and f), the reservoir samples show a higher degree of scattering without any clear lithological trends.

Fig. 10

Correlation of density, P-wave velocity, magnetic susceptibility, and thermal conductivity analyzed under dry conditions (for the number of analyzed plugs see Table 4)

Fig. 11

Correlation of permeability, porosity, and thermal conductivity analyzed under dry conditions (for number of analyzed plugs see Table 4)

Fig. 12

Selected rock properties of the reservoir samples plotted against reservoir depth. The pre-caldera lavas were classified with respect to their degree of hydrothermal alteration (weak, moderate, and strong)

Fig. 13

Cross plots of matrix permeability, porosity, P-wave velocity, magnetic susceptibility, and thermal conductivity analyzed on the reservoir samples with respect to their alteration facies in (a–d) and the chemical index of weathering (CIW, here excluding marbles and ignimbrites) in (e–h)

Table 4 Petrophysical rock properties of the outcropping andesite units and reservoir core samples of the Los Humeros geothermal field

The Teziutlán and Cuyoaco andesites predominantly exhibit low matrix porosity (< 5%) and matrix permeability (< 10^–17 m²), with these two parameters showing a very weak correlation trend (Fig. 11a). Thereby, fractures observed in some plugs result in increased permeability values of up to 10^–13 m². The vesicular Teziutlán andesite lavas form the exception and show an increased matrix porosity ranging between 10 and 30%, along with permeability values ranging between 10^–16 and 10^–13 m². Compared to the outcrop samples, the reservoir samples show an increased average matrix porosity and permeability (~ 10^–16 m2, Fig. 11b). Matrix porosities of the basaltic andesitic to rhyolitic lavas vary between 5 and 25%. The ignimbrite samples show the highest porosity (~ 33%), but a comparably low matrix permeability (5·10^–18 m²). The marbles show the opposite trend, exhibiting the highest matrix permeability (10^–14 to 10^–15 m²), but a very low matrix porosity of approximately 1%.

Thermal conductivity and thermal diffusivity of the outcrop samples varies between 0.69 and 1.87 W m⁻¹ K⁻¹ and 0.45 and 1.1·10^–6 m² s⁻¹, respectively (Figs. 10a, 18). Both parameters are positively correlated (Fig. 18) and decrease with matrix porosity (Fig. 11e). Although the reservoir core samples exhibit a higher matrix porosity and permeability as well as lower bulk density, they show a significantly increased average thermal conductivity and diffusivity compared to the outcrop samples (Fig. 11d and f, Table 4). With 0.67 W m⁻¹ K⁻¹ the ignimbrites show the lowest thermal conductivity, while some skarns and silicified lavas exhibit increased thermal conductivities of up to 3.27 W m⁻¹ K⁻¹ (sample H28-2).

Under saturated conditions, thermal conductivity and thermal diffusivity of the reservoir samples increase by about 20% and 50% in average, respectively (Table 4, Fig. 18). For samples with low porosity (< 5%), saturated thermal conductivity and thermal diffusivity are consistent with measurements performed under dry conditions, while both parameters increase for samples with an intermediate to high porosity (10 to 33%) and/or permeability (e.g., brecciated or fractured lavas and ignimbrites). Likewise, saturated thermal conductivity and thermal diffusivity of the outcrop samples is on average ~ 13 and ~ 30% higher compared to the dry measurements showing a similar trend. However, the correlation of both parameters under saturated conditions shows significant degree of scatter (Fig. 18).

Specific heat capacity shows comparatively small variations throughout the dataset and ranges between 752 and 807 J kg⁻¹ K⁻¹, whereby the ignimbrites, andesitic lavas, and marbles from the reservoir cores show the highest values. The volumetric heat capacity is controlled by the respective bulk density of the samples. Therefore, the porous ignimbrite and ignimbrite–andesite breccias feature the lowest volumetric heat capacity and the marbles feature the highest volumetric heat capacity.

P-wave and S-wave velocity of the outcrop samples range between 1510 and 6276 ms⁻¹ as well as 880 and 3730 ms⁻¹, respectively (Fig. 10c) and are positively correlated (Fig. 18). Thereby, samples that contain abundant fractures (e.g., RLM9, LH40-LH72), plot in the same range as the porous lavas (LH8-12). The sonic wave velocities are positively correlated with thermal conductivity, thermal diffusivity, and density, but decrease with increasing porosity. Similar to the observations made for thermal conductivity, both parameters show no correlation with matrix permeability. Compared to the outcrop samples, the reservoir samples show slightly lower average P-wave and S-wave velocities (Table 4). It is notable that marbles and ignimbrites, which are usually plotting opposing to each other (Figs. 10b, 11b, d, f), feature sonic wave velocities in the same order of magnitude (Fig. 10d). This phenomenon is most likely caused due to the high porosity of these ignimbrite samples and the abundance of fractures observed on the marbles, which both reduce propagation of sound waves. Saturated P-wave and S-wave velocities increased by a factor of up to 1.35 and 1.1, respectively, compared to the measurements under dry conditions. Contrary to the thermal properties, they maintain a positive linear trend (Fig. 18).

Most of the outcrop samples feature magnetic susceptibilities between 2 and 16·10^–3 SI. Only few lavas show reduced values of about 0.3 to 0.9·10^–3 SI indicating weathering/alteration (Figs. 3j, 10e). Furthermore, the results indicate mineralogical differences between the sampling locations and units. On average, the Teziutlán andesites feature higher magnetic susceptibilities (3–16·10^–3 SI) compared to the Cuyoaco andesites (2–4·10^–3 SI, Table 4). Opposite, the reservoir samples show a wider range and a high variation. As shown in Fig. 10f, the samples can be grouped in two clusters. The first one comprises ignimbrite, breccias, and lavas affected by argillic as well as propylitic alteration with a weak to strong alteration intensity. These samples exhibit magnetic susceptibilities from 1.5 up to 18.05·10^–3 SI, whereby the ignimbrites and rhyolitic lavas contain lower values than the remaining lavas. The second population contains magnetic susceptibilities below 0.45 to 0.1·10^–3 SI (below 0 in Fig. 10f) and comprises lavas and breccias with a predominantly strong hydrothermal overprint as well as marbles, skarns, and silicified lavas. In general, limestones and quartz-rich rocks are known to show a diamagnetic behavior (Hrouda et al., 2009). However, while the marbles from core H13-3 show negative values (diamagnetic), the other marbles show low, but positive values of about 0.01 to 0.63·10^–3 SI (paramagnetic).

Particle and bulk density, thermal conductivity and diffusivity, as well as volumetric heat capacity gradually increase with reservoir depth (Fig. 12; from top to base particle density: 2.2 to 2.82 g cm⁻³, bulk density: 1.62 to 2.78 g cm⁻³, thermal conductivity: 0.5 to 2.96 W m⁻¹ K⁻¹, and VHC: 1272 to 2146 J m³ K⁻¹). Likewise, sonic wave velocity increases with reservoir depth but shows a higher variation in the shallow and deep reservoir sections. Since porosity is negatively correlated with the aforementioned parameters, it gradually decreases with reservoir depth from > 30% to ~ 1% (Fig. 12c). However, low porosity lavas (< 5%) occur at several depth levels. Matrix permeability, specific heat capacity, and magnetic susceptibility indicate no correlation with reservoir depth. Particularly, matrix permeability and magnetic susceptibility show a high variability at all depth levels.

A classification with respect to alteration facies or chemical indices does not reveal any clear relationships between the different investigated parameters (Fig. 13). In general, propylitic alteration seems to cause more homogeneous rock petrophysical properties than argillic alteration (Fig. 13a–d) when compared to outcrop samples (Figs. 10 and 11). The chemical indices again show no relevant correlation with the investigated parameters (Fig. 13e–h).

Discussion

Petrography and hydrothermal alteration

The investigated Teziutlán andesite unit comprises a suite of different basaltic andesitic to andesitic, porphyritic, massive and minor porous lavas with a dark gray to medium gray color. The Cuyoaco andesite unit comprises gray to slightly reddish-violet, porphyritic, massive and often fractured lavas with andesitic to dacitic composition, in agreement with the description of the Alseseca/Cuyoaco andesites presented in Yáñez and García, (1982) and Carrasco-Núñez et al. (1997). However, the Cuyoaco andesite unit is generally classified as the “Hornblende andesite” in literature (Viggiano and Robles 1988a, 1988b; Cedillo 1997, 2000; Arellano et al. 2003). Except for sample LH51 (Fig. 3), hornblende was not identified in the studied samples.

Hydrothermal alteration observed on the reservoir samples is often restricted to fractures and varies on the cm-scale, which leads to high heterogeneity in terms of sample appearance, chemical composition, and rock properties. Consequently, each collected sample has unique features and it is challenging to correlate between individual boreholes.

In previous studies, the reservoir samples retrieved from the shallower reservoir were classified as (Teziutlán) augite andesites (Cedillo 2000; Arellano et al. 2003). However, the reservoir samples retrieved from shallow to intermediate depth levels (~ 350–1400 m bgl) consist of different, very heterogeneous and often fractured mafic, trachyandesite, andesitic and rhyolitic lavas as well as andesite/ignimbrite breccias and ignimbrites. The samples are mainly affected by argillic to propylitic alteration with alteration intensities varying from weak to strong. These ignimbrites and ignimbrite/andesite breccias belong to the caldera group and are commonly interpreted as the cap rock of the geothermal reservoir (Cedillo 2000; Arellano et al. 2003).

Samples retrieved between 1500 and ~ 2500 m bgl consist of strongly altered often fractured or brecciated basaltic to andesitic lavas. Propylitic alteration with epidote, chlorite, K-feldspar, calcite, and quartz is dominant in these samples. According to González-Partida et al. (2022) deeper samples with predominantly clinozoisite, epidote, chlorite, along with smectite/illite, magnetite, biotite, sericite, pyrite, chalcopyrite, quartz, and calcite represent high-temperature propylitic alteration (~ 250– > 350 °C). The deepest samples collected in this study comprise skarns (> 2800 m bgl). Both, skarns and marbles, contain numerous fractures mainly filled with calcite and often contain high-temperature minerals, such as garnet, diopside, apatite, and/or wollastonite. However, high-temperature alteration minerals, such as garnet, diopside, and wollastonite, were also observed at shallower depth levels as shown in fracture fillings of sample H38-4 (1950 m) in the central to northern zones of the geothermal field (also referred to as the central collapse area). This observation is well in line with findings from (i) Martı́nez-Serrano (2002), who defined the beginning of a skarn zone at about 1800 m depth in the central collapse area, (ii) Jentsch et al. (2020), who reported the highest temperature anomalies at Loma Blanca within the central collapse area, as well as (iii) Urbani et al. (2021), who discussed existing borehole-temperature profiles in the same area. However, our findings show that high-temperature secondary mineral assemblages are not limited to the central collapse area as previously described in Prol-Ledesma and Browne (1989).

Advanced silicification along with a significant transformation of the phenocrysts and groundmass was observed at different depth levels and boreholes in the reservoir. González-Partida et al. (2022) defined a deep-seated silicic or high-temperature argillic alteration zone that occurs below ~ 2500 m depth (> 350 °C) with mainly quartz but also anhydrite, pyrophyllite, dickite, and white mica. However, partial or complete silicification might also occur at shallower depth levels as observed in sample H39-2 (1650 m bgl, Fig. 5). Silicification and the generally elevated amounts of quartz at all depths could be related to deep acidic supercritical fluids with high HCL, HF, and B concentrations with pH < 4 (Prol-Ledesma 1998; Flores-Armenta et al. 2010; Izquíerdo et al. 2011). Further sources of acidic fluids might be shallow intrusions (degassing of supercritical fluids), temperature/pressure changes, boiling, or reactivation of fractures (González-Partida et al. 2022). We interpret silicification at shallower depth levels within the propylitic alteration zone as well as the local occurrence of high-temperature minerals at different depth levels as possible evidence of shallow intrusions, in agreement with the shallowest geometry of the Los Humeros post-caldera stage magmatic plumbing system (Lucci et al. 2020) and with the results from field-constrained analog modeling proposed by Urbani et al. (2020, 2021). Shallow cryptodomes (e.g., Urbani et al. 2020, 2021) and/or magma pockets (Lucci et al. 2020) could be also responsible for thermal metamorphism in limestones and lavas that locally underwent transformation to marble, hornfels, and skarn, respectively. Furthermore, the occurrence of low and high-temperature alteration minerals at the same depth or even within one sample (Martı́nez-Serrano 2002; Prol-Ledesma 1998), plus the high variability of alteration facies and intensity in the upper section of the pre-caldera group, points to the possibility of several heterogeneous heat sources. This observation is well in line with temperature profiles and modeling of shallow magma pockets as presented in Deb et al. (2021), Lucci et al. (2020), Urbani et al. (2020, 2021), and Bonini et al. (2021). For example, Deb et al. (2021) predicted the occurrence of a mafic magma pocket in the central collapse area (close to well H3) at approximately 1.5 km depth and larger mafic to trachyandesite magma pockets close to the Xalapazco crater and well H18 at about 2.5–3 km depth. Thereby, associated mineral assemblages dominated by quartz form in close vicinity of the intrusions (Heřmanská et al. 2019) often causing a sealing effect (Martı́nez-Serrano, 2002; González-Partida et al. 2022).

The temporal evolution of the different hydrothermal processes was not investigated in this study but is described for selected samples in Rochelle et al. (2021) and is discussed in Gonzalez-Partida et al. (2022). Martı́nez-Serrano (2002) and Prol-Ledesma (1998) proposed at least two different stages of hydrothermal activity due to the coexistence of high- and low-temperature minerals at shallow depths (high-temperature hydrothermal activity followed by a low-temperature stage), which was most likely caused by permeability changes in the reservoir and reactivation of faults and fracturing over time. Gonzalez-Partida et al. (2022) described preferential fluid pathways with different fluid sources which might have significantly shifted over time in some areas of the reservoir. However, fragments of altered lavas, marble, and skarn found in some andesitic lavas and in the Xáltipan ignimbrite suggest that hydrothermal alteration of the deeper reservoir, skarn metasomatism, as well as marble-contact metamorphism already existed prior to the collapse of the Los Humeros caldera and has been influenced since then by the volcanic activity of the caldera and post-caldera group. Multiple temporal and spatial alteration events could explain the high variability in terms of sample appearance, alteration intensity, alteration facies, and changes in chemical composition in our sample set.

Contrary to previous studies (Cedillo 2000; Arellano et al. 2003), a thick tuff layer, which separates the pre-caldera andesites within the reservoir into an upper and lower andesite unit, was not identified. In these previous works, the deeper andesite sections were interpreted as equivalent to the Cuyoaco andesites and described as “Hornblende andesites.” Besides a few remnants of possibly amphibole in two samples, we could not identify hornblende in the reservoir samples. However, despite a few studies and short reports from the early exploration phase of the Los Humeros geothermal field in the late 70s and early 80s (e.g., Yáñez and García 1982; Ferriz 1982; Ferriz and Mahood 1984), not much attention was paid to the pre-caldera andesites in the study area. Further detailed petrological investigations and dating are needed to better understand their temporal and spatial evolution.

Bulk geochemistry

The reservoir samples are characterized by significant compositional variations compared to the outcropping Teziutlán andesites (Additional file 1: Table S4). Calculated compositional variations for the different facies of hydrothermal alteration are shown in Fig. 14a. The SiO₂ content shows variations of (i) ± 5% in samples affected by argillic and propylitic alteration and (ii) up to + 33% in samples that underwent advanced silicification. All the studied reservoir samples are characterized by a strong enrichment in Fe₂O₃. In general, the reservoir samples are characterized by highly variable trends, from enrichment to depletion, of the other major and trace element contents (Additional file 1: Table S4).

Fig. 14

a Calculated concentration changes and standard deviation for major, trace, and rare-earth elements with respect to the different alteration facies. The concentration changes were calculated using the average chemical composition of the Teziutlán andesite unit. b Plot depth vs. Eu/Eu* after Carrasco-Núñez et al. (2017b) and plot ∑LREE/(Eu/Eu*) vs. ∑HREE showing clearly lower Eu/Eu* ratios for the rhyolitic lavas. d Plot Ti/SiO₂ vs. Eu/Eu* showing different groups for rhyolites, silicified samples, and the remaining lavas of the pre-caldera group

Different widely used chemical indices were then tested (Table 3, Figs. 5, 6, and 13) to search possible relationships between bulk compositional changes and petrophysical properties. In contrast with the existing literature (e.g., Pola et al. 2012; Lee et al. 2021) the applied chemical indices seem to not correlate with the observed degree of hydrothermal alteration. With respect to studied reservoir rocks, (i) the CIW, CIA, PIA, and AI indices are not discriminant between argillic and propylitic alteration facies, (ii) the CCPI index is lowered for argillic alteration while increased for propylitic alteration, (iii) pervasive silicification is characterized by higher CIW, CIA, PIA, and AI and lower CCPI values, and (iv) skarns are characterized by lower values for all the applied indices.

When the REE budget is investigated it is possible to observe that the outcropping Teziutlán and Cuyoaco andesite units show comparable fractionated REE patterns with LREE enrichment and gently sloping profiles with Teziutlán andesites generally slightly higher in REE contents than those of Cuyoaco andesites. To note, many samples from both Teziutlán and Cuyoaco andesite units are characterized by a significant Sm_N (chondrite-normalized) positive anomaly that can be interpreted as the response of localized interactions with carbonatic groundwaters (Johannesson and Xiaoping 1997) or with carbonatic weathering materials (Savichev and Vodyanitskii 2011), such as the underlying Mesozoic limestones.

The studied reservoir samples show comparable fractionated REE profiles to those of the Teziutlán andesites except for some samples characterized by argillic and propylitic alterations and showing a Ce_N positive anomaly (Fig. 9, Additional file 1: Table S1). Existing literature (e.g., Braun et al. 1993; Jayananda et al. 2016) have widely demonstrated that Ce anomalies are the consequence of oxidation processes of Ce + 3 to Ce + 4 and following precipitation of the latter one in CeO₂, with negative anomalies related to removal of Ce by fluids and positive anomalies related to (i) Ce precipitation from the fluid phase in oxidizing conditions (Jayananda et al. 2016) or (ii) formation of Ce-bearing secondary minerals (Fodor et al. 1994).

Few reservoir samples, such as the rhyolite H19-1 and the shallow-depth cutting from H20 and H43 wells, show a well-developed Eu negative anomaly (Figs.  9b, 14b), that is typical of differentiated felsic magmas (e.g., Lucci et al. 2018). Similar to the observation from Carrasco-Núñez et al. (2017b), the Eu anomaly combined with Ti/SiO₂ ratio allows to discriminate the rhyolites from the andesites and from skarns and silicified lavas (Fig. 14d).

A cross-correlation analysis of the REE budget (Additional file 1: Table S3) revealed linear trends for the HREE (r = 0.90–0.99) and to a lesser extent for the LREE (r = 0.77–0.99) suggesting a general immobile behavior of the HREE and a possible minor mobility of LREE during hydrothermal alteration, in agreement with the work of Whitford et al. (1988) and Morgan and Wandless (1980).

Studies on compositional changes during hydrothermal processes indicate that elemental mobility could change under different alteration conditions and are strongly dependent by reservoir conditions and original unaltered protolith (Yongliang and Yusheng 1991). In Los Azufres geothermal field (Mexico), Torres-Alvarado et al. (2010) documented basaltic material at depth of ca. 2700 m affected by propylitic alteration and silicification and characterized by relevant loss of all major elements, and also of Nb, Zr, Hf, P, Th, and REE. On the contrary, Pandarinath et al. (2008), working on shallower (< 450 m depth) rhyolites from the same geothermal field, documented a general immobile REE behavior. Moreover, Pandarinath et al. (2020) again documented inert behavior of REE in intensively and hydrothermally altered (i.e., argillic and silicic alteration) outcropping dacitic lavas from the Acoculco geothermal field (Mexico). During the last half century, it was proposed that (i) REE behave highly selective in different alteration facies (e.g., Paraspoor et al. 2009), that (ii) epidote s.l. could play a relevant role in reallocating and controlling REE concentrations (e.g., Palacios et al. 1986; Torres-Alvarado et al. 2010), (iii) and important REE loss during alteration dominated by a single SiO₂ phase, such as opaline silica and quartz (e.g., Hopf 1993), and REE mobility controlled by REE concentration in reservoir fluids and the concomitant availability of complexing ions, such as CO₃²⁻, F⁻, and PO₄³⁻, and in acidic conditions Cl⁻ and SO₄²⁻ (Hikov 2011; Yongliang and Yusheng 1991).

However, the here-developed statistical analysis showed no significant differences between the reservoir samples and the Teziutlán andesites (see Kruskal–Wallis tests in Additional file 1: Table S5), confirming the main accepted immobile behavior of REE and their capability to maintain the signature of the unaltered/unmetamorphosed magmatic precursor (e.g., Buchs et al. 2013; Rossetti et al. 2017). In this view, the REE budget variations in the reservoir samples can therefore be interpreted as compositional heterogeneities of the volcanic material. Further detailed investigations integrating the quantification of REE contents in the reservoir fluids and in primary and secondary mineral phases are therefore required to fully understand the elemental mobility during hydrothermal alteration in active geothermal settings.

To conclude, a general elemental mobility was identified for the reservoir samples making unlikely the application of classic discrimination methods (e.g., Graf 1977; MacGeehan et al. 1980; Winchester and Floyd 1977; Floyd and Winchester 1978; Cullers and Graf 1984; Rollinson 1993; McLean and Barrett 1993). Furthermore, also the Nb/La vs. Sr/Nb method proposed by Carrasco-Núñez et al. (2017b) appears unrealistic plotting silicified lavas and skarns in the field of rhyolites. However, a statistical analysis shows significant differences in REE contents between the reservoir samples and the Cuyoaco andesites (Fig. 14c; see Kruskal–Wallis tests in Additional file 1: Table S5). Furthermore, integrating the sample depth and considering that shallow samples generally show ∑HREE ~ 15, the majority of the investigated reservoir core samples can be assigned to the Teziutlán andesites.

Petrophysical properties—outcrop vs. reservoir samples

The petrophysical data reflect not only the geological heterogeneity and the impact of hydrothermal alteration observed on the reservoir samples but also minor mineralogical differences between the outcrop samples as well as the influence of fractures. Comparison of both reservoir and outcrop samples reveals significant changes in the petrophysical characteristics.

Most of the lavas collected in the outcrops are characterized by an overall low average matrix porosity (< 4%) and matrix permeability (< 10^–18 m² to 10^–16 m²) as well as an intermediate thermal conductivity (~ 1.48 W m⁻¹ K⁻¹), thermal diffusivity (~ 0.84 10^–6 m² s⁻¹), and sonic wave velocities (V_P: ~ 5400 m s⁻¹, V_S: ~ 2450 m s⁻¹). The weak correlation between matrix porosity and permeability implies that fluid flow is predominantly fracture controlled, which has been confirmed by Lelli et al. (2021). The vesicular lavas of the Teziutlán andesite unit collected from outcrops north and northeast of the Los Humeros caldera feature not only higher matrix porosities (~ 15%) and permeabilities (10^–16 m² to 10^–14 m²) but also overall lower bulk densities, thermal conductivities, and sonic wave velocities. Specific heat capacity shows no distinct differences, while magnetic susceptibility shows a high variability and reveals mineralogical differences between the sampling locations.

In contrast to the outcrop samples, the reservoir samples show a higher parameter range and scattering (e.g., porosity: < 2% to > 30%, permeability: 10^–18 m² to 10^–14 m², and thermal conductivity: 0.7 W m⁻¹ K⁻¹ to 2.7 W m⁻¹ K⁻¹). Results of the marbles and ignimbrites are well in line with measurements performed on outcrop samples collected from the Xáltipan ignimbrite and marbles from Las Minas (Weydt et al. 2022). However, the ignimbrite core samples investigated in this study only represent the unwelded to partially welded upper section of the Xáltipan ignimbrite. According to Cavazos-Álvarez et al. (2020), the Xáltipan ignimbrite shows an increased degree of welding with reservoir depth alongside a porosity reduction. Therefore, porosity of the highly welded basal sections can be smaller than < 5% (Weydt et al. 2022; Cavazos-Álvarez et al. 2020).

Taking only the lavas of the pre-caldera group into account, the reservoir samples show an increased average matrix porosity, permeability, thermal conductivity, thermal diffusivity, and heat capacity and thus enhanced reservoir quality, but reduced average sonic wave velocities and magnetic susceptibility (Table 4, Fig. 15). However, both datasets reveal a bi- to multimodal distribution for the majority of parameters investigated (Fig. 15). A distinct shift in the density functions can be observed for bulk density, porosity, permeability, magnetic susceptibility, and to a lesser extent for thermal diffusivity and thermal conductivity.

Fig. 15

Kernel probability density functions of selected properties of the pre-caldera outcrop samples (black line) and reservoir samples (gray dashed line; excluding marbles and ignimbrites). N, number of analyzed plugs; vertical lines, mean value

Results from sonic well logging (well H42) conducted within the Los Humeros geothermal field show a similar porosity distribution for the andesitic lavas than the outcrop samples (Deb et al. 2019) with two distinctive groups: massive lavas with porosities < 5% and porous/fractured lavas with porosities ranging from ~ 10 to 30%. Thereby, the average porosity of the logged interval is 9 ± 6% and the “high porosity sections” were interpreted as fracture zones that occur at different depth levels (Lorenzo-Pulido 2008; Deb et al. 2019). Deb et al. (2019) concluded that about 30% of the andesitic lavas might contain these highly porous/fractured sections, while the remaining sections are massive andesites with a porosity of 3% to 4%. In our sample set, the observed increased matrix porosity is caused by mineral dissolution and precipitation predominantly along fractures modifying the initial pore network (Figs. 4, 5, 6, 16, and 17). Results of X-ray microCT analyses performed on seven andesitic samples retrieved from well H20 show the importance of microporosity as a result of hydrothermal alteration (Cid et al. 2021). Thereby, macropores account for an average total rock volume fraction of 3.7 ± 0.9%, whereas microporosity (< 1 µm in diameter) is higher with a rock volume fraction of 6.2 ± 0.6%. The round to elongated macropores in the samples are often isolated, suggesting that micropores, which often form a network near larger grains/structures, provide an important link between larger pores. According to Cid et al. (2021) microporosity enhances the matrix permeability by approximately 80%. With simulated porosities of 6.6% up to 14.5% and permeabilities of 10^–17 m² up to 10^–15 m², the calculations are in good agreement with our results. The enhanced matrix porosity and permeability (Fig. 15, Figs. 4, 5, 6, and 17) and the change in the porosity–permeability relationship in the majority of the hydrothermally altered reservoir samples becomes important when upscaling hydraulic permeameters to reservoir scale in numerical models. Comparable to the work of Heap and Kennedy (2016), the porosity–permeability relationships of the andesitic lavas cannot be described with one linear trend (Fig. 11a and b). Farquharson et al. (2015) identified a critical porosity threshold where fluid flow changes from crack to pore dominated at about 14% in andesitic lavas. Likewise, Kushnir et al. (2016) identified a changepoint at 10.5% in basaltic andesitic lavas. For future reservoir simulations of the Los Humeros geothermal field, these changepoints in the parameter relationships need to be considered. Additionally, the reservoir samples were predominantly retrieved from the hydrothermally altered and fractured damage zones of larger faults within the geothermal field (Fig. 2). For a more precise simulation of fluid flow in a local reservoir model, it would be essential to correctly estimate the lateral extension of the fault damage zone and the hydrothermal aureoles and to distinguish between unaltered, low-porous, and hydrothermally altered and high-porous reservoir sections.

The outcrop samples show a negative correlation of thermal conductivity with porosity. Thus, an increased porosity would expect a reduced matrix thermal conductivity. Although thermal conductivity of the reservoir samples is also negatively correlated with matrix porosity, thermal conductivity is significantly increased contrary to expectations. One plausible explanation is the enrichment of calcite and quartz in the system, which can be observed most clearly in the skarns showing the highest thermal conductivities and lavas that underwent pervasive silicification. After Clauser and Huenges (1995; and references herein), calcite (~ 3.6 W m⁻¹ K⁻¹) and quartz (~ 7.7 W m⁻¹ K⁻¹), together with less abundant alteration minerals, like pyrite (~ 19.2 W m⁻¹ K⁻¹), magnetite (~ 5.1 W m⁻¹ K⁻¹), hematite (~ 11.3 W m⁻¹ K⁻¹), and garnet (~ 5.5 W m⁻¹ K⁻¹, grossular), have higher thermal conductivities compared to feldspars (~ 2.0–2.7 W m⁻¹ K⁻¹). Likewise, the presence of alteration minerals, like epidote (~ 2.5–3.1 W m⁻¹ K⁻¹) and chlorite (~ 3.1–5.3 W m⁻¹ K ¹), led to an increased bulk thermal conductivity of the rock matrix.

According to Hrouda et al. (2009) volcanic rocks are often strongly magnetic and the variation in their magnetic susceptibility can be used as a geological indicator, e.g., for mapping mineralogical heterogeneities, metamorphic zones, or to differentiate between units. Furthermore, magnetic susceptibility is a very sensitive indicator of hydrothermal alteration. The oxidation of (titano)magnetite to iron oxides like hematite can reduce the magnetic susceptibility by several orders of magnitude; thereby, the smaller the grains the faster the alteration (Hrouda et al. 2009). Magnetic susceptibility decreases in the following order: Ms-magnetite (~ 3,000·10^–3 SI) > Ms-maghemite (< 3,000–2,000 10^–3 SI) > Ms-hematite (1.3–7·10^–3 SI; Dunlop and Özdemir 2015, Hrouda et al. 2009). The outcrop samples indicate local variations of magnetic susceptibility in the study area. The older Cuyoaco andesites contain magnetic susceptibilities of ~ 2 to 4·10^–3 SI, whereas the Teziutlán andesites show a much wider range with magnetic susceptibilities of ~ 2 up to 16·10^–3 SI. Thereby, weathering/alteration observed on some samples of the Cuyoaco andesites decreases the susceptibility by about one order of magnitude. The reservoir samples can be grouped into (1) altered samples that still contain ferromagnetic minerals (10^–3–10^–2 SI); (2) altered samples exhibiting magnetic susceptibilities reduced by about one to two orders of magnitude compared to group 1 (10^–3–10^–4 SI); and (3) marbles, silicified lavas, and skarns showing a diamagnetic to paramagnetic behavior with magnetic susceptibilities that are up to three orders of magnitude smaller compared to group 1 (10^–4–10^–6 SI). The loss of magnetic signals in the hydrothermally altered rocks could be used to identify and map the hydrothermal aureoles and thus, prospective fluid pathways in the reservoir using geophysics (e.g., borehole logs, airborne magnetic surveys; Kanakiya et al. 2021).

The results of the petrophysical measurements performed on the reservoir samples are in good agreement with previous studies (Contreras et al. 1990; García-Gutiérrez and Contreras 2007). For example, matrix porosity and permeability data of mainly andesitic lavas and tuff of the Los Humeros geothermal field presented in Contreras et al. (1990) range between 2.5% and 23.3% (mean = 11.6%) and 10^–15 m² and 10^–17 m² (mean = 3.3·10^–16 m²), respectively. These data also show a high degree of scatter, similar to our observations. In general, both particle and bulk density increase with reservoir depth, and porosity decreases with reservoir depth, while matrix permeabilities show a high variability at all depths. Thus, a uniform nonpermeable layer within the pre-caldera group as described in Cedillo (2000) and Arellano et al. (2003) cannot be confirmed.

Martı́nez-Serrano (2002) and Prol-Ledesma (1998) reported self-sealing processes and the reduction of permeability due to hydrothermal alteration in the deeper part of the reservoir, while Izquíerdo et al. (2011) described the generation of millimeter-sized vugs caused by intensive leaching in silicified lavas (sample H26-4). A definitive conclusion about this is not possible at this stage, especially because of the high geological heterogeneity observed in our samples. In general terms, the precipitation of secondary phases (e.g., clay minerals, calcite, quartz, or epidote) in pores and along fractures reduces matrix porosity and permeability. However, fracturing or reactivation of fractures (thus facilitating further fluid–rock interactions) counteracts this effect. Self-sealing processes can be observed most clearly in the deepest part of the reservoir in the marble and skarn samples (porosities < 2%), which contain numerous fractures and fissures of different generations that are predominantly filled with calcite. Uniform positive density anomalies occur at about 1000 m above sea level in the central area of the Los Humeros geothermal field in density models (Cornejo 2020), which correlate with the skarn zone described in Martı́nez-Serrano (2002) confirming self-sealing processes as observed on the skarns and carbonates.

The comparison with literature data underlines the high geological variability of volcanic rocks (Table 7), which are predominantly controlled by matrix porosity and the occurrence of fractures (Mielke et al. 2015; Navelot et al. 2018). In general, hydrothermally altered andesitic lavas feature higher matrix porosities and permeabilities. Hence, alteration type, the intensity of alteration, matrix porosity, and density of fracturing control the rock parameters. The impact of hydrothermal alteration on the rock properties strongly depends on the protolith and its initial pore network (Villeneuve et al. 2019). Thus, hydrothermal alteration affects each rock unit in a different way leading to a high geological heterogeneity on the centimeter to tens-of-meter scale.

Villeneuve et al. (2019) reported that deeply buried lavas affected by propylitic alteration from different geothermal fields in New Zealand tend to have low porosity and low permeability and deform through microcracking and dilation, resulting in hydraulic active zones localized in fractured intervals. Although similar alteration types and intensities to those described in this study were observed in deep geothermal reservoirs in New Zealand (Siratovich et al. 2014, Wyering et al. 2014; Table 7), the andesitic reservoir rocks in Los Humeros show significantly higher matrix porosities and permeabilities (Table 4). Mordensky et al. (2019a, b) and Frolova et al. (2014) reported increased porosities and permeabilities in reservoir samples for argillic as well as propylitic alteration. According to Frolova et al. (2014) porosity decreases from low-temperature alteration (argillic, zeolitic; ~ 24–40%) to high-temperature alteration (propylitic, 5.6%) and increase again in secondary quartzites (12%; usually form during 300–550 °C, pH: 1–4) and quartz-feldspar metasomatites (20%). In contrast to the findings of the aforementioned studies, Heap and Kennedy (2016) reported significant smaller porosities and permeabilities in altered andesitic lavas compared to unaltered lavas collected from outcrops.

While changes in matrix porosity and permeability have been widely discussed (Frolova et al., 2014; Heap and Kennedy 2016; Mielke et al. 2016; Kushnir et al. 2016), only a few studies provide information on changes in thermal properties and the results seem to differ from study to study. Navelot et al. (2018) reported consistent thermal conductivities with increasing intensity of hydrothermal alteration, while Mielke et al. (2015) and Heap et al. (2020b) reported decreasing values. Furthermore, Heap et al. (2020b) concluded that specifically acid-sulfate alteration increases specific heat capacity and reduces thermal conductivity and diffusivity.

Navelot et al. (2018) reported higher median magnetic susceptibility values for slightly moderately altered lavas (11.97·10^–3 SI) than for fresh andesitic lavas (11.6·10^–3 SI) due to increased concentrations of iron oxides during the first phase of alteration of pyroxenes and plagioclase. However, a complete transformation of the initial minerals removes the magnetic minerals and reduces magnetic susceptibility about two orders of magnitude (0.09·10^–3 SI for advanced altered andesitic lavas). Likewise, Kanakiya et al. (2021) reported high remnant magnetization in altered lavas of the Whakaari geothermal field. This phenomenon can also be observed to some extent in our data (Fig. 15, Table 4), where some reservoir samples affected by propylitic and argillic alteration show slightly higher magnetic susceptibilities than the outcrop samples. However, more detailed investigations are required, which determine and quantify the magnetic minerals and their alteration products in our sample set in order to draw a final conclusion.

In contrast to previous studies (Pola et al. 2012; Frolova et al. 2014), the correlations between alteration intensity, chemical indices, and alteration facies show no clear trends (Figs. 8 and 13; Table 3), and thus, the indices prove not to be useful for predicting the petrophysical characteristics of the reservoir rocks. Propylitic alteration seems to cause less scattering of the parameters and deviations from the parameter trends as observed on the outcrop samples.

Data processing and implications for modeling the Los Humeros geothermal field

Rock property measurements of the outcrop and reservoir core samples provide matrix properties and thus, represent the centimeter to decimeter scale only (eventually with small-scale or single fractures). Analysis under standardized laboratory conditions is necessary to ensure the comparability of the results. Consequently, the data generated do not reflect in situ conditions, such as high reservoir temperatures, overburden pressures, confining pressures, and fluid properties at reservoir depths. Therefore, the data need to be corrected for reservoir conditions (Weydt et al., 2022) and (depending on the aim and scale of future applications) corrected to reservoir scale (upscaling). Several analytical and empirical relationships and correction functions have been identified in the past to convert the parameters to reservoir conditions (Farmer 2002; Rühaak et al. 2015; Ringrose and Bentley 2021). However, the majority of recent studies do not consider reservoir temperatures above 350 °C. Future research is needed to better estimate the in situ rock characteristics in super-hot geothermal reservoirs (including up to supercritical conditions, e.g., Kummerow et al. 2020; Heap et al. 2020a).

Since the beginning of the exploration of the Los Humeros geothermal field several conceptual geological models have been developed (Viggiano and Robles 1988a, b; Martı́nez-Serrano and Alibert 1994; Cedillo 1997, 2000; Arellano et al. 2003). Viggiano and Robles (1988a, b) classified the reservoir rocks according to the four main stratigraphic units, whereby the pre-caldera andesites were divided into augite andesites in the upper part (equivalent to the Teziutlán andesites) and hornblende andesites in the basal part (equivalent to the Alseseca and Cuyoaco andesites). Later on, Cedillo (2000) subdivided the reservoir rocks into ten units, which form mainly continuous layers in the subsurface. The incorporation of ignimbrites and tuffs forming an impermeable continuous layer between the augite and hornblende andesites led several authors to discuss whether the Los Humeros geothermal field hosts two different geothermal reservoirs or not (Izquíerdo et al. 2003, 2015; Gutiérrez-Negrín and Izquíerdo 2010). However, the revision of nine lithostratigraphic profiles (Carrasco-Núñez et al. 2017b) as well as the petrographic and petrophysical analyses presented in this study cannot confirm a sufficiently thick tuff layer, which would act as an aquitard within the andesitic reservoir. Furthermore, the hornblende andesites unit covers several hundred meters in the preliminary 3D model (Calcagno et al. 2020). Based on our findings, the reservoir samples can most likely be related to the Teziutlán andesites unit in the outcrops, which supports the findings of Carrasco-Núñez et al. (2017b) who concluded that the Cuyoaco andesites have a very limited extension and have only been identified in the bottom of a few boreholes. However, since both andesite units feature similar rock properties, it is not necessary to define separate units within a 3D numerical model.

According to Carrasco-Núñez et al. (2017b, 2018) the mafic to rhyolitic lavas identified between ~ 900 m and 1400 m bgl can be linked to the pre-caldera rhyolites (~ 270–693 ka), which has also been defined as the pre-caldera stage. However, considering the new findings presented in Lucci et al. (2020), Urbani et al. (2020, 2021) and Deb et al. (2021) some of these lavas could also possibly be related to the post-caldera stage magmatic activity. A distinct feature is their variable thickness along with no lateral continuity (Carrasco-Núñez et al. 2017b; Urbani et al., 2020, 2021; Cavazos-Álvarez et al. 2020). Based on our findings, together with additional borehole data provided by CFE, these lavas occur more frequently in the upper section of the pre-caldera group, and thus, this unit might have a much larger lateral extent as previously described (Cedillo 2000; Arellano et al. 2003; Carrasco-Núñez et al. 2017b, 2018). Due to the highly variable lateral extension, thickness, alteration intensity, hydraulic properties, and most likely fracture pattern of these lavas, it would be useful to define a separate model unit in local reservoir models with a high resolution (small grid size).

For the parametrization of reservoir models of the Los Humeros geothermal field, it has to be emphasized that the majority of the reservoir samples mainly represent the fractured and hydrothermally altered sections in close proximity to larger fractures and fault zones within the geothermal field. Data from the outcrop samples could be used to depict the unaltered and low-porous sections in the reservoir. For a large-scale regional model (with a large grid size), we suggest using data from outcrop samples in order to avoid an overestimation of matrix porosity and permeability.

Conclusions

This study presents a new dataset of petrophysical rock properties of mainly andesitic lavas from the Los Humeros geothermal field (Mexico), combined with petrographic and geochemical (XRF and ICP-MS) analyses. Hydrothermally altered borehole core samples were compared with stratigraphically equivalent outcrop samples to quantify the impact of hydrothermal alteration on the physiochemical characteristics of reservoir rocks in geothermal systems linked to active volcanic systems. The PLM and SEM observations were used to discriminate primary and secondary assemblages, alteration intensity, and facies. The bulk elemental modifications as well as selected chemical indices were discussed with the aim to highlight the possible relationship between chemical changes, alteration intensity, and rock properties.

In summary, the results showed the following:

• The investigated Cuyoaco and Teziutlán andesites comprise a variable suite of mainly unaltered, occasionally fractured, pyroxene-bearing andesitic to dacitic, as well as basaltic andesitic to andesitic lavas with a porphyritic to glomerophyric texture, respectively.

• From top to base, the reservoir samples represent shallow ignimbrites and ignimbrite/andesite breccias (~ 350–670 m bgl), followed by a variable suite of mafic to rhyolitic lavas between ~ 800 and 1400 m depth, basaltic to andesitic lavas between ~ 1500 and 2500 m, as well as skarns (~ 2900 m) and marble (pre-volcanic basement). The samples show a high geological variability prohibiting a clear correlation between the boreholes.

• Hydrothermal alteration observed on the reservoir samples predominantly occurs along cracks and fractures and is highly variable on a cm-scale. Argillic alteration was mainly observed at shallow to intermediate depths (to ~ 1400 m) followed by propylitic alteration and skarn. However, advanced silicification and high-temperature minerals were observed in samples at different depths, which could have resulted from separate, shallow intrusions.

• Significant bulk chemistry changes were identified in the reservoir samples, mainly in those affected by advanced silicification. Thus, standard discrimination methods are not applicable for our dataset. The Teziutlán and Cuyoaco andesites both show fractionated REE patterns with LREE enrichment and gently slope toward HREEs. The reservoir samples show REE budget and fractionated patterns strongly overlapping with those of Teziutlán andesites, while differences are identified with respect to the Cuyoaco andesites. The negative Eu anomaly observed in rhyolites and in shallow-depth cuttings, together with the TiO₂/SiO₂ ratio allows to discriminate felsic material in the reservoir. The general inert REE behavior during hydrothermal processes is highlighted.

• The Teziutlán and Cuyoaco andesites are predominantly characterized by low matrix porosities (< 4%) and permeabilities (< 10^–17 m²) as well as intermediate thermal conductivities, thermal diffusivities, and sonic wave velocities. Additionally, the Teziutlán andesites comprises a few porous basaltic andesitic lavas with intermediate to high matrix permeabilities. Bulk density, thermal conductivity, and sonic wave velocities are mainly controlled by matrix porosity. However, fractures significantly increase matrix permeabilities (up to 10^–13 m³) and reduce wave velocities. Specific heat capacity shows comparatively small variations, while magnetic susceptibility shows location specific trends.

• In contrast, the reservoir samples show enhanced hydraulic and thermal properties, but reduced bulk densities, sonic wave velocities, and magnetic susceptibility. In general, matrix porosity decreases with reservoir depth, while thermal conductivity and bulk density increase. Other parameters, such as matrix permeability and magnetic susceptibility, show a high degree of scatter at all depths. The significant loss of magnetic susceptibility in the reservoir samples could be helpful to identify hydrothermal aureoles as an indicator for possible fluid pathways during geophysical surveys. The correlation with alteration intensity, alteration facies, or chemical indices reveals no clear trends, which might be the result of multiple hydrothermal processes/stages.

Appendix

Additional graphs and data

See Tables 5, 6, 7

Table 5 Geographic coordinates of the outcrop samples

Table 6 Sample depth and lithological interpretation of the reservoir samples

Table 7 Overview of published rock property data with respect to hydrothermal alteration

See Figs. 16, 17, 18

Fig. 16

Microcapture photographs of selected reservoir samples showing the impact of hydrothermal alteration, brecciation, and fracturing at the cm-scale

Fig. 17

SEM images showing (a) amygdales filled with chlorite, epidote, and quartz, (b) alteration rims near larger feldspar phenocrysts, (c-d) epidote s. l. and quartz filling vugs (c) and veins (d), and pervasive silicification and relicts of phenocrysts in andesitic lavas. f metasomatic vein in fractured marble with wollastonite, garnet and apatite. Note the abundance of secondary porosity (black) in (c) to (f). Abbreviations: Adr = andradite, Ap = apatite, Cal = calcite, Ccp = chalcopyrite, Chl = chlorite, Ep = epidote, Fsp = feldspar, Grs = grossular, Kfs = K-feldspar, p = pores, Pl = plagioclase, Py = pyrite, Qz = quartz, Wol = wollastonite

Fig. 18

Cross-plots of thermal conductivity, thermal diffusivity, P-wave and S-wave velocities analyzed under dry and saturated conditions. See legend in Figs. 10 and 11