Abstract
Surface deformation due to fluid extraction can be detected by satellite-based geodetic sensors, providing important insights on subsurface geomechanical properties. In this study, we use Differential Interferometric Synthetic Aperture Radar (DInSAR) observations to measure ground deformation due to fluid extraction at the Los Humeros Geothermal Field (Puebla, Mexico). Our main goal is to reveal the pressure distribution in the reservoir and to identify reservoir compartmentalization, which can be important aspects for optimizing the production of the field. The result of the PS-InSAR (Persistent Scatterer by Synthetic Aperture Radar Interferometry) analysis shows that the subsidence at the LHGF was up to 8 mm/year between April 2003 and March 2007, which is small relative to the produced volume of . The subsidence pattern indicates that the geothermal field is controlled by sealing faults separating the reservoir into several blocks. To assess if this is the case, we relate surface movements with volume changes in the reservoir through analytical solutions for different types of nuclei of strain. We constrain our models with the movements of the PS points as target observations. Our models imply small volume changes in the reservoir, and the different nuclei of strain solutions differ only slightly. These findings suggest that the pressure within the reservoir is well supported and that reservoir recharge is taking place.
1. Introduction
A good understanding of reservoir processes and properties is crucial for optimizing subsurface operations. Uncertainty in future production and potential risks such as induced seismicity can be reduced using detailed information about the mechanical and hydraulic properties of the reservoir and its surroundings. These properties affect the ground surface movements induced by subsurface activities. Therefore, such information may be achieved by monitoring surface deformation. Relating surface deformation to the process of subsurface extraction or injection can give an indication of the spatial distribution and the amount of volume and pressure changes in the subsurface.
The deployment of DInSAR (Differential Interferometric Synthetic Aperture Radar) allows detecting small movements on the Earth’s surface (e.g., [1, 2]). It utilizes the phase difference between two SAR images to estimate displacement along the satellite line-of-sight (LOS). Time-series analysis facilitates the monitoring of gradual changes in ground movements by selecting stable point scatterers (Persistent Scatterers (PS)) in multiple interferograms, while decreasing atmospheric disturbances and improving the performance of phase unwrapping.
Numerous studies have applied the DInSAR technique and modeled the source of deformation associated with the exploitation of subsurface resources. These models most commonly relate surface deformation to subsurface extraction or injection processes through so-called influence functions that are based on analytical solutions for different nuclei of strain. For instance, Trugman et al. [3] and Samsonov et al. [4] used inflation point sources to model subsidence due to production at the world’s second largest geothermal field, Cerro Pietro. Surface subsidence at the Reykjanes geothermal field was modeled with point and ellipsoidal pressure sources [5]. Nuclei of strain solutions other than the inflation sources were also applied in several studies. As an example, Atefi Monfared and Rothenburg [6] employed Okada’s solution for expansion or contraction in one direction [7] as reservoir deformations are mostly one dimensional. Vasco et al. [8] modeled the horizontal opening of a vertically-oriented fracture due to CO2 injection using Okada’s solution.
In this paper, we present a case study of the Los Humeros Geothermal Field (LHGF). LHGF is among the largest geothermal fields in Mexico with an installed capacity of ~93.6 MW for an operational capacity of 68.6 MW and is operated by the national Mexican electrical company (Comisión Federal de Electricidad (CFE)). We show the first results of the PS-InSAR (Persistent Scatterer by Synthetic Aperture Radar Interferometry) time series over the LHGF. We analyzed 13 C-band Single-Look Complex (SLC) images acquired by the European Space Agency’s (ESA) Envisat satellite between April 2003 and March 2007 to detect surface movements due to field operations and potential discharge and recharge zones. We relate ground deformation with volume changes in the reservoir through analytical solutions for different nuclei of strain. The main goal of this study is to have a better understanding on the subsurface processes at the LHGF based on surface movements. Furthermore, we intend to reveal the pressure distribution within the reservoir and identify reservoir compartmentalization, which may help with future planning and optimization on the production of the geothermal field.
2. The Los Humeros Geothermal Field (LHGF)
The LHGF is a superhot geothermal system connected to the Los Humeros caldera complex located in the eastern sector of the Trans-Mexican Volcanic Belt (Figure 1). The basement of the LHGF is composed of granites and schists of Paleozoic age, overlain by a thick metamorphosed Mesozoic limestone succession. Volcanic activity in the area initiated during Miocene times (~10 Ma) and produced andesites that outcrop in the eastern part of the Los Humeros caldera (Figure 1) [9]. The LHGF is connected to a caldera system that has been active from 0.46 Ma until recent [9], although according to the most recent studies the Los Humeros caldera is considerably younger (0.16 Ma, [10]). The Los Humeros caldera complex was formed by at least two major rhyolitic eruptions and multiple minor to medium eruptions and lava flows. The duration of the caldera forming period was ~410 ky [9], or alternatively significantly shorter (94 ky, [10]). Recent andesitic and basaltic volcanism is poorly dated but considered to be <20-40 ky old [9].
Commercial exploitation of the LHGF started in the early 90s. Reinjection began few years after production, in 1995; however, injection rates are about an order of magnitude smaller than the produced volumes (Figure 2). Since the beginning of production, more than 50 wells have been drilled. The main production wells are located in the northern part of the LHGF. Average yearly production and injection rates during the period of the InSAR monitoring (2003–2007) are and , respectively (Figure 2). Precaldera (~10.5-1.5 Ma) andesites of low to medium permeability provide the main reservoir formation. The reservoir consists of several blocks separated by mostly normal faults (Figure 1) (e.g., [11, 12]). Two main fault systems can be distinguished, one with NE-SW to E-W striking structures (e.g., Las Papas) and a second one with NW-SE to NS orientation (e.g., Maxtaloya). These structures have major control on the geothermal system. Neotectonic deformation of the caldera floor is recorded based on field observations, associated with movements along the inner-caldera faults due to recent/active resurgence processes [11].
3. Interferometric Synthetic Aperture Radar Monitoring
We performed PS-InSAR time-series analysis to resolve ongoing ground deformation due to geothermal exploration at the Los Humeros Geothermal Field. We analyzed SLC images (descending track 212) acquired in C-band by ESA’s Envisat satellite between April 2003 and March 2007 (Table 1). We selected a subset of covering the Los Humeros caldera system for the time-series processing.
The selection of the master image follows the criteria of minimizing the perpendicular and temporal baselines; therefore, we coregistered the images to the master geometry. We used DORIS Precise orbits provided by the ESA for the orbital correction, and the topographic phase was computed and removed from the interferograms using the 1-arc-second (30 m) resolution Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM). For the time-series processing, we discarded the interferograms with no visible coherence and perpendicular baselines above 500 m. Interferograms were geocoded and imported into the MATLAB-based StaMPS (Stanford Method for Persistent Scatterers) framework [13]. The PS processing in StaMPS started with an initial selection of PS candidates based on the amplitude dispersion of the pixels [14]. The number of PS candidates was further reduced in an iterative process. We correct the interferograms for the spatially correlated and uncorrelated look angle phases, orbital errors, and master atmosphere. Spatially uncorrelated errors are proportional to the perpendicular baselines, and their phases are estimated during the PS selection and removed from the wrapped phase prior to unwrapping. All the others are estimated after unwrapping and then removed before performing unwrapping again. Even though we are using a 30 m SRTM DEM, it has been shown that the influence of DEM on the deformation rate calculation is small and inversely proportional to the wavelength [15]. In Figure 3(c), we observe that the maximum DEM error is observed outside the area of interest reaching ~0.01 rad/m, which corresponds to ~12 m of DEM error for Envisat.
The individual interferograms are generally of good quality, fairly coherent even over long time scales. On the other hand, dense vegetation especially in the northern part of the study area decreases the interferometric coherence. Most parts of the study area are located at high altitudes (above 2500 m) with strong relief. As a result, some interferograms, and consequently the mean line-of-sight (LOS) velocities, seem to be influenced by topography-related tropospheric phase delays (e.g., [1, 16]) (Figures 3(a) and 3(b), 4 and 5). We corrected the topography-related tropospheric artifacts based on a linear relationship between phase and elevation using the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [17]. We selected a region outside the geothermal field (Figures 3(a) and 5(a) black box) to estimate the tropospheric contribution and removed them from the interferograms. The correlation between elevation and unwrapped phase was identified on many, but not on all interferograms (Figure 4). The correction was performed for each interferogram individually.
The mean LOS velocities over the Los Humeros Geothermal Field vary between -2 and -8 mm/year (Figure 6), where negative values indicate movements away from the satellite. The largest subsidence was observed in the northeastern part of the field, east of the two injection wells operating during the period of the InSAR. The largest subsidence was observed in the northeastern part of the field, east of the two injection wells operating during the period of the processed InSAR. Further explanation of the estimated subsidence is addressed in Discussion and Conclusions.
(a)
(b)
(c)
4. Subsidence Modeling
Subsurface extraction or injection processes induce movements on the ground surface with a variable amplitude primarily through poro-elastic coupling. Other factors such as subsurface induced pressures compared to the in situ stress field, depth of the reservoir, properties of the reservoir rock, and over- and underburden media can also influence the observed surface movement pattern. Volume changes in the subsurface due to, for instance, the extraction of fluids result in pressure reduction leading to compaction strain, which is a function of elastic strength or compaction coefficient. Furthermore, injection processes can generate thermo-poro-elastic or thermo-poro-elasto-plastic deformation due to the temperature difference between the injected fluid and reservoir. Since cold fluids are injected into the reservoir, temperature near the injection wells drops, leading to thermal contraction. Our models only account for expansion due to the injected volumes, resulting in uplift on the surface. Therefore, there might be a biased estimation of surface movements close to the injectors due to the effect of thermal contraction that we do not take into account.
We related subsidence with fluid extraction from the geothermal field through influence functions like fast-forward models. Such models most commonly apply nuclei of strain solutions based on a center of compression (equivalent to an inflation point source), although these functions may show a mismatch with measurements (e.g., [19]). Fokker and Osinga [20] have shown that the geometry of the reservoir and the elasticity contrast with its surroundings have a major influence on the shape of the subsidence bowl. Therefore, we tested different strain nuclei to construct the influence functions.
We compared the movements of the PS (for the period of one year) as observations to the responses of a number of solution scenarios. We used nuclei of strain solutions for point sources in an elastic half-space based on a center of compression [21] (model 1, model 3) and a tensile, horizontally oriented fault (opening in vertical direction) [7] (model 2, model 4). In the case of the center of compression, we used the equation after Segall [22] to calculate the displacement vector at the surface. For a point source located at (, , ) with the associated volume change , the surface displacements , , are written as follows: where is the distance from the source to the surface observation point and is the Poisson’s ratio. Equation (1) assumes a relationship between pressure and volume change: where is the pressure change associated with the source, is the shear modulus, and is the radius of the source (, i.e., point source).
The formula for the tensile sources is more complicated. We followed Okada [23] to calculate the displacement vectors at the surface. The equation contains the same parameters as the one for the center of compression, including the coordinates describing the location of the source, the Poisson’s ratio of the elastic half-space, and the volume change associated with the tensile dislocation. We constructed our forward models using the MATLAB-based dislocation models after Nikkhoo et al. [24]. In addition to the point sources, we employed these solutions for finite rectangular sources. Finite rectangular sources are equivalent to the integration of the displacements attributed to the point sources through a rectangle [23]. We tested two configurations for the modeling: first we performed the fitting along a section, then we modeled the northern part of the field area, where a clear subsidence signal was observed. We tested numerous cross sections and found the most suitable having a sufficient number of PS (Figures 6(a) and 7). We included PS data within 130 m distance from the section as target observations for the subsidence modeling, yielding 81 PS in total. We placed three sources according to the location of the production and the injection wells. We located the sources at the main extraction/injection depth of 1650 m. The layout of the sources is illustrated in Figure 8. The easternmost and westernmost sources employ contraction due to production, while the middle one uses expansion as a result of injection. The centers of the sources are located 700 m and 800 m from each other.
The layout of the modeled area covering the northern part of the geothermal field area is shown in Figure 9. Based on the locations of the wells and the subsidence pattern, we used 7 sources in total: 6 contraction sources and 1 expansion source placed between the two injectors (Figure 9(b)). We employed nuclei of strain solutions for inflation/deflation point sources (model 3) and for tensile point sources opening in a vertical direction (model 4), equivalent to the sources used for model 1 and model 2, respectively.
We selected a Poisson’s ratio of the elastic half-space to 0.265 based on laboratory experiments on rock samples from the field. Since we simplified the subsurface to an elastic half-space, the geological units with different elastic properties are not distinguished. Considering that the Poisson’s ratios of the volcanic rocks at the LHGF are almost identical (±0.002), we found this approximation acceptable.
The sources, be they a center of compression or a tensile source, result in an extended movement signal at the surface—with the distribution of movement due to a tensile source more concentrated above it. As a result, patterns in the source, like a steep gradient in the pressure, do not propagate to the same pattern at the surface but are smoothed. Thus, as sources at different positions influence the surface movement data, inverse modeling or simultaneous parameter estimation is required.
We fitted the influence functions to the PS observations using a weighted least-squares approach. The fitting parameters were the volume changes attributed to the sources. We minimized the weighted sum of the residuals using the standard deviations of the LOS velocity estimates to construct the weighting factors. The best-fitting volume changes for the different models are listed in Table 2. We report the weighted RMS of each model to represent the goodness of the fit. In the case of the inflation/deflation sources along the section (model 1), the volume change of the westernmost source was negligible. For the tensile sources along the section (model 2), the volume changes were in the same order of magnitude. The total amount of volume changes attributed to the 7 sources is approximately two times larger for model 3 than for model 4. The RMS values in the case of all models are below 1 (Table 2), suggesting that the fits are reasonably good.
Figure 8 shows the responses derived from the strain nuclei placed at 1650 m depth along the cross section. A good fit was achieved by the tensile point source (model 2). The misfits are a bit larger for the inflation sources (model 1), and the influence function is wider and also slightly less deep. Furthermore, the strength of the westernmost source is negligible in contrast with the production data. We only show the models using point sources because our modeling results of the point and rectangular sources are almost identical. The error estimates around the InSAR observations correspond to the standard deviation of the mean velocity estimates. The standard deviation from the mean velocity was estimated from the average velocity with respect to the mean of all pixels without a selection of a reference point. The reasoning behind not selecting a reference point is because we have no external information on which area is eventually stable. The maximum standard deviation estimated in our area of interest is ~1.5 mm/year.
Figure 9 illustrates the modeling results of the northern part of the field. In contrast with the results through the section, the inflation point sources (model 3) are fitting better than the vertical tensile sources (model 4). The extent and geometry of the subsiding area only slightly differ (Figures 9(b) and 9(e)).
5. Discussion and Conclusions
We processed 13 C-band radar images acquired by the ESA’s Envisat satellite between April 2003 and March 2007, to monitor surface movements at the LHGF. The geothermal field is located at high altitudes inside the Los Humeros caldera. Therefore, it was necessary to perform corrections for artifacts due to topography-related phase delays on the interferograms.
The PS-InSAR analysis shows that the LHGF is characterized by up to 8 mm/year of subsidence during the period of the InSAR monitoring (Figure 6). We attribute this deformation to field operations given that the estimated maximum subsidence is located around the production wells. Given the concentration of deformation around the wells, the estimated source depth and the short wavelength deformation pattern, we think that the estimated deformation is only due to geothermal production and not volcano-tectonic related. The area of maximum subsidence is relatively small, located at the northern part of the geothermal field. This area appears to be isolated from the injection wells that were operational during the period of the InSAR analysis. This isolation is supported by the epicenters of the induced earthquakes [18], most of which are located west from the injectors, suggesting that the majority of the injected fluids are directed westwards (Figure 6(b)). No clear subsidence signal is observed in the southern and western part of the field, although large numbers of production wells have been drilled in these areas. This indicates a significant pressure support that might originate from deep recharge and partly from the injected volumes. Pressures cannot be supported only by the injectors, since the injected volumes are about ten times smaller than the produced volumes. Additionally, due to the low PS density of the northwestern part of the geothermal field and the lack of pressure measurements, we cannot be sure about the contribution of the injectors to the pressure support. Pressure drop in the reservoir is most likely restricted below the clearly subsiding area. As a result, the pressure distribution does not deviate much from the virgin one. As regular pressure measurements from the wells are not available to support these findings, this conclusion demonstrates the versatility of the use of surface movement data combined with inversion.
The procedure that we have employed propagates the measured surface movement directly to the driving parameters, i.e., the pressure magnitude and distribution. As such, it bypasses intermediate results like patterns of, e.g., a steeper edge of the subsidence bowl over a sealing fault. Such patterns would indeed be very hard to identify without prior knowledge about the pressure distribution due to the combination of the small magnitude of the signals and the smooth signature at the surface of a steep pressure gradient at depth.
We modeled the subsidence using influence functions based on different nuclei of strain solutions. Our model results, combined with the deformation estimation from InSAR as observations along a cross section, show that the vertical tensile point sources having deeper and steeper responses fit better with the observed subsidence pattern. However, when the entire field was used, the inflation point sources performed slightly better. Atefi Monfared and Rothenburg [6] suggested that Okada’s solution for expansion/contraction in the vertical direction is more applicable to model reservoir deformations, as they are laterally more extensive. Fokker and Osinga [20] extended this conclusion with numerical experiments, showing that the shape of the influence function depends on the geometry of the reservoir and on the elasticity contrast between the reservoir and its surroundings. Considering that the Los Humeros reservoir is laterally extensive and there is only a slight difference between the fitted responses of the nuclei of strain solutions, our results are inconclusive on this issue.
The corresponding rock compaction volume from the inversion exercise is of the order of over the period of one year. This is at least ~50 times smaller than the net production volume during the period of the InSAR of about . This ratio is much larger than the corresponding ratio of many other geothermal fields. For Cerro Pietro, for instance, the produced volume is only about 3 times larger than the modeled subsurface volume change [4]. Of course, the extracted volume cannot be directly associated with the volume change within the reservoir: the reservoir compressibility, in combination with the pressure reduction of the reservoir, results in a smaller subsurface volume change than the extracted volume. For a closed system, the ratio of bulk moduli of the rock and the water is of the order of the ratio of the produced water volume and the rock compaction volume. Since that would imply an unrealistically high rock bulk modulus of hundreds of GPa, we conclude that there must be pressure support of the porous medium originating most likely from recharge. Recharge supported by meteoric water may occur locally in the northern part of the field, but it is limited to 1-10% of the total produced fluids [25]. Minor subsidence can also be identified on the InSAR outside the production area in the north that may be attributed to these small recharge zones of the production area. Regional recharge outside the study area through deep structures may contribute to the recovery of the extracted fluid volume.
The InSAR data together with the modeling results suggest two reservoir characteristics. First, they indicate that the pressure within the reservoir is well supported suggesting that recharge is taking place. Second, they imply that the Los Humeros Geothermal Field is controlled by sealing faults that separate the reservoir into several blocks. However, additional subsurface data, for instance regular pressure measurements from the wells, are needed to improve our modeling results. This would also allow us to study fault sealing behaviour that controls reservoir compartmentalization. Still, our results make it clear that based on the subsidence pattern we have obtained a better understanding of the pressure conditions within the reservoir and potential local recharge zones. This will facilitate better quality decisions on well planning and operations. Additionally, the surface deformation pattern can contribute to our understanding on superhot geothermal systems that are of high potential for geothermal energy development.
Data Availability
The radar images we processed to arrive at the findings of this study are freely accessible from the European Space Agency via EOLI-SA upon registering to the (A)SAR On The Fly Service. The codes we used to construct our forward models are free-of-charge, and they are available from the following website: http://volcanodeformation.com/software.html. Production data from the Los Humeros Geothermal Field operated by the national Mexican electrical company (Comisión Federal de Electricidad, CFE) are confidential. One may contact the CFE for data request.
Conflicts of Interest
The authors declare that they have conflicts of interest.
Acknowledgments
The authors thank the Comisión Federal de Electricidad (CFE, Mexico) for providing their data for this research. The research leading to these results has received funding from the GEMex Project, funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 727550. More information is available on the GEMex Website: http://www.gemex-h2020.eu.