Geodetic Measurements and Numerical Models of Deformation at Coso Geothermal Field, California, USA, 2004–2016

Abstract

We measure transient deformation at Coso geothermal field using interferometric synthetic aperture radar (InSAR) data acquired between 2004 and 2016 and relative positions estimated from global positioning system (GPS) to quantify relationships between deformation and pumping. We parameterize the reservoir as a cuboidal sink and solve for best-fitting reservoir dimensions and locations before and after 2010 in accordance with sustainability efforts implemented in late 2009 at the site. Time-series analysis is performed on volume changes estimated from pairs of synthetic aperture radar (SAR) and daily GPS data. We identify decreasing pore-fluid pressure as the dominant mechanism driving the subsidence observed at Coso geothermal field. We also find a significant positive correlation between deformation and production rate.

1. Introduction

The Coso geothermal field, located near China Lake, California, is the third largest geothermal field in the United States with an installed capacity to generate ∼270 MW of electrical power. It lies within an extensional step-over between dextral faults that hosts an actively developing metamorphic core complex [1,2]. Directly beneath the field lies a partially molten magma body, located below the relatively shallow (∼5 km depth) brittle–ductile transition, which heats water migrating deep into the geothermal area following precipitation in the Sierra Nevada mountain range [3]. The buoyant hot water eventually rises via faults and fractures associated with the transtensional tectonics of the region into the reservoir comprised of highly fractured plutonic and metamorphic rocks of Mesozoic age [2,4]. Fluid temperatures within the reservoir reach ∼350 ${}^{\circ }\mathrm{C}$ at the relatively shallow depths of ∼3 km tapped by more than 80 production wells (Figure 1). Following generation of electricity and condensation, the fluids are reinjected via more than 30 injection wells around the field. The average rates of (fluid) production and reinjection between 2004 and 2016 have been $(24.30\pm 0.27)\times {10}^8$ kg/month ($924\pm 10$ kg/s) and $(11.72\pm 0.02)\times {10}^8$ kg/month ($446\pm 7$ kg/s), respectively, which, neglecting natural recharge, corresponds to a net extraction rate of ∼$12.58\times {10}^8$ kg/month (478 kg/s) [5]. In late 2009, the reservoir operators of the Coso Geothermal Plant implemented the Hay Ranch Water Project, which uses a ∼$15-\mathrm{km}$ pipeline to recharge the existing reservoir at Coso with supplemental water, thereby increasing geothermal energy production (e.g., [6,7]).

The Coso geothermal field is well known for being one of the most seismically active regions in California. Previous studies at Coso have linked such seismicity to the motion of fluids within the geothermal system and changes in local tectonic stress (e.g., [8,9,10,11,12,13]). Geodetic studies have also observed strong subsidence signals that occur in the production region near the Coso geothermal plant (e.g., [8]), leading to the suggestion that geothermal fluids are driving both the deformation and seismicity at the site.

Similar phenomena have been observed at other geothermal fields, as recently reviewed by Majer et al. [14], Brodsky and Lajoie [15], Ellsworth [16], and Mazzoldi et al. [17]. One example is Brady Hot Springs, a geothermal field in Nevada. There, Cardiff et al. [18] found evidence relating increased microseismic activity to periods of brief cessations in pumping at the site. Such behavior is explained by the adaptation of subsurface effective (normal) stress to long-term normal operations, such that “extraction of fluids inhibits fault slip by increasing the effective [normal] stress on faults; in contrast, brief pumping cessations represent times when effective stress is decreased below its long-term average, increasing the likelihood of microseismicity” [18]. Similar phenomena have been observed at the San Emidio geothermal field in Nevada as well [19].

In this study, we explore the temporal association between pumping and deformation at Coso to determine whether there is a statistical relationship between them. We also test two competing hypotheses for the geophysical processes observed at Coso: (1) Thermal contraction and/or (2) a poroelastic response to pressure changes caused by pumping. We use an updated catalog of microseismicity, pumping records from the site, as well as interferometric synthetic aperture radar (InSAR) and global positioning system (GPS) data spanning from 2004 to 2016 for our analysis.

2. Data

2.1. Interferometric Synthetic Aperture Radar (InSAR)

We analyze 38 synthetic aperture radar (SAR) images acquired between 2004 and 2016 over Coso. Of these, 34 SAR images were acquired in ascending track 349 between 2004 and 2010 by the C-band ASAR sensor aboard the Envisat satellite [20]. An additional four SAR images were acquired after 2014 in ascending track 64 by the C-band sensor aboard the Sentinel-1A satellite [21]. We use the images to create a set of interferometric pairs using an open-source InSAR processing software (GMTSAR), which utilizes generic mapping tools (GMT) to create and visualize interferometric pairs [22,23]. The topographic contribution of each interferogram is removed using a digital elevation model with 1 arcsecond posting from the Shuttle Radar Topographic Mission [24]. Unwrapping is performed using the “statistical-cost, network-flow phase-unwrapping algorithm” (Snaphu) developed by Chen and Zebker [25] to calculate the range change rate $\dot{\rho }$ in millimeters per year.

The pairs from Environmental Satellite (Envisat) are selected using a minimum spanning tree (MST) algorithm according to orbital separation (e.g., [26]). The same data set was used in a previous deformation modeling study at Coso [27]. To summarize individual interferometric pairs, we consider the field of unwrapped range change rate $\dot{\rho }(E,N)$. The mean rate between 2005 and 2010 is derived by averaging the range change rates from our 81 Envisat pairs (Figure 1), where E and N are the position coordinates easting and northing, respectively in a Universal Transverse Mercator (UTM) projection [28]. The $k^{th}$ pixel of the mean rate field $\overline{\dot{\rho }}(E,N)$ comprised of n pairs is computed as

$${\overline{\dot{\rho }}}_k(E,N)=\frac{1}{n}\sum _{i=1}^n{\dot{\rho }}_k^{\left(i\right)}(E,N).$$

(1)

We use the same averaging procedure to form a second stack, spanning 2014 to 2016, from the 10 pairs in the Sentinel-1A data set. The resulting data set of pairs from both satellites is shown in Figure 2 and is publicly available on the Geothermal Data Repository [29].

Subsidence (positive range change) as rapid as ∼30 mm/year is observed over a circular area some 3 km in radius centered on the production wells (Figure 1). This signature is consistently observed in all interferometric pairs spanning the production interval between 2004 and 2016, as shown in Figure 3. The rate and spatial extent of the deformation field are broadly consistent with InSAR data spanning 1993 to 1999 [8].

2.2. Global Positioning System (GPS)

Data acquired from continuous GPS station COSO lead to daily estimates of the three components (eastward, northward, and upward) of relative position. The station is a part of the Southern California Integrated GPS Network and is located within 5 km of the center of the identified subsidence bowl. The GPS time-series data have been analyzed using methods outlined in Blewitt et al. [32] and are publicly available [33]. We use a subset of the time series from station COSO that spans from 2004 to 2015 to analyze deformation within the deforming region. We also use data from COSO to verify the accuracy of the InSAR range change rates.

Data are also available from the campaign GPS station COSJ, located to the northwest of the deforming region. This station is a part of the MAGNET network with a publicly available time series of relative position that have been analyzed using standard procedures outlined in Blewitt et al. [32,34]. We use data from COSJ to estimate far-field deformation when verifying the accuracy of the InSAR range change rates.

Time series of relative positions from stations COSO and COSJ are shown in Figure 4 and Figure 5, respectively.

2.3. Seismic Catalog

Local and regional seismicity near the Coso geothermal field is recorded by a local borehole seismic network and analyzed by the Navy Geothermal Program Office. Kaven et al. [11] refined the initial event locations by estimating a best-fitting, one-dimensional velocity model and relocating each of the events. Resulting event hypocenters are found to have (relative) location uncertainties on the order of $300\mathrm{m}$ horizontally and $600\mathrm{m}$ vertically.

The average values of seismic velocities in this region are 5.8 and 3.5 km/s for $V_p$ and $V_s$, respectively. We use these velocities to derive an average value of Poisson’s ratio $\nu =0.21$ for our deformation analysis. Given the average $V_p$, we also estimate the density of the surrounding material to be about $2.6{\mathrm{g}/\mathrm{cm}}^3$ [35]. This is consistent with other density estimates at Coso (e.g., [2,36]).

2.4. Pumping Records

Pumping records for the monthly rates of gross injection and gross production in kilograms per month for the Coso geothermal power plant from 2005 to 2016 (Figure 6) have been published by the Division of Oil, Gas, and Geothermal Resources [5].

3. Methods

3.1. Accuracy

We start with analyzing the accuracy of unwrapped InSAR range change rates by comparison to GPS measurements. Starting with our Envisat pairs, we first convert the GPS estimates of vector position $\mathbf{u}$ relative to stable North America from both stations to range change $\mathsf{\Delta }\rho$ by taking the negative scalar product with the unit vector ${\widehat{\mathbf{s}}}^{\left(ENVI\right)}=[-0.41;-0.08;0.91]$ pointing from the pixel on the ground to the radar sensor aboard the satellite.

$${\mathsf{\Delta }\rho }^{\left(GPS\right)}=-{\widehat{\mathbf{s}}}^{\left(ENVI\right)}·{\mathbf{u}}^{\left(GPS\right)}.$$

(2)

We then linearly interpolate and extrapolate $\mathsf{\Delta }{\rho }^{(GPS,COSO)}$ and $\mathsf{\Delta }{\rho }^{(GPS,COSJ)}$ in time to find estimates corresponding to the start and end dates of the Envisat pairs. The GPS range values at individual points in time are then converted to pairs by first taking the difference between range estimates from COSO and COSJ to remove the far-field deformation and then by time corresponding to the time intervals of the Envisat pairs. The uncertainty for the range change estimated from GPS is derived from measurement uncertainty at each station. Differenced range changes are then converted to range change rates after dividing by the time interval for each pair. We find an estimated range change rate of $3.8\pm 3.3\mathrm{mm}/\mathrm{year}$ at COSO with respect to COSJ. We then difference the InSAR range change measurements for pixels corresponding to COSO with measurements for pixels corresponding to COSJ from the Envisat pairs. We use the standard error of the mean as an estimate of uncertainty for the range changes observed using InSAR. The mean range change rate at COSO with respect to COSJ as measured by Envisat is $5.8\pm 10.7\mathrm{mm}/\mathrm{year}$. We compare the difference between the GPS and InSAR data sets for unwrapped range change rates estimated at COSO with far-field deformation effects removed using a Student’s one-sample t-test (e.g., [37]). We test the null hypothesis that the means of the two data sets are equal. We find $T_{calc}=-1.3$, which is less in absolute value than the critical value $T_{\alpha /2}=2.0$ with 80 degrees of freedom. Thus, we conclude (with 95% confidence) that there is no significant difference between mean range change rate estimated from GPS and measured by Envisat.

We perform the same procedure for our Sentinel-1A (S1A) pairs. In this case, we use the Sentinel-1A unit pointing vector ${\widehat{\mathbf{s}}}^{\left(S1A\right)}=[-0.63;-0.11;0.77]$ to derive the corresponding range change rates from GPS. The estimated range change rate at COSO with respect to COSJ from GPS over the 2014–2016 time frame is $4.6\pm 2.3\mathrm{mm}/\mathrm{year}$. The mean range change rate at COSO with respect to COSJ as measured by Sentinel-1A is $3.3\pm 2.2\mathrm{mm}/\mathrm{year}$. Repeating the same test at 95% confidence, we find $T_{calc}=2.1$ to be less than the critical value $T_{\alpha /2}=2.3$ with nine degrees of freedom. Thus, we conclude that there is no significant difference between range change rate estimated from GPS and measured by Sentinel-1A.

Figure 7 and Figure 8 show histograms of the difference in range change rates between GPS and individual InSAR pairs for both the Envisat and Sentinel-1A data sets.

3.2. Estimating Volume Change of The Reservoir

A previous modeling study by Ali et al. [27] uses a poroelastic, homogeneous model to describe the roughly circular pattern of surface deformation at Coso. Their best-fitting results are shown in

Table 1. Finite element modeling results from Ali et al. [27].

Parameter Name	Best-Fitting Estimate	Uncertainty
Easting in m	428,651.3	1000
Northing in m	3,987,538.4	1000
Reservoir depth in m	2366.6	250
Reservoir half thickness in m	2053.1	250
Reservoir half length in m	3027.2	500
Young’s Modulus in GPa	25	0
Poisson’s Ratio	0.25	0
Pressure Change in MPa	$-0.3$	0
Biot’s Coefficient	1	0

.

To arrive at refined estimates of dimensional properties and volume change of the reservoir in this study, we model the deformation in each of the 91 interferometric pairs using a “cuboid” parameterization to include a single sink at a given depth in an elastic half space with uniform material properties [38]. In this model, the cuboid represents a volume element with sides of width W, length L, and height H, and an initial volume of $V_0=LWH$. The cuboid is sliced into eight equal-sized octants by three rectangular patches that are mutually orthogonal. Each rectangular patch i is a dislocation with a negative value of tensile opening and has a corresponding slip vector of $\mathbf{U}=[0,0,-U_{i,3}]$. The slip $-U_{i,3}$ on a singular patch is proportional to the ratio of its area to the total volume change:

$$\mathsf{\Delta }V=-1\times (U_{W,3}LH+U_{L,3}WH+U_{H,3}LW).$$

(3)

We favor this parameterization over other analytical deformation models which may be used to describe circular surface deformation patterns due to its flexibility in terms of model geometry and its direct interpretability in terms of the geothermal reservoir. The model parameters, including the source location, volume change, cuboid dimensions, and InSAR-related nuisance parameters (e.g., [39,40]), along with their uncertainties, are estimated using simulated annealing methods employed by the general inversion of phase technique (GIPhT) [39,41].

The temporal distribution of our InSAR data set separates our deformation modeling into two distinct time intervals in accordance with the realization of the Hay Ranch project: 2004 to mid-2010 (corresponding to Envisat pairs and the time up to the start of the project), and 2014 to 2016 (corresponding to Sentinel-1A pairs and the time after the project was initiated). Working with the same stack of Envisat pairs as Ali et al. [27], we estimate all the model parameters starting with initial estimates based on the best-fitting results from Ali et al. [27] (Table 1). We assume the Poisson’s ratio to be $\nu =0.21$, consistent with the average values of seismic velocities used to locate the seismicity. We then take the best-fitting estimates for the sink location and cuboid dimensions found from the stack inversion to be fixed and estimate the volume change for each of the 81 individual Envisat interferometric pairs. We repeat the same procedure for the Sentinel-1A pairs using their corresponding stack.

3.3. Time Series

We are interested in the temporal trend of deformation at Coso in terms of volume change of the geothermal reservoir. To analyze this trend, we perform time-series analysis using temporal adjustment on the set of volume changes $\mathsf{\Delta }{V}_{i,j}$ estimated from the individual InSAR pairs spanning the time intervals from $t_i$ to $t_j$. This procedure converts the pair-wise volume changes $\mathsf{\Delta }{V}_{i,j}$ for individual interferometric pairs into a cumulative value of volume change at each point $t_i$ and $t_j$ in time [26]. As described previously, this approach implicitly assumes that the temporal dependence and spatial dependence of the deformation field are separable functions [26,39]. Accordingly, we write the vector displacement field $\mathit{u}$ as

$$u(t,x)=f\left(t\right)G\left(x\right),$$

(4)

where $f\left(t\right)$ is a function of time t only and $G\left(x\right)$ is a function of spatial position coordinate x only. In our case we consider $G\left(x\right)$ to be the model of a cuboidal sink contracting in a half space with uniform elastic properties. For the temporal function, we consider several parameterizations.

The first assumes a constant rate during the time interval between 2004 and 2010:

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

(5)

where the initial time $t_0=2003.87$ is the start date of our data set in decimal years, and a (different) constant rate $a_2$ during the time interval between 2010 and 2016:

$$f_2\left(t_i\right)=a_2(t_i-t_1),$$

(6)

where the time $t_1=2010$ refers to the date when pumping operations were altered at Coso (e.g., [6,7,42]).

To model the full data set, we use a piecewise-linear parameterization with m breaks at times $t_m$ that form $(m-1)$ segments:

$$f_3\left(t_i\right)={\mathsf{\Sigma }}_{j=1}^{j=m}{a}_{j}D\left(t_i\right),$$

(7)

where

$$D\left(t_i\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{t}_i<t_j\hfill \\ (t_j-t_i)\hfill & \mathrm{if}{t}_j\le t_i<t_{j+1}\hfill \\ (t_{j+1}-t_j)\hfill & \mathrm{if}{t}_i\ge t_{j+1}\hfill \end{array}\right..$$

(8)

We use this function with $m=2$, $m=4$, and $m=12$.

We also consider an exponentially decaying rate parameterization

$$f_4\left(t_i\right)=a_e\left(1-exp\left(-\frac{t_i-t_0}{\tau }\right)\right),$$

(9)

where $\tau$ is a characteristic time scale found through nonlinear optimization. This trend would be consistent with drawdown as has been observed at Brady (e.g., [43,44]).

To increase the temporal sampling, we also analyze the volume change in the temporal dimension using the daily displacements derived from GPS data. We difference each daily record of displacement in our data set from station COSO with the previous day’s record, resulting in 3650 measurements of differential position, i.e., displacement. We then convert these displacements to range change using Equation (2) and divide by the time interval $\mathsf{\Delta }t=t_j-t_i$ of the paired values to arrive at the range change rate $\dot{\rho }$. We estimate the volume change and reservoir depth in a similar manner to the InSAR data set by taking the best-fitting estimates for the source location and cuboid dimensions found from the stack inversion to be fixed. We then perform a similar time-series analysis on the resulting volume change rates.

4. Results

4.1. Deformation Modeling

The results from applying the “cuboid” parameterization to mean range change rates from both the Envisat and Sentinel-1A stacks are shown in Figure 9 and Figure 10 as well as

Table 2. Results from nonlinear inversion with a “cuboid” parameterization on the stack of MST Envisat pairs spanning 2004 to 2011.

Parameter Name	Best-Fitting Estimate	Uncertainty
Easting gradient $\delta \rho /\delta x$	6.78 $\times {10}^{-8}$	2.19 $\times {10}^{-7}$
Northing gradient $\delta \rho /\delta y$	−9.93 $\times {10}^{-9}$	1.86 $\times {10}^{-7}$
Upwards gradient $\delta \rho /\delta z$	−9.97 $\times {10}^{-7}$	9.67 $\times {10}^{-7}$
Offset in cycles	−0.5	0.8
Poisson ratio $\nu$	0.21	0
Volume change rate $\dot{V}{[\mathrm{m}}^3/\mathrm{year}]$	−2.02 $\times {10}^5$	2.81 $\times {10}^4$
Centroid Easting E in m	428,580.4	531.3
Centroid Northing N in m	3,987,509.0	250.0
Centroid Depth in m	2418.2	531.3
Cuboid L in m	3.07 $\times {10}^3$	3.23 $\times {10}^2$
Cuboid W in m	2.28 $\times {10}^3$	2.25 $\times {10}^2$
Cuboid H in m	1.95 $\times {10}^3$	2.42 $\times {10}^2$

and

Table 3. Results from nonlinear inversion with a “cuboid” parameterization on the stack of Sentinel-1A pairs spanning 2014 to 2016.

Parameter Name	Best-Fitting Estimate	Uncertainty
Easting gradient $\delta \rho /\delta x$	1.64 $\times {10}^{-9}$	6.25 $\times {10}^{-8}$
Northing gradient $\delta \rho /\delta y$	−1.39 $\times {10}^{-7}$	9.38 $\times {10}^{-8}$
Upwards gradient $\delta \rho /\delta z$	6.73 $\times {10}^{-7}$	2.19 $\times {10}^{-6}$
Offset in cycles	−0.5	0.6
Poisson ratio $\nu$	0.21	0
Volume change rate $\dot{V}{[\mathrm{m}}^3/\mathrm{year}]$	−1.07 $\times {10}^6$	2.25 $\times {10}^5$
Centroid Easting E in m	428,870.4	48.4
Centroid Northing N in m	3,986,482.0	46.9
Centroid Depth in m	3098.7	195.3
Cuboid L in m	3.12 $\times {10}^3$	4.84 $\times {10}^1$
Cuboid W in m	2.34 $\times {10}^3$	4.84 $\times {10}^1$
Cuboid H in m	1.93 $\times {10}^3$	4.84 $\times {10}^1$

, respectively. We find an estimated depth of reservoir to be $2.4\pm 0.5\mathrm{km}$ for the time period covering 2004 to 2011, consistent with the results from Ali et al. [27], while the estimated depth of the reservoir from 2014 to 2016 is $3.1\pm 0.2\mathrm{km}$. We calculate the dimensionless misfit $\chi =1.7$ of the model to the Envisat data as the square root of the reduced ${\chi }^2$-test statistic ([45], p. 334). Similarly, we find the misfit of the model to the Sentinel-1A stack to be $\chi =1.2$.

We similarly estimate the reservoir depth before and after 2010 using the GPS data and source location estimates from the InSAR stacks. We find that the best-fitting estimated reservoir depth for observations before 2010 is $2.6\pm 0.5\mathrm{km}$, whereas the best-fitting estimated reservoir depth for observations after 2010 is $3.1\pm 0.6\mathrm{km}$. For these solutions, the misfit $\chi$ is 1.4 and 1.6, respectively.

4.2. Time-Series Analysis

We start by parameterizing the volume change estimates before and after 2010 by separate, single-rate temporal functions corresponding to when well operations were varied at the site (e.g., [6,7,42]). Due to a lack of data between 2011 and 2014, we also include a break at the end of the Envisat data set (3 September 2010). Interferometric pairs included in the analysis are shown in Figure 2. For the InSAR data, we parameterize all the Envisat estimates between 2004 and 2010 using a single rate, and then similarly parameterize all the Sentinel-1A estimates between 2014 and 2016 using a single rate function. We cannot estimate any volume change for the time between 2011 and 2014, when we have no InSAR coverage. The results for this parameterization (Equation (7), $m=4$) are shown in Figure 11. We calculate the dimensionless misfit to be $\chi =0.5$. We perform a two-tailed Student’s t-test to test the null hypothesis that these two estimated rates are equal (e.g., [37]). We find a significant difference between the estimated rate before 2010 and the estimated rate after 2014 with 95% confidence (

Table 4. Results from two-tailed Student’s t-test with degrees of freedom $dof$ at the 95% confidence level (e.g., [37], p. 524) for volume change rates estimated from interferometric synthetic aperture radar (InSAR) data. $H_0$: equal rates; $H_1$: rates differ.

Start of Break 1	Start of Break 2	Rate 1 $\times {10}^6$ [m3/year]	Rate 2 $\times {10}^6$ [m3/year]	$\mathit{dof}$	Test Value	Critical Value	Result
20031114	20100101	$(-1.2\pm 0.04)$	$(-1.3\pm 0.29)$	78	3.95	1.99	reject $H_0$
20031114	20141102	$(-1.2\pm 0.04)$	$(-0.80\pm 0.05)$	76	−25.16	1.99	reject $H_0$
20100101	20141102	$(-1.3\pm 0.29)$	$(-0.80\pm 0.05)$	16	−4.67	2.12	reject $H_0$

).

We also use a piecewise-linear temporal function with breaks on 1 June and 1 December of each year from 2005 to 2010 (Equation (7), $m=12$) to further examine if any seasonal trends we see in the pumping records (Figure 6) are reflected in the InSAR data before changes in well operations. We perform an F-test for model complexity to decide whether the increased complexity of the piecewise-linear model is justified (e.g., [37], pp. 536 and 627). We test the null hypothesis that the single rate parameterization and the piecewise-linear parameterization fit the data equally well at 95% confidence. We find $F_{calc}=0.94$, which is less than the critical value $F_{\alpha =0.05}=1.43$ with degrees of freedom $do{f}_1=90$ and $do{f}_2=79$. We conclude that the complexity of the piecewise-linear model is not justified with 95% confidence and continue with the single rate parameterization.

We use the same single-rate temporal functions for analyzing volume change estimated from GPS, treating the data before 2010 and after 2010 separately, as shown in Figure 12 (Equations (5) and (6)). We calculate the dimensionless misfit for the pre-2010 interval to be $\chi =2.3$. The dimensionless misfit for the post-2010 interval is $\chi =1.5$. We perform a two-tailed Student’s t-test to test the null hypothesis that the rates in the two intervals are equal (e.g., [37], p. 524). The null hypothesis is rejected with 95% confidence (

Table 5. Results from two-tailed Student’s t-test with degrees of freedom $dof$ at the 95% confidence level (e.g., [37], p. 524) for volume change rates estimated from GPS data pre-2010 and post-2010. $H_0$: no significant difference in rates; $H_1$: signficant difference in rates.

Pre-2010 Rate $\times {10}^6$ m3/year	Post-2010 Rate $\times {10}^5$ m3/year	$\mathit{dof}$	Test Value	Critical Value	Result
$(-1.2\pm 0.02)$	$(-8.10\pm 0.09)$	3638	−843.22	1.96	reject $H_0$

) We find that the rate for the post-2010 interval is significantly smaller (slower) than the pre-2010 rate.

Table 6. Estimates of subsurface volume change from multiple data sets.

Data Source	$\dot{\mathit{V}}$ [L/s]
Gross injection [5]	$427\pm 105$
Gross production [5]	$-924\pm 202$
Net production [5]	$-497\pm 228$
Cuboid sink estimated from InSAR data from 2005 to 2010	$-38\pm 1$
Cuboid sink from InSAR data from 2014 to 2016	$-25\pm 2$
Cuboid sink estimated from GPS data from 2005 to 2010	$-38\pm 1$
Cuboid sink from GPS data from 2014 to 2016	$-26\pm 3$
Average recharge [47,48]	130

compares several estimates of the subsurface volume change rate $\dot{V}$. Rates are derived from pumping records using a range for the produced brine density of ${\rho }_{prod}\in [990,1250]{\mathrm{kg}/\mathrm{m}}^3$. The injected brine density is calculated using the thermal expansion coefficient for water, a production temperature of $350{}^{\circ }$C, and an injection temperature of $160{}^{\circ }$C [46]. Results shown in Table 6 assume a density of $1000{\mathrm{kg}/\mathrm{m}}^3$ for the produced brine and $1043{\mathrm{kg}/\mathrm{m}}^3$ for the injected brine. Net production is calculated as the difference between the rates derived from gross injection and gross production.

4.3. Correlation Test

We quantify the relationship between pumping and deformation by performing correlation tests. Using the modeled values of volume change derived from temporal adjustment, we estimate the volume change at the start of each month and compare it to the monthly gross production rate. We normalize the values using a statistical Z-transform, defined as

$$\widehat{\mathbf{x}}=\frac{\mathbf{x}-\overline{\mathbf{x}}}{{\sigma }_x}$$

(10)

(e.g., [37], p. 181). The results are summarized in Figure 13. We first consider the monotonic relationship between volume change and gross production rate using Spearman’s test (e.g., [37], p. 786). We find a Spearman’s correlation coefficient $r_s=0.75$ with a corresponding p-value of $p=2.14\times {10}^{-16}$, indicating the probability of rejecting the null hypothesis of no correlation when it is actually true. This result suggests a strong monotonic relationship between volume change and gross production. We further explore the linearity of this relationship using Pearson’s test. We find a high correlation between the volume change and gross production rate, with a (Pearson’s) correlation coefficient of $R=0.75$ (e.g., [37], p. 599). We test the null hypothesis of no correlation against the alternative of non-zero correlation using a Student’s t-test (e.g., [37], p. 599). We find a p-value of $p=2.2\times {10}^{-16}$ Thus, we conclude that the positive correlation between the volume change each month and gross production rate is significant. Consequently, this result supports our hypothesis that deformation is causally associated with pumping at Coso. Performing the same tests with the rates of gross injection and net production yields p-values of $5.3\times {10}^{-7}$ and $8.0\times {10}^{-7}$, respectively.

4.4. Identifying a Driving Mechanism for the Observed Subsidence

We convert our estimates of volume change from modeling the interferometric pairs with a “cuboid” parameterization to volumetric strain rates by dividing by the initial volume $V_0$ of the single cuboid and the time span $\mathsf{\Delta }t$. The modeling results indicate that the cuboidal sink is shrinking with an average volumetric strain rate of $\overline{\dot{ϵ}}=\overline{\mathsf{\Delta }V}/\left(V_0\mathsf{\Delta }t\right)$ of the order of 90 microstrain/year or 3 picostrain/second in absolute value.

We consider two possible mechanisms: (1) Decreasing pore-fluid pressure and (2) thermal contraction of the rock matrix. Accordingly, we interpret our estimates of volumetric strain derived from our estimates of volume change as a result of temperature change and/or pressure change. Considering a decrease in pore fluid pressure as the cause of volume change, the observed volume change rate $\dot{V}$ in the cuboidal reservoir is related to a pressure change rate $\dot{P}[\mathrm{Pa}/\mathrm{year}]$ by:

$$\begin{array}{cc}\hfill {\dot{V}}^{\left(P\right)}& =\left(\frac{1}{H}\dot{P}\right)V_0\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & ={\dot{ϵ}}^{\left(P\right)}{V}_0,\hfill \end{array}$$

(12)

where $1/H{[\mathrm{Pa}}^{-1}]$ is the poroelastic expansion coefficient [49] and ${\dot{ϵ}}^{\left(P\right)}{[\mathrm{year}}^{-1}]$ is a poroelastic volumetric strain rate.

Alternatively, if we consider the volume change to be caused by thermal contraction, then the estimated rate of volume change rate in the cuboid is related to the rate of temperature change $\dot{T}[\mathrm{K}/\mathrm{year}]$ by

$$\begin{array}{cc}\hfill {\dot{V}}^{\left(T\right)}& =\left({\alpha }_T\dot{T}\right)V_0\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & ={\dot{ϵ}}^{\left(T\right)}{V}_0,\hfill \end{array}$$

(14)

where ${\alpha }_T{[\mathrm{K}}^{-1}]$ is the thermal expansion coefficient and ${\dot{ϵ}}^{\left(T\right)}{[\mathrm{year}}^{-1}]$ is a thermal volumetric strain rate.

We can also consider a combination of these two processes through a linear combination of Equations (12) and (14).

Following the procedure in Reinisch et al. [40], we define a priori confidence intervals for reasonable values of strain rate under each interpretation. Using results from well logs [50], we define the mean values for the rates of change in pressure $\dot{P}$ and temperature $\dot{T}$ with confidence intervals defined such that a rate of zero is two standard deviations away from the mean value.

$$\begin{array}{cc}\hfill \dot{P}& \in (-9.4\pm 4.7)\times {10}^5\mathrm{Pa}/\mathrm{year},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \dot{T}& \in (-2.6\pm 1.3)\mathrm{K}/\mathrm{year}.\hfill \end{array}$$

(16)

We similarly define confidence intervals for ${\alpha }_T\in (1.0\pm 0.5)\times {10}^{-5}1/\mathrm{K}$ and $1/H\in (5.0\pm 2.5)\times {10}^{-10}1/\mathrm{Pa}$ (e.g., [8,27,49,51,52,53,54]). This leads to confidence intervals for volumetric strain rate interpreted in terms of decreasing pore-fluid pressure alone, thermal contraction alone, or a combination of thermal contraction and decreasing pore-fluid pressure:

$$\begin{array}{cc}\hfill {\dot{ϵ}}^{\left(P\right)}& \in (-5.6\pm 4.2)\times {10}^{-5}{\mathrm{year}}^{-1},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\dot{ϵ}}^{\left(T\right)}& \in (-2.6\pm 2.0)\times {10}^{-5}{\mathrm{year}}^{-1},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\dot{ϵ}}^{(P+T)}& \in (-8.2\pm 4.6)\times {10}^{-5}{\mathrm{year}}^{-1}.\hfill \end{array}$$

(19)

We then compare the realized a posteriori 68% confidence interval for our estimates to the a priori 68% confidence intervals defined in Equation (17) to determine if any of the interpretations are unreasonable. Figure 14 shows the resulting comparisons for Envisat and Sentinel-1A pairs. In each case, we see no overlap between the realized 68% confidence interval for the volumetric strain rates and the 68% confidence interval for reasonable values of ${\dot{ϵ}}^{\left(T\right)}$ defined a priori. We infer that the estimated values of volumetric strain rates are too high to be attributed to thermal contraction of the rock matrix alone. In contrast, we see clear overlap between the realized 68% confidence interval for volumetric strain rates and the 68% confidence interval for reasonable values of ${\dot{ϵ}}^{\left(P\right)}$ and ${\dot{ϵ}}^{(P+T)}$ defined a priori, with the most apparent overlap with the ${\dot{ϵ}}^{\left(P\right)}$ interpretation. We infer that decreasing pore-fluid pressure is the dominant driving mechanism causing the subsidence observed at Coso.

5. Discussion

We find that the volumetric strain is most reasonably attributed to decreasing pore-fluid pressure during both time intervals. We also find a significant difference between the volume change rates estimated before 2010, between 2010 and 2011, and after 2014 from our InSAR data set.

The significant increase in volume change rate during the first six months following completion of the Hay Ranch project could be explained by the increased pumping activity at the site, as suggested by Eneva et al. ([42], their Figure 4) and apparent in Figure 6. The strong correlation between the amount of volume change and pumping rates (Figure 13) corroborates such an explanation.

For observations before 2010, we find the best-fitting cuboidal model for the reservoir to have a centroid depth of $2.4\pm 0.5$ km, consistent with previous studies (e.g., [27,42]). The results using the Sentinel-1A stack for the time interval between 2014 and 2016, however, show that the best-fitting depth for the reservoir between 2014 and 2016 is $3.1\pm 0.2$ km. The reservoir depths estimated from GPS data during both the pre-2010 and post-2010 intervals are similar to the results estimated from InSAR data. To test whether this change in depth is significant, we use a two-tailed Student’s t-test at 95% confidence. We test the null hypothesis that there is no significant difference between the centroid depth estimated for pairs before 2011 and the centroid depth estimated for pairs after 2014. We find a test statistic of $|T_{calc}=4.367|$ which is greater than the critical value $T_{\alpha /2}=1.99$ with 89 degrees of freedom. Thus, we infer a significant difference in centroid depths before 2011 and after 2014 with 95% confidence.

The results from this geodetic modeling suggest that the depth to the center of the reservoir increased after 2010. A change in reservoir depth could be explained by reservoir depletion after 2010, when the pumping rate was increased (e.g., [42]). If the fluid level in the reservoir drops, then faults that were previously saturated by fluids within the reservoir would no longer be filled with such pressurized fluids. This Coulomb effect provides one explanation for the correlation between the monthly seismicity rates and pumping rates at Coso before 2010 and the lack of correlation thereafter as observed by Reinisch [55]. An increase in reservoir depth after 2010 could also explain the decrease in maximum deformation rates observed by InSAR after 2010 (e.g., Sentinel-1A pairs), as shown in Figure 3) and found in previous studies (e.g., [42]). The magnitude of displacement at the surface decreases as the depth of the sink increases (e.g., [38]). This reasoning explains the decrease in estimated volume change rate between pairs before 2011 and pairs after 2014.

Comparing the volume change rates in Table 6, we find that the estimates of volume change rate of the reservoir from deformation modeling are an order of magnitude less than those predicted through standard density calculations using the pumping records. Under the simple assumptions of a nearly incompressible fluid in a poroelastic half space, Segall [56] interprets Skempton’s coefficient B as the “ratio of solid volume change to change in pore fluid volume.... Thus, for example, if water is uniformly withdrawn from a rock with $B=0.8$ the volumetric contraction of the rock is 80% of the volume of the extracted water.” Deng et al. [57] give the range of Skempton’s coefficient for crustal rock to be between $0.5$ and $0.9$. Given the metamorphic setting of Coso with basalt and rhyolites present (e.g., [2,7]), values for B could fall below this range as well (e.g., [58]). Considering the recharge to the system, a value of B around $0.65$ would explain the discrepancy between volume change rate estimates in Table 6.

Eneva et al. [42] have also suggested poroelastic effects to explain the subsidence at Coso. The model presented by Segall [56] for depleting a reservoir at constant rate and source location has been used to describe changes in subsidence at other geothermal fields (e.g., [59]). According to this model, if the rate of production (or net extraction) is constant, then the response in terms of subsidence would vary as a smooth function of time. To test this possibility, we perform temporal adjustment using a temporal function with exponentially decaying rate (Equation (9)) with a characteristic time scale of $\tau =18\mathrm{year}$ found through nonlinear optimization. Our cuboidal models estimated from the InSAR and GPS data sets show different reservoir depths before and after 2010, which violates the model assumptions from Segall [56]. We use the differential (day-to-day) measurements of vertical displacement corresponding to the pairs in the GPS data set. We find a misfit of $\chi =2.0$. For comparison, we perform temporal adjustment on the pair-wise vertical displacements using a piece-wise linear parameterization with a break on January 1, 2010 (Equation (7), $m=2$), after the well operations were altered. We find a misfit of $\chi =1.8$, indicating a better fit than the smooth model corresponding to a reservoir depleting at constant rate. Using an F-test, we reject the null hypothesis that the two models provide equally good fit with 95% confidence. We infer that the changes in subsidence rate, reservoir contraction, and sink depth are due to the change in injection protocol in late 2009. Apparently, the geothermal system at Coso is too complex to be explained by a model of a simple sink with constant rate of depletion and constant location.

6. Conclusions

Using methods outlined previously [40], we have advanced the characterization of the subsidence observed by InSAR at Coso geothermal field by estimating the volume change rate for each interferometric pair using a “cuboid” parameterization. Using temporal adjustment, we find a significant difference between volume change rates estimated before and after 2010 with 95% confidence. We also identify decreasing pore-fluid pressure as the dominant mechanism driving the observed deformation. We confirm a significant positive correlation between deformation and production rate.