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

: We measure transient deformation at Coso geothermal ﬁeld 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-ﬁtting 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-ﬂuid pressure as the dominant mechanism driving the subsidence observed at Coso geothermal ﬁeld. We also ﬁnd a signiﬁcant positive correlation between deformation and production rate.


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 • 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 ± 0.27) × 10 8 kg/month (924 ± 10 kg/s) and (11.72 ± 0.02) × 10 8 kg/month (446 ± 7 kg/s), respectively, which, neglecting natural recharge, corresponds to a net extraction rate of ∼12.58 × 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-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 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].  [28]. Faults from Jennings [30] and Jennings et al. [31] are shown as thin black lines. The interferogram is overlaid on a 2017 Landsat/Copernicus image from Google Earth.   Figure 3. Example interferograms from Envisat and Sentinel-1A satellites spanning the production interval. For each pair (row), the left panel shows wrapped phase in cycles, the middle panel shows range change rate in mm/year, and the right panel shows values of range change rate in mm/year taken along the profile (thick black line in the middle panel). The incidence angle (denoted as the small arrow in the northeast of the phase plot) is 21.1 • and the azimuth (from North, denoted as the large arrow in the northeast corner of the phase plot) is −11.0 • for Envisat pairs. The incidence angle is 33.9 • and the azimuth is −13.1 • for Sentinel-1A pairs. One cycle of color in wrapped phase corresponds to 28 mm of range change. Plotting conventions as in Figure 1.

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.

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 m horizontally and 600 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 ν = 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 g/cm 3 [35]. This is consistent with other density estimates at Coso (e.g., [2,36]).

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].

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 u relative to stable North America from both stations to range change ∆ρ by taking the negative scalar product with the unit vectorŝ (ENV I) = [−0.41; −0.08; 0.91] pointing from the pixel on the ground to the radar sensor aboard the satellite.
We then linearly interpolate and extrapolate ∆ρ (GPS,COSO) and ∆ρ (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 ± 3.3 mm/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 ± 10.7 mm/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 α/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ŝ (S1A) = [−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 ± 2.3 mm/year. The mean range change rate at COSO with respect to COSJ as measured by Sentinel-1A is 3.3 ± 2.2 mm/year. Repeating the same test at 95% confidence, we find T calc = 2.1 to be less than the critical value T α/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. Figures 7 and 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. Figure 7. Histogram of differences between range change rate at COSO with respect to COSJ estimated from GPS using linear interpolation and extrapolation and range change rates at COSO with respect to COSJ observed by Envisat pairs.

Figure 8.
Histogram of differences between range change rate at COSO with respect to COSJ estimated from GPS using linear interpolation and extrapolation and range change rates at COSO with respect to COSJ observed by Sentinel-1A pairs.

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. 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 = LW H. 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 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: 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 ν = 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.

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 ∆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 ∆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 u as where f (t) is a function of time t only and G(x) is a function of spatial position coordinate x only. In our case we consider G(x) 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: 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: 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: where We use this function with m = 2, m = 4, and m = 12.
We also consider an exponentially decaying rate parameterization where τ 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 ∆t = t j − t i of the paired values to arrive at the range change rateρ. 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.

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 Figures 9 and 10 as well as Tables 2 and 3, respectively. We find an estimated depth of reservoir to be 2.4 ± 0.5 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 ± 0.2 km. We calculate the dimensionless misfit χ = 1.7 of the model to the Envisat data as the square root of the reduced χ 2 -test statistic ( [45], p. 334). Similarly, we find the misfit of the model to the Sentinel-1A stack to be χ = 1.2.   [39]. The estimate of the parameter vector was based on finite element models by Ali et al. [27]. Inversion was performed using unwrapped range change rates. Results are shown in terms of unwrapped range change rate: (a) Observed range change rate, (b) modeled range change rate, (c) residual between observed and modeled, and (d) absolute value of residuals. We use the Universe Transverse Mercator (UTM) coordinate system as in Figure 1.  [39]. The estimate of the parameter vector was based on finite element models by Ali et al. [27]. Inversion was performed using unwrapped range change rates. Results are shown in terms of unwrapped range change rate: (a) Observed range change rate, (b) modeled range change rate, (c) residual between observed and modeled, (d) and absolute value of residuals. We use the UTM coordinate system as in Figure 1. 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 ± 0.5 km, whereas the best-fitting estimated reservoir depth for observations after 2010 is 3.1 ± 0.6 km. For these solutions, the misfit χ is 1.4 and 1.6, respectively.

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 χ = 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).
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 α=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.  Table 4. Results from two-tailed Student's t-test with degrees of freedom do f 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. 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 χ = 2.3. The dimensionless misfit for the post-2010 interval is χ = 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) We find that the rate for the post-2010 interval is significantly smaller (slower) than the pre-2010 rate. Table 5. Results from two-tailed Student's t-test with degrees of freedom do f 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.  Table 6 compares several estimates of the subsurface volume change rateV. Rates are derived from pumping records using a range for the produced brine density of ρ prod ∈ [990, 1250] kg/m 3 . The injected brine density is calculated using the thermal expansion coefficient for water, a production temperature of 350 • C, and an injection temperature of 160 • C [46]. Results shown in Table 6 assume a density of 1000 kg/m 3 for the produced brine and 1043 kg/m 3 for the injected brine. Net production is calculated as the difference between the rates derived from gross injection and gross production. Table 6. Estimates of subsurface volume change from multiple data sets.

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 aŝ (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 × 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 × 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 × 10 −7 and 8.0 × 10 −7 , respectively.

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 ∆t. The modeling results indicate that the cuboidal sink is shrinking with an average volumetric strain rate of¯ = ∆V/(V 0 ∆t) 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 rateV in the cuboidal reservoir is related to a pressure change rateṖ [Pa/year] by: where 1/H [Pa −1 ] is the poroelastic expansion coefficient [49] and˙ (P) [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 changė where α T [K −1 ] is the thermal expansion coefficient and˙ (T) [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Ṗ and temperatureṪ with confidence intervals defined such that a rate of zero is two standard deviations away from the mean value. P ∈ (−9.4 ± 4.7) × 10 5 Pa/year, We similarly define confidence intervals for α T ∈ (1.0 ± 0.5) × 10 −5 1/K and 1/H ∈ (5.0 ± 2.5) × 10 −10 1/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: 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˙ (T) 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˙ (P) and˙ (P+T) defined a priori, with the most apparent overlap with the˙ (P) interpretation. We infer that decreasing pore-fluid pressure is the dominant driving mechanism causing the subsidence observed at Coso. Overlain are 68% confidence intervals for reasonable values defined a priori (red) and realized values from deformation modeling (green). Also shown is the best fitting estimate of mean strain rate (black) found from temporal adjustment with a single rate temporal function.

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 ± 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 ± 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 α/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 τ = 18 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 χ = 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 χ = 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.

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.