Open AccessArticle

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

Ignacio F. Garcia-Meza

^1,2

J. Alejandro González-Ortega

^3,*

Olga Sarychikhina

³,

Eric J. Fielding

and

Sergey Samsonov

⁴

Posgrado en Ciencias de la Tierra, Centro de Investigación Científica y de Educación Superior de Ensenada, Ensenada 22860, BC, Mexico

Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91011, USA

División de Ciencias de la Tierra, Centro de Investigación Científica y de Educación Superior de Ensenada, Ensenada 22860, BC, Mexico

⁴

Canada Centre for Mapping and Earth Observation, Natural Resources Canada, Ottawa, ON K1S3L3, Canada

Author to whom correspondence should be addressed.

Remote Sens. 2024, 16(20), 3788; https://doi.org/10.3390/rs16203788

Submission received: 9 August 2024 / Revised: 20 September 2024 / Accepted: 20 September 2024 / Published: 12 October 2024

(This article belongs to the Special Issue Advanced Remote Sensing Technology in Geodesy, Surveying and Mapping)

Download

Browse Figures

Versions Notes

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.

Keywords:

3D surface displacements; Cerro Prieto; SAR sensors; interferometry; InSAR; inversion; subsidence

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 (

M_{L} \geq

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 ≤

M_{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

M_{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. 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).

dlos = d_{u} Cos (θ_{inc}) - d_{n} Sin (θ_{inc}) Cos (α_{h} - \frac{3 π}{2}) - d_{e} Sin (θ_{inc}) Sin (α_{h} - \frac{3 π}{2})

(1)

where

d_{u}

d_{n}

, and

d_{e}

are the displacement vector components: up, north, and east, respectively;

θ_{inc}

is the incidence angle;

α_{h}

the platform heading; and

α_{h} - \frac{3 π}{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)).

d = [\begin{matrix} dlosa \\ dlosd \\ dlose \\ dlosw \end{matrix}]

(2)

G = [\begin{matrix} Cos {(θ}_{Sa}) & \sin {(θ}_{Sa}) \cos (∝_{Sa}) & \sin {(θ}_{Sa}) \sin (∝_{Sa}) \\ Cos {(θ}_{Sd}) & \sin {(θ}_{Sd}) \cos (∝_{Sd}) & \sin {(θ}_{Sd}) \sin (∝_{Sd}) \\ Cos {(θ}_{Uwe}) & \sin {(θ}_{Uwe}) \cos (∝_{Uwe}) & \sin {(θ}_{Uwe}) \sin (∝_{Uwe}) \\ Cos {(θ}_{Uew}) & \sin {(θ}_{Uew}) \cos (∝_{Uew}) & \sin {(θ}_{Uew}) \sin (∝_{Uew}) \end{matrix}]

(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 (

m_{WLS}

), as outlined in (Equation (4)). This solution is described by the data (

d

) and the weight matrix (

W_{e}

), which defines the relative contribution of each single error to the total prediction error. The weight matrix is typically diagonal, and

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.

{m_{WLS}}^{est} = {{[G^{T} W_{e} G]}^{- 1} G^{T} W}_{e} d

(4)

The sensitivity vector components of each dataset used in this study are detailed in Table 2. 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), 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 (

\leq M_{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.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs16203788/s1, Figure S1: Perpendicular baseline versus time plot for Sentinel 1A ascending and descending orbital passes; Figure S2: LOS displacement maps from Sentinel 1A unwrapped interferograms; Figure S3: LOS displacement maps from UAVSAR unwrapped interferograms; Figure S4: Mean errors of LOS displacement maps from Sentinel 1A and UAVSAR; Figure S5: Continuous REGNOM sites GPS velocities in Igb14-NOAM; File S1: RADARSAT_UAVSAR_east.grd; File S2: RADARSAT_UAVSAR_north.grd; File S3: RADARSAT_UAVSAR_up.grd; File S4: SENTINEL_UAVSAR_east.grd; File S5: SENTINEL_UAVSAR_north.grd; File S6: SENTINEL_UAVSAR_up.grd.

Author Contributions

Conceptualization, I.F.G.-M. and J.A.G.-O.; methodology, I.F.G.-M. and J.A.G.-O.; formal analysis, I.F.G.-M.; investigation, I.F.G.-M.; data curation, I.F.G.-M., J.A.G.-O. and S.S.; writing—original draft preparation, I.F.G.-M.; writing—review and editing, I.F.G.-M., J.A.G.-O., O.S. and E.J.F.; visualization, I.F.G.-M., J.A.G.-O., O.S. and E.J.F.; supervision, J.A.G.-O.; project administration, J.A.G.-O. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Council of Humanities Science and Technology (CONAHCYT) and the CICESE Earth Science postgraduate internal project number C0F046.

Data Availability Statement

The authors confirmed that the data supporting the findings of this study are available within the article and Supplementary Materials files.

Acknowledgments

We thank the Assistant Editor and three anonymous reviewers for their thorough and thoughtful reviews of this manuscript. Sentinel-1 data were provided by the European Space Agency (ESA) through the Alaska Satellite Facility (ASF). UAVSAR data were provided through “https://uavsar.jpl.nasa.gov/ (accessed on 22 April 2024)”. We also thank the Mexican Institute of Water Technology (IMTA) for the precise leveling data used in this work. Some of the figures were generated using the Generic Mapping Tools [73] (GMT; Wessel et al., 2019). We would also like to thank the team who carried out part of this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We thank the UAVSAR team and NASA pilots for their help in acquiring and processing the UAVSAR data, and the Earth Surface and Interior focus area for supporting our research.

Conflicts of Interest

The authors declare no conflict of interest.

References

Fuhrmann, T.; Garthwaite, M.C. Resolving three-dimensional surface motion with InSAR: Constraints from multi-geometry data fusion. Remote Sens. 2019, 11, 241. [Google Scholar] [CrossRef]
Mehrabi, H.; Voosoghi, B.; Motagh, M.; Hanssen, R.F. Three-Dimensional Displacement Fields from InSAR through Tikhonov Regularization and Least-Squares Variance Component Estimation. J. Surv. Eng. 2019, 145, 04019011. [Google Scholar] [CrossRef]
Fialko, Y.; Simons, M.; Agnew, D. The complete (3-D) surface displacement field in the epicentral area of the 1999 Mw 7.1 Hector Mine earthquake, California, from space geodetic observations. Geophys. Res. Lett. 2001, 28, 3063–3066. [Google Scholar] [CrossRef]
Wright, T.J.; Parsons, B.E.; Lu, Z. Toward mapping surface deformation in three dimensions using InSAR. Geophys. Res. Lett. 2004, 31, L01607. [Google Scholar] [CrossRef]
González, P.J.; Fernández, J.; Camacho, A.G. Coseismic three-dimensional displacements determined Using SAR data: Theory and an application test. Pure Appl. Geophys. 2009, 166, 1403–1424. [Google Scholar] [CrossRef]
Hu, J.; Li, Z.W.; Ding, X.L.; Zhu, J.J.; Zhang, L.; Sun, Q. Resolving three-dimensional surface displacements from InSAR measurements: A review. Earth-Sci. Rev. 2014, 133, 1–17. [Google Scholar] [CrossRef]
Polcari, M.; Palano, M.; Fernández, J.; Samsonov, S.V.; Stramondo, S.; Zerbini, S. 3D displacement field retrieved by integrating Sentinel-1 InSAR and GPS data: The 2014 South Napa earthquake. Eur. J. Remote Sens. 2016, 49, 1–13. [Google Scholar] [CrossRef]
Delbridge, B.G.; Bürgmann, R.; Fielding, E.; Hensley, S.; Schulz, W.H. Three-dimensional surface deformation derived from airborne interferometric UAVSAR: Application to the Slumgullion Landslide. J. Geophys. Res. Solid Earth 2016, 121, 3951–3977. [Google Scholar] [CrossRef]
Hu, X.; Bürgmann, R.; Schulz, W.H.; Fielding, E.J. Four-dimensional surface motions of the Slumgullion landslide and quantification of hydrometeorological forcing. Nat. Commun. 2020, 11, 2792. [Google Scholar] [CrossRef]
Handwerger, A.L.; Booth, A.M.; Huang, M.; Fielding, E.J. Inferring the Subsurface Geometry and Strength of Slow-Moving Landslides Using 3-D Velocity Measurements From the NASA/JPL UAVSAR. J. Geophys. Res. Earth Surf. 2021, 126, e2020JF005898. [Google Scholar] [CrossRef]
Vasco, D.W.; Rutqvist, J.; Jeanne, P.; Samsonov, S.V.; Hartline, C. Using geodetic data in geothermal areas. Lead. Edge 2020, 39, 883–892. [Google Scholar] [CrossRef]
Yazbeck, J.; Rundle, J.B. A Fusion of Geothermal and InSAR Data with Machine Learning for Enhanced Deformation Forecasting at the Geysers. Land 2023, 12, 1977. [Google Scholar] [CrossRef]
Xu, X.; Sandwell, D.T.; Tymofyeyeva, E.; Gonzalez-Ortega, A.; Tong, X. Tectonic and Anthropogenic Deformation at the Cerro Prieto Geothermal Step-Over Revealed by Sentinel-1A InSAR. IEEE Trans. Geosci. Remote Sens. 2017, 55, 5284–5292. [Google Scholar] [CrossRef]
Carnec, C.; Fabriol, H. Monitoring and modeling land subsidence at the Cerro Prieto Geothermal Field, Baja California, Mexico, using SAR interferometry. Geophys. Res. Lett. 1999, 26, 1211–1214. [Google Scholar] [CrossRef]
Glowacka, E.; González, J.; Fabriol, H. Recent Vertical Deformation in Mexicali Valley and its Relationship with Tectonics, Seismicity, and the Exploitation of the Cerro Prieto Geothermal Field, Mexico. Pure Appl. Geophys. 1999, 156, 591–614. [Google Scholar] [CrossRef]
Hanssen, R.F. Radar Interferometry; Springer: Dordrecht, The Netherlands, 2001; Volume 2. [Google Scholar] [CrossRef]
Sarychikhina, O.; Glowacka, E.; Mellors, R.; Vidal, F.S. Land subsidence in the Cerro Prieto Geothermal Field, Baja California, Mexico, from 1994 to 2005. J. Volcanol. Geotherm. Res. 2011, 204, 76–90. [Google Scholar] [CrossRef]
Sarychikhina, O.; Glowacka, E.; Robles, B.; Nava, F.A.; Guzmán, M. Estimation of Seismic and Aseismic Deformation in Mexicali Valley, Baja California, Mexico, in the 2006–2009 Period, Using Precise Leveling, DInSAR, Geotechnical Instruments Data, and Modeling. Pure Appl. Geophys. 2015, 172, 3139–3162. [Google Scholar] [CrossRef]
Sarychikhina, O.; Glowacka, E.; Robles, B. Multi-sensor DInSAR applied to the spatiotemporal evolution analysis of ground surface deformation in Cerro Prieto basin, Baja California, Mexico, for the 1993–2014 period. Nat. Hazards 2018, 92, 225–255. [Google Scholar] [CrossRef]
Trugman, D.T.; Borsa, A.A.; Sandwell, D.T. Did stresses from the Cerro Prieto Geothermal Field influence the El Mayor-Cucapah rupture sequence? Geophys. Res. Lett. 2014, 41, 8767–8774. [Google Scholar] [CrossRef]
Samsonov, S.V.; Feng, W.; Fialko, Y. Subsidence at Cerro Prieto Geothermal Field and postseismic slip along the Indiviso fault from 2011 to 2016 RADARSAT-2 DInSAR time series analysis. Geophys. Res. Lett. 2017, 44, 2716–2724. [Google Scholar] [CrossRef]
Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef]
Lanari, R.; Casu, F.; Manzo, M.; Zeni, G.; Berardino, P.; Manunta, M.; Pepe, A. An Overview of the Small BAseline Subset Algorithm: A DInSAR Technique for Surface Deformation Analysis. Pure Appl. Geophys. 2007, 164, 637–661. [Google Scholar] [CrossRef]
López-Quiroz, P.; Doin, M.-P.; Tupin, F.; Briole, P.; Nicolas, J.-M. Time series analysis of Mexico City subsidence constrained by radar interferometry. J. Appl. Geophys. 2009, 69, 1–15. [Google Scholar] [CrossRef]
Tong, X.; Schmidt, D. Active movement of the Cascade landslide complex in Washington from a coherence-based InSAR time series method. Remote Sens. Environ. 2016, 186, 405–415. [Google Scholar] [CrossRef]
Asaka, T.; Nonaka, T.; Hashiba, H.; Iwashita, K.; Sugimura, T. ALOS/PALSAR SBAS Analysis of Surface Deformation in the Pacific Side of Chiba Prefecture, Japan. Trans. Jpn. Soc. Aeronaut. Space Sci. Aerosp. Technol. Jpn. 2018, 16, 593–598. [Google Scholar] [CrossRef]
Han, Y.; Liu, G.; Liu, J.; Yang, J.; Xie, X.; Yan, W.; Zhang, W. Monitoring and Analysis of Land Subsidence in Jiaozuo City (China) Based on SBAS-InSAR Technology. Sustainability 2023, 15, 11737. [Google Scholar] [CrossRef]
Lomnitz, C.; Mooser, F.; Allen, C.R.; Brune, J.N.; Tatcher, W. Sismicidad y tectónica de la región norte del Golfo de California, México. resultados preliminares. Geofísica Int. 1970, 10, 37–48. [Google Scholar] [CrossRef]
Elders, W.A.; Bird, D.K.; Williams, A.E.; Schiffman, P. Hydrothermal flow regime and magmatic heat source of the Cerro Prieto geothermal system, Baja California, Mexico. Geothermics 1984, 13, 27–47. [Google Scholar] [CrossRef]
Suárez-Vidal, F.; Mendoza-Borunda, R.; Nafarrete-Zamarripa, L.M.; Rámirez, J.; Glowacka, E. Shape and Dimensions of the Cerro Prieto Pull-Apart Basin, Mexicali, Baja California, Mexico, Based on the Regional Seismic Record and Surface Structures. Int. Geol. Rev. 2008, 50, 636–649. [Google Scholar] [CrossRef]
Argus, D.F.; Gordon, R.G.; Heflin, M.B.; Ma, C.; Eanes, R.J.; Willis, P.; Peltier, W.R.; Owen, S.E. The angular velocities of the plates and the velocity of Earth’s centre from space geodesy. Geophys. J. Int. 2010, 180, 913–960. [Google Scholar] [CrossRef]
DeMets, C.; Gordon, R.G.; Argus, D.F. Geologically current plate motions. Geophys. J. Int. 2010, 181, 1–80. [Google Scholar] [CrossRef]
González-Ortega, J.A.; González-García, J.J.; Sandwell, D.T. Interseismic Velocity Field and Seismic Moment Release in Northern Baja California, Mexico. Seismol. Res. Lett. 2018, 89, 526–533. [Google Scholar] [CrossRef]
Glowacka, E.; Sarychikhina, O.; Nava, F.A. Subsidence and Stress Change in the Cerro Prieto Geothermal Field, B.C., Mexico. Pure Appl. Geophys. 2005, 162, 2095–2110. [Google Scholar] [CrossRef]
Truesdell, A.H.; Thompson, J.M.; Coplen, T.B.; Nehring, N.L.; Janik, C.J. The origin of the Cerro Prieto geothermal brine. Geothermics 1981, 10, 225–238. [Google Scholar] [CrossRef]
Gutiérrez-Negrín, L.C.; Maya-González, R.; Luis Quijano-León, J. Present Situation and Perspectives of Geothermal in Mexico. In Proceedings of the Proceedings World Geothermal Congress, Melbourne, Australia, 19–25 April 2015; pp. 19–25. Available online: https://www.geothermal-energy.org/pdf/IGAstandard/WGC/2015/01002.pdf (accessed on 11 September 2024).
Gutiérrez-Negrín, L.C.; Romo-Jones, J.M.; Izquierdo-Montalvo, G.; Canchola-Féliz, I. 2021 Annual Report IEA Geothermal. 2022. Available online: https://drive.google.com/file/d/1FuRq-r8rcVlLmd3nGKUOSuMqw91cS0N3/view?usp=share_link (accessed on 6 December 2023).
Albores, A.; Reyes, A.; Brune, J.N.; Gonzalez, J.; Garcilazo, L.; Suarez, F. Seismicity studies in the region of the Cerro Prieto geothermal field. Geothermics 1980, 9, 65–77. [Google Scholar] [CrossRef]
Frez, J.; González, J.J. Crustal Structure and Seismotectonics of Northern Baja California. In The Gulf and Peninsular Province of the Californias; American Association of Petroleum Geologists: Tulsa, OK, USA, 1991. [Google Scholar] [CrossRef]
Hauksson, E.; Stock, J.; Hutton, K.; Yang, W.; Vidal-Villegas, J.A.; Kanamori, H. The 2010 Mw 7.2 El mayor-cucapah earthquake sequence, Baja California, Mexico and Southernmost California, USA: Active seismotectonics along the Mexican pacific margin. Pure Appl. Geophys. 2010, 168, 1255–1277. [Google Scholar] [CrossRef]
Wei, S.; Fielding, E.; Leprince, S.; Sladen, A.; Avouac, J.; Helmberger, D.; Hauksson, E.; Chu, R.; Simons, M.; Hudnut, K.; et al. Superficial simplicity of the 2010 El Mayor–Cucapah earthquake of Baja California in Mexico. Nat. Geosci. 2011, 4, 615–618. [Google Scholar] [CrossRef]
Fletcher, J.M.; Teran, O.J.; Rockwell, T.K.; Oskin, M.E.; Hudnut, K.W.; Mueller, K.J.; Spelz, R.M.; Akciz, S.O.; Masana, E.; Faneros, G.; et al. Assembly of a large earthquake from a complex fault system: Surface rupture kinematics of the 4 April 2010 El Mayor–Cucapah (Mexico) Mw 7.2 earthquake. Geosphere 2014, 10, 797–827. [Google Scholar] [CrossRef]
Gonzalez-Ortega, A.; Fialko, Y.; Sandwell, D.; Alejandro Nava-Pichardo, F.; Fletcher, J.; Gonzalez-Garcia, J.; Lipovsky, B.; Floyd, M.; Funning, G. El Mayor-Cucapah (Mw 7.2) earthquake: Early near-field postseismic deformation from InSAR and GPS observations. J. Geophys. Res. Solid Earth 2014, 119, 1482–1497. [Google Scholar] [CrossRef]
Majer, E.L.; McEvilly, T.V.; Albores, A.; Diaz, S.C. Seismological studies at Cerro Prieto. Geothermics 1980, 9, 79–88. [Google Scholar] [CrossRef]
Fabriol, H.; Munguía, L. Seismic activity at the Cerro Prieto Geothermal Area (Mexico) from August 1994 to December 1995, and its relationship with tectonics and fluid exploitation. Geophys. Res. Lett. 1997, 24, 1807–1810. [Google Scholar] [CrossRef]
Gonzalez, M. Two Mw 4.8 Cerro Prieto, Baja California, Mexico, Earthquakes on 1 June and 10 September 1999: Strong-Motion Observations. Bull. Seismol. Soc. Am. 2001, 91, 1456–1470. [Google Scholar] [CrossRef]
Majer, E.L.; McEvilly, T.V. Seismological studies at Cerro Prieto field: 1978–1982. In Proceedings of the Fourth Symp. on the Cerro Prieto Geothermal Field, Baja California, Mexico, 10–12 August 1982; pp. 145–151. Available online: https://www.geothermal-library.org/index.php?mode=pubs&action=view&record=1005959 (accessed on 11 September 2024).
Fabriol, H.; Glowacka, E. Seismicity and fluid reinjection at Cerro Prieto Geothermal Field: Preliminary Results. In Proceedings of the Twenty-Second Workshop on Geothermal Reservoir Engineering, Stanford University, Stanford, CA, USA, 27–29 January 1997; pp. 11–17. Available online: https://pangea.stanford.edu/ERE/pdf/IGAstandard/SGW/1997/Fabriol.pdf (accessed on 11 September 2024).
Mellors, R.J. Comparison of Four Moderate-Size Earthquakes in Southern California Using Seismology and InSAR. Bull. Seismol. Soc. Am. 2004, 94, 2004–2014. [Google Scholar] [CrossRef]
Dawson, J.; Tregoning, P. Uncertainty analysis of earthquake source parameters determined from InSAR: A simulation study. J. Geophys. Res. Solid Earth 2007, 112, B09406. [Google Scholar] [CrossRef]
Lohman, R.B.; Simons, M. Locations of selected small earthquakes in the Zagros mountains. Geochem. Geophys. Geosystems 2005, 6, Q03001. [Google Scholar] [CrossRef]
Dawson, J.; Cummins, P.; Tregoning, P.; Leonard, M. Shallow intraplate earthquakes in Western Australia observed by Interferometric Synthetic Aperture Radar. J. Geophys. Res. Solid Earth 2008, 113, B11408. [Google Scholar] [CrossRef]
Funning, G.J.; Garcia, A. A systematic study of earthquake detectability using Sentinel-1 Interferometric Wide-Swath data. Geophys. J. Int. 2018, 216, 332–349. [Google Scholar] [CrossRef]
Massonnet, D.; Feigl, K.L. Radar interferometry and its application to changes in the Earth’s surface. Rev. Geophys. 1998, 36, 441–500. [Google Scholar] [CrossRef]
Rosen, P.A.; Gurrola, E.; Sacco, G.F.; Zebker, H. The InSAR scientific computing environment. In Proceedings of the EUSAR 2012; 9th European Conference on Synthetic Aperture Radar, Nuremberg, Germany, 23–26 April 2012; pp. 730–733. Available online: https://ieeexplore.ieee.org/document/6217174 (accessed on 11 September 2024).
Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The Shuttle Radar Topography Mission. Rev. Geophys. 2007, 45, RG2004. [Google Scholar] [CrossRef]
Goldstein, R.M.; Werner, C.L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. 1998, 25, 4035–4038. [Google Scholar] [CrossRef]
Chen, C.W.; Zebker, H.A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization. J. Opt. Soc. Am. A 2001, 18, 338. [Google Scholar] [CrossRef] [PubMed]
Agram, P.S.; Jolivet, R.; Simons, M. The Generic InSAR Analysis Toolbox; User Guide; Caltech: Pasadena, CA, USA, 2012; p. 100. Available online: https://ia601802.us.archive.org/34/items/manualzilla-id-5636918/5636918.pdf (accessed on 23 September 2024).
Hensley, S.; Zebker, H.; Jones, C.; Michel, T.; Muellerschoen, R.; Chapman, B. First deformation results using the NASA/JPL UAVSAR instrument. In Proceedings of the 2009 2nd Asian-Pacific Conference on Synthetic Aperture Radar, Xi’an, China, 26–30 October 2009; pp. 1051–1055. [Google Scholar] [CrossRef]
Yu, C.; Li, Z.; Penna, N.T. Interferometric synthetic aperture radar atmospheric correction using a GPS-based iterative tropospheric decomposition model. Remote Sens. Environ. 2018, 204, 109–121. [Google Scholar] [CrossRef]
Jolivet, R.; Grandin, R.; Lasserre, C.; Doin, M.P.; Peltzer, G. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. Geophys. Res. Lett. 2011, 38, L17311. [Google Scholar] [CrossRef]
Donnellan, A.; Parker, J.; Hensley, S.; Pierce, M.; Wang, J.; Rundle, J. UAVSAR observations of triggered slip on the Imperial, Superstition Hills, and East Elmore Ranch Faults associated with the 2010 M 7.2 El Mayor-Cucapah earthquake. Geochem. Geophys. Geosystems 2014, 15, 815–829. [Google Scholar] [CrossRef]
Lin, J.; Stein, R.S. Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults. J. Geophys. Res. Solid Earth 2004, 109, B02303. [Google Scholar] [CrossRef]
Toda, S.; Stein, R.S.; Richards-Dinger, K.; Bozkurt, S.B. Forecasting the evolution of seismicity in southern California: Animations built on earthquake stress transfer. J. Geophys. Res. Solid Earth 2005, 110, B05S16. [Google Scholar] [CrossRef]
Yang, X.-M.; Davis, P.M. Deformation due to a rectangular tension crack in an elastic half-space. Bull. Seismol. Soc. Am. 1986, 76, 865–881. [Google Scholar] [CrossRef]
CICESE. Red Geodésica del Noroeste de México (REGNOM). Centro de Investigación Científica y de Educación Superior de Ensenada. 2019. Available online: http://regnom.cicese.mx/index.html (accessed on 11 December 2023).
Flores-Armenta, M. Geothermal activity and development in Mexico–keeping the production going. In Short Course on Geothermal Development and Geothermal Wells; 2012; Available online: https://api.semanticscholar.org/CorpusID:134565263 (accessed on 11 September 2024).
Romo-Jones, J.M.; Gutiérrez-Negrín, L.C.; Sánchez-Cornejo, C.; González-Alcántar, N.; García-Gutiérrez, A. 2017 Mexico Country Report. IEA Geothermal. 2018, p. 10. Available online: https://drive.google.com/file/d/1uzwpWtRXY8Pblq9lU0cWMlTmFpwxa18t/view (accessed on 11 September 2024).
IEA Geothermal. Available online: https://iea-gia.org/publications-2/annual-reports/ (accessed on 6 December 2023).
FLores-Armenta, M.; Remírez-Montes, M.; Lilibeth, M.-A. Geothermal activity and development in Mexico-keeping the production going. In Proceedings of the Short course VI Utilization of Low- and Medium-Enthalpy Geothermal Resources and Financial Aspects of Utilization, Santa Tecla, El Salvador, 23–29 March 2014; Available online: https://gogn.orkustofnun.is/unu-gtp-sc/UNU-GTP-SC-18-03.pdf (accessed on 11 September 2024).
Flores-Armenta, M.; Gutiérrez-Negrín, L.; Nieva, D. 2015 Mexico Country Report. IEA Geothermal. 2015. Available online: https://drive.google.com/file/d/1CjpI-_NQDX1a2iWlnbQH6LgM6fIV_HnV/view (accessed on 11 September 2024).
Wessel, P.; Luis, J.F.; Uieda, L.; Scharroo, R.; Wobbe, F.; Smith, W.H.F.; Tian, D. The Generic Mapping Tools Version 6. Geochem. Geophys. Geosystems 2019, 20, 5556–5564. [Google Scholar] [CrossRef]

Figure 1. Study area. (a) Location of the region described in (b). (b) Red square indicates the CPSO area described in (c). Footprints of Sentinel 1A ascending and descending SAR images are denoted by big squares (black and dark blue), and footprints of UAVSAR east and west passes are denoted by rectangles (gray). Black arrows denote the different sensors’ line of sight direction. Red lines are the main fault systems. (c) Gray contours show the subsidence displacement rate (cm/yr) from leveling measurements (2012–2015) surveyed by the Mexican Institute of Water Technology (IMTA). Recent earthquakes

M_{L}

> 5 occurred before the study period are denoted by red stars. 1.-

M_{L}

5.4, May/2006; 2.-

M_{L}

5.3, September/2009; 3.-

M_{L}

6.0, December/2009. Main tectonic faults are indicated by continuous red lines. Large-scale tectonic motion is shown by black arrows. The CPGF area is marked as dashed gray lines. Abbreviations NA = North American plate; PA = Pacific plate; EZ = exploitation zone of the CPGF, RZ = recharge zone, CPF = Cerro Prieto Fault, GF = Guerrero Fault, HF = Hidalgo Fault, LF = L Fault, IMF = Imperial Fault, MF = Morelia Fault, SF = Saltillo Fault, SF’ = Saltillo Fault continuation, CPV = Cerro Prieto Volcano. Vector data of the CPGF limits were taken from [19].

M_{L}

> 5 occurred before the study period are denoted by red stars. 1.-

M_{L}

5.4, May/2006; 2.-

M_{L}

5.3, September/2009; 3.-

M_{L}

Figure 2. Timeframe of imagery from spaceborne Sentinel 1A/RADARSAT-2 and airborne UAVSAR missions. Abbreviations: SLC = Single Look Complex, T166 = Sentinel ascending orbital pass, T173 = Sentinel descending orbital pass, MF1 = RADARSAT-2 ascending orbital pass, MF4N = RADARSAT-2 descending orbital pass, 08514S3 = Segment #3 of the flight line 08514, 26515S2 = Segment #2 of the flight line 26515. Black dots are the UAVSAR acquisition times. Gray boxes indicate temporal matching between sensors used for 3D decomposition. Dates format: YYYY/MM/DD (e.g., 1 February 2012).

Figure 3. Method flowchart for InSAR data processing. ECMWF global model is the European Centre for Medium-Range Weather Forecasts. ¹ Jackknife test for uncertainty estimations. This workflow was elaborated based on the implemented processing software (² and ³).

Figure 4. Synthetic data of the components of the 3D displacement vector of the CPGF [17]. (a) Synthetic data. (b) Calculated 3D surface displacement data. In (a,b), colors denote the vertical displacement and red color vectors represent the horizontal displacements (east-north). (c) Differences between synthetic and calculated 3D surface displacement data. (d) Flowchart for 3D inversion code validation by using synthetic data. LOSdisp = Line of Sight displacement, WLSS = Weighted Least Squares Solution.

Figure 5. Maps of average LOS displacement rate (mm/yr). (a,b) RADARSAT-2 ascending and descending orbital passes, respectively. Stable reference point used in (a,b) is located to the northwest of the Cerro Prieto basin out of the map’s data frame [21]. (c,d) UAVSAR east and west flight segments, respectively. Maps cover 2.8 years. Areas with a coherence value below 0.27 are masked. The red flag in (c,d) shows the location of the reference point. The color palette corresponds to the LOS displacement rate. Black arrows denote the sensors’ line of sight direction. Main faults are denoted by continuous red lines. Abbreviations Ifg = interferogram, LOS = line of sight, CPF = Cerro Prieto Fault, IMF = Imperial Fault, MF = Morelia Fault, SF = Saltillo, SF’ = Saltillo Fault continuation [30], and CPV = Cerro Prieto Volcano.

Figure 6. Maps of average LOS displacement rate (mm/yr). (a,b) Sentinel 1A ascending and descending orbital passes, respectively. (c,d) UAVSAR east and west flight segments, respectively. Maps (a–d) cover one year. Areas with a coherence value below 0.2 are masked. In (a), the orange squares mark the location of specific points in the exploitation (EZ) and recharge (RZ) zones. The red flag shows the location of the reference point. Notation as in Figure 5.

Figure 7. Maps of displacement vector components, derived from the RADARSAT-2 and UAVSAR datasets’ combination for the February/2012–November/2014 period. (a) Vertical displacement rate. Negative values indicate subsidence. (b) North displacement rate. Negative values indicate southward movement. (c) East displacement rate. Negative values indicate westward movement. Areas of low coherence (<0.27) are masked out. In (a), the orange squares mark the location of specific points in the exploitation (EZ) and recharge (RZ) zones. Black dots a and b indicate the location of MBIG and NVLX GPS sites, respectively. SAR stands for Synthetic Aperture Radar. ADEW stands for Ascending/Descending/East/North, which refers to the flight direction combination of the different SAR geometries. Notation as in Figure 5.

Figure 8. Maps of displacement vector components, derived from the Sentinel 1A and UAVSAR datasets’ combination for the April/2015–April/2016 period. (a) Vertical displacement rate. Negative values indicate subsidence. (b) North displacement rate. Negative values indicate southward movement. (c) East displacement rate. Negative values indicate westward movement. Areas of low coherence (<0.2) and errors (<20 mm/yr) are masked out. In (a), the orange squares mark the location of the exploitation (EZ) and recharge (RZ) zones. The green inverted triangle marks the Ejido Nuevo León location and black dots a and b indicate the location of the MBIG and NVLX GPS sites, respectively. SAR stands for Synthetic Aperture Radar. ADEW stands for Ascending/Descending/East/North, which refers to the flight direction combination of the different SAR geometries. Notation as in Figure 5.

Figure 9. Contour maps of vertical displacement rate (cm/yr) from (a) leveling measurements (2012–2015) and (b) 3D displacement vector decomposition (2012–2014); (c) contour map of the difference (residual) between leveling measurements (IMTA) and InSAR vertical displacement rate. Blue triangles are benchmarks used for interpolation and contouring. In (a,b), the contours are every 1 cm, and in (c), are every 0.4 cm. In (a,b), RP is the reference area centered at the “10037” benchmark location and is represented by a black tringle. Tectonic faults are shown as continuous red lines. The red flag in (b,c) shows the location of the InSAR reference point. InSAR stands for Interferometric Aperture Radar. Notation as in Figure 5.

Figure 10. (a) Total RMSE and difference histogram of vertical displacement rate (cm/yr); (b) correlation coefficient between leveling data (2012–2015) and vertical InSAR (February/2012–Novemeber/2014).

Figure 11. (a) The vertical displacement rate (mm/yr) obtained here vs. other works [13,19,67] in the exploitation and recharge zones in the CPSO. See Figure 7a and Figure 8a for zones’ location. (b) Production and injection wells (number of wells) vs. total electricity generated in the CPGF. In (a), continuous and dashed lines represent a location in the exploitation and recharge zones, respectively. Texts in parentheses denote the works of other authors and are associated with a color for clarity.

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 ⁴	$α_{h}$ (°)	$θ_{i n c}$ (°)	N	M
Sentinel 1A ¹	T166	asc	3 April 2015–9 April 2016	347.39	37.2	26	49
Sentinel 1A ¹	T173	dsc	3 April 2015–9 April 2016	192.7	34.85	26	49
UAVSAR ²	08514S3	east	1 February 2012–1 April 2016	84.8	50.4	12	20
UAVSAR ²	26515S2	west	1 February 2012–1 April 2016	−94.8	43.85	12	20
RADARSAT-2 ³	MF1	asc	13 September 2011–24 July 2016	349.4	38.4	53	434
RADARSAT-2 ³	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.

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	$θ_{i n c}$ (°) ¹	$α_{h}$ (°)	d_u	d_n	d_e
Sentinel 1A	35	349	0.82	0.11	0.56
UAVSAR	45 ²	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].

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

Location	Lat	Long	North ²	East ²	Up ²	RMSE ² 3D	North ³	East ³	Up ³	RMSE ³ 3D
MBIG ¹	32.4089	−115.1960	−23.3	7.9	−96.7	6.8	−6.1	7.0	−90.0	6.5
NVLX ¹	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.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Garcia-Meza, I.F.; González-Ortega, J.A.; Sarychikhina, O.; Fielding, E.J.; Samsonov, S. 3D Surface Velocity Field Inferred from SAR Interferometry: Cerro Prieto Step-Over, Mexico, Case Study. Remote Sens. 2024, 16, 3788. https://doi.org/10.3390/rs16203788

AMA Style

Garcia-Meza IF, González-Ortega JA, Sarychikhina O, Fielding EJ, Samsonov S. 3D Surface Velocity Field Inferred from SAR Interferometry: Cerro Prieto Step-Over, Mexico, Case Study. Remote Sensing. 2024; 16(20):3788. https://doi.org/10.3390/rs16203788

Chicago/Turabian Style

Garcia-Meza, Ignacio F., J. Alejandro González-Ortega, Olga Sarychikhina, Eric J. Fielding, and Sergey Samsonov. 2024. "3D Surface Velocity Field Inferred from SAR Interferometry: Cerro Prieto Step-Over, Mexico, Case Study" Remote Sensing 16, no. 20: 3788. https://doi.org/10.3390/rs16203788

APA Style

Garcia-Meza, I. F., González-Ortega, J. A., Sarychikhina, O., Fielding, E. J., & Samsonov, S. (2024). 3D Surface Velocity Field Inferred from SAR Interferometry: Cerro Prieto Step-Over, Mexico, Case Study. Remote Sensing, 16(20), 3788. https://doi.org/10.3390/rs16203788

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu