Relatively stable pressure effects and time-increasing thermal contraction control Heber geothermal field deformation

Domain: earth_science
Published: 2024-06-17
Editorial outcome: minor_revision

Abstract
Due to geological complexities and observational gaps, it is challenging to identify the governing physical processes of geothermal field deformation including ground subsidence and earthquakes. In the west and east regions of the Heber Geothermal Field (HGF), decade-long subsidence was occurring despite injection of heat-depleted brines, along with transient reversals between uplift and subsidence. These observed phenomena contradict current knowledge that injection leads to surface uplift. Here we show that high-yield production wells at the HGF center siphon fluid from surrounding regions, which can cause subsidence at low-rate injection locations. Moreover, the thermal contraction effect by cooling increases with time and eventually overwhelms the pressure effects of pressure fluctuation and poroelastic responses, which keep relatively stable during geothermal operations. The observed subsidence anomalies result from the siphoning effect and thermal contraction. We further demonstrate that thermal contraction dominates long-term trends of surface displacement and seismicity growth, while pressure effects drive near-instantaneous changes.

Introduction
Technologies to utilize natural geothermal systems as a source of renewable, carbon-free energy have existed for many decades. However, it is well-documented that geothermal energy production can cause significant geohazards like earthquakes and ground subsidence, which threaten the sustainability of long-term operations and expose nearby populations to potentially damaging impacts1–4. Various physical mechanisms have been suggested for the geohazards, including reservoir pressure depletion associated with fluid loss5,6, poroelastic contraction7,8, thermoelastic contraction9–11, and aseismic fault slip12,13. Such mechanisms can be active concurrently, but it is challenging to quantify their relative strengths and spatiotemporal evolution features, especially for large-scale fields with complex operation history.

The Heber Geothermal Field (HGF) is a representative field located in the Imperial Valley, southern California, USA, where geothermal resources are abundant (Fig. 1a). Many large energy-producing fields, like the Salton Sea Geothermal Field, have been setup since the early 1980s. All these fields operate in regimes with high crustal strain rates and heat flow owing to the tectonic transition from strike-slip to ridge-transform style tectonics8,12. The geohazards associated with geothermal operations in the valley have always been a major concern in both scientific and industrial communities14,15.

Hydrothermal flow at the HGF is controlled by a matrix permeability structure consisting of a shallow sedimentary reservoir and three faults16 (Fig. 1b). The upper surface of the reservoir is bound by clay caprock; and two northwest-trending dextral strike-slip faults (plate boundary and subsidiary faults) control its lateral extent. The boundary fault, also known as the Dixieland fault, is part of an active, blind extension of the Superstition Hills fault17–20. In addition, there is a northeast-trending normal-slip fault, with a dip angle of ~77° below 1.68 km, which is a primary pathway for deep, high-temperature brine upwelling into the shallow reservoir. We thus refer to it as the feeder fault. With the effect of regional groundwater flow from natural recharge, the overall shape of the isothermal contour of 182 °C within the reservoir resembles a lopsided mushroom above the depth of ~2 km16,21.

The HGF resource was initially estimated to support a capacity of up to 500 MW of electrical power generation. One double flash power plant and one binary power plant were completed in 1985; however, due to low production rates, the binary plant designed with a supercritical cycle and a single turbine-generator was discontinued in 1987 and subsequently decommissioned22. In 1992, the binary plant was transferred into a subcritical cycle incorporating 6 modules with 33 MWe net output. At present, the field has one double flash and three binary plants, with a capacity of 89 MWe23. A total of 69 vertical and deviated wells were drilled during the development history (Fig. 2a), mainly in the years of 1985, 1987, 1993, and 2006–2009 (Supplementary Table 1). Of these, there are 28 injection and 31 extraction wells as well as 10 wells that were first operated for ~2 years of heat extraction, and then used for reinjecting heat-depleted brines since 1993.

Along with the complex operation history, injection, and production fluid temperatures vary with location and time (Fig. 2b). In the west and east regions, where only fluid injection occurs, the average temperatures of fluid injected into their respective wells increased by ~5 °C after 2005 and decreased by ~27 °C after mid-2006, respectively. The temperature history in the central region is more complicated because all forms of operations are active: extraction, injection, and extraction later converted to injection. The temperatures of fluid extracted from 11 deviated extraction wells near the feeder fault were ~188 °C in 1985, declined sharply by ~55 °C in 1986, and then increased to ~150 °C in 1988; from 1988 to 2017, the average fluid temperature decreased by ~21 °C. In 1993, 20 vertical extraction wells began large-scale heat production with fluid temperatures of ~174 °C, which decreased to ~161 °C by 2017. At the deviated wells near the feeder fault, temperature decreased at approximately 0.72 °C/yr, which is ~33% higher than at the vertical wells (0.54 °C/yr). Such irregular temperature variations indicate that thermal effects within the reservoir are inhomogeneous in space and time.

To resolve the physical mechanisms of HGF deformation and quantify the spatiotemporal evolution of their strengths, we compile and analyze more than two decades of geodetic and seismic observations. We construct a realistic 3D thermo-hydro-mechanical model of the reservoir and bounding faults to simulate the full operation history of the HGF over 32 years, with results constrained by the surface displacement observations and reported geothermal brine temperatures (see Data Availability). Here we show that the pressure effects of pressure fluctuation and poroelastic responses keep relatively stable during the operation history, but the strength of thermal contraction associated with the reinjection of cool edgeothermal fluids increases gradually with time and eventually overwhelms the pressure effects. Long-term trends of surface displacement and seismicity growth at the HGF are dominated by thermal contraction, while the pressure effects drive their short-term changes.

Results
Decade-long surface displacement anomalies
Vertical and horizontal ground displacements at the HGF have been measured by leveling and interferometric synthetic aperture radar (InSAR) techniques with relatively high temporal and spatial resolution. Leveling surveys with a dense network of benchmarks have been conducted approximately annually since 199423 (Supplementary Fig. 1). The displacement time series of each benchmark are used to characterize the vertical displacement trends of the HGF and also three subregions: west, central and east. Also, we use the time series to generate vertical displacement rate maps of four time periods: January 1994 to December 2004, February 2006 to November 2010, November 2010 to December 2014, and December 2014 to January 2018 (Supplementary Fig. 2), with the Kriging method24. InSAR monitoring is implemented with the Envisat satellite images (May 2006 to September 2010) and the Sentinel-1 satellite images (April 2015 to December 2017, June 2018 to August 2019) based on the SqueeSAR approach25. More details on geodetic data processing are provided in the Methods section.

The leveling displacement maps and time series show that the central region of the HGF experienced long-term subsidence with the maximum rates increasing from ~1.5 cm/yr to ~4 cm/yr in 2005 (Fig. 3). In contrast, both west and east regions are characterized by reversals of vertical displacements around 2005. The west region was subsiding at rates up to 1.5 cm/yr from 1994 to 2005 and then transitioned to uplift with rates up to 2 cm/yr. The east region shows a reversal in displacement trends: initially uplifting at rates up to 0.5 cm/yr and then subsiding with rates up to 2 cm/yr. The displacement rate and trend changes that initiated around 2005 and lasted for more than 10 years are referred to as decade-long transients. The leveling-observed transients are confirmed by InSAR measurements, but the leveling-measured subsidence at the HGF center is significantly larger than the InSAR result (Fig. 3b), which is most likely due to differences in the resolution and error sources between the two independent methods. Moreover, the InSAR-derived horizontal displacement maps of 2006-2010 and 2018-2019 reveal the central region was moving eastward with rates up to 2 cm/yr (Supplementary Fig. 3).

To investigate potential relationships between surface displacement and geothermal operations, we compare vertical displacement trends of the west, central, and east regions of the HGF with the respective injection and extraction volume changes of geothermal fluid (Fig. 3a). In the west region, with only fluid injection (Fig. 2a and Supplementary Fig. 4), the vertical displacement is characterized by subsidence before 2005 and uplift after 2005 with temporal variations which are similar to the changes of injection rates. It is common that fluid injection causes surface uplift due to pore pressure increases within the target reservoir26,27. However, the ground of the west region was subsiding from 1994 to 2005 even with a net increase of injected fluid, which is an unexpected phenomenon. The same displacement anomaly of ground subsiding with a net increase of injected fluid is also seen in the east region after 2005. In the central region with net fluid extraction since early 1985, subsidence was relatively constant until the post-2005 transient of acceleration, which corresponds to the prominent increase of production rates after 2005. The comparisons above indicate that reservoir pressure depletion associated with fluid loss is not the sole driver of ground deformation.

Seismicity rate growth
Geothermal field deformation refers to not only surface displacement but also induced earthquakes. Low levels of seismicity at the HGF were first detected by Southern California Seismic Network (SCSN) in 1988, three years after the first power plant came online. Since then, there have been multiple episodes of transient and rapid increases in seismicity rates superimposed on a steady background rate (Fig. 4). The most notable one started in 2005 and lasted for ~2 years. To verify the temporal growth of seismicity from the SCSN catalog, we further use the method of template matching to improve the temporal resolution of the regional earthquake catalog and account for variations in detection levels due to network changes (see Methods for more details). From 2001 to 2018, 487 earthquakes are detected within the HGF, including twelve 2 < M < 3 and one M > 3 events, about eight-fold larger than the count of the SCSN catalog. The template matched catalog confirms the increasing trend of seismicity with time. In addition, we find two interesting phenomena: (1) rapid increases of earthquake within the template-matched catalog occurred at the time with prominent changes of both injection and production rates; (2) the temporal growth of seismicity rates is similar to the vertical displacement trend of the HGF shown by the normalized leveling time-series, suggesting that the physical mechanisms driving surface displacement also control the occurrence of earthquakes.

Fault slip modeling
The HGF is located within an active plate boundary zone, which implies a possibility of ground displacement associated with fault slip. More specifically, the strike-slip faults, which bound the HGF laterally could represent active, blind extensions of the Superstition Hills fault to the north (Fig. 1). As the surface displacement at the HGF consists of two parts—long-term displacement initiated before 2005 and transient displacement after 2005 (Fig. 3)—we thus conduct three scenarios of fault slip modeling with available displacement observations. First, the long-term vertical displacement measured by leveling from 1994 to 2004 is inverted with slip on the normal-slip feeder fault. Second, the vertical transient displacement from 2006 to 2010, obtained through removing the long-term displacement trend, is also inverted with the feeder fault. Third, the horizontal displacement measured by InSAR from 2006 to 2010 is modelled with the strike-slip plate boundary (Dixieland) fault. More details on inversions for fault slip refer to Methods.

Our inversions show that the magnitudes of normal slip on the feeder fault are required to be 0.4 to 1 m above the depth of 12 km to fit the long-term vertical displacement (Supplementary Fig. 6), which indicates slip rates up to 4 to 10 cm/yr. In contrast, the slip rates of the northwest extension of the nearby plate boundary fault are less than 10 mm/yr17–19, much smaller than the inversion result. On the other hand, such large cumulative slip would cause catastrophic, irreversible casing deformation of the wells intersecting with the feeder fault; but there were no related reports on such damage, and production rates have remained relatively constant. Therefore, it is implausible that the feeder fault was slipping at the prominent rate. For the transient vertical displacement characterized by not only subsidence but also uplift, major displacement features cannot be fit by normal slip on the feeder fault (Supplementary Fig. 7). To match the InSAR-measured horizontal displacement from 2006 to 2010, the magnitudes of dextral strike-slip on the plate boundary fault need to be up to 30 to 82 cm below the depth of 8 km (Supplementary Fig. 8), which suggests 6 to 16 cm/yr of slip rates, much higher than the reported slip rates and the relative plate motion rate (only ~4 cm/yr)28,29. Consequently, the dominant reason for the vertical and horizontal displacement is unlikely related to fault slip.

Thermo-hydro-mechanical (THM) modeling
To resolve the deformation mechanisms, we further investigate the effects of geothermal operations through simulating the injection and extraction history of the HGF with a 3D coupled thermo-hydro-mechanical (THM) model (Supplementary Fig. 9). Geological structures and geothermal gradients of the geothermal system are determined according to previous studies16,21,30–32 (Fig. 1b) and seismic velocity profiles of the study region from the Southern California Earthquake Center Community Velocity Model-Harvard (CVM-H)33 (Supplementary Fig. 10). Variations in rock properties with depth are relatively well resolved in the CVM-H but there is little-to-no geophysical evidence of lateral variations on scales smaller than the HGF. Thus, to honor the natural complexity of the system while keeping the problem numerically tractable, we simplify it with five layers: caprock at the surface, two reservoir layers, a basal layer, and basement rock. Additionally, we embed the two NW-trending strike-slip faults and the NE-trending feeder fault in our model as individual formations with a unique set of thermo-hydro-mechanical properties.

The model layers and faults are assumed to behave as linear thermo-poro-elastic media. The density, Young’s modulus, and Poisson’s ratio of each layer are calculated based on the local seismic velocity structure. The corresponding mechanical properties of the three faults are set with typical values34–36: 2500 kg/m3, 10 GPa, and 0.3. For the thermophysical properties, we set the specific heat and thermal expansion coefficients to be 1000 J/(kg·°C) and 10−5 °C−1, respectively, which are typical values for rocks37. The thermal conductivities of the two reservoir layers, basal formation, and feeder fault are set to be two-fold larger than the values of the other formations (i.e., caprock, basement, and the two strike-slip faults). Another critical aspect of a THM reservoir model is the permeability structure. According to published permeability models of the HGF16,30, the two reservoir layers are characterized with relatively high permeability, and the feeder fault serves as the dominant upward channel of deep hot geothermal brine. Moreover, given that the strike-slip faults are part of an active plate boundary system, it is plausible that they have damage zones with enhanced permeability38. We thus prescribe finite widths and relatively high permeabilities to the faults.

The mechanical properties of the five layers are well constrained by the seismic velocity profiles, but all other properties are assigned based on qualitative models or prior data from general earth materials, and thus have some uncertainties. Due to the intercoupling of dynamic and nonlinear physical processes within the geothermal system, inverting for model parameters is computationally infeasible. Therefore, we calibrate the uncertain parameters through an iterative procedure with a subset of surface displacement observations and fluid temperature records. As hydrogeological properties of the reservoir have a major influence on surface displacement39,40, we first calibrate its porosity, permeability, Biot coefficient, and three thermophysical parameters (specific heat, thermal conductivity and expansion coefficient) successively. Second, we calibrate the same parameters of the three faults as above, including their thicknesses. Third, the heat flow rate within the feeder fault was estimated mainly based on the initial temperature distribution prior to production30–32. In our starting model, only the feeder fault is assumed to be the upwelling conduit of deep heat flux. To identify the replenishment mechanism of heat flow at the HGF, we conduct additional simulations under two scenarios: (I) one upwelling channel of the feeder fault with different rates of heat flux; (II) two channels of the feeder and subsidiary strike-slip faults with different flow rates. Lastly, we test the influences of different fluid properties, including density, viscosity, compressibility, thermal conductivity, and specific heat, which were initially assigned to be 1000 kg/m3, 2.5 × 10−4 Pa·s, 4.2 × 10−10 Pa−1, 0.68 W/(m·°C), and 4200 J/(kg·°C), respectively.

For each parameter combination, we simulate the operation of the HGF from Apr 1985 through Dec 2010; the simulation results are compared with the vertical displacement time series (Jan 1994-Nov 2010) at 64 leveling benchmarks, InSAR displacement maps in both vertical and horizontal directions (May 2006 to Sep 2010), and temperature records of 16 injection and 16 extraction wells. The benchmarks and wells are selected to give approximately uniform coverage within the HGF (Supplementary Fig. 1 and 4). The optimal parameter values are identified based on the misfits of simulation results to displacement and temperature observations. As the observation period of the leveling data is four-fold longer than the InSAR observation period, we set their weight ratio to be 4:1 to calculate the misfits of surface displacement. To minimize potential uncertainties of the optimal solution associated with the choices of test parameter values, we conduct two rounds of model calibration (more details refer to Methods).

After calibration, only the values of three reservoir parameters are changed relative to initial settings (Supplementary Tables 2-14). First, the permeability of the lower reservoir is decreased from 20 mD (2 × 10−14 m2) to 10 mD. Although the upper and lower reservoirs are characterized with the same porosity of 18%, the upper reservoir (20 mD) is more permeable than the lower reservoir. The porosity and permeability of the feeder and subsidiary faults are equal to the values of the upper reservoir. We test whether the reservoir layers and the three faults have different permeabilities in the horizontal and vertical directions, but simulation results reveal that they deform like isotropic porous media. Second, the thermal expansion coefficients and conductivities of the two reservoir layers are increased from 1.0 × 10−5 °C−1 to 1.3 × 10−5 °C−1 and from 2 W/(m·°C) to 3 W/(m·°C), respectively, higher than the corresponding parameter values of the other formations. The reservoir is thus characterized with faster heat conduction and a more prominent deformation response to temperature changes. On the other hand, our simulations show that the fits to displacement and temperature observations reduce with increasing rates of heat flow upwelling through the feeder fault, and that the root-mean-squared (RMS) misfits of displacement become larger when the subsidiary strike-slip fault is included as another upward channel of heat flux (Supplementary Table 15). Therefore, the feeder fault is the major heat flow channel of the HGF, with replenishment rates unlikely exceeding ~20 kg/s. Moreover, the simulations with different fluid properties show that the model misfits are lowest with the initial parameters settings (Supplementary Table 16).

The model calibration also provides insights into the sensitivities of simulation results to changes in the rock and fluid properties. We quantify the sensitivities with the variation of RMS misfits to displacement and temperature observations (Fig. 5). Overall, the simulated surface displacements are more sensitive to the parameter changes than the temperature. Among all test parameters, the influences of hydraulic properties and fault thickness are largest, indicating that surface displacement and heat transfer of the HGF are mainly associated with fluid flow. Beyond those, the variation of RMS misfits to displacement observations are up to ~16-17% with test thermal expansion coefficients and specific heat of the reservoir, which is an indicator of a non-negligible thermoelastic effect. The temperature is most sensitive to the permeability changes of the feeder fault that serves as the replenishment channel of heat flux. For the fluid properties, the density, viscosity, compressibility, and specific heat all have a prominent influence on the simulation results for surface displacement.

Our parameter calibration improves the model fit to the displacement and temperature observations. The weighted RMS misfits of displacement are reduced from 2.89 cm to 2.45 cm; the RMS misfits of temperature are decreased from 21.55 °C to 20.3 °C (Supplementary Tables 4-16), respectively, with major improvement from the permeability optimization of the lower reservoir and the optimization of the thermal expansion coefficients and conductivities of the two reservoir layers. We further use the calibrated model to simulate the operation of the HGF from Apr 1985 to May 2017, ~6.5 years longer than the time span of model calibration, with available production data. Major features of observed surface displacement in both vertical and horizontal directions can be replicated with our model (Fig. 6 and Supplementary Fig. 13). After 32 years of operation, the largest surface subsidence at the HGF center is up to 50 cm. In addition, there is also significant horizontal displacement towards the center with the maximum magnitude up to ~30 cm. However, the prominent leveling subsidence at the HGF center after 2010 and the observations in the peripheral region of the HGF are not reproduced well by our model (Fig. 3b). One possible reason for the difference between the simulated and observed displacements is that the model parameters calibrated with available observations represent the rock properties in average and cannot reflect laterally inhomogeneous hydromechanical and thermophysical properties and their variations in the depth to basement rock of the study region41. The effect of inhomogeneous rock properties increases with the duration of geothermal simulation, as shown by the decreasing goodness-of-fit between the model predictions and the leveling time series at the HGF center (Supplementary Fig. 14). Another possible reason is that the displacement observations in the peripheral region of the HGF, where no geothermal operations occur, may be related to other physical processes like natural recharge of regional groundwater. In addition, it is worth noting that the leveling-observed subsidence at the HGF center is larger than the InSAR monitoring results (Supplementary Figs. 2 and 3).

Our calibrated model can also predict the reported temperature variations for most wells (Supplementary Fig. 16). For the six deviated extraction wells intersecting the feeder fault (Fig. 2a), there are systematic biases between the model predictions and observations, which could indicate a mismatch between actual temperature distribution in 1985 and the initial temperature conditions used in the model (Supplementary Fig. 11). In addition, our model cannot reproduce ~30 °C temperature decreases of 3 selected vertical extraction wells (identifiers 2591216, 2591217, and 2591221) after 2000, which are likely due to inhomogeneous hydraulic properties of the reservoir. Strong heat convection associated with fluid flow may occur in some localized regions due to higher porosity and permeability than the calibrated values.

Reservoir pressure and temperature evolution
Our simulations show that the movement of fluid within the reservoir is controlled to a first order by flow from injection wells towards extraction wells (Fig. 7a). Pore pressure variations are widespread within the reservoir, but temperature decreases are only seen in localized regions around injection wells. To inspect the spatiotemporal changes of reservoir pressure and temperature, we inspect the model results at eight probe points across the model domain located at the bottom of the upper reservoir (1.5 km depth). Pressure time series at these points show that dramatic fluctuations mainly occurred within four periods: 1985-1988, 1992-1994, 2003-2006, and 2012-2015 (Fig. 7c). The maximum pressure fluctuation – pressure decrease – occurred near the feeder and subsidiary faults from 1985 to 1988, dropping approximately 3 MPa from the initial pressure of 14.7 MPa to ~11.7 MPa (points P5 and P6). In contrast, pressure within the plate boundary fault (point P3) was decreasing slowly before 2006 by ~1.5 MPa relative to its initial value, and then started increasing but at a lower rate.

In the western region of the HGF, although the net injection rates are around 0.6 × 109 kg per month on average from 1994 to 2005, pressure decreased by 0.5-1.0 MPa (points P1-P3). From the Darcy velocity field (Fig. 7a), we can deduce that the pressure decrease is attributable to the loss of fluid that was flowing toward the center. High production rates around 3.5 × 109 kg/mo at the HGF center created a low-pressure zone siphoning fluid from surrounding regions. After 2005 the average injection rates increased up to ~1.6 × 109 kg/mo, and the pressure rose gradually, indicating that the injection volume became larger than the loss. For the central region, probe points P4-P6 show different trends: pressures declined continuously by ~2 MPa among the vertical extraction wells (point P4) but declined abruptly near the feeder and subsidiary faults (points P5-P6) during the first three years of operation. The rapid decline is likely due to exhaustive exploitation, which is probably the major reason why the operator decommissioned the binary plant in 1987. To boost power generation capacity, 9 production wells were closed in 1988 or 1989 and then converted to injection around 1993. In the east region with average injection rates of ~2.1 × 109 kg/mo, the pressure at locations within and ~1.5 km far from the well cluster (points P7 and P8) fluctuates around average values of ~16 MPa and ~14.9 MPa, respectively, across the whole time period.

In contrast to the sharp changes of pore pressure, fluid temperature varies smoothly over time (Fig. 7c), indicating insensitivity to abrupt changes in production rates and slower heat transfer than pressure diffusion. In addition, the sizes of localized temperature anomalies scale with injection rates. For instance, around probe points P5 and P7, the areas with temperature decrease over 20 °C eventually merged due to dense injection wells and large injection volumes (Supplementary Fig. 17). These two findings reveal that heat transfer within the reservoir is mainly through the convection of fluid flow driven by geothermal operations, which is supported by model sensitivity tests showing that well temperatures are more sensitive to changes in hydraulic properties rather than the thermophysical properties.

Quantification of inhomogeneous thermo-poro-elastic effects
As the calibrated model can simulate major features of observed surface displacements and fluid temperatures, we further use it to quantify the relative influence and spatiotemporal evolution of pore pressure changes, poro- and thermo-elastic effects. The geothermal operation history of the HGF is simulated under three end-member scenarios: (I) the Biot coefficients of the two reservoirs and the three faults are set to be zero; (II) the thermal expansion coefficients of the two reservoir layers and the three faults are assigned to be zero; and (III) both kinds of parameters are set to be zero. The effect of pressure changes can be nullified through comparing the simulate results of model I with those of the calibrated (reference) model. The thermoelastic effect can be isolated through comparison of the reference model with model II. The poroelastic effect of contraction and expansion can be obtained directly from model III.

The pressure changes within the HGF reservoir cause both subsidence and uplift, mainly in the central and northeast regions (Fig. 8a), respectively, with magnitudes varying slightly during the operation history of the HGF (Supplementary Fig. 21). By May 2017, the maximum subsidence and uplift are ~8 cm and ~4 cm, respectively. The poroelastic effects of contraction and expansion produce similar spatiotemporal features of vertical displacements as the pressure changes, but with magnitudes larger by ~20%. In contrast, the thermoelastic effect only causes ground subsidence; we thus call it thermal contraction. The subsidence is initially localized in the northeast region but gradually extends to the central and west regions (Supplementary Fig. 23). Around 2000, the subsidence in the central region became greater than the other regions. The effect of thermal contraction has been increasing steadily with increasing steadily with time since the operation began. After only ~2 years of operation, vertical displacement at the northeast corner caused by thermal contraction exceeds the displacement associated with the joint effect of pressure fluctuation and poroelastic responses (Fig. 8b). By 2017, the thermoelastic subsidence was more than twice as large as the displacements caused by the two pressure effects.

To further compare the differences of pressure effects and thermal contraction, we inspect their temporal variations at four specific locations (Fig. 8c, 8d). In contrast to the instantaneous response and proportional relation of poroelastic responses to pressure fluctuation, thermoelastic effect varies nonlinearly and asynchronously with temperature. For example, at point P1 the ground experienced thermoelastic subsidence with temperature increase during some time periods; at point P5, the thermal contraction effect increased ~13 years earlier than temperature decrease. These phenomena also demonstrate that the thermal effect can disturb broad regions beyond localized areas around injection wells with temperature decrease42,43. The strength of thermal contraction at a specific location not only depends on the amplitude of temperature decrease and its distance off the temperature decrease front but also the injection volume of heat-depleted brines. We can see that thermoelastic subsidence at point P6 with ~3 °C of temperature decrease is comparable to the value at point P7 with 53.5 °C of temperature decrease. Moreover, these probe points reveal that abrupt changes of vertical displacements are associated with reservoir pressure changes, which respond quickly to geothermal operations.

Physical mechanisms of HGF deformation
Although the three mechanisms of pressure fluctuation, poroelastic responses and thermal contraction all contribute to the deformation of the HGF, their effects vary significantly in space and time. In the west region, the ground subsidence anomaly from 1995 to 2005 is jointly controlled by the three mechanisms with comparable strengths (point P1 in Fig. 8c, d). The pressure decrease was associated with the siphoning effect of the HGF center with a pressure drop of up to 3 MPa relative to the initial pressure (Fig. 8a). After 2005 with the injection volume exceeding the loss, the reservoir pressure begun to increase, which is the reason for the decade-long displacement transient of turning to surface uplift. For the central region, the transient acceleration in subsidence after 2005 is mainly caused by thermoelastic contraction, which is also responsible for the displacement anomaly of turning to subsidence around 2005 in the east region (Fig. 8c, d).

In addition, the simulated displacement time series at 8 probe points show that the short-term changes and long-term trend of vertical surface displacements are associated with the two pressure effects and thermal contraction (Fig. 8 and Supplementary Fig. 24), respectively. To identify the physical mechanisms controlling deformation of the whole study region, we further normalize the simulated vertical displacements at 8 probe points caused by the pressure and thermal effects. Leveling-observed surface displacement trend and temporal growth of seismicity within the HGF is highly correlated with the normalized vertical displacements linked to thermal contraction (Fig. 4c), with the correlation coefficients up to 0.85 and 0.97, respectively. In contrast, the coefficients of the leveling displacements and seismicity growth relative to the normalized displacements associated with the pressure effects are only 0.09 and 0.23, respectively. We thus conclude that thermal contraction dominates the long-term trend of HGF deformation including surface displacement and seismicity, while the pressure effects drive short-term changes.

Methods
Leveling data analysis
There are 158 leveling stations including a reference benchmark (A-33) at the HGF (Supplementary Fig. 1). By the end of 2018 these benchmarks were surveyed 25 times: January 1994, March 1995, January 1996, December 1996, December 1998, March 2000, January 2001, December 2001, December 2002, December 2003, December 2004, February 2006, January 2007, January 2008, December 2008, October 2009, November 2010, December 2011, November 2012, December 2013, December 2014, December 2015, December 2016, January 2018, and December 2018. However, only 103 stations were included during every survey. The observation errors of leveling range from 0.2 cm to 0.5 cm. To characterize vertical displacement trends of the HGF, we divide the study area into three regions (west, central, and east) to analyze the leveling observations. The displacement time series of the 103 stations are first normalized by their root-mean-square values with the following expression,

n(i) = u(i) / sqrt( (1/N) * sum_{i=1}^{N} u(i)^2 ) (1)

where u(i) and n(i) are the measured and normalized displacement, respectively, at each time node; and N is the total number of surveys. We then average the normalized time series of the stations within each sub-region to represent its displacement trend. For the whole study area with both surface uplift and subsidence, absolute values of all normalized time series are averaged to represent the overall trend of surface displacement.

The time series of leveling displacements are also used to generate vertical displacement rate maps of four time periods: 1994/01-2004/12, 2006/02-2010/11, 2010/11-2014/12, 2014/12-2018/01, which are divided based on the observed displacement rate changes (Supplementary Fig. 1) and time periods of InSAR data (Supplementary Fig. 3). We first calculate the displacement rates of the stations with at least four surveys during each time period (left column in Supplementary Fig. 2), and then derive the maps (right column in Supplementary Fig. 2) through interpolating the discrete points with the Kriging method24.

InSAR data processing
Surface displacement of the study area is also measured with SAR images of the Envisat and Sentinel-1 satellites. All the data are processed with the SqueeSAR approach25, which utilizes permanent and distributed scatterers to measure the surface displacement in the line-of-sight (LOS) direction to the satellites. The LOS displacement rate maps of ascending and descending tracks of each time period are further decomposed into vertical and east horizontal components using the approach proposed by Eneva et al.23 (see Supplementary Fig. 3). Although most scatterers are aligned along roads and canals, their spatial coverage is dense enough to detect displacement features of the HGF. We also adopt the Kriging method24 to generate complete displacement images through interpolating the scatterers for comparison with the leveling rate maps.

Earthquake template matching
Template matching is an important observational tool for characterizing repetitive seismic sources46. Multi-station template matching often involves cross-correlating the waveforms that were recorded from a cataloged earthquake across a network against years of continuously recorded data to detect similar, smaller magnitude earthquakes. Here, we use all earthquakes of the SCSN catalog occurred in the Heber area from 2001 to 2018 as templates. Parameters used in template creation and event detection are consistent with those in our previous work47. Templates are 20 second in length and begin 5 second prior to P-phase arrivals on vertical channels and 5 second prior to S-phase arrivals on horizontal channels. All seismic waveforms are filtered between 5-15 Hz, and these templates are then cross-correlated against a five-station network of seismometers, namely MONP from network AZ, and DRE, EMS, ERP, and WES from network CI. Newly detected events have to exceed a threshold of 15 times the hourly median absolute deviation of the network normalized cross-correlation coefficients (NNCCC) in addition to the NNCCC exceeding a coefficient of 0.2.

Fault slip inversion
First, the normal-slip feeder fault and strike-slip plate boundary fault are set to be 20 km wide and discretized with 0.46 km (long) × 2.0 km (wide) and 1.46 km (long) × 2.0 km (wide) rectangular dislocation patches, respectively. Second, we conduct linear slip inversions with the bounded variable least squares algorithm48. The Green’s functions relating observations to slip on the patches are calculated with the elastic half-space dislocation model49. To add virtual observations and avoid singularities in the slip distribution, we apply the modified Laplacian operator50, which can smooth slip on boundary patches of a finite fault reasonably, during each inversion with the smoothing factor set to be 0.1.

Reservoir model setup—3D geometric model and tetrahedral meshing
We construct a 3D model to represent the HGF. The spatial scale of the model domain is 16 km long (E-W), 14 km wide (N-S), and 8 km deep, which is two-fold larger than the scope of the HGF reservoir (Supplementary Fig. 9). The geographical coordinates of the model center are coincident with well 2591224 (−115.535°E, 32.712°N). In the depth direction, based on the geological model (Fig. 1b) and seismic velocity model of the study area (Supplementary Fig. 10), we divide the model into five layers: caprock (0-0.55 km deep), upper reservoir (0.55-1.65 km deep), lower reservoir (1.65-3.15 km deep), basal layer (3.15-5.3 km deep), and basement (5.3-8 km deep). In addition, all operation wells are added into the model based on their trajectories of open-hole sections, which are mainly located within the upper reservoir. Lastly, the two NW-trending strike-slip and one NE-trending normal fault are also included in our model with the thicknesses initially set to be 100 m and 300 m, respectively. For the strike-slip faults, due to unknown surface traces, we simplify them as two lines based on the southeast extension of the Superstition Hills fault and their locations presented by James et al.16. We test the influence of different model sizes on simulation results and find that larger model sizes have no obvious improvement in fitting the leveling, InSAR and temperature observations (Supplementary Table 3).

The 3D model is discretized with tetrahedral elements for finite element analysis. According to the sizes and roles (significance) of wells, faults and layers, we classify them into four categories: (1) wells, (2) two strike-slip faults, (3) caprock, feeder fault and reservoirs, (4) basal layer and basement. The minimum and maximum sizes of tetrahedral elements used to discretize wells are 3 m and 50 m, respectively. For the two strike-slip faults, the minimum and maximum sizes are 5 m and 130 m, respectively. For the third category of model domain, the minimum and maximum sizes are 7 m and 270 m, respectively. For the fourth category of model domain, the minimum and maximum sizes are 9 m and 410 m, respectively. After discretization, there are 1166303 domain elements, 77696 boundary elements, and 5012 edge elements in total. We also test other mesh schemes and find that model simulation with the preferred mesh scheme can obtain an acceptable balance between resolving large stress gradients and computational efficiency.

Reservoir model setup—Parameters setting for mechanical, hydraulic and thermophysical properties
The mechanical parameters of Young’ modulus E and Poisson’s ratio v of the five layers are calculated based on the velocity and density profiles of the study area (Supplementary Fig. 10). For the faults in the model, their values of density, Young’s modulus, and Poisson’s ratio are initially assumed to be 2500 kg/m3, 10 GPa, and 0.3, respectively, which are widely-used values34–36. The Biot coefficient α, corresponding to the effective stress coefficient for bulk deformation, is an important parameter characterizing the influence of pore pressure changes on the deformation of rock skeleton. We calculate the parameter values with two equations: α = (1 + v + 2φ(1-2v)) / (3(1-v))51; α = 1 - K/Ks52. φ is the porosity of rock. K is the drained bulk modulus equal to E / [3(1-2v)]. Ks is the unjacketed bulk modulus, often called the solid-grain modulus, and here is assumed to be 100 GPa. Extreme values of α approaching 1 and 0 corresponds to unconsolidated material (or highly damaged rock) and a porosity-free aggregate53, respectively. There is some difference between the calculation results with the two equations (Supplementary Table 2). We thus test the influence of different Biot coefficients of the reservoirs and faults on the simulation results during subsequent model calibration.

For the hydraulic parameters of porosity and permeability, we first set their starting values and then conduct model calibration with a subset of available surface displacement observations and fluid temperature records. Guided by previously published 2D permeability models16,30, the two reservoir layers, the feeder and subsidiary faults are set with relatively high values in contrast to the other formations (Supplementary Table 2). The porosity and permeability of caprock (predominantly shales) are set to be 3% and 0.28 mD54. For the two reservoirs, we initially set the porosity and permeability to be 18% and 20 mD, respectively, based on two published studies21,30. The porosity and permeability of the subsidiary and feeder faults are initially assigned to be 20% and 20 mD, respectively, close to the values of the reservoirs. The hydraulic properties of the other formations are assigned based on the corresponding lithology. For the thermophysical properties of rocks within different formations, we set their values according to Robertson (1988)37.

Reservoir model setup—Initial and boundary conditions
We set two kinds of initial conditions for the model: (1) temperature distribution; (2) stress equilibrium between gravity and hydrostatic pressure. For the initial temperature, we refer to the published temperature models21,30–32 and assume its distribution to be axially symmetric (Supplementary Fig. 11). In the depth direction, the temperatures at the surface and the depth of 10 km below the model center are set to be 22 °C and 370 °C, respectively. The gradients of temperature increase (ΔTd) in the depth direction below the center are 55.2 °C/km and 11 °C/km, respectively, for the upper two layers. For the other three layers at the center and the model margin (10 km away from the center), the ΔTd values are assigned to be 25 °C/km. In the horizontal direction, the temperature is assumed to decrease linearly with the radius far from the center at different depths. For the initial stress equilibrium, we exert both gravity and hydrostatic pressure on the model. However, there are minor residual stresses that cannot be balanced with this method; we thus use the simulated stress at t=0 to balance the residual stress.

The boundary conditions consist of mechanical, hydraulic and thermal conditions. The mechanical boundary conditions include: (1) free upper surface, (2) roller boundaries for the bottom and side surfaces, and (3) fixed open-hole sections for all wells. The hydraulic boundary conditions comprise: (1) zero pressure at the upper surface, (2) an upwelling heat fluid source on the bottom of the feeder fault with the mass flux initially assigned to 10 kg/s30–32, (3) hydrostatic equilibrium at the bottom and side surfaces, and (4) flux discontinuity for injection and production wells. The thermal boundary conditions include: (1) thermal insulation on the top surface, (2) a heat source on the bottom of the feeder fault with the temperature of heat flux assigned to be 320 °C at the depth of 8 km, (3) heat exchange on the bottom and side surface (open boundary), and (4) heat flux boundaries for the open-hole sections of injection wells. It is worth noting that there is no conductive heat flux across the extraction wells.

Thermo-poro-elasticity formulation
The theory of thermo-poro-elasticity is derived from momentum, mass and energy conservation26,52,55. The quasi-static equilibrium of the solid matrix of porous medium (with shear modulus G, Poisson’s ratio v, permeability k and porosity φ) can be described with the Navier-Stokes equation,

G * nabla^2 * u_i + (G / (1-2v)) * (partial^2 u_k / (partial x_i partial x_k)) = alpha * (partial p / partial x_i) + beta * K * (partial T / partial x_i) - F_i (2)

where u_i is displacement component in direction x_i of a Cartesian coordinate system. The quantities p and F_i are pore pressure and body force, respectively. The Biot-Willis coefficient α corresponds to the effective stress coefficient for bulk deformation, controlling the magnitude of poroelastic strain induced by fluid injection. β is the volumetric expansion coefficient of the rock. K is the drained bulk modulus of porous media.

The mass conservation equation describing laminar Darcy-type flow (with density ρ_f, viscosity η, and bulk modulus K_f) within porous rock is represented as

S_epsilon * (partial p / partial t) - (k / eta) * nabla^2 * p = Q_m / rho_f (3)

where Q_m is the mass source/sink term. The expression of the constrained specific storage coefficient S_epsilon can be formulated as26,

S_epsilon = alpha * ( (1/K) - (1/Ks) ) + phi * ( (1/Kf) - (1/Ks) ) (4)

where 1/Ks and 1/Kf are unjacketed bulk and pore compressibility, respectively. This formulation shows that the storage coefficient depends on not only the elastic moduli of the porous medium but also the fluid properties. In this study the porous media is assumed to be isotropic and microscopically homogeneous, which indicates 1/Ks equal to 1/Ks/phi.

Heat transfer in porous media through conduction of fluid-solid mixture and convection of fluid flow obeys the law of energy conservation. Assuming local thermal equilibrium between the fluid (with specific heat C_f and thermal conductivity κ_f) and rock (with density ρ_s, specific heat C_s and thermal conductivity κ_s), and neglecting the dissipation of mechanical energy due to deformation of the solid, the energy conservation can be written as

(phi * rho_f * C_f + (1-phi) * rho_s * C_s) * (partial T / partial t) + rho_f * C_f * nabla * (q * T) - (phi * kappa_f + (1-phi) * kappa_s) * nabla^2 * T = Q_e (5)

where T is the temperature of the fluid-solid mixture, and Q_e is an external energy source. q is the Darcy flux, which can be expressed as (rho_f * k / eta) * nabla * p without consideration of the gravitational potential and has been involved in the mass conservation equation.

There are three kinds of two-way coupling within the theory of thermo-poro-elasticity (Supplementary Fig. 12). However, the effects of convective heat transfer, poro- and thermo-elastic responses of skeletal rock deformation associated with pressure and temperature changes are expected to dominate over other effects like the influence of rock deformation on the temperature and pressure fields56–58. For computational efficiency, only the three major effects are included in our modeling. We employ COMSOL Multiphysics® (version 6.0, https://www.comsol.com/) to solve the governing equations. The average time of simulating 9400 days (from 1985/4/1 to 2010/12/30) of fluid injection and extraction is about 30 hours with our workstation (CPU: Intel Xeon Gold 6240C 2.60 GHz; RAM: 384 GB DDR4 ECC REG 2933 MHz).

Model parameter calibration
An iterative calibration strategy involving dozens of steps is adopted to optimize the model parameters assigned based on qualitative models or prior data from general earth materials (Supplementary Table 2). In each step only one kind of parameter is calibrated, and the optimal parameter derived from preceding step is inherited in next step. First, the two hydraulic parameters (porosity and permeability), Biot coefficient, and three thermophysical parameters (specific heat, thermal conductivity and expansion coefficient) of the two reservoir layers are calibrated successively. Second, we test different values for the thickness, elastic, hydraulic, and thermophysical parameters of the three faults. Third, we test whether both of the feeder fault and subsidiary fault serve as the upwelling channel of heat flux. Fourth, we test the influence of different properties of water on simulation results. After above simulation tests, we conduct a second round of model calibration to further optimize the parameters with the most prominent influence on simulated displacement and temperature. As part of the calibration process, we also perform sensitivity tests of surface displacement and well temperature observations on different model parameters. More details are listed below.

1st round calibration—Step I: Parameter calibration for the reservoir layers
First, the porosities of the two reservoir layers are calibrated simultaneously. We adopt the grid searching method to pinpoint the optimum values based on the fit of simulation results to leveling, InSAR and temperature observations. The search step is set to be 2%. When the porosities of the two reservoirs are both equal to 18%, the weighted RMS misfit of simulated surface displacement to the leveling and InSAR observations is least (Supplementary Table 4). In contrast, the RMS misfits to well temperature observations are least with the porosities of the upper and lower reservoirs equal to 20% and 18%, respectively. However, the weighted RMS misfit of displacement is increased by 11% (3.25 cm to 2.89 cm) with in the two groups of porosity, three-fold larger than the improvement of well temperature (21.55 °C to 20.71 °C). We thus select the optimal porosities of the two reservoirs to be 18%, which is the same as the initial model setting (Supplementary Table 2). Second, we calibrate the horizontal and vertical permeabilities of the upper reservoir also with the grid searching method. The search step is set to be 10 mD. The RMS misfits of displacement are least with the two permeabilities equal to 20 mD (Supplementary Table 5). In contrast, the RMS misfits of temperature are least when the horizontal and vertical permeabilities are equal to 20 mD and 30 mD, respectively. However, the RMS misfit of displacement is increased by 7% with in the two groups of porosity, larger than the improvement (5%) of well temperature. We thus select the optimal permeability of the upper reservoir to be 20 mD. Third, the permeabilities of the lower reservoir are optimized with the same method as above. The misfits to displacement and temperature observations are smallest when the horizontal and vertical permeabilities both equal to 10 mD (Supplementary Table 6). Fourth, we test the influence of different Biot coefficients of the two reservoirs on the simulation results. The RMS misfits of displacement vary little with the Biot values varying from 0.52 to 1 (Supplementary Table 7). The RMS misfits of temperature keep constant with different Biot values. We thus retain the initial setting of the reservoir Biot coefficients. Lastly, we calibrate three thermophysical parameters: expansion coefficient, conductivity, and specific heat. The misfits of displacement are least with the three parameters equal to 1.3 × 10−5 °C−1, 3 W/(m·°C) and 1000 J/(kg·°C) (Supplementary Table 8), respectively. Simulated temperature shows no or little variations with different values of the thermophysical parameters.

Above simulations show that the displacement and temperature vary little with different Biot coefficients of the reservoir. Considering that the scale of the reservoir is much larger than the three faults, the influence of different Biot coefficients of the faults is also expected to be little. We thus do not include the parameter in subsequent model calibration.

1st round calibration—Step II: Parameter calibration for the three faults
We first test different values for the thickness of the feeder fault. Both the displacement and temperature misfits are smallest when the values are equal to 300 m (Supplementary Table 9). Second, we calibrate the hydraulic parameters of porosity and permeability. The displacement and temperature misfits are lowest with the two parameters equal to 20% and 20 mD (Supplementary Tables 9 and 10), respectively. Third, the two elastic parameters, Young’s modulus and Poisson’s ratio, are initially assigned with typical values for faults34–36 and then optimized, but the simulation results are insensitive to changes in these parameters. Fourth, we test the influences of three thermophysical parameters: expansion coefficient, conductivity, and specific heat. The displacement and temperature misfits vary little with different thermal expansion coefficients, and are smallest with the initial settings. As the two elastic parameters and thermal expansion coefficient of the feeder fault have little influences on simulation results, the effects of the three parameters of the two strike-slip (plate boundary and subsidiary) faults with thickness less than the feeder fault are also expected to be little. We thus only calibrate the other five parameters for the strike-slip faults: thickness, two hydraulic parameters, thermal conductivity, and specific heat (Supplementary Tables 11-14). Simulation results show that the displacement and temperature misfits are smallest with the initial settings.

1st round calibration—Step III: Simulation tests with different scenarios of heat flux upwelling
We test two scenarios: (I) only the feeder fault serving as the upwelling conduit of heat fluid; (II) both the feeder and subsidiary faults serving as the upwelling conduit. The rates of heat flux are also tested with a wide range from 1 kg/s to 100 kg/s. Under scenario I, the RMS misfits of displacement and temperature are smallest with flow rates less than 20 kg/s (Supplementary Table 15). Under scenario II, the minimum RMS misfit of displacement is 2.54 cm, ~4% larger than the least value of scenario I, indicating little heat flux upwelling through the subsidiary fault. The rate of heat flow within the feeder fault is probably less than 20 kg/s.

1st round calibration—Step IV: Simulation tests with different properties of water
We calibrate five parameters including density, viscosity, compressibility, thermal conductivity and specific heat for water. The RMS misfits of displacement and temperature vary with different values of density, viscosity, compressibility and specific heat but the parameter of thermal conductivity has no influence on simulation results (Supplementary Table 16). Our calibration shows that the initial parameters settings for the water properties are suitable for simulating geothermal operations of the HGF.

2nd round calibration
Due to the possibility that the optimal parameter values obtained from one step of the iterative calibration may change in later steps, we thus conduct a second round of model calibration to further optimize the parameters with prominent influences on simulated displacement and temperature. We first recalibrate the hydraulic parameters and thermal expansion coefficients of the reservoirs (Supplementary Tables 4-6 and 8). Simulation results shows that both displacement and temperature misfits are lowest with the optimal parameter values determined from the first round of calibration, with no additional reductions in the misfits. As a consequence, we stop the optimization of model parameters.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Jiang et al. hypothesize that Thermal contraction drives long-term deformation of the Herber geothermal field in southern California. To test this hypothesis, they employ a thermo-poro-mechanical model and investigate surface deformation, injection/production, and hydrogeological datasets.\nThe hypothesis is plausible, and the overall approach is suitable. However, I found several shortcomings in this manuscript that I detailed below.\n1- The manuscript is poorly structured, and observations, results, and interpretations are mixed and confusing. I read the manuscript several times to understand what the authors tried to do and learn the paper's logic. Authors need to heavily restructure the manuscript. Also, in every figure, they need to introduce all symbols, signs, and colors.\n2- Throughout the manuscript, errors and uncertainties are disregarded. Without knowledge of observation errors and model uncertainties, it is not possible to distinguish between competing mechanisms (e.g., thermal contraction, pore pressure, and poroelastic stresses) that drive surface deformation.\n3- The implemented model is not adequate for studying geometrical sites. As the authors provided in Figure 1b, a geothermal site is a zone of heat anomaly and, thus, different from the surrounding area in terms of mechanical properties and heat content. As is shown in similar cases (e.g., Im and Avouac et al. 2021 Natue ), the mechanical model needs to be set up to account for property contrast within the reservoir and surrounding rocks. The model setting used here implements basic layered rheology, which I believe is inappropriate; thus, results are unreliable. This might also be why modeled vertical deformation deviates from observed in most places, particularly following 2010.\n4- I could not find a model misfit plot in the manuscript that shows how the model and observations agree. Figure 2 and supplemental figures that show modeled and observed vertical land motion suggest that agreement is not good, particularly following 2010. Need to provide maps of model misfit distributions for different episodes and discuss the possible mismatches and causes. Also, I am unsure why vertical deformation is normalized (e.g., figure 2). For simplicity, authors can show observed and modeled vertical deformation and misfits along two profiles (south-north and east-west) for different time steps.", "reviewer_2": "Review of the manuscript “Thermal contraction gradually dominates Heber geothermal field deformation”, by Jiang et al.\nThis is an interesting study that shows the importance of thermal effects on the subsurface response to geothermal energy operations. The study focuses on the commonly overlooked thermal effects related to cooling around injection wells. Despite thermal effects usually receive little attention, the findings in this paper are known and the Discussion could benefit from putting in context the findings of this study with published work. In particular, the seismicity rate increase in the long term in geothermal systems has been shown to be due to cooling by Parisio et al. (2019); cooling effects on induced seismicity have also been shown at The Geysers (Jeanne et al., 2015); and the effect of cooling in the far field, i.e., in areas not affected by cooling, has the potential to reactivate distant faults in the long term (Kivi et al., 2022). Overall, I think this is an interesting study that could be accepted for publication provided that the following comments are addressed.\nMajor comments:\n- Some relevant previous studies highlighting the importance of thermal effects in the long term are missing\n- The Abstract should emphasize more the “long-term” effect of cooling\n- Lines 145-150: it should be explained why there is subsidence from the beginning in the west region, where injection occurs and therefore uplift is expected. If subsidence is due to cooling, I would expect subsidence to be delayed for a few years until a relatively big volume of rock is cooled down and, thus, contracted\n- A major limitation of the numerical simulations is the fact that fluid properties are not a function of pressure and temperature. According to the explanation of lines 282-288, fluid properties are constant. This assumption yields significant errors in terms of pore pressure distribution in regions with cooling because water viscosity increases at lower temperatures, affecting pore pressure gradients\n- Lines 391-398: the reasons of the four observations should be explained\nMinor comments:\n- I recommend to avoid mixing units, e.g., ºC and ºF (better use ºC in the whole manuscript)\n- Figure 1d: the site has a tendency to subsidence, but this figure shows “normalized uplift”, which is confusing. It may be better to show “normalized subsidence” and also include the thermoelastic effect with negative values\n- Line 231-232: please, state here which are the “three thermophysical parameters”\n- Line 257 (and other lines in the manuscript): instead of “parameters”, the term “rock properties” could be used\n- Typos: line 53: “refer it to” -> “refer to it”; line 113: “benchmarks, and has been” -> “benchmarks, has been”\n- Caption Figure 4, b) I’d say that it should say “Temperature distribution at the same time” (without point)\n- Caption Fig. 5: indicate where the location of the 8 probe points can be found\nReferences\nJeanne, P., Rutqvist, J., Dobson, P. F., Garcia, J., Walters, M., Hartline, C. and Borgia, A. (2015). Geomechanical simulation of the stress tensor rotation caused by injection of cold water in a deep geothermal reservoir. Journal of Geophysical Research: Solid Earth, 120(12), 8422-8438.\nKivi, I. R., Pujades, E., Rutqvist, J. and Vilarrasa, V. (2022). Cooling-induced reactivation of distant faults during long-term geothermal energy production in hot sedimentary aquifers. Scientific Reports, 12(1), 2065.\nParisio, F., Vilarrasa, V., Wang, W., Kolditz, O. and Nagel, T. (2019). The risks of long-term re- injection in supercritical geothermal systems. Nature Communications, 10(1), 4391.", "reviewer_3": "Please see attached pdf file.\n[Editorial Note: This PDF is displayed across the next four pages]\nReview for Nature Communication of\nThermal contraction gradually dominates Heber geothermal field deformation\nBy G. Jiang et al.\nGeneral comments.\nThis manuscript present a very interesting case study of geothermal long-term monitoring combining deformation, hydrological and seismicity data and their interpretation using thermo-hydro-mechanical modeling, in relation with the harnessing of the Heber reservoir and geothermal field in California, US. The presented data set in simply impressive by its time duration, demonstrating brilliantly that long term observations is required to understand the dynamics of a reservoir. The modelling performed seems apparently fine, although it is hard to judge without trying the codes. The number of tests performed is also impressive and the supplementary materials is an encyclopedia of tests in support of the results.\nHowever, the manuscript suffers from issues on the content\n• The title indicates temperature effect as the main cause. However, this is not clearly shown. The demonstration of the results is quite poor, possibly mainly because shown results are not properly explained. For example, at places it is hard to understand the link between the figure caption and what is shown on the figures. In addition, although we may guess it, it is not clear what is the general aim of the manuscript, except demonstrating a powerful method with just a case study. As a consequence, there is a lack of generalisation. What does this study bring for other geothermal fields? What are the lessons learned?\n• There is often a mis-use of the word “deformation”, when referring to leveling and InSAR observations, which produce displacement data. This is general in the whole manuscript and also the Supplementary Materials. I have mentioned some of them, but please check thoroughly.\n• I do not understand the ground reason of your area separations between west, central, east zones. This seems to be very artificial to me… The demonstration can be done without these artificial boundaries.\n• You rise the point of global movements of the plate boundary, which is known to be pretty active, as your study addresses observations over several 10s years. I am not convinced that the your displacement inversion can assess whether or not the Dixieland fault can be the reason of the observed changes, the cover of leveling network and InSAR coverage is simply much to small. Indeed, Dixieland fault takes place in a much larger regional networks of faults such as the Imperial and Weinert-El Centro faults, which seem to be at least as active than this one. To convince the reader, you need to refer to other studies and/or include additional data to exclude that there were no global movements during the several years. Actually, if there were, how would that affect your study? You need to reconsider your demonstration. In addition, in the supplementary material, your inversion procedure lead me to conclude you may invert the wrong component (vertical), which you also corrected for signals possibly associated to what you want to invert for…\n• Table S4 to S16. Impressive number of tests performed. However, it seems that some minima are not found as yet by lack of search in the model space. For example table S4.I. Looking at the column 16% for the upper reservoir, we see that the lower reservoir misfit decreases when the the lower reservoir porosity decreases. May be if you would take 14% for the porosity of the lower reservoir, the misfit would be even smaller? How are you sure you found the minimal value of misfit for all parameter tested? This rises the question of the validity of the final result.\nAnd the form\n• There are many problems in the figures and their captions. Some figures are not clearly explained enough. Every single mention in the figure should be explained in the caption. Many are missing. Some figures are not showing what is said in the caption.\n• When comparing figures of the same area for observations and interpolation or inversion, I would strongly suggest to use the same axis. Otherwise, it is challenging to find correspondence.\nAs a whole, this is an interesting case study, however the main message is not brought clearly enough, because the demonstration suffers from lack of clarity and some clear evidence that you reached the best model ultimately. In addition, whether the conclusion found can be generalized to other cases and how is not brought.\nI suggest that the above mentioned issues and the specific comments below need to be solved before any decision can be taken on the manuscript.\nSpecific comments.\nLine 49-52. The permeability is found along faults, which are described. However, it is not clear enough which is which with respect to the figure 1. In the text you do not mention the feeder fault, is it one on them?\nLine 49-51. You describe the matrix permeability from 0.55 km until ~2 km (?) but you do not indicate what is above or below. This is lacking also here.\nLine 51. “The normal fault” is very vague. Which one is it? How do you know it is a normal fault? What does it bring to the study?\nFigure 1a. Please indicate what “F.” means. In the caption, HGF meaning is missing. It should be given at the first occurrence in figures. HGF is messing in caption 1a. Same applies for SCSN: what is it? What are yellow circles? What are the beach balls? The plate boundary is certainly not so linear, is it? It would be more realistic to represent it with its real path, as feeder and subsidiary faults.\nFigure 1d. The red line is usually called “cumulative seismicity” However, it is also often shown with the energy released. This information may not be available, as seismic events in geothermal system are small and estimating their energy sometimes challenging, but quantification would be nice here; it would be nice to indicate what is the variability of energies of the earthquakes shown. In the leveling time series, why the error bar increases with time and there is no error on the first point? The cyan color is hard to follow. May be plot it on top of the leveling points?\nLine 70-71. I do not get it. The resource of a geothermal field cannot be designed… The resource is what we try to access, so we design the capacity of the power plant to match the resource potential. This sentence needs to be rephrasing. Or am I missing something here?\nLine 72. There is a tendency now to try not to “exploit” natural resources, but rather harness. Could you check in all the manuscript?\nLine 85-86. Possibly you need a reference here.\nLine 86. It would be important to give information on the seismic network (number and type of stations) deployed since the beginning of the exploration and over the years. This would give us an impression of if there is more earthquake due to the better recording capability of indeed because of more earthquakes.\nLine 89. I would suggest to change “indicating” with “suggesting”.\nLine 94-96. Actually in this paragraph, there is nothing new, this is known. Good to know for HGF, but so what?\nLine 114, 139, 235, and other locations possibly, e.g., supplementary figure captions. A space missing between time and series.\nFigure 2. The reason(s) of the separation of west, central and east regions is far for clear. From what is explained, there are no structural structures that justify this separation. Why not north, central and south?\nFigure 2e. I do not understand why you need to normalize the vertical displacements. Note the term “vertical deformation” is inappropriate, as what is observed with both leveling and InSAR are displacements, and actually you indicate also “per year”, which makes it technically velocities.\nLine 148. “The same scenario”. This is confusing. Which same scenario? In addition, if this is the same scenario between two areas separated the artificial boundaries (west and east), it means that there may be no need to separate them…\nFigure 3. The figure shows a series of temperature points. However, there is not clarity... for instance in Fig. 3a, I do not see 11 injection wells values. Are there all mixed in the figure? Are they averaged? Or only one is shown? Similar issues for the other sub-figures.\nLines 169-171. These interesting observations are hard to believe, but is I generally trust data. This paragraph is just describing data, but do not suggest any potential reasons for the temperature changes. Possibly this is not the place, however, after displacement observations, you suggest a reason for the observed changes. Therefore I would expect the same for temperature observations.\nAt line 178-180, there seems to be indeed a reason proposed: convection. Is that the case? If yes, I would make is even clearer.\nLine 183. Indeed an active fault may move. It would be also interesting to relate the natural seismic activity along the entire Diexiland fault with time, to make the potential link between the local seismicity and the more regional one.\nLine 186. Two major issues here:\n1. I would be very cautious to invert local data for global interpretation.\n2. Inversion of leveling data does not account for horizontal displacements. However, due to the local geodynamics setting strike-slip faults (e.g., the Dixieland fault?) mostly horizontal displacement would be expected.\nLine 197. How did you measure horizontal displacements? The explanation resides in the InSAR/leveling integration. At places, this is not clear enough on the temporal evolution of the measurements and how you inverted each period of time. How to combine both of them is sometimes vague.\nLine 201. Minor point: I would suggest you indicate the code used here also (although is it given in sup info)\nLine 223. “strike-slip” it is not clear which fault is which mechanism. You could introduce them better in figure 1.\nLine 228. Why inverting for layer parameter is not feasible? Does it not depend on which ones? You need more explanations here.\nFigure 4. What are P1, P2, P8 etc. Please report extensively all elements within the figure caption.\nFigure 5. This figure is very unclear, mainly because the sub-figure are too small, and axis are really tiny. In addition, it is not possible to understand what are the different curves of different colors.\nLine 358. Figure 15 does not exist.\nLine 491. The initials of Jean-Philippe Avouac are J. P.\nSupplementary Materials:\nPlease check the use of “deformation”.\nIn § Fault slip inversion: I have an issue on what you want to invert for and the procedure. You remove the geothermal signal. Why and how do you do this? Then you subtract the pre-2005 trends. I really do not understand why. This may be the signal you want to invert for… You may need a reference here to make sure they are the whole area signal. After removal, what is left??\nin the second step and (1), you use “corrected vertical-component data”. Here 2 things: 1. what are data corrected for? The previous trends? 2. why do you invert only the vertical data? Inversion procedure with Okada model is clearly possible with all components, so I do not understand why only vertical. In addition, as the plate boundary is a strike-slip process, there may be little vertical signal but large horizontal one. So this inversion results may be completely wrong, as you invert the wrong component in which you removed some potentially valid signal… You need to review here seriously the approach.\nIn the thermo-poro-electicity formulation, in the last paragraph you indicate primary effects and secondary effects. You could list them, and justify how to classify them as primary or secondary.\nFigure caption of S1, S3. Check the use of “deformation” here.\nFigure caption of S1. The normalization method for uplift is unclear. Why do you need to normalize? Is the normalized value inverted for??\nFigure S2. In the caption, please add the reference for the kriging here too. In the figure, you should get the same axis limits for both the data and the interpolated info.\nFigure S5. The sub figure are too small, making this figure useless. You also need to tell readers when… time is not indicated or clear enough.\nFigure S6. There are several place where the term “excess” is used for vertical displacement. I am not convince with excess it is. I did not get where this excess is coming from. Need more explanation here and before possibly.\nFigure S7. Where is A-33? Show it on the figures.\nFigure S16, S17, S18. The number of sub figure is large. You could gain place by removing the axis for the inner figure, as they are the same for all. This would make the figure much easier to read. In addition, on the all 3 pages for each figure, I count 96 (4*8*3) points. Not clear where the 3 otehrs are…"}, "author_response": "Response to Reviewer #1\nJiang et al. hypothesize that Thermal contraction drives long-term deformation of the Herber geothermal field in southern California. To test this hypothesis, they employ a thermo-poro-mechanical model and investigate surface deformation, injection/production, and hydrogeological datasets. The hypothesis is plausible, and the overall approach is suitable. However, I found several shortcomings in this manuscript that I detailed below.\nAns: Thank you very much for your efficient work of reviewing our MS and also for your constructive comments.\nAs you know that deformation of geothermal fields can be caused by reservoir pressure depletion, poroelastic and thermoelastic contraction, but the spatial and temporal evolution of their strengths is still in doubt. One major reason for this knowledge gap is due to the rarity of long-term seismic and geodetic observations combined with the extreme difficulty of setting up a realistic 3D thermo-hydro-geological model with limited subsurface geological and geophysical data (like reservoir geometry, formation lithology and fault distribution) needed to simulate decades of geothermal operations.\nIn this study, we quantify spatiotemporal evolution of the relative strengths of the three physical mechanisms for the Heber Geothermal Field (HGF) based on a 3D thermo-hydro-geological model, which is calibrated by multiple datasets, including geodetic observations over 20 years, precise well trajectories and operational records of fluid rates and temperatures. Our most important findings are:\n1. The effects of pressure flunctation and poroelastic response keep relatively stable in the operation history of the HGF. In contrast, the strength of thermal contraction effect increases with time: after ~30 years of operations, the thermal effect is 2-3 times larger than the two pressure effects.\n2. It is common that fluid injection causes surface uplift due to pore pressure increases within the target reservoir; however, the ground in west an east regions of the HGF was subsiding during some time periods with net increases of injected fluid. We show that the deformation anomaly is related to the effect of the HGF center siphoning fluid from surrounding regions and ever-increasing effect of thermal contraction.\n3. The long-term trends of surface deformation and seismicity growth at the HGF are controlled by the thermal effect. In contrast, the two pressure effects only drive instantaneous changes.\nWe have thoroughly revised the MS based on the comments from you and the other two reviewers to make these points as clear as possible. Following are our point-by-point response.\n1-1 The manuscript is poorly structured, and observations, results, and interpretations are mixed and confusing. I read the manuscript several times to understand what the authors tried to do and learn the paper's logic. Authors need to heavily restructure the manuscript.\nAns: We are really sorry that the original MS was not structured with a straightforward logic. During this revision, we rewrite the main text into three sections: (1) Introduction (without subheading), (2) Results with 7 subheaded sections, and (3) Discussion, based on the formatting instructions of Nature Communications.\nIn the Introduction section, we begin with the background of our study, and then present the operation history of the HGF, and a brief summary of our major results and conclusions of our study.\nFor the Results section, we formulate it with the following logic: (1) presenting the observation results of surface displacement and seismicity with two sections; (2) presenting the results of fault slip modeling (to test whether the observed surface displacement is due to potential slip on the two strike-slip and one normal-slip faults); (3) presenting the results of Thermo-Hydro-Mechanical (THM) modeling (including model setup, model calibration, and model fit); (4) presenting the simulation results of reservoir pressure and temperature; (5) presenting the quantification results of spatio-temporal evolution of Thermo-Poro-Elastic effects; (6) presenting the physical mechanisms of the HGF deformation.\nIn the Discussion section, we compare our major findings with previous studies, and discuss their implications for the safety of geothermal operation and mitigation of geohazards.\n1-2 Also, in every figure, they need to introduce all symbols, signs, and colors.\nAns: Thank you for your suggestion. For this revision, we revised the figures to match the revisions to the main text, and we have made sure the meaning of each symbol is introduced in the captions.\nAns: In our study, we used leveling and InSAR observations of surface displacement as well as well temperature records to constrain our 3D THM model. The observation errors of leveling data range from 0.2 cm to 0.5 cm. The uncertainties of InSAR-based velocities are less than 2 mm/yr depending on the satellite and observation periods. In the figure below (Response Figure 1), you can see that the errors are much less the displacement magnitudes of the HGF, especially in the main deforming region where signals are nearly an order of magnitude larger than the uncertainty. More importantly, as we show in the main text, the deformation trends measured by leveling are consistent with the InSAR results, which confirms that the surface displacement observations are robust.\nResponse Figure 1: Uncertainty in InSAR based ground veloicities compared to the size of the velocity estimate. Top: East rates. Bottom: Vertical rates. ENV=Envisat; SNT1a=Sentinel-1 early time period (2015/05-2017/12); SNT1b=Sentinel-1 later time period (2018/06-2019/08).\nFor the temperature records, we do not have any information on their observation errors because they come from a public repository provided by the local regulatory body. But, you can see the records of different wells in each cluster are consistent with each other (Response Figure 2).\nResponse Figure 2: Operation temperature changes of the Heber Geothermal Field.\nUnfortunately, it is computationally impractical to estimate the uncertainties of a FEM model, especially of a 3D model, and we did not find any related references offering a way to roughly estimate these. Alternatively, to improve the robustness of our results, we use the observations of surface deformation and well temperatures to calibrate our 3D model. In addition, we also test the influences of different parameters. Our calibrated model can fit the spatiotemporal variations of observed surface displacement and well temperatures, which indicates that we can use the model to quantify the spatiotemporal evolution of thermo-poro-elastic effects.\n3-1 The implemented model is not adequate for studying geometrical sites. As the authors provided in Figure 1b, a geothermal site is a zone of heat anomaly and, thus, different from the surrounding area in terms of mechanical properties and heat content. As is shown in similar cases (e.g., Im and Avouac et al. 2021 Natue ), the mechanical model needs to be set up to account for property contrast within the reservoir and surrounding rocks. The model setting used here implements basic layered rheology, which I believe is inappropriate; thus, results are unreliable.\nAns: Yes, a geothermal site is a region of heat anomaly; but it does not indicate that the hydrogeological properties of the heat region are different from the surrounding rocks. The heat anomaly is mainly associated with the location of deep source of heat fluid. For the HGF, current knowledges on the lateral variations of rock properties are limited. To honor the natural complexity of the system while keeping the problem numerically tractable, we simplify the study region with five layers: caprock at the surface, two reservoir layers, a basal layer, and basement rock. Additionally, we embed the two NW-trending strike-slip faults and the NE-trending feeder fault in our model as individual formations with a unique set of thermo-hydro-mechanical properties. The NE-trending feeder fault is the primary pathway for deep, high-temperature brine upwelling into the shallow reservoir.\nIn our 3D model, the parameter settings mainly reflect the rock properties within the geothermal field rather than the surrounding region. The major reasons for such setting are: (1) obvious surface deformation only occurred within the field; (2) we only focus on the deformation mechanisms of the field; (3) only the deformation observations and temperature records within the field are used to calibrate our model; (4) we know little about the rock properties of the surrounding region.\n3-2 This might also be why modeled vertical deformation deviates from observed in most places, particularly following 2010.\nAns: During our model calibration, we simulate the operation of the HGF from Apr 1985 through Dec 2010; the simulation results are compared with the vertical displacement timeseries (Jan 1994 - Nov 2010) at 64 leveling benchmarks, which are are selected to give approximately uniform coverage within the HGF. The leveling observations after 2010 are not included for model calibration, and only used to validate our model.\nWe count the number of leveling stations with bad misfit of the model predictions to the observations. The total number is 30 (Figure S14). Most of these stations are located in the peripheral region of the HGF. In addition, there are 3 stations (1-A-5, 1-A-4, and C-4) located at the HGF center. The difference between the simulated and observed displacements in the peripheral regions is likely due to laterally inhomogeneous hydromechanical and thermophysical properties, and/or unmodeled variations in the depth to basement rock of the study region. The effect of inhomogeneous rock properties increases with the duration of geothermal simulation. For the difference between the simulated and observed displacements at the HGF center, there is another reason that the leveling-observed subsidence is obviously larger than the InSAR monitoring results (Figures S2 and S3). Both of the two observation results are included during model calibration.\nOverall, our model can replicate major features of observed surface displacement in both vertical and horizontal directions (Figs. 6 and S13).\n4-1 I could not find a model misfit plot in the manuscript that shows how the model and observations agree. Figure 2 and supplemental figures that show modeled and observed vertical land motion suggest that agreement is not good, particularly following 2010. Need to provide maps of model misfit distributions for different episodes and discuss the possible mismatches and causes.\nAns: Thank you for your suggestion. We now plot a figure (Figure 6 in the main text) to show model predictions to surface deformation observations and model misfits.\n4-2 Also, I am unsure why vertical deformation is normalized (e.g., figure 2).\nAns: The displacement time series of the leveling stations are normalized for two purposes primarily related to visualization: (1) to represent the deformation trends of the west, central and east zones of the HGF, and (2) to represent temporal variation of deformation of the whole region. By doing this, we can compare vertical deformation trends of the west, central and east regions with the respective net injection volumes of geothermal fluid of every month. In addition, we can also compare the deformation trend of the whole region with temoporal changes of observed seismicity in the Heber geothermal field. To be clear, our misfit calculations are based on actual observations.\n4-3 For simplicity, authors can show observed and modeled vertical deformation and misfits along two profiles (south-north and east-west) for different time steps.\nAns: Thank you for your suggestion. We plot the displacement profiles during this revison. Given that Figure 6 in the main text has shown the model predictions to surface deformation observation. We put the figure of displacement profiles in supplementary materias (Figure S15).\nResponse to Reviewer #2\nThis is an interesting study that shows the importance of thermal effects on the subsurface response to geothermal energy operations. The study focuses on the commonly overlooked thermal effects related to cooling around injection wells. Despite thermal effects usually receive little attention, the findings in this paper are known and the Discussion could benefit from putting in context the findings of this study with published work. In particular, the seismicity rate increase in the long term in geothermal systems has been shown to be due to cooling by Parisio et al. (2019); cooling effects on induced seismicity have also been shown at The Geysers (Jeanne et al., 2015); and the effect of cooling in the far field, i.e., in areas not affected by cooling, has the potential to reactivate distant faults in the long term (Kivi et al., 2022). Overall, I think this is an interesting study that could be accepted for publication provided that the following comments are addressed.\nAns: Thank you very much for your careful review and constructive comments.\nDuring this revision, we read through the three papers you recommended and other literatures like Simone et al. (2013; http://dx.doi.org/10.1016/j.pce.2013.01.001). The major findings of these studies are summarized below:\n(1) Simone et al. (2013): thermal effects can induce a significant perturbation on the stress in the intact rock affected by the temperature drop; this perturbation is likely to trigger induced seismicity in the surroundings of critically oriented fractures near the injection well.\n(2) Jeanne et al. (2015): thermal effects can cause stress tensor orientation changes; the vertical stress reduction is significant and propagates far below the injection well.\n(3) Parisio et al. (2019): thermal effects dominate the geomechanical response of geothermal reservoirs.\n(4) Kivi et al. (2022): thermal stresses are transmitted much ahead of the cooled region and are likely to destabilize faults located far away from the reinjection well.\nIn contrast, our study focuses on not only the thermal effect but also the pressure effects of pressure fluctuation and poroelastic response. We try to resolve the physical mechanisms of observed surface displacement, especially of two abnormal phenomena: in the west and east regions of the Heber Geothermal Field, where only injection operations occur, the ground was subsiding from 1994 to 2005 and after 2005, respectively. In addition, we also strive to quantify the spatiotemporal evolution of relative strengths of three physical mechanisms: reservoir pressure depletion, poroelastic effects, and thermoelastic contraction.\nAlthough we now know that deformation of geothermal fields can be caused by the pressure and thermal effects, the spatial and temporal evolution of their strengths has not been well understood. One major reason for this knowledge gap is due to the rarity of long-term geodetic observations combined with the extreme difficulty of setting up a realistic 3D thermos-hydro-geological model with limited subsurface geological and geophysical data (like reservoir geometry, formation lithology and fault distribution) needed to simulate decades of geothermal operations. With exception to Jeanne et al. (2015), the published studies above were conducted with 2D conceptual models. While they represent important advances for the reasons mentioned above, we feel that they are not sufficient to characterize the spatiotemporal variation of different physical mechanisms. In the discussion section of the revised manuscript, we now cite the papers on the thermal effect and compared our results with their findings.\nFollowing are our point-by-point responses to the more specific comments:\nMajor comments:\n1. Some relevant previous studies highlighting the importance of thermal effects in the long term are missing\nAns: Thank you for your suggestion. We cited the following papers on the thermal effect:\n(1) De Simone, S., Vilarrasa, V., Carrera, J., Alcolea, A. & Meier, P. Thermal coupling may control mechanical stability of geothermal reservoirs during cold water injection. Phys. Chem. Earth 64, 117–126 (2013).\n(2) Jeanne, P., Rutqvist, J., Dobson, P. F., Garcia, J., Walters, M., Hartline, C. and Borgia, A. (2015). Geomechanical simulation of the stress tensor rotation caused by injection of cold water in a deep geothermal reservoir. Journal of Geophysical Research: Solid Earth, 120(12), 8422-8438.\n(3) Kivi, I. R., Pujades, E., Rutqvist, J. and Vilarrasa, V. (2022). Cooling-induced reactivation of distant faults during long-term geothermal energy production in hot sedimentary aquifers. Scientific Reports, 12(1), 2065.\n(4) Parisio, F., Vilarrasa, V., Wang, W., Kolditz, O. and Nagel, T. (2019). The risks of long-term re-injection in supercritical geothermal systems. Nature Communications, 10(1), 4391.\nWe think this highlights the importance of thermal effects, as suggested, and also helps clarify the novelty of our study.\n2. The Abstract should emphasize more the “long-term” effect of cooling\nAns: Thank you for your suggestion. In the revised abstract, which has a limit no more than 150 words, we emphasized the effects from the following two aspects:\n(1) thermal contraction increases gradually with time and eventually overwhelms the pressure effects of pressure fluctuation and poroelastic response, which keep relatively stable during geothermal operations;\n(2) thermal contraction dominates long-term trends of surface displacement and seismicity growth, while the pressure effects drive near-instantaneous changes.\n3. Lines 145-150: it should be explained why there is subsidence from the beginning in the west region, where injection occurs and therefore uplift is expected. If subsidence is due to cooling, I would expect subsidence to be delayed for a few years until a relatively big volume of rock is cooled down and, thus, contracted\nAns: The reason why ground in the west region was subsiding from 1994 to 2005 with a net increase of injected fluid is that the average injection rate (0.6×109 kg/month) in the west region is much less than the production rates in the central region (~3.5 × 109 kg/month). The high production rates at the HGF center caused a pressure drop of up to 3 MPa relative to the initial pressure of 14.7 MPa (Response Figure 3); this created a low-pressure center that siphons fluid from surrounding regions. The fluid volume siphoned by the HGF center from the west region are more than the injection volume before 2005.\nResponse Figure 3: Simulated pressure changes at a time node during 1994 and 2005.\n4. A major limitation of the numerical simulations is the fact that fluid properties are not a function of pressure and temperature. According to the explanation of lines 282-288, fluid properties are constant. This assumption yields significant errors in terms of pore pressure distribution in regions with cooling because water viscosity increases at lower temperatures, affecting pore pressure gradients\nAns: During this revision, we discussed this problem with Josh Taron from USGS, and add him to be our coauthor. As the functions of fluid properties changing with pressure and temperature are not available, we instead test the influences of different values of water properties on simulation results of surface displacement and reservior temperature. From table S16, you can see that the influence is little.\nTable S16. Misfits of model predictions with different properties of water.\nP.S. The parameter values in our tests are set based on the reservoir temperature (injection: 75-85°C; production: 150-175°C) and pressure (~12-18 MPa).\n5. Lines 391-398: the reasons of the four observations should be explained\nAns: Thank you for your suggestion. The reasons were presented in the last section of the original MS (Spatiotemporal evolution of thermo-poro-elastic effects). We would like to illustrate the reasons for the interesting phenomena here.\n(1) The first phenonemon that ground was subsiding during some time periods with pressure increase, especially at points P5 and P6, is due to the effect of thermal contraction.\n(2) The two phenomena that surface subsidence at point P3 with ~1.3 MPa of pressure decrease and 1.6°C of temperature decrease is comparable to the subsidence at point P7 with ~1.8 MPa of pressure increase and 53.5°C of temperature decrease, and that surface subsidence at point P6 is up to 34 cm with only 3°C of temperature decrease and ~1.5 MPa of pressure increase, are also caused by the cooling effect, whose strength depends on temperature decrease and varies with locations.\n(3) The last phenomenon that temperature decrease at point P5 started from 2009, ~13 years later than the time of ground subsidence occurred, is caused by the fact that the thermal contraction effect can disturb broad ranges over the temperature decrease front.\nResponse Figure 4: Temporal evolution of pressure, temperature and their effects at four probe points.\nOur original intention of introducing these phenoma is to show the features of the thermal effect:\n(1) thermoelastic effect varies nonlinearly and asynchronously with temperature;\n(2) thermal contraction effect can disturb broad ranges over the temperature decrease front.\nThese two features cause prominent thermal effect at some locations with little temperature decrease.\nDuring this revision, we realized that the original describtion is hard for readers to follow, and rewrote the paragraph. You can find it in the last paragraph of the section\n“Quantification of inhomogeneous Thermo-Poro-Elastic effects”.\nMinor comments:\n1. I recommend to avoid mixing units, e.g., ºC and ºF (better use ºC in the whole manuscript)\nAns: Thank you for your suggestion. We made units consistent, replacing the unit of ºF with ºC where necessary.\n2. Figure 1d: the site has a tendency to subsidence, but this figure shows “normalized uplift”, which is confusing. It may be better to show “normalized subsidence” and also include the thermoelastic effect with negative values\nAns: Thank you for your suggestion. The lable “normalized uplift” is replaced with “normalized Uz”.\nThe vertical displacements (Uz) associated with the thermoelastic effect are shown with negative values in the new Figure 4c.\n3. Line 231-232: please, state here which are the “three thermophysical parameters”\nAns: Thank you for your suggestion. The three thermophysical parameters are specific heat, thermal conductivity and thermal expansion coefficient. We added them into the corresponding sentence in the revised MS.\n4. Line 257 (and other lines in the manuscript): instead of “parameters”, the term “rock properties” could be used\nAns: Thank you for your suggestion. The word “parameters” was replaced with “rock and water properties” in the revised MS. In addition, we eliminated the use of ambigious words.\n5. Typos: line 53: “refer it to” -> “refer to it”; line 113: “benchmarks, and has been” -> “benchmarks, has been”\nAns: Thank you for your careful review. We corrected these errors during this revision.\n6. Caption Figure 4, b) I’d say that it should say “Temperature distribution at the same time” (without point)\nAns: Thank you for your suggestion. We revised the sentence in the figure caption.\n7. Caption Fig. 5: indicate where the location of the 8 probe points can be found\nAns: Thank you for your suggestion. During this revision, we replaced the original Figure 5 with new Figures 7 and 8, where we plot the location of the 8 probe points.\nResponse to Reviewer #3\nThis manuscript present a very interesting case study of geothermal long-term monitoring combining deformation, hydrological and seismicity data and their interpretation using thermohydro-mechanical modeling, in relation with the harnessing of the Heber reservoir and geothermal field in California, US. The presented data set in simply impressive by its time duration, demonstrating brilliantly that long term observations is required to understand the dynamics of a reservoir. The modelling performed seems apparently fine, although it is hard to judge without trying the codes. The number of tests performed is also impressive and the supplementary materials is an encyclopedia of tests in support of the results.\nAns: We appreciate your careful review and insightful comments, which indeed help improve the quality of our MS. We also appreciate the recognition that long-term observations and sophisticated numerical models are needed to understand geothermal hazards. For some additional context:\nWe started this project in 2017, and spent great efforts to collect and analyze Leveling and InSAR observations of more than 20 years. We found that the west, central and east regions of the Heber Geothermal Field (HGF) experienced different deformation trends. In doing so we uncovered a striking and previously unrecognized phenomena of deformation trend reversals. Further, we found that the west and east regions were subsiding during the time periods with a net increase of injected fluid. In addition, we found that the growth of observed seismicity is consistent with surface displacement trend through normalization analysis of leveling time-series.\nIt is also time-consuming and challenging work to set up a realistic 3D hydrogeological model for the HGF and to calibrate its parameters with multiple datasets including surface deformation observations, precise well trajectories and operational records of fluid rates and temperatures. The main purposes of our study are: (1) to quantify the spatiotemporal evolution of relative strengths of three physical mechanisms: reservoir pressure depletion, poroelastic effects, and thermoelastic contraction; (2) and to resolve the physical mechanisms of complex surface displacement trends observed in different locations and time periods.\nOur numerical simulations are implemented with the software COMSOL Multiphysics (version 6.0).\nThe model files are available from Zenodo (https://zenodo.org/records/10460410). We revised the MS based on the constructive and comprehensive comments from you and the other two reviwers. Following are our point-by-point responses.\nHowever, the manuscript suffers from issues on the content\n1-1. The title indicates temperature effect as the main cause. However, this is not clearly shown The demonstration of the results is quite poor, possibly mainly because shown results are not properly explained.\nAns: Thank you for this point. The meaning of the original title is that the strength of thermal effect increases gradually with time and eventually become the dominant cause. But, to avoid misleading the readers, we revised the title into “Relatively stable pressure effects and time-increasing thermal contraction control Heber geothermal field deformation”, based on our major findings:\n1. The effects of pressure flunctation and poroelastic response keep relatively stable in the operation history of the HGF. In contrast, the strength of thermal contraction effect increases with time: after ~30 years of operations, the thermal effect is 2-3 times larger than the two pressure effects.\n2. The long-term trends of surface deformation and seismicity growth at the HGF are controlled by the thermal effect. In contrast, the two pressure effects only drive instantaneous changes.\n3. It is common that fluid injection causes surface uplift due to pore pressure increases within the target reservoir; however, the ground in west an east regions of the HGF was subsiding during some time periods with net increases of injected fluid. Here we show that the deformation anomaly is related to the effect of the HGF center siphoning fluid from surrounding regions and ever-increasing effect of thermal contraction.\nWe apologize that the original MS was not structured with a straightforward logic. For this revision we significantly overhauled the structure of the text and revised the figures. We re-wrote the main text into three sections: (1) Introduction (without subheading); (2) Results with 7 subheaded sections; (3) Discussion, based on the formatting instructions of Nature Communications.\nIn the Introduction section, we begin with the background of our study, and then present the operation history of the HGF, and a brief summary of our major results and conclusions of our study.\nFor the Results section, we formulate it with the following logic: (1) presenting the observation results of surface displacement and seismicity with two sections; (2) presenting the results of fault slip modeling (to test whether the observed surface displacement is due to potential slip on the two strike-slip and one normal-slip faults); (3) presenting the results of Thermo-Hydro-Mechanical (THM) modeling (including model setup, model calibration, and model fit); (4) presenting the simulation results of reservoir pressure and temperature; (5) presenting the quantification results of spatio-temporal evolution of Thermo-Poro-Elastic effects; (6) presenting the physical mechanisms of the HGF deformation.\nIn the Discussion section, we compare our major findings with previous studies, and discuss their implications for the safety of geothermal operation and mitigation of geohazards.\n1-2. For example, at places it is hard to understand the link between the figure caption and what is shown on the figures.\nAns: Again, we apologize that our original figures were not presented clearly. For this revision, we reploted the figures to support/illustrate the major points of the main text. The captions were also rewritten to give clear explanations to each symbol in the figures.\n1-3. In addition, although we may guess it, it is not clear what is the general aim of the manuscript, except demonstrating a powerful method with just a case study. As a consequence, there is a lack of generalisation. What does this study bring for other geothermal fields? What are the lessons learned?\nAns: Thank you very much for these points. It is clear to us now that the manuscript needed to be restructured to make the motivation clear. Our study has two major aims: (1) to quantify the spatiotemporal evolution of relative strengths of three physical mechanisms: reservoir pressure depletion, poroelastic effects, and thermoelastic contraction; (2) and to resolve the physical mechanisms of complex surface deformation trends observed in different locations and time periods.\nWith the revised manuscript, we improved the generalisation of our MS by comparing our major findings with previous studies. In addition, we add Discussion section to introduce the implications of our method and major findings. To reiterate:\n(1) As you know, although the pressure and thermal effects have been recognized as the causal mechanisms of geothermal field deformation, the spatial and temporal evolution of their strengths was not well understood. One major reason for this knowledge gap is due to the rarity of long-term seismic and geodetic observations combined with the extreme difficulty of setting up a realistic 3D thermo-hydro-geological model with limited subsurface geological and geophysical data (like reservoir geometry, formation lithology and fault distribution) needed to simulate decades of geothermal operations. Our framework developed for resolving the spatiotemporal evolution of the HGF deformation mechanisms through integrating multiple geodetic, geophysical and geological data is applicable for other complex geotherm systems.\n(2) Benefited from the 3D realistic THM model of the HGF, we find distinct features of pressure and thermal effects. First, the pressure effects respond much faster to changes in geothermal operations than the thermal effect. Second, different from the instantaneous response and proportional relation of poroelastic effects to pressure fluctuation, the effect of thermal contraction varies nonlinearly and asynchronously with temperature. Third, the pressure effects keep relatively stable during geothermal operations whereas the thermal effect increases gradually with time and eventually overwhelms the pressure effects. Lastly, thermal contraction dominates long-term growth of seismicity, while the pressure effects drive near-instantaneous changes.\nThese major findings advance current understanding of deformation mechanisms associated with geothermal energy production. In general the pressure effects are dominant at the early stage of geothermal operations, whereas the thermal effect increases with time and eventually dominates.\nIn addition, these findings also provide theoretical implications for mitigating geohazards at geothermal fields and for improving the sustainability and safety of geothermal operations.\n2. There is often a mis-use of the word “deformation”, when referring to leveling and InSAR observations, which produce displacement data. This is general in the whole manuscript and also the Supplementary Materials. I have mentioned some of them, but please check thoroughly.\nAns: Thank you for your careful review. We agree, and where appropriate we replaced “deformation” with “displacement” when describing the leveling and InSAR observation results.\n3. I do not understand the ground reason of your area separations between west, central, east zones. This seems to be very artificial to me… The demonstration can be done without these artificial boundaries.\nAns: We divide the study area into the three zones based on two main reasons: (1) the west, central and east regions of the Heber geothermal fields experienced different displacement trends (for more details please refer to our response to your general comments); and (2) there are differences in geothermal operations across the field, with only injection operations in the west and east zones, but both injection and extration operations occurred in the central zones.\n4. You rise the point of global movements of the plate boundary, which is known to be pretty active, as your study addresses observations over several 10s years. I am not convinced that the your displacement inversion can assess whether or not the Dixieland fault can be the reason of the observed changes, the cover of leveling network and InSAR coverage is simply much to small. Indeed, Dixieland fault takes place in a much larger regional networks of faults such as the Imperial and Weinert-El Centro faults, which seem to be at least as active than this one. To convince the reader, you need to refer to other studies and/or include additional data to exclude that there were no global movements during the several years. Actually, if there were, how would that affect your study? You need to reconsider your demonstration. In addition, in the supplementary material, your inversion procedure lead me to conclude you may invert the wrong component (vertical), which you also corrected for signals possibly associated to what you want to invert for…\nAns: Thank you for these important comments. In our study region the leveling benchmarks are distributed within a region of ~10km*10km, but the coverage area of InSAR data is much larger, over ~20km*20km. These observations, especially InSAR, have good coverage of the southeast extension of the plate boundary (Dixieland) fault. The surface displacements associated with potential fault slip should be distributed along the faults. However, we only observed surface displacements at the HGF, which are predominatly subsidence and uplift. In the other regions, we did not observe any deformation signals, which can represent the “global movements”. In other words, the observed geodetic anomalies occurred within a localized zone within the HGF on a far-smaller scale than the Dixieland fault.\nTo convince the reader, we further conduct three scenarios of fault slip modeling to available displacement observations. First, the long-term vertical displacement measured by leveling from 1994 to 2004 is inverted with normal slip on the feeder fault. Second, the vertical transient displacement from 2006 to 2010, obtained through removing the long-term displacement trends, is inverted with normal slip on the feeder fault. Third, the horizontal displacement measured by InSAR from 2006 to 2010 is modelled with the strike-slip Dixieland fault. These modeling can prove whether the long-term and transient displacement is associated with fault slip.\nOur inversions show that neither the long-term vertical displacement nor the transient displacement observation can be reasonably fit by normal slip on the feeder fault. For the InSAR-measured horizontal displacements, the magnitudes of dextral strike-slip on the plate boundary fault need to be up to 30 to 82 cm below the depth of 8 km, which suggests 6 to 16 cm/yr of slip rates, much higher than the reported slip rates and the relative plate motion rate (only ~4 cm/yr).\n5. Table S4 to S16. Impressive number of tests performed. However, it seems that some minima are not found as yet by lack of search in the model space. For example table S4.I. Looking at the column 16% for the upper reservoir, we see that the lower reservoir misfit decreases when the the lower reservoir porosity decreases. May be if you would take 14% for the porosity of the lower reservoir, the misfit would be even smaller? How are you sure you found the minimal value of misfit for all parameter tested? This rises the question of the validity of the final result.\nAns: We understand your concern. During this revsion, we conduct more simulations with smaller porosities of the two reservoir layers during the first round of model calibration (see the following table). You can see that the RMS misfit to surface displacement observations is 4.01 cm with the porosity of the lower reservoir equal to 14%\nTable S4. Misfits of model predictions with different reservoir porosities to displacement and temperature observations. Orange and green values correspond to the first and second rounds of model calibration, respectively.\nWe realized that there would be some potential uncertainties of the optimal solution associated with the choose of test parameter values. We thus conduct two rounds of model calibration to minimize potential uncertainties.\nAnd the form\n• There are many problems in the figures and their captions. Some figures are not clearly explained enough. Every single mention in the figure should be explained in the caption. Many are missing. Some figures are not showing what is said in the caption.\nAns: During the revision, we reploted the figures in the main text and rewrote the corresponding captions.\n• When comparing figures of the same area for observations and interpolation or inversion, I would strongly suggest to use the same axis. Otherwise, it is challenging to find correspondence.\nAns: Thank you for your suggestion. During the process of reploting the displacement maps, we used the same spatial range.\nAs a whole, this is an interesting case study, however the main message is not brought clearly enough, because the demonstration suffers from lack of clarity and some clear evidence that you reached the best model ultimately. In addition, whether the conclusion found can be generalized to other cases and how is not brought. I suggest that the above mentioned issues and the specific comments below need to be solved before any decision can be taken on the manuscript.\nAns: We are very appreciative for your careful review. As we mentioned before, the MS has been revised based on the constructive and comprehensive comments from you and the other two reviewers to address these issues of clarity and supporting evidence. We believe the revised manuscript allows the readers to clearly understand the motivation, the data, the method, and the main results with clear supporting evidence and robust conclusions.\nSpecific comments.\nLine 49-52. The permeability is found along faults, which are described. However, it is not clear enough which is which with respect to the figure 1. In the text you do not mention the feeder fault, is it one on them?\nAns: The feeder fault is the normal-slip fault. We introduced it in Line 53 of the original MS. During this revision, we add the symbols for strike-slip faults and the normal-slip feeder fault into Figure 1a.\nLine 49-51. You describe the matrix permeability from 0.55 km until ~2 km (?) but you do not indicate what is above or below. This is lacking also here.\nAns: The matrix permeability structure consists of a shallow sedimentary reservoir and three faults. Above and below the reservoir are capping clays and basement, respectively (Response Figure 5).\nResponse Figure 5: Temporal evolution of pressure, temperature and their effects at four probe points. Simplified geologic and hydrothermal profile of the HGF (modified from James et al., 1987).\nLine 51. “The normal fault” is very vague. Which one is it? How do you know it is a normal fault? What does it bring to the study?\nAns: The normal fault is the feeder fault. We label it in Figure 1a. James et al. (1987) presented this fault in their figure 1.\nResponse Figure 6\nFigure 1a. Please indicate what “F.” means. In the caption, HGF meaning is missing. It should be given at the first occurrence in figures. HGF is messing in caption 1a. Same applies for SCSN: what is it? What are yellow circles? What are the beach balls? The plate boundary is certainly not so linear, is it? It would be more realistic to represent it with its real path, as feeder and subsidiary faults.\nAns: Thank you for your careful review. “F.” indicates “Fault”. HGF is Heber geothermal field; we defined it in the main text. SCSN is Southern California Earthquake Data Center; we add the full name into the caption of Figure 1a. The yellow circles are SCSN catalog. Beach balls are focal mechanisms of five M>2 earthquakes. For the plate boundary (Dixieland) fault, there is no relative information on its surface trace in our study region. So we simplify it as a line based on its northwest extension and Figure 1 in James et al. (1987).\nFigure 1d. The red line is usually called “cumulative seismicity” However, it is also often shown with the energy released. This information may not be available, as seismic events in geothermal system are small and estimating their energy sometimes challenging, but quantification would be nice here; it would be nice to indicate what is the variability of energies of the earthquakes shown. In the leveling time series, why the error bar increases with time and there is no error on the first point? The cyan color is hard to follow. May be plot it on top of the leveling points?\nAns: Thank you for your suggestion. We revised the label of left axis into “cumulative seismicity”. We agree that having some estimates of seismic energy release would be interesting, but as you mentioned this topic is deserving of its own study, especially for geothermal fields. We are planning to write a follow up study focusing on the induced seismicity, and we prefer to leave the energy calculations in that study to keep the scope of the prevent this study tractable to understand.\nThe “error bar” are derived from the flunctuation of leveling time series (see Response Figure 7). From the following figure, you can see that, for the first point, the displacement of each benchmark is equal to 0. So the “error” is 0. The “error bar” only reflect the flunctuation of leveling time series rather than the observation error. The observation errors of leveling data range from 0.2 cm to 0.5 cm.\nResponse Figure 7: Spatial distribution of leveling stations and time series of observed vertical displacement.\nThe cyan color is changed to blue color. We plot the curve on top of the leveling points.\nLine 70-71. I do not get it. The resource of a geothermal field cannot be designed… The resource is what we try to access, so we design the capacity of the power plant to match the resource potential. This sentence needs to be rephrasing. Or am I missing something here?\nAns: Thank you for your careful review. We revised the sentence to be “The HGF resource was initially estimated to support a capacity of up to 500 MW of electrical power generation.”\nLine 72. There is a tendency now to try not to “exploit” natural resources, but rather harness. Could you check in all the manuscript?\nAns: Thank you for your suggestion. We deleted the word “exploit” in the revised MS. We also checked the MS and confirm that the word “exploit” is not used in the other places.\nLine 85-86. Possibly you need a reference here.\nAns: Thank you for your suggestion. We rewrote the paragraph during the revesion.\nLine 86. It would be important to give information on the seismic network (number and type of stations) deployed since the beginning of the exploration and over the years. This would give us an impression of if there is more earthquake due to the better recording capability of indeed because of more earthquakes.\nAns: This is a good point. The stations used are from the Southern California Seismic Network (SCSN; see https://www.fdsn.org/networks/detail/CI/), and for this revision, we compared the network catalog with a matched-filter catalog based on cross-correlation of template events with continuous data. The methodology, which we now describe in the Methods section, adjusts for network coverage as it scans through continuous data. The two catalogs show a similarly increasing trend of seismicity, although the MF catalog provides better temporal resolution, which confirms the trends in seismicity are real rather than artifacts of network coverage. In addition, the template matched catalog shows that rapid increase of earthquake within the template matched catalog occurred at the time with prominent changes of injection and production rates.\nLine 89. I would suggest to change “indicating” with “suggesting”.\nAns: OK. Thank you!\nLine 94-96. Actually in this paragraph, there is nothing new, this is known. Good to know for HGF, but so what?\nAns: During this revision, we deleted the paragraph and put the relative contents on surface displacement and seismicity into the first two subsections of “Results”, respectively.\nFrom the observations of surface displacement and seismicity, we found an interesting phenomenon that the temporal growth of seismicity rates is similar to the overall trend of ground deformation shown by the normalized leveling time-series, suggesting that physical mechanisms of surface deformation also control the occurrence of earthquakes at the HGF.\nLine 114, 139, 235, and other locations possibly, e.g., supplementary figure captions. A space missing between time and series.\nAns: Thank you for your careful review! We changed the word “timeseries” into “time series” in the revised MS.\nFigure 2. The reason(s) of the separation of west, central and east regions is far for clear. From what is explained, there are no structural structures that justify this separation. Why not north, central and south?\nAns: As we introduce previously, the study area is divided into the three zones based on two main reasons: (1) the west, central and east regions of the Heber geothermal fields experienced different displacement trends, and (2) there are only injection operations in the west and east zones, but both injection and extration operations occurred in the central zones.\nFigure 2e. I do not understand why you need to normalize the vertical displacements. Note the term “vertical deformation” is inappropriate, as what is observed with both leveling and InSAR are displacements, and actually you indicate also “per year”, which makes it technically velocities.\nAns: The displacement time series of the leveling stations are normalized for two purposes primarily related to visualization of a complex history of displacements in space and time: (1) we needed a way to represent the displacement trends of the west, central and east regions of the Heber geothermal field, which show different displacement trends; and (2) to represent temporal variation of displacement of the whole region. By doing this we can: (1) compare vertical displacement trends of the west, central and east regions with the respective net injection volumes of geothermal fluid of every month; and (2) compare the displacement trend of the whole region with temoporal changes of observed seismicity at the Heber geothermal field. Again, we stress that our misfit calculations are based on the true displacements rather than the normalizaed ones.\nIn addition, to compare the observation results of leveling and InSAR, the InSAR LOS displacement rate maps of ascending and descending tracks are decomposed into vertical and east components. The leveling-measured displacement rate maps are calcuated based on the cumulative displacements of benchmarks and the corresponding observation durations.\nBased on your suggestions, we replaced the word “vertical deformation” is replaced with “vertical displacement”.\nLine 148. “The same scenario”. This is confusing. Which same scenario? In addition, if this is the same scenario between two areas separated the artificial boundaries (west and east), it means that there may be no need to separate them…\nAns: “The same scenario” indicates that ground was susiding with fluid injection. We revised the words into “The same phenomenon”. Although this phenomenon has been observed in both west and east zones, it ocurred in different time periods: 1994-2005 (west) and 2005 to now (east). More importantly, the deformation anomalies are caused by different mechansims: the siphonic effect of the HGF center (west) and thermoelastic contraction (east).\nThe reasons for dividing the study area into 3 zones refer to our above responses to your comment on Figure 2e.\nFigure 3. The figure shows a series of temperature points. However, there is not clarity... for instance in Fig. 3a, I do not see 11 injection wells values. Are there all mixed in the figure? Are they averaged? Or only one is shown? Similar issues for the other sub-figures.\nAns: Thank you for your comments. The temperature observations of injection and extraction wells are mixed in Figure 3. We clarify it in the caption.\nLines 169-171. These interesting observations are hard to believe, but is I generally trust data. This paragraph is just describing data, but do not suggest any potential reasons for the temperature changes. Possibly this is not the place, however, after displacement observations, you suggest a reason for the observed changes. Therefore I would expect the same for temperature observations.\nAns: In our assessment, the temperature changes of the 11 deviated wells are probably associated with geothermal operations. Consequently, for this revision, we put the contents on temperature records into the section of Introduction as a background of the HGF.\nAt line 178-180, there seems to be indeed a reason proposed: convection. Is that the case? If yes, I would make is even clearer.\nAns: Yes, you are right. During this revision, we explain the reason for the temperature changes in the subsection of Results (Reservoir pressure and temperature evolution).\nLine 183. Indeed an active fault may move. It would be also interesting to relate the natural seismic activity along the entire Diexiland fault with time, to make the potential link between the local seismicity and the more regional one.\nAns: It is indeed a very interesting topic, but our MS mainly focuses on the spatio-temporal evolution of deformation mechanisms. The topic requires a detailed investigation including hypocenter relocation of seismicity and analysis of its spatiotemporal migration, which is beyond the scope of this study. We will prepare another paper to investigate the influence of geothermal operations on the stress accumulation and release of active faults.\nLine 186. Two major issues here:\n1. I would be very cautious to invert local data for global interpretation.\nAns: As you know, the coverage of InSAR data is very large with respect to the signals we are modeling. Our purpose of inverting the leveling and InSAR observations is to convince the reader that the local deformation signals are not associated with potential fault slip, rather than for global interpretation. We have performed additional tests for this revision.\n2. Inversion of leveling data does not account for horizontal displacements. However, due to the local geodynamics setting strike-slip faults (e.g., the Dixieland fault?) mostly horizontal displacement would be expected.\nAns: Yes, we agree with you. For that reason the vertical and horizontal deformation observations are modelled with normal-slip feeder fault and strike-slip Dixieland fault, respectively.\nLine 197. How did you measure horizontal displacements? The explanation resides in the InSAR/leveling integration. At places, this is not clear enough on the temporal evolution of the measurements and how you inverted each period of time. How to combine both of them is sometimes vague.\nAns: Horizontal displacements are obtained through standard methods of decomposing line-of-sight InSAR maps given the presence of both ascending and descending track coverage in the same time; this gives vertical and east horizontal components. We are fortunate to have good ascending and descending coverage for this study, and the details of the processing are available in Eneva et al. (2019; GRC).\nFor this revisison, we replotted a figure to show the temporal evolution of vertical surface displacements observed by InSAR and leveling in the main text. For the horizontal deformation, we present its temporal evolution in Figure S3.\nTo convince the reader that the local deformation signals are not associated with potential fault slip, we conduct three scenarios of fault slip modeling to available deformation observations. First, the long-term vertical displacement measured by leveling from 1994 to 2004 is inverted with slip on the normal-slip feeder fault. Second, the vertical transient displacement from 2006 to 2010, obtained through removing the long-term displacement trend, is also inverted with the feeder fault. Third, the horizontal displacement measured by InSAR from 2006 to 2010 is modelled with the strike-slip plate boundary (Dixieland) fault.\nLine 201. Minor point: I would suggest you indicate the code used here also (although is it given in sup info)\nAns: Our numerical simulations are implemented with the software COMSOL Multiphysics (version 6.0). The model files are available from Zenodo (https://zenodo.org/records/10460410). We explain it in Code Availability of the revised MS based on the instructions of Nature Communications.\nLine 223. “strike-slip” it is not clear which fault is which mechanism. You could introduce them better in figure 1.\nAns: The strike-slip faults include the plate boundary (Dixieland) fault and the subsidiary fault. We added two strike-slip symbols in Figure 1a.\nLine 228. Why inverting for layer parameter is not feasible? Does it not depend on which ones? You need more explanations here.\nAns: The major reason is that physical processes of geothermal operations are intercoupling, dynamic and nonlinear. Further, with such a broad parameter space, achieving good results would require massive supercomputing power, which we do not have resources to use.\nFigure 4. What are P1, P2, P8 etc. Please report extensively all elements within the figure caption.\nAns: Thank you for your suggestion. We add specifications to these probe points into the figure caption.\nFigure 5. This figure is very unclear, mainly because the sub-figure are too small, and axis are really tiny. In addition, it is not possible to understand what are the different curves of different colors.\nAns: Our goal for this figure is to show the distinct features of pressure effects and thermal effect. During this revision, we realized that it is unneccesary to show the simulated displacements and pressure/temperature changes of all probe points in the main text, and deleted the figure. Instead, we plot a new figure to show spatiotemporal evolution of thermo-poro-elastic effects. In the new figure, we only present the pressure effects and thermal effect at 4 probe points with good representativeness. The results of the other 4 probe points are presented in the supplementary materials.\nLine 358. Figure 15 does not exist.\nAns: Sorry! It should be Figure S15.\nLine 491. The initials of Jean-Philippe Avouac are J. P.\nAns: Thank you! We corrected the mistake in the revised MS.\nSupplementary Materials:\nPlease check the use of “deformation”.\nAns: Thank you! We checked the words “deformation” and “timeseries” during this revision.\nIn § Fault slip inversion: I have an issue on what you want to invert for and the procedure. You remove the geothermal signal. Why and how do you do this? Then you subtract the pre-2005 trends. I really do not understand why. This may be the signal you want to invert for… You may need a reference here to make sure they are the whole area signal. After removal, what is left??\nAns: We are sorry for ambigious introduction. During this revison, we rewrote the text on details of slip inversions. Here we would like to clarify some important points. First, the deformation signals are only observed at the HGF. Second, the surface displacement at the HGF consists of two parts: long-term component initiated before 2005 and transient component after 2005 (Fig. 3). Therefore, there are two end-member possibilities: (1) the long-term displacement associated with fault slip, and (2) the transient displacement associated with fault slip.\nWe thus conduct three scenarios of fault slip modeling to available deformation observations. First, the long-term vertical displacement measured by leveling from 1994 to 2004 is inverted with slip on the normal-slip feeder fault. Second, the vertical transient displacement from 2006 to 2010, obtained through removing the long-term displacement trend, is also inverted with the feeder fault. Third, the horizontal displacement measured by InSAR from 2006 to 2010 is modelled with the strike-slip plate boundary (Dixieland) fault.\nin the second step and (1), you use “corrected vertical-component data”. Here 2 things: 1. what are data corrected for? The previous trends?\nAns: Yes, the data are corrected with pre-2005 (long-term) displacement trends to isolate the transient displacement signals in the vertical direction.\n2. why do you invert only the vertical data? Inversion procedure with Okada model is clearly possible with all components, so I do not understand why only vertical. In addition, as the plate boundary is a strike-slip process, there may be little vertical signal but large horizontal one. So this inversion results may be completely wrong, as you invert the wrong component in which you removed some potentially valid signal… You need to review here seriously the approach.\nAns: You are correct that the Okada model provides 3D displacments; but, since our data come from different measurement techniques with different capabilities, this is a necessary limitation in some of the tests depending on the time-period considered. Prior to 2005 we have only vertical displacements, for example. To be very clear, the vertical data is inverted with the normal-slip feeder fault. In addition, we also inverted the horizontal InSAR data with the strike-slip Dixieland fault.\nIn the thermo-poro-electicity formulation, in the last paragraph you indicate primary effects and secondary effects. You could list them, and justify how to classify them as primary or secondary.\nAns: Thank you for your suggestions. The three effects include: (1) influence of rock deformation on rock porosity and permeability; (2) thermal influence on fluid density and viscosity; (3) frictional heating due to rock deformation. Palciauskas and Domenico (1982) and McTigue (1986) have shown they are secondary effects.\nResponse Figure 8: Flowchart of two-way thermo-hydro-mechanical coupling. The arrows labeled with blue words indicate the effects included in our simulations. This can be found in the supplemental material\nFigure caption of S1, S3. Check the use of “deformation” here.\nAns: Thank you. Fixed.\nFigure caption of S1. The normalization method for uplift is unclear. Why do you need to normalize? Is the normalized value inverted for??\nAns: As we’ve discussed already, we only use normalized displacements for visualization and general analyses; we invert the true displacements. Further, to characterize vertical deformation trends of the HGF, we divide the study area into three regions (west, central and east) to analyze the leveling observations. The displacement time series of the 103 stations are first normalized by their root-mean-square values with the following expression,\n, 2 𝑁𝑁 𝑢𝑢(𝑖𝑖) where and are the measured a𝑛𝑛n(d𝑖𝑖 )no=rm𝑢𝑢a(𝑖𝑖li)zed ∑de𝑖𝑖=fo1rm𝑁𝑁−a1tion, respectively, at each time node; and N is the total number of surveys. We then average the normalized timeseries of the stations within each 𝑢𝑢(𝑖𝑖) 𝑛𝑛(𝑖𝑖) sub-region to represent its deformation trend. For the whole study area with both surface uplift and subsidence, absolute values of all normalized time series are averaged to represent the deformation trend.\nFigure S2. In the caption, please add the reference for the kriging here too. In the figure, you should get the same axis limits for both the data and the interpolated info.\nAns: Thank you for your suggestions. We replotted Figure S2 during this revision and added a reference for the Kriging method.\nFigure S5. The sub figure are too small, making this figure useless. You also need to tell readers when… time is not indicated or clear enough.\nAns: Thank you for these comments. We replotted production rates and temperature variation of all wells of the HGF with 3 figures.\nFigure S6. There are several place where the term “excess” is used for vertical displacement. I am not convince with excess it is. I did not get where this excess is coming from. Need more explanation here and before possibly.\nAns: The “excess” vertical-component displacement of 2006–2010 is obtained from removing the pre-2005 (long-term) trend. We explained it in the main text and methods for this revision.\nFigure S7. Where is A-33? Show it on the figures.\nAns: The station A-33 is shown in Figure S1.\nFigure S16, S17, S18. The number of sub figure is large. You could gain place by removing the axis for the inner figure, as they are the same for all. This would make the figure much easier to read. In addition, on the all 3 pages for each figure, I count 96 (4*8*3) points. Not clear where the 3 otehrs are…\nAns: Thank you for your suggestions and careful review. There are indeed 3 points that we did not present their simulation results. For this revision, we replotted Figures S16-S18."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have done an excellent work revising the manuscript and all comments are addressed satisfactorily.\nI do not have any further comments and recommend publishing the manuscript.", "reviewer_2": "This study presents an interesting case of the thermo-poro-mechanical response of the subsurface to fluid injection and production in a geothermal field. I have several concerns that are described in an attached file becuase I include an equation for clarity of my explanations.\n[Editorial Note: This PDF is displayed across the next three pages]\nReviewer #2 (Remarks on code availability):\nThe readme file just states that they provide a model that can be used with COMSOL. COMSOL is a comercial software whose license is very expensive, so few people will be able to open it.\nReviewer #2 (Remarks to the Author):\nReview of the manuscript “Relatively stable pressure effects and time-increasing thermal contraction control Heber geothermal field deformation”, by G. Jiang et al.\nThe authors have made significant changes to the manuscript in response to the reviewers’ comments. However, the changes have only partially addressed the comments and new comments have arisen as a result of the modifications to the manuscript.\n1) The authors justify the use of constant water properties because “the functions of fluid properties changing with pressure and temperature are not available”. These functions are well known and can be easily implemented in COMSOL. It is weird that the misfits shown in Table S16 as a function of water properties do not change for viscosity. The minimum and maximum tested viscosity values change by a factor of 2. Changing viscosity has a similar effect to changing permeability (in the same proportion). The misfit changes when changing reservoir permeability from 10 to 20 mD. Thus, changing viscosity by a factor of 2 should also have an effect on the misfit.\n2) The upwards heat flow through the Feeder Fault is a variable density-driven flow in which the less dense hot water from depth migrates upwards through the Feeder Fault because of buoyancy. If water density is assumed constant, how can upwards heat flow be reliably modeled?\n3) Figure 1 would benefit from including the direction of the maximum and minimum horizontal stresses to better understand the fault stability analysis. Apart from this, the arrows indicating the slip direction in the strike-slip faults indicate a left-lateral movement of the fault, which is the contrary to the regional movement – the San Andreas fault is a right-lateral fault. To induce a left-lateral movement of the Dixieland Fault, the maximum horizontal stress should have an azimuth around 90°, which would be perpendicular to the Feeder Fault, making it very stable (because it is a normal fault). Please, double check the orientation of the arrows indicating the movement of the strike-slip faults because they seem to have the wrong orientation. The magnitudes of the stresses used in the model should also be provided in the manuscript.\n4) Figure 4c: after 2005, pore pressure effects have a good fitting with the leveling-observed displacement trend, so the model seems to be overestimating thermal effects. Why is there this discrepancy in the ground displacement?\n5) Also, the observed trend in the east region, which presents an uplift followed by subsidence (Figure 3a), is not reproduced by the numerical model, which yields a continuous subsidence from the beginning of operations (Figure 8 and Supplementary Fig. 24).\n6) Transient and rapid increases in seismicity rate are likely linked to fault reactivation and, thus, fault slip. As indicated in the text, “The most notable one started in 2005”, coinciding with the time of the reversals in the ground deformation trend. The thermo-hydro-mechanical elastic model shows that the increase in the injection rates in 2005 has an effect on the change in the trends of ground deformation, but fault reactivation and its consequent slip are likely playing a non-negligible role as well. Furthermore, the reactivated fault may continue slipping aseismically once it is reactivated. Note that surface displacement is a reflection of all the deformation that occurs below the ground and fault slip also affects ground deformation.\n7) Fault slip modeling: it is argued that the necessary slip rate on the Feeder Fault to explain the observed ground deformation would be in the order of 4 to 10 cm/yr. This is equivalent to 0.001 to 0.003 I.tm/s – a plausible value for aseismic slip of the fault. It is argued that the slip rate is larger than the regional slip, but it should be borne in mind that injection/extraction operations induce ground deformation and may destabilize faults (see the seismicity around the faults, indicating that they are slipping). An aseismic slip rate of 0.008 I.tm/s has been inferred at Brawley, California, as a result of geothermal operations (Wei et al., 2015). This slip rate is in the same order of magnitude of the one estimated in the Heber geothermal field. Thus, the hypothesis of fault slip as responsible, at least partially, of the observed ground deformation cannot be disregarded.\n8) Another indication that faults reactivate is the observation of abrupt changes in temperature, ∼30 °C drop after 2000. Fault slip usually enhances permeability because of dilatancy. Thus, fault reactivation and the consequent permeability enhancement could lead to a sudden increase in flow rate, affecting fluid temperature.\n9) Specific heat of rock: typical values range from 500 to 600 J/(kg·°C), lower than the assumed value by roughly a factor of 2.\n10) Figures 7a and 7b: it would be helpful to add dimensions to the axes.\n11) Figure 7: we do the abrupt pressure drops, e.g., around 1994 in P4, occur in the model?\n12) Lines 362-365: it is said that “heat transfer within the reservoir is mainly through the convection of fluid flow”. It should be indicated that it is “forced convection” because natural convection cannot be simulated in this model because fluid properties are considered constant.\n13) Quantification of thermo-poro-elastic effects: I disagree with the way the individual effects are calculated. There are three models: (I) Biot’s coefficients equal to 0; (II) thermal expansion coefficients equal to 0; and (III) both Biot’s and thermal expansion coefficients equal to 0. Strain are proportional to effective stress changes, Δ ∆ ∆ ′, which are a function of total stress changes, , pore pressure changes, , and temperature changes, , as Δ ′ = Δ – Δ + Δ , where is Biot’s coefficient, is the thermal expansion coefficient and is the bulk modulus. Thus, model III gives the poroelastic effect, i.e., the deformation as a result of the stress changes caused by pore pressure changes (and not the effect of pore pressure changes as indicated in lines 375-376). The effect of pore pressure changes can be obtained by subtracting to the calibrated model the deformation of model II and model III. And thermal effects can be obtained by subtracting to the calibrated model the deformation of model I and model III.\n14) Discussion, lines 442-444: it is claimed that an improvement in the understanding of the spatial and temporal evolution of the strength of pore pressure and thermal effects is achieved in this study. However, the manuscript does not provide such better understanding. It is shown that there is a competition between pore pressure and thermal effects, with pore pressure presenting a rapid response and thermal effects becoming dominant in the long term, but this was already known. And it is shown that thermal contraction affects not only the cooled region, but also further away, which was already known. There is no quantification or clue how these processes can be quantified or estimated in other geothermal sites. In general, the observed coupled thermo-hydro-mechanical processes are already known. The novelty and interest of this study resides in the application to a geothermal site with data covering decades.\n15) The two-way coupling approach is partially neglected. In particular, rock deformation is assumed to have negligible effects on temperature and pressure fields. While this may be true in porous media, it may not be like that in fractured media, like the Heber field, which has several faults. Deformation may significantly enhance fracture/fault permeability because it scales with the cube of the aperture, affecting both the pressure and temperature field.\n16) Calibration: it is a bit surprising that the initial guess (with few exceptions of slight changes in model properties) of model properties gives the best fitting. Calibrating a 3D THM model is a challenging task.\nMinor comments:\n1) Line 341: you refer to “injection volume” but the units are given in kg\n2) Line 561: the assumed values of fault properties are justified as “widely-used values”, but the two provided references are papers of the first author of this study. References to papers by scientists that are not authoring this manuscript should be provided\n3) Line 599: the mechanical boundary condition of “fixed injection and production wells” is not very clear. What is fixed? According to the figures in the Results, the surface displacement at the wells is not zero.\n4) The units on both sides of Equation (3) are different: units are 1/s in the left-hand side, while they are m3/s in the right-hand side.\n5) Line 700: please, provide references to support “typical values of faults”. Faults can take wide range of values, depending on the brittleness of the rock and accumulated slip.", "reviewer_3": "Review of revised manuscript „Relatively stable pressure effects and time increasing thermal contraction control Heber geothermal deformation\nBy G. Jiang et al.\nSubmitted to Nature Communications\nGeneral comments\nThis version of the manuscript is much easier and smoother to read than the first one.\nThe arguments are logically exposed, and the demonstration is much clearer.\nHowever, the initial part of the argument is still a bit weak, especially concerning some terms (displacement, levelling, ...) actually many things around Figure 3, which prevent full trust in the results. Once this step passed, and admitting the observations, the demonstration is better conducted. Therefore, I strongly suggest to revise once more the manuscript to make it smooth.\nSpecific comments\nLine 18. The term “siphonic effect” is a bit awkward here, and should be explained.\nLine 19 or 20. “by cooling“ could be added after “Thermal contraction “. It seems obvious, but this is clearer by adding the reason.\nFigure 3 is still confusing. Fig 3a. Uz is not defined. Is it position or velocity. I do not understand why you normalize (line 104). Real values of displacement should be giving more information for possible inversion of the source location. In the west region, you show injection rate in blue. However, there are no wells in the west region, according to map Fg 3b. Where are the wells or where does this curve come from? I would add an indication in the figure 3a west, central and east of the three periods indicated below (e.g., using P1, P2 and P3, which would correspond to the periods indicated in fig 3b. I do not understand the grey zone, as it seems there is also a decade-long displacement transient from 1995 until 2005 (e.g., fig. 3a.). The term “transient” should therefore been defined. In Fig. 3b the caption indicates maps obtained from levelling and InSAR, but the top of each subfigure indicates only levelling. I understand the space is rather small to indicate both, but this is confusing.\nLine 123. There is confusing reference to displacement, although you refer mostly to rate of displacement, which is velocity (for instance in this line 123, where 1.5 cm/year is indeed a velocity).\nLine 124. There is inconsistency here with the main conclusion of the paper. You have both subsidence AND uplift with an increase of injection rate. Therefore, this observation contradicts line 15 of the abstracts. This needs further explanation and clarification.\nLine 128-129. I do not see where in Fig 3b, levelling-measured subsidence is large than InSAR results. If this is indicated in supplementary info, you should indicate this here. Now on the result itself, I do not understand why the levelling and InSAR height changes are different. You should give a reason for this difference.\nLine 135. I do not follow exactly what you mean. “... vertical displacement ... declining...” This is not very clear.\nLine 138-139. I do not see this in figure 3. The west region indeed subsides before 2005, but I do see a reverse, so an increase only after 2005. However, before 2005, I do no see an increase of the injection rate. Then this sentence does not make sense. In addition, abnormal is a weird word for scientific expression. I would replace with “unexpected” or something like this. However, I do not see that in figure 3.\nLine 139-140. The sentence “the same displacement anomaly...” is unclear. Which anomaly? And it is in another region and another time... so quite confusing!\nFigure 4 and line 168. Is the “thermal contraction” the same as thermoelastic effect (in the figure)?\nLine 171. It would be nice to precise if there any relationship between the observed horizontal displacement rate (2 cm/year) and the possibility you invoke here.\nLine 320. Is it necessary to mention “Strong”?\nLine 344 and 346. What is the unit “kg/mo”?"}, "author_response": "Response to Reviewer #2\nThe authors have made significant changes to the manuscript in response to the reviewers’ comments. However, the changes have only partially addressed the comments and new comments have arisen as a result of the modifications to the manuscript.\nAns: We are really sorry that the revised MS did not completely resolve your concerns. During this revision, we have tried to address your comments as thoroughly as possible. Following are our point-by-point responses. We sincerely appreciate your careful review and critical comments.\n1-1. The authors justify the use of constant water properties because “the functions of fluid properties changing with pressure and temperature are not available”. These functions are well known and can be easily implemented in COMSOL.\nAns: We only found the equation of water density changing with pressure and temperature (e.g., Wagner and Pruß, 2002),\n, 𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 where is the reference densit𝜌𝜌y =wi𝜌𝜌th0 +the𝑑𝑑 𝑑𝑑te(m𝑃𝑃p−er𝑃𝑃a0tu)r+e o𝑑𝑑f𝑑𝑑 (𝑇𝑇 a−nd𝑇𝑇 0t)he pressure of . However, this equatio𝜌𝜌n0 is unable to be implemented in COMSOL due to that𝑇𝑇 t0he two items ( and 𝑃𝑃)0 are not constants.\n𝑑𝑑𝑑𝑑 𝑑𝑑𝑑𝑑 For other properties of viscosity, thermal conductivity, and heat capac𝑑𝑑i𝑑𝑑ty, w𝑑𝑑e𝑑𝑑 did not find the corresponding functions.\nDuring this revision, we learned about that the material library of COMSOL provides the changing water properties with pressure and temperature, which is obtained through interpolation. We plot figures (attached below) to show their changes with the pressure (350 to 460 K) and temperature (14 to 17 MPa) variations of the Heber geothermal reservoir. From the following figures, you can see that:\n(1) the water properties are characterized with obvious changes within the given temperature range;\n(2) the water properties are characterized with little changes within the given pressure range;\n(3) the density, thermal conductivity, and heat capacity change slightly with the pressure and temperature variation;\n(4) the viscosity has obvious changes with temperature but sharp temperature decrease is only localized near injection wells.\nThe material library from COMSOL was too expensive to purchase, but to resolve your concern, we used a trial version and re-ran our model with the changing water properties. The updated simulation results from 0 d to 9400 d show that the misfits to leveling, InSAR and temperature observations are 2.86 cm, 2.67 cm and 21.6°C, respectively. In contrast, the misfits of the calibrated model are 2.39 cm, 2.68 cm and 20.3°C, respectively. You can see that the model misfits become larger with the variable water properties. We tested the influence of different water properties on simulation results with results in Supplementary Table S16.\nReference: Wagner, Wolfgang, and Andreas Pruß. “The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use.” Journal of Physical and Chemical Reference Data 31, no. 2 (2002): 387-535.\n1-2. It is weird that the misfits shown in Table S16 as a function of water properties do not change for viscosity. The minimum and maximum tested viscosity values change by a factor of 2. Changing viscosity has a similar effect to changing permeability (in the same proportion). The misfit changes when changing reservoir permeability from 10 to 20 mD. Thus, changing viscosity by a factor of 2 should also have an effect on the misfit.\nAns: During this revision, we checked the model settings and found that the value of water viscosity in our model was fixed to be the initial value of 2.5E-4 Pa*s. Thank you very much for careful review! We corrected this error and re-ran our simulations with the test water properties. Table S16 (we attach it below) lists the updated results for the model misfits, which are still smallest with the intial values of water properties.\n2. The upwards heat flow through the Feeder Fault is a variable density-driven flow in which the less dense hot water from depth migrates upwards through the Feeder Fault because of buoyancy. If water density is assumed constant, how can upwards heat flow be reliably modeled?\nAns: According to Lippman and Bodvarsson’ studies (attached below), exploration wells verified the presence of significant hydrothermal flow through the feeder fault zone, with temperatures of nearly 245°C and rates around 14.6 kg/s (see the following conceptual figure from Lippman and Bodvarsson, 1985). To identify the replenishment mechanism of heat flow at the HGF, we conducted additional simulations under two scenarios: (I) one upwelling channel of the feeder fault with different rates of heat flux set as a boundary condition at the base of the fault; (II) two channels of the feeder and subsidiary strike-slip faults with different flow rates.\n1. Lippman, M.J., and Bodvarsson, G.S. (1983). A modeling study of the natural state of the Heber Geothermal Field, California. Geothermal Resources Council Trans., Vol. 7, p. 441-447.\n2. Lippman, M.J., and Bodvarsson, G.S. (1983). The generating capacity of the Heber Geothermal Field, California. Proceedings Ninth Workshop Geothemal Reservoir Engineering, p. 157-166.\n3. Lippman, M.J., and Bodvarsson, G.S. (1985). The Heber geothermal field, California: natural state and exploitation modeling studies. Journal of Geophysical Research, 90 (B1), 745-758, doi:10.1029/JB090iB01p00745.\n3-1. Figure 1 would benefit from including the direction of the maximum and minimum horizontal stresses to better understand the fault stability analysis.\nAns: Thank you very much for your suggestion! We add the arrows into Figure 1a (we attached it below) to show the orientation of regional background principal stress.\nThe stress orientation is from the SCEC Community Stress Model (https://www.scec.org/research/csm), which compiles multiple sources of information.\n3-2. Apart from this, the arrows indicating the slip direction in the strike-slip faults indicate a left-lateral movement of the fault, which is the contrary to the regional movement – the San Andreas fault is a right-lateral fault. To induce a left-lateral movement of the Dixieland Fault, the maximum horizontal stress should have an azimuth around 90°, which would be perpendicular to the Feeder Fault, making it very stable (because it is a normal fault). Please, double check the orientation of the arrows indicating the movement of the strike-slip faults because they seem to have the wrong orientation.\nAns: Thank you very much for your careful review! Yes, we plotted the arrows with wrong orientations. The two strike-slip faults in our study region are indeed right-lateral. We corrected the error in Figure 1a.\n3-3. The magnitudes of the stresses used in the model should also be provided in the manuscript.\nAns: In our model, we consider the field to be in stress equilibrium at time zero, and did not consider the regional background stress. Alternatively, the mechanical boundary conditions include: (1) free upper surface, (2) roller boundaries for the bottom and side surfaces, and (3) fixed open-hole sections of injection and production wells. In addition, we used the simulated stress at t=0 to make the model achieve initial stress balance, and therefore the model simulates the impacts of geothermal operations on the state of stress and fluid pressure and temperature.\n4. Figure 4c: after 2005, pore pressure effects have a good fitting with the leveling-observed displacement trend, so the model seems to be overestimating thermal effects. Why is there this discrepancy in the ground displacement?\nAns: Leveling-observed surface displacement trend and temporal growth of seismicity within the HGF is highly correlated with the normalized vertical displacements linked to thermal contraction (Fig. 4c), with the correlation coefficients up to 0.85 and 0.97, respectively. In contrast, the coefficients of the leveling displacements and seismicity growth relative to the normalized displacements associated with the pressure effects are only 0.09 and 0.23, respectively.\nDuring this revision, we further calculate the correlation coefficients for the time period after 2005. The coefficient between the leveling observations after 2005 and the constribution of pore pressure effects is -0.02. In contrast, the coefficient between the leveling observations after 2005 and the constribution of thermal effect is 0.89.\n5. Also, the observed trend in the east region, which presents an uplift followed by subsidence (Figure 3a), is not reproduced by the numerical model, which yields a continuous subsidence from the beginning of operations (Figure 8 and Supplementary Fig. 24).\nAns: Yes, our model cannot simulate the ground uplift of the east region before 2005; but our model does not overfit the subsidence after 2005 (see the following two figures: one shows the locations of leveling benchmarks; another shows the observed and simulated leveling time series). Therefore, we do not think that our model overestimates the thermal effect. It is hard for us to present potential reasons for the bad fit before 2005 aside from the obvious sources of uncertainty related to reservoir structure. However, based on our modeling we gained the following insights: (1) the observed surface displacements far outside the HGF scope are harder to be fit than the observations within the scope with geothermal operations; and (2) the model paramaters calibrated with available observations represent the rock properties in average. For the peripheral region of the HGF without geothermal operations, the displacement observations are prone to be disturbed by other physical processes like natural recharge of regional groundwater. We note this directly in the text.\n6. Transient and rapid increases in seismicity rate are likely linked to fault reactivation and, thus, fault slip. As indicated in the text, “The most notable one started in 2005”, coinciding with the time of the reversals in the ground deformation trend. The thermohydro-mechanical elastic model shows that the increase in the injection rates in 2005 has an effect on the change in the trends of ground deformation, but fault reactivation and its consequent slip are likely playing a non-negligible role as well. Furthermore, the reactivated fault may continue slipping aseismically once it is reactivated. Note that surface displacement is a reflection of all the deformation that occurs below the ground and fault slip also affects ground deformation.\nAns: Thank you for your comments! But there is no robust evidence to support the hypothesis that aseismic fault slip plays an important role in surface deformation, and we carefully tested this possibility in our analyses. In addition to your note, we also note that the transient displacement seen in the vertical direction is dominated by uplift, and our tests show that such uplift cannot be fit by reasonable slip on the normal-slip feeder fault (see the supplementary Fig. 7).\n7. Fault slip modeling: it is argued that the necessary slip rate on the Feeder Fault to explain the observed ground deformation would be in the order of 4 to 10 cm/yr. This is equivalent to 0.001 to 0.003 μm/s – a plausible value for aseismic slip of the fault. It is argued that the slip rate is larger than the regional slip, but it should be borne in mind that injection/extraction operations induce ground deformation and may destabilize faults (see the seismicity around the faults, indicating that they are slipping). An aseismic slip rate of 0.008 μm/s has been inferred at Brawley, California, as a result of geothermal operations (Wei et al., 2015). This slip rate is in the same order of magnitude of the one estimated in the Heber geothermal field. Thus, the hypothesis of fault slip as responsible, at least partially, of the observed ground deformation cannot be disregarded.\nAns: This is a good comment, but we do not find supporting evidence for this. Because the surface displacement at the HGF consists of two parts – long-term displacement initiated before 2005 and transient displacement after 2005 – we conducted three scenarios of fault slip modeling of the available geodetic observations. First, the long-term vertical displacement measured by leveling from 1994 to 2004 is inverted with slip on the normal-slip feeder fault. Second, the vertical transient displacement from 2006 to 2010, obtained through removing the long-term displacement trend, is also inverted with the feeder fault. Third, the horizontal displacement measured by InSAR from 2006 to 2010 is modelled with the strike-slip plate boundary (Dixieland) fault.\nOur inversions show that: (1) vertical transient displacement cannot be fit by slip on the feeder fault; (2) although the long-term vertical displacement can be fit by slip on the feeder fault, the slip rates are required to be 4~10 cm/yr. If such slip exist, the cumulative slip on the feeder fault is up to 1.2~3 m during 30 years of geothermal operation. Further, as we presented in Figure 2 (attached below), there are several wells intersecting with the feeder fault which are the main production wells for the geothermal plant. If long-term aseismic slip exists, such large cumulative slip would cause catastrophic, irreversible casing deformation of the wells intersecting with the feeder fault; but there were no related reports on such damage and production rates have remained relatively constant, implying that not major repairs to the production wells in response to fault slip has been needed. We note this in the text.\n8. Another indication that faults reactivate is the observation of abrupt changes in temperature, 30 °C drop after 2000. Fault slip usually enhances permeability because of dilatancy. Thus, fault reactivation and the consequent permeability enhancement could lead to a sudden increase in flow rate, affecting fluid temperature.\nAns: We present the temperature changes of all wells in Figure 2 above. We believe you are referring to the abrupt changes of the wells in the east region around 2006. All the abrupt changes can be predicted by our model based on injection and production alone (see Figure 3 below). In other words, we show that the temperature variations are associated with geothermal operations with no need to invoke fault permeability changes.\n9. Specific heat of rock: typical values range from 500 to 600 J/(kg·°C), lower than the assumed value by roughly a factor of 2.\nAns: According to Table 12 from Robertson E. C. (1988), the value of specific heat varies from 800 J/(kg·°C) to 1200 J/(kg·°C) with the temperature increasing from 50°C to 200°C. We note this in the text.\nRobertson E. C. (1988), Thermal properties of rocks, USGS Open-File Report 88-441.\n10. Figures 7a and 7b: it would be helpful to add dimensions to the axes.\nAns: Thank you for your suggestion! We revised Figures 7a and 7b.\n11. Figure 7: we do the abrupt pressure drops, e.g., around 1994 in P4, occur in the model?\nAns: The abrupt pressure changes are totally caused by geothermal operations. Again, we are using the injection and production values reported by the geothermal operator as inputs to the model.\n12. Lines 362-365: it is said that “heat transfer within the reservoir is mainly through the convection of fluid flow”. It should be indicated that it is “forced convection” because natural convection cannot be simulated in this model because fluid properties are considered constant.\nAns: OK, we revised the sentence. Thank you for your suggestion!\n13. Quantification of thermo-poro-elastic effects: I disagree with the way the individual effects are calculated. There are three models: (I) Biot’s coefficients equal to 0; (II) thermal expansion coefficients equal to 0; and (III) both Biot’s and thermal expansion coefficients equal to 0. Strain are proportional to effective stress changes, Δ ′, which are a function of total stress changes, Δ , pore pressure changes, Δ , and temperature changes, Δ , as Δ ′ = Δ – Δ + Δ , where is Biot’s coefficient, is the thermal expansion coefficient and is the bulk modulus. Thus, model III gives the poroelastic effect, i.e., the deformation as a result of the stress changes caused by pore pressure changes (and not the effect of pore pressure changes as indicated in lines 375-376). The effect of pore pressure changes can be obtained by subtracting to the calibrated model the deformation of model II and model III. And thermal effects can be obtained by subtracting to the calibrated model the deformation of model I and model III.\nAns: Thank you very much for your comments! During this revision, we checked a classic reference on the topic, “Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology” written by H. F. Wang, where we found the equation of total strain of rock skeleton (see the attached screenshot for the expression).\nFrom equation 2.36 from Wang’s book, we can see that the total strain consists of two parts: poroelastic strain and free strain. The free strain is associated with pressure changes. Therefore, the equation of total strain of our reference model can be simplified as:\nThe corresponding equations of models I, II, and III can be simplified as:\n𝓔𝓔𝑑𝑑𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 𝓔𝓔𝑑𝑑𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑇𝑇𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃 +𝓔𝓔𝑑𝑑𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 +𝓔𝓔𝑑𝑑ℎ𝑃𝑃𝑃𝑃𝑒𝑒𝑇𝑇𝑇𝑇\nI:\nII:\n𝓔𝓔𝑑𝑑𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 𝓔𝓔𝑑𝑑𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑇𝑇𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃 +𝓔𝓔𝑑𝑑ℎ𝑃𝑃𝑃𝑃𝑒𝑒𝑇𝑇𝑇𝑇\nIII:\n𝓔𝓔𝑑𝑑𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 𝓔𝓔𝑑𝑑𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑇𝑇𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃 +𝓔𝓔𝑑𝑑𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃\n𝓔𝓔𝑑𝑑𝑇𝑇𝑇𝑇𝑇𝑇𝑇𝑇 = 𝓔𝓔𝑑𝑑𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑇𝑇𝑇𝑇𝑃𝑃𝑇𝑇𝑃𝑃𝑃𝑃\nWe agree with you that “model III gives the poroelastic effect”; but the effect of pore pressure changes can be obtained by comparison of the reference model with model I. The thermal effects can be isolated through comparison of the reference model with model II. In our MS, the poroelastic contribution was mistaken as the pressure contribution; and the pressure contribution was mistaken as the poroelastic contribution. We corrected the error during this revision.\n14. Discussion, lines 442-444: it is claimed that an improvement in the understanding of the spatial and temporal evolution of the strength of pore pressure and thermal effects is achieved in this study. However, the manuscript does not provide such better understanding. It is shown that there is a competition between pore pressure and thermal effects, with pore pressure presenting a rapid response and thermal effects becoming dominant in the long term, but this was already known. And it is shown that thermal contraction affects not only the cooled region, but also further away, which was already known. There is no quantification or clue how these processes can be quantified or estimated in other geothermal sites. In general, the observed coupled thermo-hydromechanical processes are already known. The novelty and interest of this study resides in the application to a geothermal site with data covering decades.\nAns: Respectfully, we disagree with these characterizations of the paper. We are not aware of previous studies that indicate that all of this ‘was already known’ in a generalized way. There are theoretical models for all of these phenomena, of course, but we have not claimed to invent new theory. In practice there are geodetic observations from almost every active geothermal field, or course, but this is the first study to have decades of time-varying surface displacement changes with advanced modeling capable of replicating the complex signals seen in time in space. Fundamentally, our prior understanding of deformation and seismicity at geothermal fields may be limited if not biased by the fact that most sites do not have such wealth of information. In other words, while the individual phenomena considered here have each been proposed to explain limited geodetic observations in the past, it is clear from this study that multiple phenomena act in concert. This finding would not have been possible in practice without long-term geodetic observations, and thus we will need far more cases where decades of geodetic data are available to establish general patterns. With that said, the modeling approach taken can be adopted at any geothermal field, which we note in the text.\nBelow we provide a little more discussion on what makes this study unique:\nAs you know that deformation of geothermal fields can be caused by reservoir pressure depletion, poroelastic and thermoelastic contraction. Such mechanisms can be active concurrently, but it is challenging to quantify their relative strengths and spatiotemporal evolution features, especially for large-scale fields with complex operation history. One major reason is due to the rarity of long-term seismic and geodetic observations combined with the extreme difficulty of setting up a realistic 3D thermo-hydro-geological model with limited subsurface geological and geophysical data (like reservoir geometry, formation lithology and fault distribution) needed to simulate decades of geothermal operations.\nIn this study, we quantify spatiotemporal evolution of the relative strengths of the three physical mechanisms for the Heber Geothermal Field (HGF) based on a 3D thermo-hydro-geological model, which is calibrated by multiple datasets, including geodetic observations over 20 years, precise well trajectories and operational records of fluid rates and temperatures. Our most important findings are:\n1. The effects of pressure flunctation and poroelastic response keep relatively stable in the operation history of the HGF. In contrast, the strength of thermal contraction effect increases with time: after ~30 years of operations, the thermal effect is 2-3 times larger than the two pressure effects.\n2. It is common that fluid injection causes surface uplift due to pore pressure increases within the target reservoir; however, the ground in west an east regions of the HGF was subsiding during some time periods with net increases of injected fluid. We show that the deformation anomaly is related to the effect of the HGF center siphoning fluid from surrounding regions and ever-increasing effect of thermal contraction.\n3. The long-term trends of surface deformation and seismicity growth at the HGF are controlled by the thermal effect. In contrast, the two pressure effects only drive instantaneous changes.\nThe finding that pressure and thermal effects dominate the early- and late-stage deformation of geothermal fields, respectively, helps to reconcile the controversy of differing deformation mechanisms reported for the same geothermal field, like at the Coso field (Fialko and Simons, 2000; Im et al., 2021) and the Brady Hot Springs field (Reinisch et al., 2020; Ali et al., 2016). In addition, the framework, developed for resolving the spatiotemporal evolution of the HGF deformation mechanisms through integrating multiple geodetic, geophysical and geological data, is applicable to other complex geothermal systems.\n15. The two-way coupling approach is partially neglected. In particular, rock deformation is assumed to have negligible effects on temperature and pressure fields. While this may be true in porous media, it may not be like that in fractured media, like the Heber field, which has several faults. Deformation may significantly enhance fracture/fault permeability because it scales with the cube of the aperture, affecting both the pressure and temperature field.\nAns: Thank you for this comment. We agree, in principle, and it is true that the cubic flow law is a well known empirical relationship between fracture aperture and hydraulic conductivity. However, incorporating this level of detail in our model would dramatically increase both the computational burden and the parameter space, and it is unclear that the added complexity would yield a ‘better’ result. It is also extremely challenging and uncertain to try and characterize the fracture system without significant relevant data, which we do not have. If we gained access to fluid pressure and temperature observations from within the reservoir, we could test for various flow regimes and their effect on the deformation field; unfortunately these are proprietary and likely never to be shared publicly.\nDuring our model calibration, we test the influences of different properties (mechanical, hydrogeological and thermophysical parameters) of the geothermal reservoir and faults on simulation results. These simulations reflect to some extent the influence of rock deformation on rock properties. In addition, our model calibration shows that the influences of hydraulic properties and fault thickness are largest, indicating that surface displacement and heat transfer of the HGF are mainly associated with fluid flow, and that the influence of rock properties is relatively minor. In our study, we focus on the dominant deformation mechanisms of the Heber geothermal field. Therefore, the coupling between rock deformation on temperature and pressure is not considered in our study.\n16. Calibration: it is a bit surprising that the initial guess (with few exceptions of slight changes in model properties) of model properties gives the best fitting. Calibrating a 3D THM model is a challenging task.\nAns: Yes, we totally agree! This work is based on more than two years of calibration work including: (1) designing and adjusting the calibration strategies; (2) adjusting the initial parameter settings based on prior studies on the Heber geothermal field, and the properties of water and rock. One can think of the starting parameters as the most informed guess using all the available constraints we could find, and the calibration procedure as testing variations from that set of parameters.\nMinor comments:\n1) Line 341: you refer to “injection volume” but the units are given in kg\nAns: Thank you very much for your careful review! We revised the sentence to be “In the western region of the HGF, although the net injection rates are around 0.6×109 kg per month on average from 1994 to 2005”.\n2) Line 561: the assumed values of fault properties are justified as “widely-used values”, but the two provided references are papers of the first author of this study. References to papers by scientists that are not authoring this manuscript should be provided\nAns: Thank you very much for your suggestion. We add three papers of Chang and Segall (2016), Zbinden et al. (2020), and De Simone et al. (2013) into our MS as relevant references:\n1. Chang, K. W., and P. Segall (2016), Injection-induced seismicity on basement faults including poroelastic stressing, J. Geophys. Res. Solid Earth, 121, 2708–2726, doi:10.1002/2015JB012561.\n2. Zbinden, D., Rinaldi, A. P., Diehl, T., & Wiemer, S. (2020). Hydromechanical modeling of fault reactivation in the St. Gallen deep geothermal project (Switzerland): Poroelasticity or hydraulic connection? Geophysical Research Letters, 47, e2019GL085201.\nhttps://doi.org/10.1029/2019GL085201\n3. De Simone, S., V. Vilarrasa, J. Carrera, A. Alcolea, and P. Meier (2013), Thermal coupling may control mechanical stability of geothermal reservoirs during cold water injection, Phys. Chem. Earth, 64, 117–126.\n3) Line 599: the mechanical boundary condition of “fixed injection and production wells” is not very clear. What is fixed? According to the figures in the Results, the surface displacement at the wells is not zero.\nAns: Thank you very much for your careful review! The open-hole sections of all wells are fixed to the material layers in our model. This implies there would be some deformation of the open hole section as the surrounding material deforms. We revised the sentence in MS.\n4) The units on both sides of Equation (3) are different: units are 1/s in the left-hand side, while they are m3/s in the right-hand side.\nAns: Thank you for your careful review! Again we checked the equations in H.F. Wang’s book (relevant section attached below). In equations 4.65 and 4.66 from that book, the term Q presents the volume of fluid per unit bulk volume per unit time. In our equation, Qm presents the corresponding mass source. Our equation is consistent with the equation 4.66 from Wang (2000).\n5) Line 700: please, provide references to support “typical values of faults”. Faults can take wide range of values, depending on the brittleness of the rock and accumulated slip.\nAns: This is a great point. We took the approach of assuming that because the faults are either high-conductivity flow zones (feeder fault) or plate boundary faults, that they should have reduced elastic moduli in the damage zone. You are correct there is significant variations possible. We based the initial values for the surrounding rock on the seismic velocity model converted to elastic parameters, and then reduced the moduli of the fault zones accordingly.\nResponse to Reviewer #3\nThis version of the manuscript is much easier and smoother to read than the first one. The arguments are logically exposed, and the demonstration is much clearer. However, the initial part of the argument is still a bit weak, especially concerning some terms (displacement, levelling, …) actually many things around Figure 3, which prevent full trust in the results. Once this step passed, and admitting the observations, the demonstration is better conducted. Therefore, I strongly suggest to revise once more the manuscript to make it smooth.\nAns: Thank you again for reading through our manuscript carefully. We have made further revisions based on the comments from you and Reviewer #2. Following are our point-by-point responses to your specific comments.\n1. Line 18. The term “siphonic effect” is a bit awkward here, and should be explained.\nAns: Thank you for your suggestion. We should have refered to this as ‘siphoning’ rather than ‘siphonic’. To be clear, the siphoning effect refers the phenomenon introduced in Line 14: high-yield production wells at the HGF center siphon fluid from surrounding regions. During this revision we are unable to elaborate in the abstract because of tight word limits (150 words). But, in the text we use the definite article “the” to indicate the phenomenon described in Line 14, and fix the phrasing.\n2. Line 19 or 20. “by cooling“ could be added after “Thermal contraction “. It seems obvious, but this is clearer by adding the reason.\nAns: Thank you for your suggestion. We added the words “by cooling” into Line 17, where the words “thermal contraction” first appear.\n3-1. Figure 3 is still confusing. Fig 3a. Uz is not defined. Is it position or velocity.\nAns: Thank you for your careful review! “Uz” is vertical displacement. We add the explanation into the caption.\n3-2. I do not understand why you normalize (line 104). Real values of displacement should be giving more information for possible inversion of the source location.\nAns: The displacement time series of the leveling stations are normalized for two purposes primarily related to visualization: (1) to represent the deformation trends of the west, central and east zones of the HGF, which have similar temporal characteristics but at different scales, and (2) to represent temporal variation of deformation of the whole region. By doing this, we can compare vertical deformation trends of the west, central and east regions with the respective net injection volumes of geothermal fluid of every month (Figure 3a). In addition, we can also compare the deformation trend of the whole region with temoporal changes of observed seismicity in the Heber geothermal field (Figure 4).\nTo be clear, the true values of displacement were used for model calibration, and misfits in both the hydro-geo-mechanical simulations and the fault slip inversions were calculated with actual observations and model predictions. We simplified the text to eliminate any confusion about how the normalized observations are used.\n3-3. In the west region, you show injection rate in blue. However, there are no wells in the west region, according to map Fg 3b. Where are the wells or where does this curve come from?\nAns: Thank you for your comment! As we introduced in Line 65, there are 28 injection and 31 extraction wells as well as 10 wells that were first operated for ~2 years of heat extraction, and then used for reinjecting heat-depleted brines since 1993. From Figure 1a and Figure 2 (we attach it below), you can see that there are only extraction wells in the west and east regions. The panels in Figure 3b are characterized with limited space; if we plot all wells in one panel, there will be too much clutter and overlap (as shown in Figure 1a). Consquently, we plot the injection wells, extraction wells, and extraction-injection wells in three panels of Figure 3b, respectively. See response to your comment #3-6 for Figure 3b.\n3-4. I would add an indication in the figure 3a west, central and east of the three periods indicated below (e.g., using P1, P2 and P3, which would correspond to the periods indicated in fig 3b.\nAns: We are wondering if there is a misunderstanding of Figure 3a. The three panels in Figure 3a show the vertical displacement trends of the west, central, and east zones of the Heber geothermal field.\n3-5. I do not understand the grey zone, as it seems there is also a decade-long displacement transient from 1995 until 2005 (e.g., fig. 3a.). The term “transient” should therefore been defined.\nAns: Thank you for your suggestion! We defined the “transients” in the MS as: The displacement rate and trend changes that initiated around 2005 and lasted for more than 10 years are referred to as “decade-long transients”.\n3-6. In Fig. 3b the caption indicates maps obtained from levelling and InSAR, but the top of each subfigure indicates only levelling. I understand the space is rather small to indicate both, but this is confusing.\nAns: In Figure 3b (we attach it below), there are 3 panels. The first panel shows the vertical displacement rate map derived from interpolation of leveling observations from January 1994 to December 2004. The last two panels show the vertical displacement rate map derived from interpolation of leveling and InSAR observations from 2006 to 2010. The purpose to show the InSAR map is to make a comparison with the leveling results during a comparable time-period.\n4. Line 123. There is confusing reference to displacement, although you refer mostly to rate of displacement, which is velocity (for instance in this line 123, where 1.5 cm/year is indeed a velocity).\nAns: Thank you for your comments, but we think that “displacement rate” are also common words and are widely used without any potential ambiguities.\n5. Line 124. There is inconsistency here with the main conclusion of the paper. You have both subsidence AND uplift with an increase of injection rate. Therefore, this observation contradicts line 15 of the abstracts. This needs further explanation and clarification.\nAns: Thank you for your comment! We revised the sentence of Line 15 to be “Here we show that high-yield production wells at the HGF center siphons fluid from surrounding regions, which can cause subsidence at low-rate injection locations.”.\n6. Line 128-129. I do not see where in Fig 3b, levelling-measured subsidence is large than InSAR results. If this is indicated in supplementary info, you should indicate this here. Now on the result itself, I do not understand why the levelling and InSAR height changes are different. You should give a reason for this difference.\nAns: From the last two panels of Figure 3b (we attach it above), you can see that the leveling-measured subsidence at the HGF center is larger than the InSAR result. These are, of course, different geodetic methods subject to different influences and sources of error and different resolutions in both time and space. Hence, one possible reason relates to velocities for both datasets being derived from linear regression for the period indicated. As the leveling shows, the rate of change in the vertical direction is non-linear, which suggests there would be some bias between the two vertical results, since there are different covariance structures and different number of data. We added a brief explanation into the main text.\nReferences: Eneva, M., Adams, D., Hsiao, V., Falorni, G. and Locatelli, R., 2019. Surface Deformation at the Heber Geothermal Field in Southern California. In Proceedings of the 44th Workshop on Geothermal Reservoir Engineering, SGP-TR-214. Stanford University, Stanford, California, February (pp. 11-13).\n7. Line 135. I do not follow exactly what you mean. “… vertical displacement … declining…” This is not very clear.\nAns: In the first panel of Figure 3a (we attach it below), you can see that the vertical displacement is characterized by subsidence before 2005 and uplift after 2005. The temporal variation trend (red error bars) is similar to the changes of injection rates (blue curve). During this revision, we rewrote the sentence to be “the vertical displacement is characterized by subsidence before 2005 and uplift after 2005 with the temporal variations which are similar to the changes of injection rates”.\n8. Line 138-139. I do not see this in figure 3. The west region indeed subsides before 2005, but I do see a reverse, so an increase only after 2005. However, before 2005, I do no see an increase of the injection rate. Then this sentence does not make sense. In addition, abnormal is a weird word for scientific expression. I would replace with “unexpected” or something like this. However, I do not see that in figure 3.\nAns: From 1994 to 2005, both the vertical displacement and injection rates are characterized by a declining trend. As you know, fluid injection is expected to cause surface uplift. Consequently, it is an unexpected phenomenon that the ground of the west region was subsiding from 1994 to 2005 even with a net increase of injected fluid. During this revision, we replaced the word “abnormal” with “unexpected”.\n9. Line 139-140. The sentence “the same displacement anomaly…” is unclear. Which anomaly? And it is in another region and another time... so quite confusing!\nAns: In the third panel of Figure 3a (we attach it below), you can see that the ground of the east region was subsiding from 2005 to 2017 with a net increase of injected fluid. This is the displacement anomaly that we want to highlight. During this revision, we replaced the sentence to be “In the east region, the ground was also subsiding from 2005 to 2017 with a net increase of injected fluid”.\n10. Figure 4 and line 168. Is the “thermal contraction” the same as thermoelastic effect (in the figure)?\nAns: Yes, they are the same. During this revision, we replaced the words “thermoelastic effect” in Figure 4 with “thermal contraction”.\n11. Line 171. It would be nice to precise if there any relationship between the observed horizontal displacement rate (2 cm/year) and the possibility you invoke here.\nAns: According to our modeling, for a distributed slip model to even come close to matching the InSAR-measured horizontal displacement from 2006 to 2010, the magnitudes of dextral strike-slip on the plate boundary fault need to be up to 30 to 82 cm below the depth of 8 km (Supplementary Fig. 8), which suggests 6 to 16 cm/yr of slip rates. If that were true, the slip rates would be much higher than the reported slip rates and the relative plate motion rate (only ~4 cm/yr), which cannot be reconciled. Consequently, the dominant reason for the horizontal displacement is unlikely related to fault slip.\n12. Line 320. Is it necessary to mention “Strong”?\nAns: We know that heat convection associated with fluid flow occurs everywhere within the geothermal reservoir. But the convection is stronger in some localized regions with higher higher porosity and permeability, which can cause prominent temperature decrease. So we think it is necessary to keep the word “strong”.\n13. Line 344 and 346. What is the unit “kg/mo”?\nAns: “kg” is the abbreviation for \"kilogram\"; \"mo\" is the most common abbreviation for \"month\".\nREVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nReview of the manuscript “Thermal contraction gradually dominates Heber geothermal field deformation”, by Jiang et al.\nThe authors have corrected some errors that they had in figures and in the description of the methods (see, for example, response to comment 13 of reviewer 2). The results have not changed, so I hope that the error was just in the description and that the analysis was performed correctly.\nThere are some aspects in which the authors and I continue to disagree, like the origin of the needed slip rate of the Feeder Fault to explain the observed deformation (I have explained, giving reference values from other sites, that a slip rate of 4-10 cm/yr can be achieved with aseismic slip of the fault). I do not think that such disagreement in the interpretation of the causes of the ground deformation should prevent this manuscript from being published, I think it is normal that we have different views, especially in complex sites like the studied one. Apart from this, I think there is already a quite good understanding of thermal effects and their temporal evolution with respect to poromechanical effects. I recall a figure from the report on induced seismicity of the National Academy of Sciences (2013) showing production and injection rates at The Geysers, together with monitored induced seismicity. Low-magnitude seismicity correlates with injection volumes, i.e., the poromechanical response in immediate and mainly local, affecting fractures in the vicinity of wells. Apart from this, there is an increase in the frequency of magnitude 3 earthquakes after 20 years from the start of water injection, and of magnitude 4 earthquakes after 30 years. This shows that the thermomechanical response takes time and that it extends further away than the cooled region, destabilizing distant faults. Similar observations have been shown through numerical simulations of geothermal systems. This is why I was saying that the explained processes at HGF are already known.\nResponse to Reviewer #2\nThe authors have corrected some errors that they had in figures and in the description of the methods (see, for example, response to comment 13 of reviewer 2). The results have not changed, so I hope that the error was just in the description and that the analysis was performed correctly.\nAns: Thank you very much for your comment. During last revision, we carefully checked our model and simulation results based on the comments from you and Reviwer #3. We found that the values of water properties in our model were fixed to their initial values and that the strengths of poroelastic effect was mistaken as the pressure effect. We corrected these errors during last revision; they did not change the conclusions of the paper.\nThere are some aspects in which the authors and I continue to disagree, like the origin of the needed slip rate of the Feeder Fault to explain the observed deformation (I have explained, giving reference values from other sites, that a slip rate of 4-10 cm/yr can be achieved with aseismic slip of the fault). I do not think that such disagreement in the interpretation of the causes of the ground deformation should prevent this manuscript from being published, I think it is normal that we have different views, especially in complex sites like the studied one.\nAns: We totally understand your concern, and appreciate your perspective. However, if the slip rates on the Feeder fault are indeed 4~10 cm/yr, the cumulative slip would be 1.2~3 m during 30 years of geothermal operation. If this kind of long-term aseismic slip exists, such large cumulative slip would cause catastrophic, irreversible casing/borehole deformation of the wells intersecting with the feeder fault; but there were no related reports on such damage and production rates have remained relatively constant, implying that major repairs to the production wells (or re-drilling) in response to fault slip has not been needed.\nApart from this, I think there is already a quite good understanding of thermal effects and their temporal evolution with respect to poromechanical effects. I recall a figure from the report on induced seismicity of the National Academy of Sciences (2013) showing production and injection rates at The Geysers, together with monitored induced seismicity. Low-magnitude seismicity correlates with injection volumes, i.e., the poromechanical response in immediate and mainly local, affecting fractures in the vicinity of wells. Apart from this, there is an increase in the frequency of magnitude 3 earthquakes after 20 years from the start of water injection, and of magnitude 4 earthquakes after 30 years. This shows that the thermomechanical response takes time and that it extends further away than the cooled region, destabilizing distant faults. Similar observations have been shown through numerical simulations of geothermal systems. This is why I was saying that the explained processes at HGF are already known.\nAns: Thank you very much for your comment. During this revision, we found the report on induced seismicity of the National Academy of Sciences (2013) from the following link:\nhttps://www.resolutionmineeis.us/sites/default/files/references/nas-2013.pdf, and the figure that you mentioned in page 65 of the report (also attached below).\nFrom this figure, we see that: (1) the number of EQs with magnitudes of M≥1.5 (not only low-magnitude seismicity) correlates with injection volumes; (2) the annual number of EQs with magnitudes of M≥3.0 is very low before 1985 and is around 25 after 1985. However, neither of these phenomena reveal one of our findings that thermal contraction dominates long-term trends of surface displacement and seismicity growth, while pressure effects drive near-instantaneous changes.\nIn addition, our study also has other important findings like: (1) the pressure effects are relatively stable during the operation history of the HGF whereas the thermal effect increases gradually with time and eventually overwhelms the pressure effects; (2) the subsidence anomalies in west an east regions of the HGF during some time periods with net increases of injected fluid are related to the effect of the HGF center siphoning fluid from surrounding regions and ever-increasing effect of thermal contraction, respectively. Moreover, the framework, developed for resolving the spatiotemporal evolution of the HGF deformation mechanisms through integrating multiple geodetic, geophysical and geological data, is applicable to other complex geothermal systems."}]

Decision letter