Spatio–Temporal Analysis of Deformation at San Emidio Geothermal Field, Nevada, USA Between 1992 and 2010

Abstract

Although subsidence has been observed at the San Emidio geothermal field in Nevada using interferometric synthetic aperture radar since the early 1990s, the spatial extent and temporal evolution of the subsidence have not heretofore been quantified. Furthermore, the weather conditions and geographic location of San Emidio negatively affect interferometric image quality, causing low correlation amongst pairs. To address this, we introduce a new method for selecting pairs in areas of low correlation and small deformation signal using a minimum spanning tree method with a measure of image quality as the weighting criterion. We validate our pair selection approach by comparing our data products to SqueeSAR${}^{\mathrm{TM}}$ data products from a previous study at San Emidio. We also develop a deformation model which characterizes the spatial extent of subsidence at San Emidio in terms of volume change of the reservoir. After applying this deformation model to our data set of interferometric pairs, we examine the temporal relationship of the observed deformation with production and injection operations associated with geothermal power production.

1. Introduction

1.1. Interferometric Synthetic Aperture Radar

Since the late 1980s, interferometric synthetic aperture radar (InSAR) has proved itself to be a powerful geodetic technique capable of measuring surface deformation from a variety of sources. The principle of synthetic aperture radar dictates that image resolution depends on the length of the antenna and bandwidth of the system (e.g., [1,2]), allowing for image resolution on the order of a meter. The two-dimensional structure of radar data along with the pair-wise nature of interferometry allow InSAR to capture both the spatial and temporal extent of deformation with uncertainty on the order of millimeters to centimeters.

The availability of radar systems with diverse wavelengths and cadences increases the applicability of InSAR to a wide variety of geophysical studies, including monitoring geothermal processes. In particular, differential InSAR (DInSAR) has been widely used to monitor surface deformation at many geothermal sites in the Western U.S. Several studies have used observed subsidence to develop subsurface reservoir models (e.g., [3,4,5]). In addition, data sets of multiple InSAR pairs have been used to determine temporal variations in subsidence at geothermal fields (e.g., [6,7]).

As noted by Barbour et al. [8], traditional DInSAR methods may fall short in agricultural areas. In such cases, a common alternative method is to make use of both persistent and distributed scatterer processing to improve data quality and enhance signal detection in InSAR pairs (SqueeSAR${}^{\mathrm{TM}}$, [9]). While effective, this method can be computationally expensive. Recently, Xu et al. [10] analyzed Sentinel-1A data at the Cerro Prieto geothermal field, California, U.S. using a coherence-based small baseline subset method with common-point stacking to make atmospheric corrections.

1.2. San Emidio Geothermal Field, Nevada

San Emidio geothermal field is located $\sim 100$ km north of Reno, Nevada in the Basin and Range Province. The field resides in the San Emidio desert with the Fox Range to its west and the Northern Lake Range to its east [11]. The faulting regime at San Emidio is predominantly N-striking and W-dipping and consists of two major faults: the Northern Lake Range fault and the San Emidio fault [11,12]. A stress analysis has been performed, applying a “normal faulting stress regime to the San Emidio area faults, with a minimum horizontal stress direction oriented 115 [degrees]. This is consistent with the $s_{\mathrm{hmin}}$ determined through inversion of fault data by Rhodes [11]. Under these stress conditions, north–northeast striking, steeply dipping fault segments have the highest dilation tendency, while north-northeast striking 60 degrees dipping fault segments have the highest tendency to slip. Interestingly, the San Emidio geothermal field lies in an area of primarily north striking faults, which have moderate dilation tendency and moderate to low slip tendency” [13].

In the south–east portion of the region, near the operating power plant, the San Emidio fault is intersected by a NNE-striking normal fault. There is also an ENE-striking sinistral-normal fault located in the northeast portion of the region that connects two strands of the Northern Lake Range fault system [11,12]. The San Emidio fault (SEF) lies to the west of and in the footwall of the Lake Range fault (LRF). Southeast of the operating field is the topographic expression of a right step over in the LRF. The same seems to occur with the SEF, as seen in multiple geophysical data sets. The permeability created by minor dilation and high fault density within this step over is thought to facilitate deep fluid circulation at San Emidio [12]. Earlier studies demonstrated spatial correlation between surface deformation, faulting, and wellfield production and injection [12].

In 1987, a binary power plant began operating under the supervision of Empire Geothermal Power LLC in the SW region of San Emidio with a capacity of 3.6 MW [14]. The plant used shallow production wells between 30 to 100 m depth until the early 2000s. To address the issue of cooling in the shallow reservoir, new production wells were installed at depths between 500 and 700 m along the San Emidio fault and the original production wells were converted to injection wells [14]. In 2008, the ownership of the power plant changed to U.S. Geothermal Inc., who commissioned a new plant with an operating capacity of 14.7 MW in 2012 [14]. Ormat Technologies acquired U.S. Geothermal Inc. in 2018 and is now the current operator of the plant. Locations of the wells, power plant, and faults are shown on the maps in Figure 1.

The geographic location of San Emidio and its resulting weather conditions present challenges to measuring the deformation observed at San Emidio by interferometric synthetic aperture radar (InSAR) due to poor correlation. In addition, the signature of the deformation signal at San Emidio is small. A previous study by Eneva et al. [12] employs the SqueeSAR method to improve data quality and enhance signal detection in InSAR pairs at San Emidio.

In this study, we explore improving pair selection using a minimum spanning tree (MST) algorithm with a weighting criterion defined to avoid pairs with poor quality while simultaneously selecting pairs with clearer deformation signal. We establish a deformation model that spatially characterizes the subsidence signal at San Emidio. Using time-series analysis on a data set of pairs spanning from 1992 to 2010, we estimate how deformation changes over time at the field.

2. Data and Methods

2.1. Data

We select the same sets of synthetic aperture radar (SAR) scenes as used by Eneva et al. [12] to facilitate comparison between the resulting interferometric (InSAR) data products. This includes 98 epochs from the C-band Envisat (ENVI) satellite [16]: 53 from descending track 27 and 45 from ascending track 120. Also included are 38 epochs from the descending track 27 of the European Remote-Sensing (ERS) satellite missions ERS-1 and ERS-2 [17]. These epochs are listed in

Table A1. List of epochs in data set.

Epoch (YYYY–MM–DD)	Satellite	Track	Frame
1992–05–03	ERS	T27	2799
1992–06–07	ERS	T27	2799
1992–07–12	ERS	T27	2799
1992–08–16	ERS	T27	2799
1993–01–03	ERS	T27	2799
1993–041–8	ERS	T27	2799
1993–06–27	ERS	T27	2799
1993–09–05	ERS	T27	2799
1993–12–19	ERS	T27	2799
1995–04–11	ERS	T27	2799
1995–05–16	ERS	T27	2799
1995–06–20	ERS	T27	2799
1995–08–29	ERS	T27	2799
1995–11–07	ERS	T27	2799
1995–11–08	ERS	T27	2799
1995–12–12	ERS	T27	2799
1996–01–16	ERS	T27	2799
1996–01–17	ERS	T27	2799
1996–04–30	ERS	T27	2799
1996–05–01	ERS	T27	2799
1996–08–14	ERS	T27	2799
1996–10–23	ERS	T27	2799
1996–11–27	ERS	T27	2799
1997–03–12	ERS	T27	2799
1997–07–30	ERS	T27	2799
1997–10–08	ERS	T27	2799
1998–05–06	ERS	T27	2799
1998–06–10	ERS	T27	2799
1999–01–06	ERS	T27	2799
1999–08–04	ERS	T27	2799
1999–09–08	ERS	T27	2799
2000–04–05	ERS	T27	2799
2000–05–10	ERS	T27	2799
2000–06–14	ERS	T27	2799
2000–07–19	ERS	T27	2799
2000–08–23	ERS	T27	2799
2000–09–27	ERS	T27	2799
2000–12–06	ERS	T27	2799
2001–01–10	ERS	T27	2799
2003–10–29	ENVI	T120	801
2003–12–03	ENVI	T120	801
2004–06–23	ENVI	T27	2799
2004–06–30	ENVI	T120	801
2004–09–01	ENVI	T27	2799
2004–09–08	ENVI	T120	801
2004–10–06	ENVI	T27	2799
2004–10–13	ENVI	T120	801
2004–11–10	ENVI	T27	2799
2004–11–17	ENVI	T120	801
2004–12–22	ENVI	T120	801
2005–01–19	ENVI	T27	2799
2005–02–23	ENVI	T27	2799
2005–03–02	ENVI	T120	801
2005–03–30	ENVI	T27	2799
2005–05–04	ENVI	T27	2799
2005–05–11	ENVI	T120	801
2005–06–08	ENVI	T27	2799
2005–06–15	ENVI	T120	801
2005–07–13	ENVI	T27	2799
2005–07–20	ENVI	T120	801
2005–08–17	ENVI	T27	2799
2005–09–21	ENVI	T27	2799
2005–10–26	ENVI	T27	2799
2005–11–30	ENVI	T27	2799
2005–12–07	ENVI	T120	801
2006–01–11	ENVI	T120	801
2006–02–08	ENVI	T27	2799
2006–02–15	ENVI	T120	801
2006–03–15	ENVI	T27	2799
2006–04–19	ENVI	T27	2799
2006–04–26	ENVI	T120	801
2006–05–24	ENVI	T27	2799
2006–05–31	ENVI	T120	801
2006–06–28	ENVI	T27	2799
2006–08–02	ENVI	T27	2799
2006–09–06	ENVI	T27	2799
2006–10–11	ENVI	T27	2799
2006–10–18	ENVI	T120	801
2006–11–15	ENVI	T27	2799
2006–11–22	ENVI	T120	801
2006–12–20	ENVI	T27	2799
2006–12–27	ENVI	T120	801
2007–01–24	ENVI	T27	2799
2007–02–28	ENVI	T27	2799
2007–03–07	ENVI	T120	801
2007–04–04	ENVI	T27	2799
2007–05–09	ENVI	T27	2799
2007–06–13	ENVI	T27	2799
2007–07–18	ENVI	T27	2799
2007–08–22	ENVI	T27	2799
2007–10–03	ENVI	T120	801
2007–10–31	ENVI	T27	2799
2007–11–07	ENVI	T120	801
2007–12–12	ENVI	T120	801
2008–01–09	ENVI	T27	2799
2008–01–16	ENVI	T120	801
2008–02–13	ENVI	T27	2799
2008–02–20	ENVI	T120	801
2008–03–19	ENVI	T27	2799
2008–03–26	ENVI	T120	801
2008–04–23	ENVI	T27	2799
2008–04–30	ENVI	T120	801
2008–05–28	ENVI	T27	2799
2008–06–04	ENVI	T120	801
2008–07–02	ENVI	T27	2799
2008–07–09	ENVI	T120	801
2008–08–06	ENVI	T27	2799
2008–08–13	ENVI	T120	801
2008–09–10	ENVI	T27	2799
2008–09–17	ENVI	T120	801
2008–10–15	ENVI	T27	2799
2009–03–04	ENVI	T27	2799
2009–03–11	ENVI	T120	801
2009–04–08	ENVI	T27	2799
2009–04–15	ENVI	T120	801
2009–05–13	ENVI	T27	2799
2009–05–20	ENVI	T120	801
2009–06–17	ENVI	T27	2799
2009–06–24	ENVI	T120	801
2009–07–22	ENVI	T27	2799
2009–07–29	ENVI	T120	801
2009–08–26	ENVI	T27	2799
2009–09–02	ENVI	T120	801
2009–09–30	ENVI	T27	2799
2009–10–07	ENVI	T120	801
2009–11–11	ENVI	T120	801
2009–12–09	ENVI	T27	2799
2009–12–16	ENVI	T120	801
2010–01–13	ENVI	T27	2799
2010–01–20	ENVI	T120	801
2010–02–17	ENVI	T27	2799
2010–02–24	ENVI	T120	801
2010–03–24	ENVI	T27	2799
2010–03–31	ENVI	T120	801
2010–04–28	ENVI	T27	2799
2010–05–05	ENVI	T120	801
2010–06–09	ENVI	T120	801

.

Interferograms are calculated by combining pairs of compatible images using GMT5SAR [18,19], an interferometric processing software package that utilizes generic mapping tools (GMT) to create and visualize interferometric pairs. GMT5SAR allows for removal of noise by applying an adaptive (“modified Goldstein”) filter that depends on coherence [19,20,21]. An example of one such interferogram spanning 10th November 2004 to 10th September 2008 from track 27 of the Envisat satellite is shown in Figure 1. Additional examples are shown in Figure A1, Figure A2, Figure A3, Figure A4, Figure A5 and Figure A6.

Table 1. Unit vector $\widehat{\mathbf{s}}=[s_E,s_N,s_U]$ pointing from a point on the ground to the sensor aboard the satellite and satellite incidence angle $\theta$ for European Remote-Sensing (ERS) track 27 and Envisat tracks 27 and 120.

Satellite	Track	$\widehat{\mathbf{s}}$	$\mathit{\theta }$
ERS	T27	$[0.35;-0.08;0.93]$	${21}^{\circ }$
ENVI	T27	$[0.35;-0.08;0.93]$	${21}^{\circ }$
ENVI	T120	$[-0.37;-0.08;0.92]$	${22}^{\circ }$

shows the unit vector $\widehat{\mathbf{s}}$ pointing from the ground to the sensor aboard the satellite and incidence angle for each satellite.

2.2. Methods

2.2.1. Selecting and Weighting Pairs

Given each of the three satellite track’s set of epochs (at points in time), we selected a data set of pairs (spanning intervals in time). To do so, we considered image quality as measured by the amount of phase noise within an interferometric pair (e.g., [16]). Factors that contribute to quality include: (a) the orbital separation (also known as the perpendicular component $B_{\perp }$ of the baseline), (b) the time span of the pair (also known as the temporal baseline), and (c) the difference in Doppler centroid frequencies (e.g., [22,23,24]). One popular approach is to select minimum spanning tree (MST) pairs using orbital separation as a weighting criterion (e.g., [25]). This method falls short, however, when imaging mountainous areas such as northern Nevada, which receive considerable precipitation in winter months.

To address these complications, we weighted each pair by a measure of its quality [26]. We used an empirical measure of quality $\gamma (i,j)$ for a pair spanning from time $t_i$ to time $t_j$ as defined by the product of three individual components [2,22,27]:

$$\gamma (i,j)={\gamma }_g(i,j)·{\gamma }_t(i,j)·{\gamma }_D(i,j),$$

(1)

where subscripts g, t, and D represent geometric, temporal, and Doppler contributions, respectively. We defined these individual components of $\gamma (i,j)$ in a manner similar to that of Refice et al. [22] and incorporated additional complexity in the geometric and temporal components.

The geometric component is defined as:

$${\gamma }_g(i,j)=\left(1-\frac{|b_{\perp }(i,j)|}{w_b{b}_{crit}}\right)·\left[\mathcal{H}\right(w_b{b}_{crit}-|b_{\perp }(i,j)\left|\right)],$$

(2)

where $b_{crit}$ represents the critical baseline (about 1100 m for Envisat and ERS [22,28]), $w_b$ is a specified weighting value for the critical baseline (set to be $0.5$ to favor pairs with orbital baselines less than 550 m), and $\mathcal{H}$ represents the Heaviside “step” function. The Heaviside function was used to filter pairs with orbital separation greater in absolute value than $w_b{b}_{crit}$.

The temporal component is defined using an exponential decay function:

$$\begin{array}{cc}\hfill {\gamma }_t(i,j)=\exp \left(\frac{-|\Delta t|}{b_{temp}}\right)·\left[\min \right[& \min [\mathcal{H}(m_{t_i}-3.5),1-\mathcal{H}(m_{t_i}-11.5)],\hfill \\ & \min [\mathcal{H}(m_{t_j}-3.5),1-\mathcal{H}(m_{t_j}-11.5)]],\hfill \end{array}$$

(3)

where $\Delta t$ is the time span (in absolute value) of the pair in days, $b_{temp}$ is a specified temporal decay constant (set to 80 days following Refice et al. [22] to favor pairs with shorter time spans), and $m_{t_i}$ is the numerical value of the month of $t_i$. In this case, we used the Heaviside function to filter pairs with epochs occurring during winter months (mid-November through mid-March).

Finally, the Doppler component of $\gamma (i,j)$ is defined as:

$${\gamma }_D(i,j)=\left(1-\frac{|f_{Dc}\left(i\right)-f_{Dc}\left(j\right)|}{B_a}\right)·\left[\mathcal{H}\right(B_a-|f_{Dc}\left(i\right)-f_{Dc}\left(j\right)\left|\right)],$$

(4)

where $f_{Dc}$ is the Doppler centroid frequency [Hz] for each image and $B_a$ is the azimuthal bandwidth [Hz].

We normalized $\gamma$ to be between zero and one and set a MST weighting criterion of

$$w(i,j)=1-{\gamma }_{norm}(i,j).$$

(5)

The resulting subset of pairs included in

Table A2. Data set of MST pairs selected using MST methods and $\gamma (i,j)$ as a weighting criterion. First and second epochs of each pair are given in YYYY–MM–DD format, perpendicular baseline $b_{\perp }$ and temporal baseline $\Delta t$ are given in meters and days, respectively.

Epoch 1	Epoch 2	${\mathit{b}}_{\mathbf{\perp }}$	$\mathbf{\Delta }\mathit{t}$
1992–06–07	1992–07–12	−507.8	35
1992–06–07	1993–04–18	151.1	315
1992–07–12	1992–08–16	−279.5	35
1992–08–16	1993–01–03	−329.8	140
1993–04–18	1993–09–05	−188.8	140
1993–04–18	1993–12–19	407.6	245
1993–06–27	1993–09–05	502.7	70
1993–06–27	1995–04–11	13.1	653
1995–04–11	1995–05–16	−256.8	35
1995–05–16	1995–06–20	71.9	35
1995–06–20	1995–08–29	223.5	70
1995–08–29	1995–11–08	401.6	71
1995–08–29	1996–01–16	−305.5	140
1995–08–29	1996–08–14	−42.1	351
1995–11–07	1995–11–08	−197.8	1
1995–11–07	1995–12–12	99.9	35
1995–11–07	1996–05–01	173.4	176
1996–01–16	1996–01–17	−268.2	1
1996–04–30	1996–05–01	−90.4	1
1996–05–01	1996–10–23	57.8	175
1996–08–14	1997–07–30	148.9	350
1996–10–23	1996–11–27	257.5	35
1997–03–12	1997–07–30	−126.5	140
1997–07–30	1997–10–08	445.7	70
1997–10–08	1998–05–06	57.3	210
1998–05–06	1998–06–10	297.2	35
1998–05–06	1999–08–04	−447	455
1998–06–10	1999–01–06	−2326.6	210
1999–08–04	2000–04–05	−182.5	245
1999–09–08	2000–04–05	410	210
2000–04–05	2000–06–14	0.5	70
2000–05–10	2000–08–23	−434.2	105
2000–06–14	2000–07–19	−359.2	35
2000–06–14	2000–09–27	350.5	105
2000–08–23	2000–09–27	−282.1	35
2003–10–29	2004–09–08	231.9	315
2003–12–03	2004–09–08	315.6	280
2004–06–23	2004–10–06	326.7	105
2004–06–23	2004–11–10	−385.3	140
2004–09–08	2004–11–17	401.8	70
2004–10–13	2005–05–11	456.9	210
2004–11–10	2005–03–30	−206.1	140
2004–11–10	2008–09–10	48.8	1400
2004–11–17	2004–12–22	268.3	35
2004–11–17	2005–06–15	−272.4	210
2004–12–22	2005–03–02	−115.1	70
2005–01–19	2005–02–23	−38.4	35
2005–05–11	2005–07–20	201.5	70
2005–06–15	2005–07–20	−497.3	35
2005–07–20	2006–01–11	132.7	175
2005–11–30	2006–04–19	227.7	140
2006–02–15	2006–05–31	232.1	105
2006–03–15	2006–04–19	−253.9	35
2006–04–26	2007–10–03	527.7	525
2006–06–28	2006–08–02	332.7	35
2006–11–15	2007–04–04	361.7	140
2006–11–22	2006–12–27	−149.4	35
2006–11–22	2007–11–07	−475.2	350
2006–11–22	2008–08–13	−393	630
2006–12–20	2007–02–28	−19.3	70
2006–12–27	2007–03–07	−296.5	70
2007–01–24	2007–04–04	133	70
2007–05–09	2007–07–18	−7.2	70
2007–06–13	2007–07–18	7.2	35
2007–06–13	2007–08–22	234	70
2007–10–03	2007–11–07	199.3	35
2007–10–03	2007–12–12	−186.1	70
2008–01–09	2008–02–13	−208.5	35
2008–01–16	2008–03–26	121.2	70
2008–02–13	2008–04–23	186.3	70
2008–02–20	2008–04–30	148.3	70
2008–03–19	2008–04–23	−248.7	35
2008–03–26	2008–04–30	−322.1	35
2008–04–30	2008–07–09	−19.5	70
2008–05–28	2008–07–02	118	35
2008–05–28	2008–08–06	253.3	70
2008–06–04	2008–07–09	−60.1	35
2008–07–09	2008–08–13	226.6	35
2008–07–09	2008–09–17	−144.6	70
2008–09–10	2009–04–08	309.8	210
2008–09–17	2009–04–15	120.9	210
2008–10–15	2009–05–13	−30.2	210
2009–03–04	2009–06–17	−16	105
2009–03–11	2009–04–15	−200.1	35
2009–04–08	2009–05–13	250.2	35
2009–04–15	2009–05–20	57	35
2009–05–13	2009–06–17	271.2	35
2009–05–13	2009–07–22	−186.6	70
2009–05–20	2009–06–24	−86.6	35
2009–06–17	2009–08–26	−59.2	70
2009–06–24	2009–07–29	235.4	35
2009–07–22	2009–09–30	−57.6	70
2009–07–29	2009–09–02	88.4	35
2009–09–02	2009–11–11	−210.6	70
2009–09–30	2009–12–09	−339.8	70
2009–09–30	2010–04–28	196.3	210
2009–10–07	2009–11–11	154.6	35
2009–10–07	2009–12–16	−171.8	70
2009–11–11	2010–05–05	−119.7	175
2010–01–13	2010–02–17	−165.8	35
2010–01–20	2010–02–24	−149.6	35
2010–02–17	2010–04–28	161.8	70
2010–02–24	2010–03–31	227.8	35
2010–02–24	2010–05–05	−87.8	70
2010–03–24	2010–04–28	−238.7	35
2010–05–05	2010–06–09	129.4	35

optimizes the quality of the interferograms. This data set includes three trees in a minimum spanning forest (MSF) comprised of data from both the ERS and Envisat satellites, as graphed in Figure 2. These pairs are publicly available on the Geothermal Data Repository [29].

2.2.2. Deformation Modeling

We used simulated annealing to derive a model that best describes the deformation. This nonlinear approach to deformation modeling was implemented in an open-source software package named the general inversion of phase technique (GIPhT) [26,30]. For a given deformation model and a corresponding set of initial estimates and bounds for the parameters, GIPhT inverts data values from pixels of an observed interferogram using simulated annealing and outputs refined estimates of the model parameters.

We considered a simple model for the deformation observed at San Emidio. To describe the reservoir, we assumed a cuboidal model with a single sink in an elastic half space with uniform material properties. The cuboid represents a volume element with sides of width W, length L, and height H, giving an initial volume of $V_0=LWH$. This volume element is defined by three rectangular patches that are mutually orthogonal. We estimated seven model parameters, including source location (easting, northing, and elevation), cuboid dimensions (length, width, and height), and volume change $\Delta V$. We assumed a uniform Poisson’s ratio of $\nu =1/4$. To start, we chose initial estimates of the location and dimension parameters based on results from Eneva et al. [12] (their Figure 2b). In addition to the seven parameters corresponding to the deformation model, we also estimated nuisance parameters corresponding to contributions from atmospheric effects, orbital errors, and an initial offset per epoch (e.g., [30]).

We started by deriving best-fitting estimates for the dimensions and location of the cuboidal sink using the ENVI T27 stack of (unwrapped) range change rates derived by Eneva et al. [12] using SqueeSAR. To avoid complications arising from signals due to sources other than those due to production at the site (e.g., precipitation, mining), we focused on modeling only the deformation observed around the production wells as outlined in Figure 1. Once best-fitting estimates of the cuboidal reservoir’s location and dimensions were found, they were used in modeling the individual MST pairs to arrive at estimates of reservoir volume change.

2.2.3. Time-Series Analysis

To determine any temporal trends in the deformation, we performed time-series analysis on the estimated volume changes derived from the spatial deformation modeling. To handle the pair-wise nature of these volume changes, we used a graph-theoretic approach to the procedure known as temporal adjustment, which converts a series of individual, pair-wise volume changes $\Delta V_{i,j}$ into cumulative volume change at individual points $t_i$ and $t_j$ in time [26,31]. We considered several different parameterizations for the temporal function $f\left(t\right)$. For simplicity, we started with a single-rate parameterization:

$$f_1\left(t_i\right)=a_1(t_i-t_0),$$

(6)

where $t_0=1992.43$ is the start date of our InSAR data set in decimal years.

We also considered a piecewise-linear parameterization with m breaks at times $t_k$ that form $(m-1)$ intervals:

$$\begin{array}{cc}\hfill f_2\left(t_i\right)& ={\Sigma }_{k=1}^{k=m}{a}_{k}D\left(t_i\right),\mathrm{where}\hfill \\ \hfill D\left(t_i\right)& =\left\{\begin{array}{cc}0\hfill & \mathrm{if}{t}_i<t_k\hfill \\ (t_i-t_k)\hfill & \mathrm{if}{t}_k\le t_i<t_{k+1}\hfill \\ (t_{k+1}-t_k)\hfill & \mathrm{if}{t}_i\ge t_{k+1},\hfill \end{array}\right.\hfill \end{array}$$

(7)

where $a_k$ is a parameter to be estimated.

3. Results

3.1. Analysis of Data Quality

To determine the quality of our data selection, we compared the deformation measured from our MST data sets to those from Eneva et al. [12] derived by TRE ALTAMIRA using the SqueeSAR procedure [26]. We worked with the phase gradient. For each of the three sets of data, we calculated the eastward component of the phase gradient directly from wrapped phase using quadtree resampling [32]. As described by Ali and Feigl [32], this quantity is one component of the deformation gradient [33]. We used quadtree resampling [32] to smooth interferometric phases and derive the east component of the phase gradient fields to form a stack [34]. Both of these functionalities were incorporated into GIPhT. We also converted the SqueeSAR results from Eneva et al. [12] (their Figure 2b,d) into wrapped phase and performed the same quadtree resampling and gradient procedure. Finally, we compared the fields derived from each of the three SqueeSAR stacks to the corresponding gradients of pairs selected by MST. Dividing by the time interval $\Delta t$, we found a velocity gradient that can be thought of as a strain rate expressed in microstrain/yr or picostrain/s:

$$\dot{\psi }=\frac{\Delta \varphi }{\Delta X\Delta t},$$

(8)

where the quotient of the difference in phase $\Delta \varphi$ divided by the difference in easting position $\Delta X$ defines the eastward component of the phase gradient. The resulting velocity gradient data sets are summarized in

Table 2. Mean for SqueeSAR (SQR) and minimum spanning tree (MST) stack gradients for each satellite track.

Data Set	SQR Mean	MST Mean
	[picostrain/s]	[picostrain/s]
ERS T27	$-0.8\times {10}^{-3}$	$-0.5\times {10}^{-3}$
ENVI T27	$2.4\times {10}^{-3}$	$2.4\times {10}^{-3}$
ENVI T120	$-1.7\times {10}^{-3}$	$-0.8\times {10}^{-3}$

. Additional time-series analyses of the MST gradient data sets are shown in Figure A7, Figure A8 and Figure A9.

The mean and sample standard deviation for each differenced set of gradients are shown in

Table 3. Distribution parameters for differenced SqueeSAR and MST stack gradients for each satellite track.

Data Set	Sample Mean	Sample std. Deviation
	[picostrain/s]	[picostrain/s]
ERS T27	$2.9\times {10}^{-4}$	0.34
ENVI T27	$-4.2\times {10}^{-5}$	0.28
ENVI T120	$8.6\times {10}^{-4}$	0.34

in terms of strain rate. We see that none of the three sets of differences was significantly different than zero. To determine if the means of the SqueeSAR stack gradients were different from the means of the MST stack gradients, we performed a Student’s T-test. The results are summarized in

Table 4. Results of Student’s T-test for equal means (e.g., [35], p. 521). $H_0$: means are equal; $H_1$: means are significantly different; p-value: probability of rejecting $H_0$ when it is true.

Data Set	p-Value
ERS T27	0.94
ENVI T27	0.97
ENVI T120	0.82

and Figure 3. All three data sets failed to reject the null hypothesis of equal means at 95% confidence.

3.2. Deformation Modeling

The best-fitting estimates of the cuboidal reservoir model parameters derived from modeling the ENVI T27 stack of (unwrapped) range change rates from Eneva et al. [12] are shown in

Table 5. Best-fitting estimates of a single cuboidal sink model after simulated annealing on unwrapped range change rates from ENVI T27 SqueeSAR stack.

Parameter Name	Best-Fitting Estimate	Uncertainty
Centroid Easting in m	296,109	375
Centroid Northing in m	4,472,750	380
Centroid Depth in m	500	75
Cuboid Length in m	2000	500
Cuboid Width in m	200	50
Cuboid Thickness in m	500	100
Volume Change Rate (${\mathrm{m}}^3/\mathrm{yr}$)	$-4.6\times {10}^3$	$0.5\times {10}^3$

. Results are shown in Figure 4 and Figure 5. We defined the cost of the inversion as the L1 misfit of the unwrapped range change rate [6] and found it to be 0.7 mm/yr. These best fitting estimates of location and dimensions of the cuboidal reservoir were used to model all individual pairs in our MST data set in terms of wrapped phase change, resulting in estimates of reservoir volume change which were then used for subsequent time-series analysis.

3.3. Time-Series Analysis

We worked with volume change rates derived from modeling deformation in terms of wrapped phase change from the individual MST pairs using the cuboidal sink model. Starting with a single-rate parameterization (Equation (6)), we found a best-fitting estimate of volume change rate to be $(-1.1\pm 0.1)\times {10}^4{\mathrm{m}}^3/\mathrm{yr}$. We defined a dimensionless misfit $\chi$ of the model to the data as the square root of the reduced ${\chi }^2$ statistic ([36], p. 334). We found $\chi =1.6$. Results are shown in Figure 6.

We also explored the possibility that the rate changes over time (Equation (7)). We tried a three-segment piecewise-linear parameterization with a break during the gap in the data set. We found the best-fitting estimates of volume change rates for the two intervals of data coverage to be $(-1.3\pm 0.1)\times {10}^4{\mathrm{m}}^3/\mathrm{yr}$ for the time interval from 1992 to 2001 and $(-0.9\pm 0.1)\times {10}^4{\mathrm{m}}^3/\mathrm{yr}$ for the time interval from 2003 to 2010. We found a corresponding misfit of $\chi =1.5$.

To determine if the increased complexity of the piecewise-linear parameterization is justified, we used an F-test for model complexity (e.g., [35], p. 627). The results are shown in

Table 6. Results of F-tests for model complexity (e.g., [35], p. 627). $H_0$: more complex parameterization does not provide a significantly improved fit to the data over the single-rate parameterization; $H_1$: more complex parameterization provides a significantly improved fit to the data over the single-rate parameterization; $dof$: denotes degrees of freedom.

${\mathit{dof}}_1$	${\mathit{dof}}_2$	Test Value	Critical Value	Result
105	104	−10.44	1.38	fail to reject $H_0$

. We found that the added complexity of the piecewise-linear parameterization is not justified at 95% confidence. We concluded that the temporal trend of the volume change of the modeled reservoir is best explained as a constant rate.

4. Discussion

When we compared our MST data sets (selected using quality $\gamma$ as a weighting criterion) to the corresponding data sets from SqueeSAR, we found no significant difference in mean strain rates estimated from the stacks of gradients at 95% confidence in all three InSAR data sets. Thus, our modification to the MST method of pair selection introduced by Refice et al. [22] allows MST methods to be applicable to InSAR deformation analysis in areas of low correlation and small signal. This provides another viable option for InSAR time-series analysis which may be less time-intensive and computationally expensive than existing methods.

We found that the deformation at the San Emidio geothermal field is well explained by parameterizing the reservoir as a cuboidal sink (Figure 5). This sink aligns with faults in the area and spatially encompasses the majority of production and injection wells at the site. We found a best-fitting depth of 500 m, which is consistent with lower bounds of likely reservoir depths based on cross-section analyses [14,37].

When modeling the observed deformation as a function of time, we found that the best-fitting parameterization is a constant rate of volume change. This is in agreement with the line-of-sight time-series analysis by Eneva et al. [12].

To explore the possibility that the observed deformation is related to pumping at the site, we compared the cumulative values of volume change derived from temporal adjustment with records from Ormat of cumulative monthly gross production (Figure A10 and Figure A11) as reported to the State of Nevada. We normalized the estimated cumulative volume change values (Figure A11) and the observed cumulative gross production values using the statistical Z-transform (e.g., [35]). We then tested the correlation using Pearson’s test (e.g., [35], p. 599) with the null hypothesis that there is no correlation between pumping and deformation. We found a statistically strong correlation between normalized cumulative volume change and normalized cumulative gross production with a correlation coefficient of $R=>0.99$ (Figure 7). The corresponding p-value, or the probability of rejecting the null hypothesis of no correlation when it is true, is less than ${10}^{-100}$. Thus, we inferred a significant correlation between cumulative deformation and pumping.

5. Conclusions

We have developed a new method to select a good set of interferometric pairs in areas of poor interferometric image correlation using a minimum spanning tree algorithm with a seasonally- and spatially-weighted measure of quality as the weighting criterion. This selection procedure yields a temporally averaged, phase gradient rate field that is equivalent in mean to that produced by the SqueeSAR procedure, thus providing a viable alternative to selecting an optimal data set when working in areas of poor correlation with a small deformation signature.

Using nonlinear inversion, we determined that the deformation at San Emidio is well explained by parameterizing the reservoir as a cuboidal sink. Temporal analysis of the volume change rates estimated from individual interferometric pairs using this deformation model suggests a constant rate of volume change of $(-1.1\pm 0.1)\times {10}^4{\mathrm{m}}^3/\mathrm{yr}$ between 1992 and 2010, consistent with the time-series analysis in terms of line-of-sight displacement performed by Eneva et al. [12]. When comparing cumulative volume change to cumulative gross production, we found a strong, positive correlation. This suggests that deformation at San Emidio was influenced by geothermal production at the site between 1998 and 2010.