3D Surface Velocity Field Inferred from SAR Interferometry: Cerro Prieto Step-Over, Mexico, Case Study

Abstract

The Cerro Prieto basin, a tectonically active pull-apart basin, hosts significant geothermal resources currently being exploited in the Cerro Prieto Geothermal Field (CPGF). Consequently, natural tectonic processes and anthropogenic activities contribute to three-dimensional surface displacements in this pull-apart basin. Here, we obtained the Cerro Prieto Step-Over 3D surface velocity field (3DSVF) by accomplishing a weighted least square algorithm inversion from geometrically quasi-orthogonal airborne UAVSAR and RADARSAT-2, Sentinel 1A satellite Synthetic Aperture-Radar (SAR) imagery collected from 2012 to 2016. The 3DSVF results show a vertical rate of 150 mm/yr and 40 mm/yr for the horizontal rate, where for the first time, the north component displacement is achieved by using only the Interferometric SAR time series in the CPGF. Data integration and validation between the 3DSVF and ground-based measurements such as continuous GPS time series and precise leveling data were achieved. Correlating the findings with recent geothermal energy production revealed a subsidence rate slowdown that aligns with the CPGF’s annual vapor production.

1. Introduction

Ground deformation is the superficial expression of several geophysical processes like seismicity, volcanic eruptions, landslides, and subsidence, among others. Precise estimates of the 3D displacement field and its space–time evolution are crucial to comprehend such processes that have originated in the Earth’s subsurface interior. The Interferometric Synthetic Aperture Radar (InSAR) has become a valuable tool for mapping precise ground deformation with high spatial resolution in the radar line of sight (LOS) direction. Nonetheless, accurately estimating the 3D displacement vector from satellite InSAR data is challenging. The acquisition geometries of SAR imaging from space, specifically the look angle and near polar orbits, result in InSAR data having a higher sensitivity to the vertical displacement component compared to the horizontal ones (easting and northing), and the north component exhibiting the highest estimation error, ranging from centimeters to decimeters, even when combining multiple spatial LOS geometries [1,2]. On the other hand, approaches for approximating three-dimensional surface displacement are the following: (a) Combine geodetic data (GPS + InSAR) taking advantage of the GPS three-dimensional precision. However, GPS networks are not as dense as wanted, providing disadvantages to this integration method. (b) Azimuthal offsets or along-track interferometry are used to constrain the northward motion displacement, which is capable of achieving a high-density spatial solution from satellite Synthetic Aperture Radar (SAR) data. However, its applicability is restricted to surface deformation events, earthquakes, or rapid landslides, where horizontal displacements exceed decimeter-scale magnitudes (>10–15 cm), because the along-track interferometry is much less sensitive. (c) Use two different incidence angles like the right–left looking LOS direction, (d) the combination of independent track images, and (e) multi-sensor combination with a hybrid InSAR phase and SAR amplitude pixel offset tracking time-series analysis, which is suitable for measuring the decimeter- to meter-scale displacements but is less precise than conventional InSAR [3,4,5,6,7,8,9,10].

In addition, 3D surface deformation based on InSAR data for geothermal activities is essential to investigate the complex realistic three-dimensional structural, coupled fluid flow and geomechanics modeling that could help to understand the potential and distribution of geothermal resources as well as the surface movements and their relation to volume changes within the reservoir [11,12]. Thus, the significance of the Cerro Prieto Step-Over (CPSO) 3D surface velocity field (3DSVF) lies in its role as a weighted input for estimating changes in the local stress field due to surface deformation, which has been poorly resolved, temporary and spatially [13].

There have been several studies about the CPSO surface velocity field. Nonetheless, they have been limited in obtaining one- or two-dimensional solutions by carrying out geodetic-leveling campaigns and implementing InSAR from polar orbit acquisitions only. Additionally, despite the Global Positioning System (GPS) enabling the determination of three-dimensional displacement, there were only scarce GPS observations in the area with few continuous GPS stations installed. It was first observed by precise leveling and InSAR data that the maximum vertical displacement rate in the CPSO was 120 mm/yr during the 1993–1998 period, related to subsidence [14,15,16]. During the 2005–2016 period, the maximum vertical displacement rates from 140 mm/yr to 160 mm/yr and up to ~180 mm/yr were observed, and ~40 mm/yr for the horizontal rate, with the north component neglected [17,18,19,20]. Space–time variations of the surface deformation rate were significantly associated with the development and operations of the Cerro Prieto Geothermal Field (CPGF).

In this study, we use InSAR data to achieve the 3DSVF model of the CPSO. Interferometric pairs were formed by using the SAR images from (1) spaceborne C-band Copernicus Sentinel-1 satellites, operated by the European Space Agency (ESA), (2) airborne L-band NASA UAVSAR (Uninhabited Aerial Vehicle Synthetic Aperture Radar), operated by the Jet Propulsion Laboratory (JPL), and (3) spaceborne C-band RADARSAT-2, owned and operated by MDA (formerly MacDonald Dettwiler and Associates) [21]. The SAR images from spaceborne missions were acquired in ascending and descending orbits, whereas the SAR images from UAVSAR were acquired in the East and West direction flight lines. The combined spaceborne and airborne data provided independent and semi-orthogonal geometric imagery (Figure 1b). The decomposition of the three-dimensional displacement vector (up, north, and east) was performed by using an inversion scheme formalism based on a weighted least-square solution and employing a subset of independent LOS displacement rate maps derived using a Small Baseline Subset (SBAS) approach. The SBAS approach utilizes differential interferograms to properly constrain the spatial decorrelation. The acquisition pairs are distinguished for having a small orbital separation [22]. The benefits of using SBAS in surface deformation analyses have been proved and enhanced during the last +20 years. It has been applied for subsidence analysis on large urban areas, accurately monitoring slow deformation rates and notable landslides, where spatial and temporal variations were retrieved by enhancing images’ coherence. In addition, the SBAS approach is very convenient for quantifying surface deformation at a millimeter scale, providing space–time features of the observed deformation, since this innovative technique has the capability of producing InSAR time series, which is helpful for analyzing specific signal patterns that represent the investigated geophysical process [23,24,25,26,27].

In the following analysis, we took advantage of the high spatial resolution and moderate temporal coverage of space and airborne SAR images. Using a singular inversion scheme formalism and the obliquity and independence among SAR geometries, we obtained the CPSO three-dimensional surface velocity model caused by tectonic and anthropogenic influence. This is a valuable data source that can help us better understand the subsidence process affecting the region’s infrastructure and farmlands. InSAR, three-dimensional results were compared with precise leveling data (2012–2015) surveyed by the Mexican Institute of Water Technology (IMTA), as well as GPS observations from continuous sites (cGPS) of the Northwestern Mexico Geodetic Network (REGNOM), and interpreted in terms of the recent CPGF activity.

2. Study Area

The CPSO is in the Mexicali Valley (MV) in the northeastern part of Baja California, Mexico, which extends from the international US–Mexico border to the Colorado River Delta in the Gulf of California and from the Colorado River to the Cucapah mountain ranges. The MV is characterized by smooth topography with an average height from ~10 to 20 m above mean sea level (mamsl), where the Cerro Prieto Volcano (CPV) is the most representative topographic feature with 225 mamsl. The CPSO is a complex area delimited by the two main strike-slip faults of the region, Cerro Prieto (CPF) to the southwest and Imperial (IMF) to the northeast, alongside the Morelia and Saltillo dipping faults to the north and south, respectively [28,29,30] (Figure 1). In addition to the tectonic deformation resulting from the relative motion between the North American and Pacific Plates with a velocity of 48–52 mm/yr and a tectonic strain rate of ~1.5 μstrain/yr [31,32,33], anthropogenic activity plays a predominant role, constituting the majority of the total strain rate (~95%). This is primarily attributed to fluid extraction in the CPGF, contributing to a high subsidence rate in the CPSO and resulting in a significant Coulomb stress change rate (0.15–1.7 bar yr⁻¹) sufficient to trigger moderate seismicity in the area [15,20,34].

The extensional tectonics, thick sedimentary layer (sandstones interbedded with shales), high heat flow, active faults, and abundant groundwater resources filtered from the Colorado River have created the ideal conditions for geothermal reservoirs’ formation, like the CPGF in the CPSO, which hosts geothermal fluids at ~2400 m depth. The CPGF is the third largest geothermal field in the world and the most important in Mexico with an installed and operating capacity of 520 MW, producing 2510.6 GWh of gross electricity during 2021 [35,36,37].

The MV is a seismically active zone and absorbs almost all of the plate boundary deformation at this latitude, with the most relevant earthquakes (${\mathrm{M}}_{\mathrm{L}}\ge$ 6) occurring along the Imperial and Cerro Prieto faults. The dominant focal mechanism of these earthquakes is the right-lateral strike-slip [38,39]. Recently, the 2010 Mw 7.2 El Mayor-Cucapah Earthquake Sequence produced the multi-faults rupturing within the Sierra Cucapah mountain range and the MV [40,41,42,43]. Moreover, the area between CPF’s northern tip and the southern end of the IMF, known as the Mexicali Seismic Zone, is dominated by swarms, mainly related to dipping faults bounding the CPSO interior, with depth intervals from 2 to 6 km [39,44,45,46].

Some studies suggest that part of the seismicity in the CPSO is induced by anthropogenic activity related to the CPGF exploitation [34,47,48]. According to the Northwestern Mexico Seismic Network (RESNOM), 138 events occurred during the study period (February/2012–April/2016), with magnitudes of 1.1 ≤ ${\mathrm{M}}_{\mathrm{L}}$ ≤ 4.2 in the CPSO, and taking into account that the InSAR technique is only able to retrieve surface deformation caused by shallow events of ${\mathrm{M}}_{\mathrm{W}}$ ≥ 4.5 [49,50,51,52,53], we considered the detected surface deformation mainly as aseismic and anthropogenic in this study.

3. Materials and Methods

The SAR Single Look Complex (SLC) images acquired over the study area by Sentinel-1A and UAVSAR were collected and processed. Sentinel-1A and UAVSAR data were retrieved from the NASA online repository hosted by the University of Alaska (Alaska Satellite Facility, ASF) and the NASA/JPL UAVSAR website “https://uavsar.jpl.nasa.gov/cgi-bin/data.pl (accessed on 22 April 2024)” data products, respectively. SAR SLC images from the ascending (T166) and descending (T173) tracks of Sentinel-1A, as well as from the east-oriented (EO, 08514S3) and west-oriented (WO, 26515S2) flight segments of UAVSAR, were employed. Additionally, LOS displacement maps derived from RADARSAT-2 images captured during ascending (MF1) and descending (MF4N) orbital passes were incorporated [21]. The basic parameters of SAR and datasets used in this study are presented in

Table 1. Basic parameters of SAR datasets used in this study. α_h: heading; θ_inc: incidence angle; N: number of acquisitions (SLC); M: number of interferograms per subset; asc = ascending orbital pass; dsc = descending orbital pass.

Sensor	Subset	Orbit	Time Span 4	${\mathsf{\alpha }}_{\mathbf{h}}$ (°)	${\mathbf{\theta }}_{\mathbf{i}\mathbf{n}\mathbf{c}}$ (°)	N	M
Sentinel 1A 1	T166	asc	3 April 2015–9 April 2016	347.39	37.2	26	49
Sentinel 1A 1	T173	dsc	3 April 2015–9 April 2016	192.7	34.85	26	49
UAVSAR 2	08514S3	east	1 February 2012–1 April 2016	84.8	50.4	12	20
UAVSAR 2	26515S2	west	1 February 2012–1 April 2016	−94.8	43.85	12	20
RADARSAT-2 3	MF1	asc	13 September 2011–24 July 2016	349.4	38.4	53	434
RADARSAT-2 3	MF4N	dsc	1 October 2011–11 August 2016	−170.3	44	58	344

¹ Right side looking. ² Left side looking/time acquisition is variable. ³ Left and right side looking. ⁴ Time span represents the whole period of time processed per data source, but not all was included in final analysis.

. The interferometric processing and 3D displacement vector decomposition were carried out for the spanned time intervals by spaceborne (ascending and descending orbital passes) and airborne (EO and WO flight segments) data (Figure 2). Two subsets were used, RADARSAT-2+UAVSAR and Sentinel-1A+UAVSAR, covering different periods for later comparison between leveling and cGPS data. Thus, the RADARSAT-2+UAVSAR subset covers the period from 1 February 2012 to 18 November 2014, while the Sentinel-1A+UAVSAR subset covers the period from 3 April 2015 to 1 April 2016.

In this study, the interferometric processing includes two steps. The first step comprises the unwrapped differential interferograms’ generation from SLC SAR images using the conventional two-pass differential InSAR approach [54]. Interferometric pairs from Sentinel 1A and UAVSAR images were generated by using the open-source software InSAR Scientific Computing Environment (ISCE) version 2 [55] (Figure 3). The 1 Arc-second (30 m) global digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) [56] was employed for co-registration, removing the topographic contribution from the interferometric phase, and geocoding. During the interferometric processing of Sentinel-1A images, the precise orbit determination (POD) ephemeris and calibration auxiliary data supplied by the ESA were utilized. The UAVSAR processing from raw data to a co-registered stack of SLC images was performed by the UAVSAR team, and included motion compensation and topographic phase correction. We started with the UAVSAR co-registered SLC stack to make interferograms and unwrap them. For each dataset, the SAR images were aligned or co-registered to the chosen reference scene in radar coordinates. The interferogram processing was completed at approximately 30 m pixel size, achieved through looks or pixel averaging in range and azimuth directions. Differential interferograms (without phase contributions of the Earth’s curvature and topography) were filtered using the Goldstein–Werner power-spectral filter [57] to enhance the signal-to-noise ratio and improve the phase unwrapping performance. The phase unwrapping was carried out using the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU) [58]. Finally, unwrapped differential interferograms were geocoded to the WGS84 geographic coordinate system. Interferometric processing of the RADARSAT-2 pairs was performed with the GAMMA software with an equivalent workflow [21].

The second step of the interferometric processing consisted of a time series analysis employing a Small Baseline Subset (SBAS) approach, which involves combining and inverting a collection of unwrapped interferograms to determine the surface average LOS displacement rate and LOS displacement time series [22,59]. To avoid the spatio-temporal decorrelation in Sentinel 1A data, only the interferometric pairs with a short perpendicular spatial baseline (<200 m) and relatively short temporal baseline were used during SBAS processing. For each scene, two subsequent scenes were considered for interferometric pair formation (Figure S1), “double redundancy” consideration. For UAVSAR images’ acquisition, NASA’s aircraft uses a Precision GPS-controlled Autopilot that typically achieves a repeat-track narrow tube of 5 m radius to achieve short perpendicular baselines of less than 10 m [8,60]. However, only two subsequent scenes were considered for interferometric pair formation to prevent temporal decorrelation. In addition, the unwrapped differential interferograms were visually examined to identify and exclude from the SBAS processing those interferograms displaying obvious effects, such as temporal decorrelation, phase unwrapping errors, and atmospheric contribution (Figures S2 and S3). The SBAS approach was performed using the Generic InSAR Analysis Toolbox (GIAnT), Caltech, Pasadena, CA, USA [59] (Figure 3).

The interferograms were referenced to the same coherent pixel (reference point) considered stable. The reference point was established in the CPV (−115.3 W/32.4 N), situated outside the boundaries of the CPSO area. During SBAS processing, several steps were undertaken, including stratified tropospheric delay correction, the removal of phase ramps, correction for DEM errors, the estimation of temporal coherence to assess inversion quality, and ultimately, the determination of the average LOS surface displacement rate. Additionally, the processing involved the generation of LOS displacement time-series and the determination of associated errors by using the Jackknife test resampling technique for each pixel and epoch. The ECMWF weather models’ atmospheric correction was specifically performed on Sentinel-1A data. This effect is mainly correlated with the steep topography, such as mountain ranges [61]. The stratified component of the tropospheric phase delay is mapped from global atmospheric models that provide estimates of the air temperature, atmospheric pressure, and humidity as a function of elevation (topography), and then is removed from the interferograms [62].

We developed an inversion scheme algorithm based on the three-dimensional vector decomposition equation projected on the range component (dlos) in the LOS direction [16] (Hanssen, 2001).

$$\mathrm{dlos}={\mathrm{d}}_{\mathrm{u}}\mathrm{Cos}\left({\mathsf{\theta }}_{\mathrm{inc}}\right)-{ \mathrm{d}}_{\mathrm{n}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{inc}}\right)\mathrm{Cos}\left({\mathsf{\alpha }}_{\mathrm{h}}-{\displaystyle \frac{3\mathsf{\pi }}{2}}\right)-{ \mathrm{d}}_{\mathrm{e}}\mathrm{Sin}\left({\mathsf{\theta }}_{\mathrm{inc}}\right)\mathrm{Sin}({\mathsf{\alpha }}_{\mathrm{h}} - {\displaystyle \frac{3\mathsf{\pi }}{2}})$$

(1)

where ${\mathrm{d}}_{\mathrm{u}}$, ${ \mathrm{d}}_{\mathrm{n}}$, and ${ \mathrm{d}}_{\mathrm{e}}$ are the displacement vector components: up, north, and east, respectively; ${\mathsf{\theta }}_{\mathrm{inc}}$ is the incidence angle; ${\mathsf{\alpha }}_{\mathrm{h}}$ the platform heading; and ${\mathsf{\alpha }}_{\mathrm{h}}-{\displaystyle \frac{3\mathsf{\pi }}{2}}$ is the look azimuth, positive clockwise from north [6]. The incidence angle was considered variable during the 3D inversion process; as an example, the UAVSAR incidence angle changes up to 20° over the study area. Equation (1) can be represented in its matrix form per each SAR geometry and the three-dimensional displacement components, respectively (Equations (2) and (3)).

$$\mathrm{d}=\left[\begin{array}{c}\mathrm{dlosa}\\ \mathrm{dlosd}\\ \mathrm{dlose}\\ \mathrm{dlosw}\end{array}\right]$$

(2)

$$\mathrm{G}=\left[\begin{array}{ccc}\mathrm{Cos}{(\mathsf{\theta }}_{\mathrm{Sa}})& \sin {(\mathsf{\theta }}_{\mathrm{Sa}}\left)\cos \right({\text{∝}}_{\mathrm{Sa}})& \sin {(\mathsf{\theta }}_{\mathrm{Sa}}\left)\sin \right({\text{∝}}_{\mathrm{Sa}})\\ \mathrm{Cos}{(\mathsf{\theta }}_{\mathrm{Sd}})& \sin {(\mathsf{\theta }}_{\mathrm{Sd}}\left)\cos \right({\text{∝}}_{\mathrm{Sd}})& \sin {(\mathsf{\theta }}_{\mathrm{Sd}}\left)\sin \right({\text{∝}}_{\mathrm{Sd}})\\ \mathrm{Cos}{(\mathsf{\theta }}_{\mathrm{Uwe}})& \sin {(\mathsf{\theta }}_{\mathrm{Uwe}}\left)\cos \right({\text{∝}}_{\mathrm{Uwe}})& \sin {(\mathsf{\theta }}_{\mathrm{Uwe}})\sin \left({\text{∝}}_{\mathrm{Uwe}}\right)\\ \mathrm{Cos}{(\mathsf{\theta }}_{\mathrm{Uew}})& \sin {(\mathsf{\theta }}_{\mathrm{Uew}}\left)\cos \right({\text{∝}}_{\mathrm{Uew}})& \sin {(\mathsf{\theta }}_{\mathrm{Uew}})\sin \left({\text{∝}}_{\mathrm{Uew}}\right)\end{array}\right]$$

(3)

where S: Sentinel, U: UAVSAR, a: ascending, d: descending, e: east, and w: west.

The 3D inversion process, performed pixel by pixel, was addressed using a Weighted Least Square Solution (${\mathrm{m}}_{\mathrm{WLS}}$), as outlined in (Equation (4)). This solution is described by the data ($\mathrm{d}$) and the weight matrix (${\mathrm{W}}_{\mathrm{e}}$), which defines the relative contribution of each single error to the total prediction error. The weight matrix is typically diagonal, and $\mathrm{G}$ represents the data kernel (Equation (3)). The data kernel is constructed with the directional cosines, taking into account the variable incidence angle and heading of each SAR geometry.

$${{\mathrm{m}}_{\mathrm{WLS}}}^{\mathrm{est}}={{[{\mathrm{G}}^{\mathrm{T}}{\mathrm{W}}_{\mathrm{e}}\mathrm{G}]}^{-1}{\mathrm{G}}^{\mathrm{T}}\mathrm{W}}_{\mathrm{e}}\mathrm{d}$$

(4)

The sensitivity vector components of each dataset used in this study are detailed in

Table 2. Sensitivity vector components of each dataset. α_h: heading; θ_inc: incidence angle; d_u, d_n, and d_e are the values of the sensitivity vector.

Sensor	${\mathbf{\theta }}_{\mathbf{i}\mathbf{n}\mathbf{c}}$ (°) 1	${\mathsf{\alpha }}_{\mathbf{h}}$ (°)	du	dn	de
Sentinel 1A	35	349	0.82	0.11	0.56
UAVSAR	45 2	85	0.77	0.65	0.06

¹ It is the average value of the incidence angle. ² UAVSAR incidence angle varies in a 40° range [63].

. The 3D inversion of the UAVSAR and RADARSAT-2 data used the same Equations (3) and (4), with the substitution of RADARSAT-2 ascending and descending tracks for Sentinel-1 ascending and descending tracks.

The coding algorithm was developed using the Matlab R2016a (9.0.0.341360) computing tool and underwent validation through synthetic surface displacement data (Figure 4). These synthetic data were derived from an anthropogenic deformation model created by [17] and implemented in the Coulomb 3.1 software [64,65] (Figure 4a). The deformation model consists of seven rectangular tensional cracks embedded in an elastic half-space [66], modeling the compaction effect of reservoirs or aquifers within the sedimentary strata delimited by faults. An accurate description of the thee-dimensional synthetic displacement field of the CPSO area was achieved. Based on Equation (1), the four LOS displacement models (ascending/descending/east/west), geometrically oblique and independent among them, were calculated by using the SAR sensors’ geometry values (see Table 1) and the 3D synthetic data as input. The calculated synthetic LOS displacement values served as input to calculate the new 3D surface displacement field (Figure 4b) by employing the inversion algorithm code based on Equation (4). Subtracting the 3D synthetic data from the 3D calculated displacement field, a discrepancy technically around zero (~10⁻¹⁵) is obtained (Figure 4c). Various tests were conducted, where Gaussian random errors were incorporated to synthetic data, by repeating the process previously described. Therefore, the difference between synthetic and calculated 3D surface displacement data is consistent with the added random errors (Figure 4d), which validates the algorithm code functioning and makes it feasible for real data (InSAR LOS displacement) application aimed to retrieve the three-dimensional surface displacement field of the CPSO.

4. Results

4.1. SBAS Approach

Figure 5 and Figure 6 present the LOS displacement rate maps derived from the SBAS approach. The maps show a LOS displacement rate of up to ~120 mm/yr, with a maximum on the north-eastern side of the CPSO. This area corresponds to the recharge zone (RZ) as proposed [15,34], while in the extraction zone (EZ) of the CPGF, the maximum LOS displacement rate reaches ~70 mm/yr (Figure 6a). The error estimate of the LOS displacement rate was determined for each dataset. The UAVSAR data showed higher estimation errors, primarily due to the limited availability of UAVSAR images, leading to significant temporal separation between images within its datasets and due to the four times larger radar wavelength that reduces sensitivity to displacements. For UAVSAR interferometric data covering one year, the average error of coherent pixels is ±10 mm/yr, with certain areas displaying higher error estimates exceeding ±50 mm/yr (Figure S4a,b). For UAVSAR interferometric data covering 2.8 years, the maximum estimated error reaches ±35 mm/yr. The Sentinel 1A maximum estimated error reaches ±5 mm/yr (Figure S4c,d).

4.2. Inversion Problem

The 3D velocity solution was obtained pixel by pixel, within an array of 818 by 649 elements and with a 30 m spatial resolution in every coherent pixel over the CPSO area.

4.2.1. 3DSVF (2012–2015 Period)

The 3D displacement vector decomposition for 2012–2015 was carried out using a combination of RADARSAT-2 and UAVSAR interferometric data. The interferometric data used, derived from the SBAS approach, are depicted in Figure 5. The results of the 3D displacement vector decomposition are presented in Figure 7.

It has been observed that error estimates within the UAVSAR datasets significantly influence the north displacement component, leading to unrealistic displacements beyond the expected range where the UAVSAR geometry exhibits high sensitivity (Table 2). Consequently, pixels with low temporal coherence (<0.27) were masked out, as depicted in Figure 7. This approach enabled a focused analysis of how errors in the UAVSAR dataset specifically impact the north and up surface velocity maps, in contrast to the east component map (Figure 7c).

The vertical displacement outlines a clear subsidence footprint bounded by the CPF and IMF, with a maximum rate of 130 mm/yr (Figure 7a), located in the northeastern CPSO area (RZ) [15,34]. The maximum vertical displacement in the CPGF extraction zone reaches 135 mm/yr. The subsidence area is delimited by tectonic faults, which constitute the limits of the Cerro Prieto basin and the CPSO.

The horizontal displacement (north and east) components indicate the movement of the CPSO boundary zones toward its center, and the area of maximum subsidence, as expected for a volume decrease at the roughly 2.7 km depth of the geothermal reservoir along the axis of the CPGF. Consequently, the areas located to the north–northwest of the Cerro Prieto basin demonstrate southward and eastward movement, with a maximum rate of approximately 40 mm/yr for both components. Similarly, the areas located to the south and east of the Cerro Prieto basin exhibit northward and westward movement, displaying similar maximum magnitudes. The horizontal displacement magnitude decreases toward the subsidence bowl’s central zone.

4.2.2. 3DSVF (2015–2016 Period)

The 3D displacement vector decomposition for 2015–2016 was carried out using a combination of Sentinel-1A and UAVSAR interferometric data. The interferometric data used, derived from the SBAS approach, are depicted in Figure 6. The results of 3D displacement vector decomposition are presented in Figure 8. In Figure 8, the pixels exhibiting the errors of estimates of more than 20 mm/yr were masked out and excluded from analysis. In this time period, the intervals between UAVSAR acquisitions were larger, so more area is masked out.

For 2015–2016, the derived rate maps of displacement vector components present similar patterns to those observed for 2012–2015. The maps illustrating the displacement field reveal the affected area associated with land subsidence, expanding in a NE-SW direction from the southern tip of the IMF to the northern tip of the CPF. The maximum rate of vertical displacement exceeds 150 mm/yr around the recharge zone in the period from April 2015 to April 2016, while a second high-velocity area, ~120 mm/yr, is located to the east of the Ejido Nuevo León in the vicinity of the CPGF (Figure 8a). Regarding the horizontal north and east displacement components (Figure 8b,c), the maps show an average surface displacement rate of up to 40 mm/yr. The convergence effect of displacements is evident toward the basin’s center, corresponding to the subsidence phenomenon, an area delimited by the MV main fault system.

4.3. 3DSVF Validaton

The vertical displacement rate obtained by the RADARSAT-2+UAVSAR subsets’ combination was compared with precise leveling data surveyed by the Mexican Institute of Water Technology (IMTA) during 2012–2015. The contour maps of leveling data (Figure 9a) and vertical displacement InSAR data (Figure 9b) were generated using the Kriging interpolation algorithm. We used the same pixel array dimension by considering the distribution of 47 out of 66 benchmarks coincident between both subsets in the study area. We highlight that the InSAR reference point is in a stable place different from that of leveling data (10037 benchmark: RP in Figure 9a) because UAVSAR flight segments do not cover the area where leveling RP is located. Nonetheless, the comparison indicates a reasonably good agreement in both magnitude and trend between the results of these two independent techniques, as illustrated in Figure 9. The total root-mean-square error (RMSE) of 0.96 cm/yr and a correlation coefficient of 0.91 are obtained between these two distinct measurements (Figure 10). However, notable disparities, reaching up to 2 cm/yr, are noticeable in the southeastern part of the study area, in proximity to the southern tip of the Saltillo Fault (Figure 9c). The smallest differences between datasets are observed in the CPGF operation zone and the northeast side of the CPSO in the recharge zone, with rates ranging from 0.1 to 0.4 cm/yr (Figure 9c).

Additionally, we conducted a comparison of the InSAR 3DSVF derived from both subsets’ combinations, RADARSAT-2+UAVSAR and Sentinel 1A+UAVSAR (

Table 3. InSAR-derived three-dimensional surface displacement rate (mm/yr) at GPS station locations and its three-dimensional RMSE.

Location	Lat	Long	North 2	East 2	Up 2	RMSE 2 3D	North 3	East 3	Up 3	RMSE 3 3D
MBIG 1	32.4089	−115.1960	−23.3	7.9	−96.7	6.8	−6.1	7.0	−90.0	6.5
NVLX 1	32.3935	−115.1832	8.3	−9.5	−83.7	7.8	25.8	−8.3	−76.3	1.2

¹ Continuous GPS stations of the Geodetic Network of the Northwestern Mexico, located in the CPSO. ² Sentinel 1A-UAVSAR subset. ³ RADARSAT-2-UAVSAR subset.

), with GPS velocity data (Figure S5) from REGNOM continuous geodetic stations [67], referenced to a GPS site located in the Cerro Prieto Volcano (CPIG-GPS) [33], the same reference point pixel location used in the InSAR time series analysis obtained here.

We used two cGPS stations located inside the CPSO. The total RMSE < 8 mm/yr is obtained during overlapping time periods, for the three-dimensional surface displacement rate obtained with InSAR (Table 3). These results demonstrate a robust agreement between both surface displacement techniques, in line with other studies [8].

5. Discussion

The SBAS approach was performed on multiple SAR datasets to obtain LOS surface displacement rates from independent and semi-orthogonal geometric imagery. Two In-SAR SBAS dataset combinations were used to derive the 3DSVF of the CPSO for two distinct periods.

This study demonstrates the feasibility of using InSAR to estimate 3D surface deformation, as evidenced by the obtained 3DSVF. The study results are tested and validated by precise leveling and GPS observations during overlapping periods. For the 2012–2015 period, a comparison of vertical displacement rates between InSAR and leveling measurements yields a total RMSE of 0.96 cm/yr (Figure 10), similar to other studies in the area [19]. We can also notice an overestimated bias in the vertical displacement rate from InSAR (mu = 0.65 cm/yr), possibly due to the different reference point used in leveling data, where the Cerro Prieto Volcano vertical displacement is −0.5 cm/yr [33]. Our analysis comparison for vertical displacement and 3DSVF generally shows high consistency between InSAR results, leveling measurements, and GPS observations during both periods.

The subsidence process is mainly associated with vertical displacements. Nonetheless, it also includes horizontal motion (east/north) that converges in a concentric pattern of water table drawdown at the center of the Cerro Prieto basin in a NE-SW orientation. Horizontal displacements are relatively small in the maximum subsidence area, significant on its periphery, and delimited by the main active faults. This is a tectonic fault controlling system described by a transtensional regime, due to dextral strike-slip faults where seismic activity ($\le {\mathrm{M}}_{\mathrm{L}}$6) is present, along with smaller earthquakes and swarms in the oblique normal faults, the latter being characterized to the north and south periphery sections of the CPSO [30]. Faults in the CPSO behave as a geological limit that could be also described by the high spatial resolution of InSAR data (Figure 7 and Figure 8). We compared our estimated vertical displacement rate for different periods with estimates from previous works using pixels selected in the recharge and exploitation zones (Figure 11a). For the period of April/2015–April/2016, we obtained vertical displacement rates of ~150 mm/yr and ~100 mm/yr in the recharge and exploitation zones (CPGF), similar to [13] for the period of October/2014–July/2016. Conversely, for the February/2012–November/2014 period, we found ~135 mm/yr and ~98 mm/yr in the recharge and exploitation zones, respectively, which differed significantly compared to the one estimated by the 2D LOS displacement decomposition for the November/2012–March/2014 period [19], where the vertical displacement rate in the recharge zone is ~160 mm/yr, although the maximum rate in the area is up to 180 mm/yr, while in the exploitation zone, the vertical displacement rate is ~80 mm/yr. Continuous GPS time series from a site located in Ejido Nuevo León at the eastern limits of the CPGF exploitation zone shows a vertical displacement rate of 103 mm/yr for the period of October/2015–April/2015, since the MBIG site started working in October 2015 (Figure S5a), equivalent to the one estimated here. Therefore, this variability suggests the apparent influence of the geometric components of the SAR sensors, the selected InSAR reference point during the temporal analysis, and the 2D vs. 3D LOS approaches, leading to a vertical displacement average rate of ~155 ± 5 (160 and 150) mm/yr and 90 ± 10 (100 and 80) mm/yr for the recharge and exploitation zones, respectively, between both periods analyzed here.

The subsidence rate observed in this study suggests a deceleration during the period from April 2015 to April 2016 compared to the period from February 2012 to November 2014, as reported in previous works [19]. Specifically, this comparison shows an approximate slowing down rate of about 6% (10 mm/yr ↓) in the recharge zone. However, if we compare our results between these two periods, we can observe a subsidence rate increment of about 13% (19 mm/yr ↑). In contrast, there is a clear rate increment in the exploitation zone of around 19% (18.8 mm/yr ↑) when comparing our results to those of other authors [17], as while comparing our results, the subsidence rate increment is only about 9% (9.8 mm/yr ↑). Although this provides solid insights about the CPSO surface deformation, subsidence rate variations could not be directly linked to the production of geothermal electricity in the CPGF, but are strongly related to steam extraction and water reinjection data, which unfortunately are not accessible to the public. Only the number of wells in operation for extraction and injection is released to the public. Ref. [36] observed a reduction of 20% in the average production rate per well in 2013 compared to 2008. A reduction in electricity generation is observed between 2011 and 2017, with a total production of 4547 GWh and 3554 GWh for those years, respectively [68,69] (Figure 11b). Additionally, the IEA Geothermal annual reports [70] indicate a decrease in production wells from 160 to 147 during the study period, while the number of injection wells increased from 18 to 29 (Figure 11b). It is important to note that a comprehensive analysis of CPGF energy production [71,72] is beyond the scope of this work.

An examination of the vertical displacement rates derived from the 3D approximation using both dataset combinations, RADARSAT-2+UAVSAR and Sentinel 1A+UAVSAR, in comparison to the 2D approximation by other authors, suggests that omitting the north component in surface deformation studies, influenced by both tectonic and anthropogenic activities, may result in over- or underestimations of vertical displacement rates. This might be explained by the weighting contribution of every displacement component during the vector decomposition process, propagating such a contribution to the up and east components when the north one is neglected in the analysis, which is particularly noteworthy when the velocity vectors of this north component exhibit large magnitudes (>40 mm/yr). In addition, it is important to remark that SAR sensors’ geometries, the time span coverage among datasets, and InSAR velocity calculations play an important role in the variability of the displacement components’ (east/north/up) estimations, describing the main limitations of this study. However, the findings of this study indicate that accurately considering all three-dimensional displacement components enables a more precise estimation of the stress changes in the neighboring faults.

6. Conclusions

The 3DSVF for the CPSO was exclusively derived from InSAR data. The combination of spaceborne and airborne InSAR data, sourced independently and from semi-orthogonal geometries, proved essential for accurately estimating the north–south displacement component. This comprehensive approach addresses the limitations associated with re-lying solely on spaceborne data and enhances the overall effectiveness of the 3DSVF estimation in capturing diverse displacement components.

Time series analyses of InSAR data in the line of sight (LOS) direction, along with the 3DSVF results in the up, north, and east directions, reveal that the Imperial and Cerro Prieto faults serve as the physical boundaries of the CPSO. The deformed area is predominantly situated between the traces of the faults.

Subsidence rates, reaching up to approximately 150 mm/yr, align with the average rate estimated over the last decade. The variability in the surface displacement rates for equivalent periods is influenced by the geometric properties of different Synthetic Aperture Radar (SAR) systems, such as the azimuth and incidence angle, as well as the location of the InSAR reference point. Omitting the north–south displacement component during the decomposition of the LOS displacement vector can lead to both the overestimation and underestimation of the other components (vertical and east) of the displacement vector. Therefore, this exclusion affects the accuracy of vertical and eastward displacement components, particularly when the north–south component is significant. This observation underscores the significance of considering the 3DSVF in seismic hazard assessments, especially in relation to subsidence and 3D surface deformation phenomena within the CPSO.