Spatio-temporal evolution of aseismic ground deformation in the Mexicali Valley ( Baja California , Mexico ) from 1993 to 2010 , using differential SAR interferometry

Ground deformation in Mexicali Valley, Baja California, Mexico, the southern part of the MexicaliImperial valley, is influenced by active tectonics and human activity, mainly that of geothermal fluid extraction in the Cerro Prieto Geothermal Field. Significant ground deformation, mainly subsidence (∼ 18 cm yr), and related ground fissures cause severe damage to local infrastructure. The technique of Differential Synthetic Aperture Radar Interferometry (DInSAR) has been demonstrated to be a very effective remote sensing tool for accurately measuring the spatial and temporal evolution of ground displacements over broad areas. In present study ERS-1/2 SAR and ENVISAT ASAR images acquired between 1993 and 2010 were used to perform a historical analysis of aseismic ground deformation in Mexicali Valley, in an attempt to evaluate its spatio-temporal evolution and improve the understanding of its dynamic. For this purpose, the conventional 2-pass DInSAR was used to generate interferograms which were used in stacking procedure to produce maps of annual aseismic ground deformation rates for different periods. Differential interferograms that included strong co-seismic deformation signals were not included in the stacking and analysis. The changes in the ground deformation pattern and rate were identified. The main changes occur between 2000 and 2005 and include increasing deformation rate in the recharge zone and decreasing deformation rate in the western part of the CPGF production zone. We suggested that these changes are mainly caused by production development in the Cerro Prieto Geothermal Field.


Introduction
Any process involving large-scale extraction of groundwater, including geothermal energy development, can be accompanied by ground deformation.Ground deformation related to geothermal energy production has been observed and studied in many geothermal fields all over the world (e.g.Mossop and Segall, 1997;Carnec and Fabriol, 1999;Fialko and Simons, 2000;Hole et al., 2007;Allis et al., 2009;Keiding et al., 2009).Ground deformation associated with geothermal energy production can be reduced or controled through monitoring combined with aquifer management (Lowe, 2012).
Differential synthetic aperture radar interferometry (DIn-SAR) is a space-based geodetic technique widely applied for large-scale monitoring (spatially continuous) of ground deformation at a low cost (Massonnet and Rabaute, 1993;Zebker et al., 1994;Massonnet and Feigl, 1998), and particularly for subsidence regardless of its cause (Massonnet et al., 1997;Fielding et al., 1998;Avallone et al., 1999;Carnec and Fabriol, 1999;Wright and Stow, 1999;Carnec and Delacourt, 2000).Compared with ground-based geodetic techniques, such as precise leveling and GPS surveys, which provide accurate and precise information but only at a number of measuring points in a deforming surface, DInSAR allows detecting, measuring and monitoring ground deformation over large areas (∼ 10 4 km 2 ) with fine spatial resolution (∼ 10 2 m 2 ) and high accuracy (∼ 1 cm) (Gabriel et al., 1989;Bürgmann et al., 2000;Hanssen, 2001).The systematic cov-erage and frequent observations (e.g.approximately monthly repeat-cycle for ERS and ASAR) provided by SAR sensors allows the extraction of information on the spatial and temporal evolution of ground deformation at more frequent intervals than traditional techniques (Massonnet and Feigl, 1998;Baer et al., 2002;Tomás et al., 2005;Herrera et al., 2007).The availability since 1992 of a large archive of satellite acquisitions allows an analysis of historical deformation.
In the present study ERS-1/2 SAR and ENVISAT ASAR images acquired between 1993 and 2010 were used to perform an analysis of historical ground deformation caused by geothermal energy development in the Cerro Prieto Geothermal Field (CPGF) in Mexicali Valley, in an attempt to evaluate its spatio-temporal evolution and improve the understanding of its dynamic.

Study area
Ground deformation in Mexicali Valley, Baja California, Mexico, the southern part of the Mexicali-Imperial Valley, is influenced by active tectonics and human activity.Mexicali Valley is located in a complex tectonic environment, at the tectonic boundary between the Pacific and North American plates.This area is characterized by high tectonic seismicity, heat flow and surface deformation, related to the tectonic regime of the area.Besides the tectonic deformation, extraction of fluids in the CPGF which is the oldest and largest Mexican geothermal field in operation and the second largest in the world produces deformation of large magnitude (Glowacka et al., 1999).Significant ground deformation (mainly subsidence) and related ground fissures cause severe damage to local infrastructure such as roads, irrigation canals and other facilities.A detailed plan-view map of the study area is presented in Fig. 1.
Ground deformation in the Mexicali Valley has also been measured and monitored using the DInSAR technique (Carnec and Fabriol, 1999;Hanssen, 2001;Sarychikhina et al., 2011Sarychikhina et al., , 2015)).Carnec and Fabriol (1999) and Hanssen (2001) used ERS1/2 images acquired during 1993-97 and 1995-97, respectively.Significant local subsidence has been detected.Due to the high degree of temporal and spatial decorrelation, reliable data was obtained only over the mainly desert area of CPGF.Although agricultural activity in the area limited the investigation, interferometric monitoring revealed that the ground deformation is associated with the geothermal fluid withdrawal and was in agreement with the leveling data.More recent studies by Sarychikhina et al. (2011Sarychikhina et al. ( , 2014Sarychikhina et al. ( , 2015) ) and Sarychikhina and Glowacka (2015) used Envisat ASAR images acquired during 2003-10 to evaluate the spatial pattern and rate of land subsidence.The Differential Interferograms Stacking (DIS) procedure was applied to estimate the annual ground deformation rate for different periods.The results were compared to 1994-97 and 1997-2006  Here we extend these studies by using aditional SAR images (ERS 1/2 SAR) and the DIS procedure to obtain ground deformation rate for different periods from 1993 to 2009.

Methodology for deformation rate estimation
The conventional 2-pass DInSAR method using ERS 1/2 SAR and ENVISAT ASAR images was used to produce interferograms.The DIS procedure was used to estimate ground deformation rates in the study area for different time periods.Here, the DIS procedure, which consisted in averaging of multiple differential interferograms, has been applied in order to overcome limitations of conventional DInSAR, such as low coherence over long temporal separations and phase distortions introduced by atmospheric effects.
Data processing was performed using the GAMMA commercial software, developed by Wegmüller and Werner (1997) by first performing two-pass interferometry and unwrapping of differential interferograms, and then the stacking of multiple differential interferograms and geocoding.The selection of interferograms for stacking was based on their quality in terms of high coherence and low noise level.
Unfortunately, it was impossible to form the annual stacks of coherent short temporal baseline (35-140 days) differential interferograms with reasonable coherence level from ERS 1/2 dataset because of important gaps in data acquisition and geometrical decorrelation of some short time period differential interferograms.For this reason, the pairs from ERS 1/2 were grouped into two time periods (stacks): 1993-97 and 1998-2000.
Because here we attempt to present results of only aseismic deformation the differential interferograms that presented strong co-seismic deformation signals were not included in the stacking and analysis.It was impossible to generate an aseismic LOS displacement rate map for 2008 due to a high level of moderate sized seismicity during this year.

Historical deformation analysis
Figure 2 illustrates the maps of annual LOS displacement rate for different time periods obtained using the DIS method.All maps of LOS displacement rate show a roughly NE-SW oriented elliptical-shaped pattern of deformation with two bowls exhibiting high deformation rates.The similarity between maps of LOS displacement rate from ascend-ing and descending tracks (for 2007), along with the known vertical displacement from ground-based measurements including leveling surveys, suggest that the observed LOS displacements may be interpreted as reflecting mostly vertical surface displacement.However, Fig. 3, which presents the annual LOS displacement rates along the profile (A-A' in Fig. 1) across the subsidence zone, shows important differences (∼ 2 cm of magnitude) between ascending and descending datasets in the central part of the CPGF.This indicates the influence of a horizontal displacement component, probably related to horizontal displacement at and beyond the margin of the subsidence bowl, as proposed by Segall (1989).
The negative rate values in Fig. 2 indicate an increase in the LOS range that corresponds to ground subsidence.The subsiding area boundaries appear to correlate with faults, as can be seen in Fig. 2 where they are superimposed onto the maps of LOS displacement rate.The first region of subsidence, with larger negative-valued LOS displacement rate for all analysed periods (Fig. 3), is located in the area between the eastern limits of the CPGF and the Saltillo fault, which was proposed as a recharge zone in previous studies (Glowacka et al., 1999(Glowacka et al., , 2005;;Sarychikhina, 2003).The maximum LOS displacement rate in this region increased from about −12.5 to −17 cm yr −1 during 1993-97, and 2005-07 and 2009, respectively.The largest negative-valued LOS displacement rate (−17 cm yr −1 ) corresponds to about 18 cm yr −1 of ground subsidence.The area of this subsidence bowl is also evidently increasing over time.This increasing LOS displacement (subsidence) rate and areal extend of subsidence is probably related to geothermal fluid production causing continued fluid-pressure declines in the area of extraction in the CPGF that extend to the recharge zone.The 1994-97 leveling survey data zone (Glowacka et al., 1999(Glowacka et al., , 2005;;Sarychikhina et al., 2011) projected to the LOS direction shows a maximum LOS displacement rate of about −8 cm yr −1 in the recharge which is not in accord with our results.This discordance is probably due to limited coverage of the leveling network.
The second region of high subsidence rate is located in the CPGF production zone.Because the ascending and descending pass datasets present important differences in this zone, we compared the data only from the same pass.
For the descending pass data, the maximum LOS displacement rate occurring in the eastern part of the CPGF production zone is about −10 cm yr −1 for all compared periods (1993-97, 1998-2000, 2005-07).In the western and central parts of the CPGF production zone the LOS displacement rate is lower in 2005-07 than in 1993-97 and 1998-2000.It could be related to decreased production in the western part of CPGF (CP I).The LOS displacement rates for individual years during 2005-07 period are similar.
For ascending pass data, we found important differences in LOS displacement rate and pattern between the 2004 and the 2007, 2009 periods.However, due to the fact that the 2004 LOS displacement rate was estimated using only 2 differential interferograms, we think that this difference is due to possible errors (e.g.atmospheric artifacts) in the 2004 dataset.There is no significant difference between the 2007 and 2009 LOS displacement rates.

Conclusions
A comparison of the annual aseismic LOS displacement rates for 1993-97, 1998-2000, 2004-07, and 2009 periods obtained using multiple Differential Interferograms Stacking (DIS) method reveals the important changes in spatial pattern and rate of ground deformation in the study area.The main changes occur between 2000 and 2005 and in- clude increasing deformation rate in the recharge zone and decreasing deformation rate in the western part of the CPGF production zone.The observed changes may be caused by changes in production in the CPGF, an increased production in the eastern part of the field due to the newest power plant (CP IV), which started operating in 2000, and decreased production in the western part of CPGF (CP I).Since 2000 there was not considerable changes in the total CPGF production, so based on the results of this study, we can state that it took a maximum 5 years for ground deformation to stabilize after the perturbation caused by changes in the CPGF production in 2000.
During 1993-2000 and 2005-09 periods, the spatial pattern and rate of deformation could be considered constant (or without significant changes).

Figure 1 .
Figure 1.Plan-view map of the study area.The black thick dotted line frames the limits of the Cerro Prieto Geothermal Field.The borders of the evaporation pond and the Cerro Prieto volcano (CPV) are also shown (thin gray lines).A striped orange rectangle indicates the extraction area of CPIV which started the operation since 2000.Solid red lines are well-known surface traces of tectonic faults.CPF = Cerro Prieto Fault, GF = Guerrero Fault, IF = Imperial Fault, MF = Morelia Fault, SF = Saltillo fault.Dotted red lines are proposed surface fault traces based on mapped fissure zones from González et al. (1998),Glowacka et al. (2006Glowacka et al. ( , 2010a, c), c),Lira (2006) andSuárez-Vidal et al. (2007, 2008).SF' is continuation of Saltillo fault as proposed bySuárez-Vidal et al. (2008).Black square marks the location of the reference, considered stable, region.Black thick line corresponds to the profile A-A' illustrated in Fig.3.Modified fromGlowacka et al. (2010c).

Figure 2 .
Figure 2. Maps of LOS displacement rate (cm yr −1 ) showing the time evolution of the ground deformation field in and around the Cerro Prieto Geothermal Field in Mexicali Valley.Maps were obtained using the multiple Differential Interferograms Stacking (DIS) procedure.The numbers in the right lower corner indicate the number of differential interferograms used in stack.Negative displacement rates indicate an increase in the LOS range (ground subsidence).Black square marks the location of the reference, considered stable, region.Tectonic faults (solid and dotted red lines), CPGF limits (black thick dotted line), evaporation pond and CPV (thin gray lines) are superimposed onto the maps.The maps are in UTM coordinates, the marks are every 5 km.

Figure 3 .
Figure 3.Comparison between LOS displacement rates derived from multiple Differential Interferograms Stacking (DIS) procedure for different periods along the profile A-A' (location shown in Fig. 1).Solid red lines indicate location of tectonic faults which cross the profile.Faults notation as in Fig. 1.
leveling data, revealing significant changes in ground deformation pattern and rate.