On the role of thermal stress and fluid pressure in triggering seismic and aseismic faulting at the Brawley Geothermal Field, California.

Surface deformation and earthquake swarms are correlated in space and time with operations at the Brawley geothermal field in southern California. The seismicity culminated in 2012, about 2 years after the onset of geothermal activity, with a M5.4 earthquake. These earthquakes occurred at a > 5km depth, much larger than the ~1km reach of the geothermal wells, raising questions about the triggering mechanism. Surface deformation shows that aseismic slip on a normal fault intersecting the geothermal reservoir preceded the swarm and possibly triggered it. In this study, we resort to geomechanical modeling to investigate how the sequence of aseismic and seismic slip unfolded. The modeling accounts for thermo- and poro-elastic stress changes induced by the geothermal operations and allows for inelastic deformation and faulting of the reservoir and surrounding me-dium. The simulation successfully reproduces the flow rates and well-head pressures reported by the operator as well as the measured surface subsidence. By varying the model parameters, we show that the surface subsidence is due to thermal contraction and normal faulting. The fault reactivation is driven by pressure changes and thermal unclamping. The pressure-driven reactivation is rapid and influences a larger area, while the temperature-driven reactivation is more gradual and more localized near the injection wells. In our simulation, aseismic normal faulting driven by the geothermal operation leads to elastic stress release via yielding and faulting within the reservoir volume and, conversely, to stress build-up beneath the reservoir, where the 2012 swarm developed. Such a stress transfer provides a plausible explanation for the 2012 Brawley swarm. Our study shows how a geothermal operation can, in principle, contribute to seismic hazard mitigation through the aseismic release of tectonic stresses within a geothermal field but points to the difficulty of mitigating the hazard posed by stress transfers in the surrounding area.


Introduction
The hazard posed by induced earthquakes is a significant issue of relevance to geothermal energy production, oil and gas extraction, and CO 2 storage [e.g., Ellsworth, 2013;Zang et al., 2014;Bujize et al., 2019;Kim et al., 2018]. Although this has long been known, our understanding of the mechanics of induced seismicity remains limited. As a result, it is rarely possible to describe quantitatively and unambiguously how a particular earthquake might have been induced. In this study, we focus on a particularly outstanding example of a swarm of seismicity in Southern California that culminated with an Mw 5.4 event in 2012. This swarm is thought to have been induced by the operations at the Brawley geothermal field (Fig. 1) [e.g., Wei et al., 2015].
Although the Brawley swarm occurred within the geothermal field, it was not immediately recognized as induced event because it occurred in an area of natural seismicity and most of the events, including the larger M>5 earthquakes, have an hypocentral depth of more than 6 km, which is much larger than the <2km depth reach of the wells (Fig. 1b). However, the fact that the swarm occurred not long after an abrupt increase in production and injection rates is a hint that it was indeed induced. Seismicity beneath the nearby Salton Sea geothermal field shows a similar correlation with operational activity and has been inferred to be induced [Brodsky and Lajoie, 2013]. It is not straightforward to diagnose if earthquakes are induced based on spatial correlations. Seismicity can indeed be induced at possibly large distances from the operating wells due to poroelastic stress changes [e.g., Segall et al., 1994;Zhai et al., 2019;Goebel et al., 2017;Goebel and Brodsky, 2018] or to fault-controlled hydraulic connectivity [e.g., Shapiro and Dinske, 2009].
In a previous study of the 2012 Brawley swarm, Wei et al. [2015] discovered that aseismic slip on a normal fault was activated first and demonstrated that this aseismic deformation contributed to a static Coulomb stress increase at the location of the swarm. They suggested that the aseismic fault activation was triggered by a pore pressure increase due to the injection. This scenario was corroborated based on observations of surface subsidence (Fig. 1c&d), inversion for these data for fault slip, and static stress transfer calculation. A similar scenario involving aseismic slip and seismicity induced by fluid injections has been shown to explain observations from controlled experiments [e.g., Gulglielmi et al., 2015;Lengliné et al., 2017]. In their study of the Brawley swarm, however, Wei et al. [2015] didn't analyze whether their interpretation was plausible from a geomechanical point of view. They didn't verify that the aseismic normal faulting derived from a kinematic inversion of surface deformation data could quantitatively be reproduced in a physics-based model of the effect of the geothermal injection and production. They also didn't account for the possible role of pressure variations and thermal stresses in explaining surface deformation and induced seismicity.
Thermal and pressure effects have been shown to be significant factors of surface deformation in geothermal fields [e.g., Im et al., 2017;Fialko and Simons, 2000;Reinisch et al., 2020a] and can contribute to fault reactivation [Gan and Elsworth, 2014] and seismicity rate changes [Im et al., 2021]. It is therefore probable that both thermal and pressure effects contributed to the surface deformation at Brawley, and to the triggering of the 2012 swarms. Surface subsidence at Brawley considerably accelerated when the injection pressure was increased (Fig. 1c&d). This observation shows that the subsidence was not driven by pressure depletion, which is often considered a primary subsidence mechanism in geothermal fields [e.g., Reinisch et al., 2020b;Barbour et al., 2016]. Instead, it is more probable that it resulted from thermal contraction, as inferred at the Coso geothermal field [e.g., Im et al., 2021], or from aseismic slip [Wei et al., 2015;Gan and Elsworth, 2014]. The aim of this study is, therefore, to relate the geothermal field operations with surface deformation and seismicity based on geomechanical simulations accounting for thermal and poroelastic effects.
Hereafter, we first briefly present the setting of the Brawley geothermal field. We next describe our simulation method and results. Finally, we discuss how our simulations shed new light on how the 2012 Brawley swarm was probably triggered by the geothermal field operations.

2012 Setting of the Brawley Geothermal field
The Brawley geothermal field is located in the Imperial Valley, south California (Fig. 1a). The Imperial Valley is a seismically active zone along the boundary between the Pacific and the North American plates (Fig. 1a inset). The Brawley area is a zone of strain transfer between the  [Hauksson et al., 2012] distributed by the Southern California Earthquake Data Center. The blue line denotes the surface rupture trace of the normal fault activated by the Mw 4.4 earthquake of August 26, 2012. Blue and black triangles represent locations of injection and production wells, respectively. Yellow and red stars denote the relocated epicenters, with associated focal mechanisms, of the two Mw>5 strike-slip events (yellow) and the two Mw>4 normal events (red) [Wei et al., 2015]. The San Andreas (SAF), San Jacinto (SJF), and Imperial Faults (IF) faults are shown in the regional map (inset). (b) EW cross-section showing the depth distribution of seismicity along the strike-slip fault that hosted the 2012 Brawley swarm (outlined by the parallelogram in panel a). We used QTM catalog [Ross et al., 2019] [Wei et al., 2015]). (d) Field production and injection history (https://maps.conservation.ca.gov/doggr) and surface subsidence combined from InSAR and leveling [Wei et al., 2015]. Imperial Fault to the south and the San Andreas Fault to the North, which, along with the San Jacinto fault, accommodates most of the relative plate motion [e.g., Brothers et al., 2009]. Transtension there is associated with high heat flow and hydrothermal activity [Younker et al., 1982], making the valley area a good candidate for geothermal field development. Accordingly, one of the largest geothermal fields in the US, the Salton Sea geothermal field, is located north of the Brawley geothermal area. However, active seismicity has been a source of considerable concern in developing geothermal resources in the area.
The Brawley field was first developed in the 1970s but abandoned after a few years due to scaling and corrosion problems. The development resumed in 2006, and test wells were drilled in 2007. The field came online in 2009 [Matlick and Jayne, 2008]. Public records report flow rates at a total of 42 wellbores, including 19 producers and 23 injectors (https://maps.conservation.ca.gov/doggr). Based on the well completion reports, the open hole injection interval (blue vertical lines in Fig. 1b) ranges between a depth of ~630m (averaged top) and ~ 1100m (averaged bottom). The initial flow rate was rather low, at ~0.3 Megaton/month, but rapidly increased to ~1.7 Megaton/month by mid-2010 (Fig. 1d). The flow rate increase is due to the pressure change at a few injection wellbores (Fig. 3c). Injection and production rates are generally comparable, so the operations should not have resulted in any significant groundwater flow beyond the extent of the geothermal field.
The relocated instrumental seismicity (Hauksson-Yang-Shearer (HYS) catalog [Hauksson et al., 2012]) outlines a network of NE striking left-lateral and less prominent NNW striking right-lateral strike-slip faults (Fig. 1a). The seismicity has been active throughout the period between 1981 and 2019 ( Fig. 1e), but the larger events (Mw>5) occurred within four years from the beginning of the geothermal operation during the 2012 swarm. The seismicity associated with this swarm forms a well-defined NE striking lineament, outlined by the parallelogram in Figs. 1a, consistent with the NE strike of the right-lateral focal plane of the largest magnitude event in the swarm (M5.4) [Wei et al., 2013]. The time evolution of seismicity within this zone shows a slight increase at the beginning of the exploratory drilling in 2007 and a further and more prominent increase when power production started in 2009 (Fig. 1d&f). This temporal correlation supports the view that the Brawley swarm was triggered by the geothermal operation but is not very compelling and does not demonstrate causality. The correlation in space is even less obvious due to the gap between the reservoir (average injector bottom hole depth < 1.1km) and the depth range of the 2012 swarm. Most of the events occurred at a depth of 8 to 10 km, and few events occurred at a depth shallower than 6km. This gap is clear in the QTM (Quake Template Matching) catalog, which reports seismicity only for the period during 2008-2017 but has a significantly smaller detection threshold estimated to M~0.3 ( Fig. 1b) [Ross et al., 2019].
Leveling and InSAR data showed a clear signal of subsidence which started before the 2012 swarms (Figs. 1c &d) and was amplified during the swarm leading to normal faulting at the surface [Eneva et al., 2012;Wei et al., 2015]. The asymmetric pattern seen in the geodetic data prior to 2012 suggests that the fault started slipping aseismically before the onset of the swarm [Wei et al., 2015]. The normal fault also did produce a shallow Mw 4.4 event a few hours after the largest event in the 2012 swarm (Fig. 1a). This event was probably the main cause of the surface ruptures that were observed in the field. The subsidence rate significantly increased when the flow rate increased in 2010 (Fig. 1d), suggesting that the aseismic faulting resulted from the mechanical impact of the geothermal operation. We note that the seismicity is mostly confined to the footwall of the normal fault (Fig. 1b).

Simulation set-up
We use Tough-FLAC [Taron et al., 2009;Rutqvist et al., 2002] to simulate the effect of the Brawley geothermal operation. This simulator accounts for the coupling between fluid flow, thermal effects, and deformation, whether elastic or inelastic. So, stress changes due to Green elements represent the host rock (elastic and impermeable), light brown elements represent the fault zone (can fail, permeable, cohesionless), and light blue elements represent the reservoir rock (can fail, permeable). The blocks that host injection and production wellbores are represented in blue and red, respectively. The location of well-hosting blocks at 850m depth is shown in the bottom left. Blue blocks with an arrow represent the three injectors with the boxcar pressure step (Fig. 3c). (b): Initial effective stresses once steady-state is reached during the model initialization stage. The principal stresses are parallel to the model axis with σ z ~ σ y > σ x . We assume a value of σ x low enough so that the shallower 2km upper portion of the normal fault is critically stressed.The detailed input parameters are presented in Table 1. temperature changes and poro-elastic effects are accounted for. We consider a 50km × 40km × 30km domain and assume symmetry with respect to the vertical x-z plane. Simulation is thus conducted in half of the domain (Fig. 2a). The reservoir is represented by a 5.5km × 3km × 1.5km volume embedded at the top of the domain (light blue in fig. 2a). The normal fault is represented by a 20m thick fault zone striking parallel to the y-direction and dipping 55 • westward. The fault cuts across the reservoir and host rock (light brown in Fig. 2a).
The reservoir and fault zone elements are assumed to fail according to Mohr-Coulomb criteria. The appropriateness of this criterion, which is insensitive to the intermediate principal stress, to assess the possibility of faulting and triggering seismicity is well established [King et al., 1994]. The friction coefficient is set to 0.6 for both fault and reservoir, and a cohesion of 10 MPa is assigned to the reservoir elements. The fault zone is assumed cohesionless. Friction is constant once the condition for failure is reached, meaning neither hardening nor weakening is associated with inelastic strain. The host rock ( Fig. 2a light green) is fully elastic. All elements are assigned a volumetric thermal contraction coefficient 4.5 × 10 − 5 /K, a bulk modulus of 10GPa, and a Poisson's ratio of 0.2. These thermo-elastic properties are typical of quartz-rich rock ranges 2.5 -6 × 10 − 5 /K at 50 -200 • C [Cooper and Simons, 1977;Somerton, 1992]. We also conducted 'no-thermal-stress' and 'no-fault' simulations to assess the influence of fault slip and thermal stressing in our reference model. In the no-thermal-stress case, the thermal expansion coefficient is set to zero. In the no-fault case, the fault zone elements inside and outside of the reservoir have properties identical to the reservoir and host rock properties, respectively. Inertial effects are not solved for in these simulations, so that deformation is treated as quasistatic.
The stress field accounts for gravity and tectonic loading. Gravitational body forces are applied to all elements with an effective density (rock densitywater density) of 1400kg/m 3 . The principal horizontal stress directions are assumed parallel to the Ox and Oy axis with σ y > σ x in accordance with regional constraints on the stress field [Zoback, 1980;Hardebek & Hauksson, 2001]. We assume vertical stress (the z-axis in our model) and maximum horizontal stress (the y-axis in our model) are nearly equal in magnitude so that both normal fault and strike-slip fault are critically stressed as observed in the field. Our stress tensor is consistent with the normal fault that produced aseismic slip at the surface and with left-lateral strike-slip motion on the N70E striking fault observed during the M5.4 maximum event in the 2012 swarm. Horizontal stresses are applied at the side of the domain with a depth-dependent stress gradient. For the yy normal stress, the gradient is 93% of the gravitational stress gradient. For xx normal stress, the stress gradient is 33% at depth < 2km and 61% at depth >2km of gravitational stress so that the fault is critically stressed at shallow depth. We initialize the model by running simulations with account for the boundary conditions and gravity until a steady-state is reached. Therefore, the stresses after initialization (Fig. 2b) are not identical to the input stresses. All simulation results shown in this work represent the result after this initialization stage.
For the computational efficiency and stability, we assume uniform initial fluid pressure and temperature over the entire domain of 15 MPa and 180 • C, respectively. We are using an effective rock density (rock densitywater density), so that the effective stresses are consistent with an hydrostatic initial condition. The assumed initial temperature The pressure data show that several injection wells (colored) are highly pressurized during ~2010-~2014. The blue dashed line shows the assumed 'boxcar' pressure history at the three main injection wells (Fig. 2a). The injection pressure at the other five wells is kept constant.
(180 • C) is based on the reservoir temperature measurement from a previous study [Matlick and Jayne, 2008]. The initial uniform pressure represents the hydrostatic condition (pressure gradient is equalized with fluid gravity), and therefore, the fluid flow is mainly controlled by wellbore pressure change. The applied pressure (15MPa) is approximately a hydrostatic pressure at the bottom of the reservoir. We applied 15 MPa of constant boundary normal stresses to cancel out the pressure in effective stress. We consider that these simplifications are acceptable since the temperature change occurs mostly within the reservoir, the pressure dependence of water hydraulic properties (i.e., density and viscosity) are insignificant, and fluid influx from outside of the reservoir is insignificant since production and injection rates are nearly equal (Fig. 1d).
We set eight injectors and seven producers within the reservoir between 640 m and 1070 m depth. We selected only a subset of the actual wells to simplify the simulations. A Peaceman well-block pressure model [Peaceman, 1983] is employed with a virtual wellbore radius of 10 cm and a skin factor of -5 located within the well block shown in Fig. 2a. The negative skin factor reflects hydraulic and/or thermal stimulation as typically observed in geothermal injectors [e.g., Sanyal, 1987]. Initially, 0.3 MPa over-pressure and 0.1 MPa under-pressure are applied to injection and production well, respectively. From year 1.5 until year 5.5 (for four years), the injection over-pressure of 3 wellbores is increased to 1.9 MPa (Fig. 3c). This boxcar shape pressure step matches approximately the operational data (Fig. 3c). Constant reservoir and fault permeabilities are applied at 13 md and 39 md, respectively. The host rock permeability is set low at 0.3 µd. The permeabilities and well skin factor were adjusted by trial-and-error to reproduce the flow rate data.
Our model is a simplification of reality in particular because we assume a uniform initial pressure and temperature. A more realistic simulation would require the pressure and temperature gradients. Also, we assumed simultaneous injector pressure change, but the operation data show a much more complex pressure history (Fig. 3c). Despite these simplifications, the model is satisfying with regard to a few key points: (1) the pressure change (timing and magnitude) reflects the operation data (Fig. 3c); (2) the flow rate, which determines the thermal stressing rate, corresponds to operation data (Fig. 3a); (3) the simulated surface subsidence (magnitude and shape) matches the leveling observations (Fig. 5). Therefore, we believe the simulation provides a reliable firstorder estimation of the deformation and stress change induced by the geothermal operations.

Flow rate and temperature
Our simulation yields flow rates in relatively good agreement with the field data (Fig. 3a). The increased flow rate from mid-2010 to mid-2014, during which the 2012 swarm occurred, is due to the boxcar pressure increase (Fig. 3c). We note, however, that the simulated flow rates do not reproduce the monthly fluctuations (Fig. 3a) because of the simplified pressure history. It results that the pressure-driven stress Fig. 4. Temperature distribution at the end of the simulation (~2019) at a depth of 900m. The three zones of maximum temperature change correspond to the wells with increased pressure and, therefore, the higher injection rate during the 2010-2014 period (Fig. 1a&3c). Triangles denote locations of production wells.  Fig. 3). In reality, the well-head pressure fluctuation must have induced pressure-driven stress variations throughout the period of production, which are not captured by our simplified model.
Production temperature shows almost no change over the ten years of operation (Fig. 3b) because the temperature changes are most significant only near the injection wells (Fig. 4). Therefore, the thermal stress change is expected mostly near the injection area. In effect, the temperature at the end of the simulation is lowest around the three highpressure (i.e., high-flowrate) injection wells (Fig. 4), and the temperature drop does not appear to have reached the production wells (Fig. 4  triangles).

Fault slip and surface subsidence
The simulated fault slip and surface subsidence are in agreement with the field observations (Figs. 5 and 6). Comparison between the observed subsidence pattern (Fig. 5a) and our simulation (Fig. 5b-d) show that the subsidence is driven by both aseismic fault slip and thermal contraction. When fault slip is not considered (Fig. 5c), the surface above the reservoir shows a pressure-driven uplift with thermally driven subsidence localized close to the wells. In this case, fault elements outside of the reservoir are assumed fully elastic, and therefore normal fault slip is not allowed. Hence, the no-fault-slip simulation results in only pure volumetric change from thermal contraction and pressure change. In the no-fault-slip case, the gradual subsidence driven by thermal contraction is more compensated by uplift driven by pressure change between mid-2010 and mid-2014 (Fig. 6c). By contrast, the observations show a sustained subsidence throughout that period ( fig.  6c blue dashed line). As a result, the amplitude of the field subsidence by 2019 is much larger than predicted by the no-fault simulation. These results show that the subsidence observed in the field is primarily driven by aseismic normal fault reactivation and a lesser degree by thermal contraction. The impact of pressure-driven volume change seems negligible.
In our reference simulation, the subsidence rate increase at mid-2010 is primarily due to the high pore pressure (Fig. 6c red bold line). This accelerated subsidence is also predicted by the no-thermal-stress simulation ( Fig. 6c red dashed line), but it does not occur in the no-fault simulation ( Fig. 6c red dotted line). Hence, in our simulation, the strong subsidence in 2010 primarily results from pressure-driven fault reactivation. The strong initial subsidence in 2011 is followed by slower continuous subsidence at ~2cm/year ( Fig. 6c red bold line) until mid-2014. This gradual subsidence is not observed in the no-thermal-stress case ( Fig. 6c red dashed line), implying that this continuous subsidence is a thermally driven process. The no-fault case (Fig. 6 dotted line) also shows thermally driven subsidence but with a much lower rate than our reference simulation, implying that the thermally driven subsidence is not only a result of the thermal contraction itself but also a result of fault slip [e.g., Gan and Elsworth, 2014].
Our simulation shows fast subsidence at mid-2010 due to aseismic motion on the normal fault, followed by less rapid but continuous thermally driven subsidence until the end of the high flow rate period in 2014. The temporal evolution of subsidence predicted by our simulation is roughly similar to the observed one ( Fig. 6c blue dashed line). The data do, however, show a more gradual increase of subsidence in 2010. The difference may be due to our simplified scenario, which assumes simultaneous abrupt pressure increases at the three main injection wellbores. The activation of the normal fault can be due, in principle, either to the local pore pressure increase or to the reduction of the normal stress driven by thermal contraction. When thermal stress is not considered (Fig. 5d), the pore pressure increase leads to early rapid subsidence due to the fault activation. However, the predicted surface pattern misses the intense subsidence in the area around wells, which is clear in the field observation (Fig. 5a). So, both the spatial pattern and the gradual increase of subsidence after the onset of fast injection indicate a significant thermal effect similar to that produced in our reference simulation.
The activated portion of the normal fault extends well beyond the reservoir, but the slip distribution peaks within the reservoir, where cooling is most significant (Fig. 7a). Comparison between the simulation result ( Fig. 7b solid lines) and the no-thermal-stress case (Fig. 7b dashed  lines) shows that the fraction of slip driven by thermal contraction is  mainly localized at a depth less than ~ 2km. This localized area corresponds to the area of temperature decrease (Fig. 7d). The slip at greater depth must be driven by pore pressure diffusion within the fault zone and by the static stress transfer due to fault slip at a shallower depth.
In any case, our reference simulation, in which normal fault slip is driven by pore pressure and thermal contraction, predicts a subsidence pattern in better qualitative and quantitative agreement with the observations than the models ignoring fault slip or thermal effects . We conclude that the observations require a combination of reservoir volume change and normal fault slip, driven by pore pressure increase and thermal contraction.

Stress change on the deep strike-slip fault activated during the 2012 swarm
We estimated stress variations on a fault that coincides with the location and geometry (vertical, NE striking) delineated by the 2012 swarm. This fault orientation is also consistent with the right-lateral fault plane derived from the modeling of the Mw 5.4 mainshock of August 26 [Wei et al., 2015]. We find that our reference model implies significant stress changes at depths well below the injection and production wells (Fig. 8). Given that the thermal stresses are mostly localized near the injection area (Fig. 7), this is primarily due to the reactivation of the normal fault. According to our simulation, diffusion of the pressure within the fault zone also plays a role. Our reference model leads to a pore pressure increase of ~ 0.1MPa by 2012 at 6-8 km depth (Fig. 7c). Our simulation shows that slip on the normal fault results in both normal stress reduction and shear stress increase on the vertical strike-slip fault. These effects contribute to a Coulomb stress increase in a large area in the footwall of the normal fault at a depth between ~7km and 15km (Figs. 8&9). This area of significant Coulomb stress increase corresponds well to the location and depth of the Brawley swarm (Fig. 9). The model shows a zone of increased seismicity at a depth shallower than 1km, which does not coincide with any detected seismicity. It is probable that at these shallow depths, rocks are dominantly rate-strengthening so that earthquakes can't nucleate [Blanpied et al., 1991]. This would be consistent with the observation of shallow aseismic fault motion on a number of different faults in the Imperial Valley area (Donnellan et al., 2014).
The Coulomb stress increase is initially rapid when the normal fault slip is activated by mid-2010, then slowly increases until injection pressures are reduced in mid-2014. The first M>4 event occurs in a few months after the significant Coulomb stress increase in 2010, but the Brawley swarm in 2012 somewhat lags behind the Coulomb stress increase predicted by our model (Fig. 8c). The rapid increase in 2010, however, should be considered with caution, given our simplified boxcar representation of the injection pressure history. In any case, our simulation suggests that the 2012 Brawley swarm was triggered by Coulomb stress changes of the order of only 0.01MPa. It is interesting to note that the seismicity settled down a few months after the onset of the Brawley swarm, even though the flow rate remained high for a few more years (Fig. 8c).

Thermally driven aseismic slip
Our study provides support for Wei et al.'s [2015] view that the 2012 Brawley swarm resulted from the Coulomb stress increase on a deep-seated strike-slip fault driven by aseismic slip on a shallow normal fault within the Brawley geothermal field. The kinematic model of Wei et al. [2015] indeed closely resembles the prediction of our reference simulation. Our simulations provide new insight into the mechanism that drove aseismic normal faulting and on the thermal effects that were not analyzed in the previous study. In our simulations, slip on the normal fault is driven both by unclamping due to pore pressure increase and thermal contraction within the geothermal field.
The pattern and amount of slip needed to match the observations of surface deformation are reproduced assuming standard mechanical properties of rock. Our simulation shows that the subsidence observed at the surface is mainly due to fault slip. The slip is primarily driven by pressurization at injectors. However, a significant fraction (~ 40% in our simulation) of maximum subsidence at the time of the 2012 Brawley swarm (Fig. 6c) is likely of thermal origin due to the effect of thermal contraction and fault-slip induced by thermal unclamping. The role of faulting on surface deformation is strongly supported by asymmetric shape (Fig. 1c) and surface rupture, which were not observed in the nearby Salton Sea geothermal field [Barbour et al., 2016]. Fault reactivation can thus add to the effect of pressure depletion [e.g., Barbour et al., 2016] and thermal contraction [e.g., Reinisch et al., 2020a] that are generally observed at geothermal fields. All three mechanisms are actually mechanically coupled and probably contribute jointly in general to the observed surface deformation [Im et al., 2021]. Note, however, that pressure depletion considered alone would tend to stabilize faults within the reservoir and put the overburden in horizontal compression. That would not explain the observation of normal faulting at the surface and of induced seismicity. The reasons for the different behaviors at the Salton Sea and Brawley geothermal fields are unclear, and we suspect that the role of thermal stresses at the Salton Sea field may have been underestimated.
We did not analyze here why slip on the normal fault was mostly aseismic. We note that the low initial normal stress at shallow depth and the further decrease of effective normal stress due to thermal contraction and pore pressurization should favor aseismic slip according to both theory and in-situ experimental results [Guglielmi et al., 2015;Cappa et al., 2019].

Possible assistance from direct pressurization
Our model doesn't involve a direct hydraulic connection between the geothermal field and the swarm. This mechanism has been documented from the spatio-temporal evolution of seismicity observed in a number of examples of seismicity induced by fluid injection [e.g., Fischer et al., 2008;Shapiro et al., 1997]. The direct effect of pore pressure diffusion is estimated to be very small in our simulation since we assumed very low permeability of 3 × 10 − 19 m 2 beneath the reservoir based on the previous study [County of Imperial Planning Department,1979;Matlick and Jayne, 2008;openei.org] except within the preexisting normal fault. In our simulation, the high permeability of normal fault at a deeper area However, the choice of geometric and hydraulic parameters in our simulation is somewhat arbitrary, and we have neglected the effect of fault reactivation on permeability [e.g., Guglielmi et al., 2015;Im et al., 2018;Cappa and Rutqvist, 2011]. Assuming a larger fault-zone permeability or a permeability increase with shear strain (for both the normal and strike-slip faults) would have resulted in different pressure distribution. The high permeability zone associated with the normal fault could also be much thicker, and therefore, a wider area of the strike-slip fault could have been directly pressurized. Abundant seismicity connecting the reservoir and swarm area was detected in December 2010 (Fig. 10) when the significant subsidence initiated (Fig. 10 inset). This observation may imply some stimulation below the reservoir. A shallower weak swarm was also detected roughly two days before the main 2010 swarm (Fig. 10b). These events might have been driven by shallow poro-thermo elastic stress change from the geothermal operation or by a possible vertical hydraulic connection between the reservoir and the swarm area. If the strike-slip fault was stimulated, the increased pressure at the injectors might have reached the depth of the swarm. Alternatively, such stimulation may redistribute the natural pressure distribution if the initial pressure was far from hydrostatic. The pressure change at depth can contribute to an additional Coulomb stress increase.

A swarm triggered by small stress perturbations
Our model is used to validate a possible scenario to connect the geothermal field with the observations (seismicity and subsidence). We acknowledge that the comparison of the modeling results with the observation is more qualitative than quantitative. Coulomb stress variations are estimated only to first order because the details of the faults geometries are not known. We did not carry out the sensitivity tests needed to assess uncertainties due to the computation efforts. However, we believe that our best-fit simulation provides a reliable first-order estimation since it successfully reproduced flow rate, wellbore pressure, and surface subsidence. The pattern and amplitude of subsidence and Coulomb stress changes in our simulation, which we interpret to have triggered the swarm in 2012, are similar to the observation and those estimated by Wei et al. [2015] based on the normal fault slip distribution estimated from surface deformation. These stress changes are estimated to be of the order of 0.01MPa. This value is much smaller than the typical 1-10MPa stress drop during earthquakes. So, these events classify as 'triggered' according to the terminology suggested by McGarr et al. [2002].
The geothermal operation probably contributed to initiating the swarm but was not a significant source of the elastic strain released by the swarm. The bulk of the released strain was probably of tectonic origin. The Brawley area is indeed a zone of continuous strain build-up and seismicity (Fig. 1e). Geothermal operation temporarily boosted the seismicity rate. However, the boost disappeared a few months after the onset of the Brawley swarm, even though the geothermal flow rate and injection pressures remain at the same level until mid-2014 and Coulomb stress kept increasing (Fig. 8c). This observation suggests that fault reactivation during the swarm released most of the shear stress that was initially available to drive earthquakes at depth.
The stress changes that triggered the 2012 swarm are about ten times smaller than the Coulomb stress changes estimated to drive aftershocks following large earthquakes [King et al., 1994]. Swarms may be sensitive to smaller stress changes than regular earthquakes. This will be expected if they occur in a zone of high pore pressure, as is commonly assumed [e.g., Thomas et al., 2012]. According to the rate-and-state model of earthquake nucleation, the seismicity rate is multiplied by exp(Δτ/a(σ-P)) [Dieterich, 1994], where Δτ is the Coulomb stress change, a is friction rate parameter, σ is normal stress, and P is pressure.
Accordingly, a small Coulomb stress change can result in some significant triggering if 'a' is small or the pore pressure is high. However, a high pore pressure should, in principle, favor aseismic slip, so it is unclear that a large pore pressure is a right explanation for the high sensitivity of the seismicity in the swarm area.

Implication for mitigation strategy of induced earthquakes
According to our simulations, the geothermal operations at Brawley released elastic shear stress, which had accumulated as a result of tectonic loading, via a combination of aseismic inelastic deformation and triggered seismicity. A similar strain release occurred at the Coso geothermal field and can explain the lack of aftershocks following the 2019 Ridgecrest earthquakes [Im et al., 2021]. To mitigate induced seismicity, one may consider releasing the elastic shear stress that can potentially drive large seismic events. However, the 2012 Brawley swarm is a notable example that shows that relatively large seismic events can be triggered by only a small stress change in the underburden, showing that it would probably be very difficult to control the impact of a release of stress in one area on the seismicity in the surrounding area. A Coulomb stress release in an area results necessarily in a Coulomb stress increase in another area where earthquakes can thus be triggered. Stress changes in geothermal reservoirs are inevitable since heat extraction itself and the associated thermal stresses imply stress redistribution within and around the geothermal reservoir. Assessing the presence of preexisting faults around the reservoir, their initial state of stress, and the expected stress changes are therefore required to analyze the potential earthquake risk associated with the development of a geothermal field. The field-scale modeling developed in this study can then be used to assess seismic hazards and explore possible mitigation strategies.

Conclusion
This study demonstrates that pore pressure changes and thermal unclamping due to the geothermal operations at Brawley probably drove aseismic motion on a normal fault intersecting the geothermal reservoir. Thus, the study provides a geomechanical foundation to the scenario initially proposed by Wei et al. [2015]. We have indeed shown that the 2012 swarm was probably triggered by stress build-up in the footwall of the normal fault intersecting the geothermal field with the possible assistance from direct pressurization due to fault-enabled pressure diffusion beneath the reservoir. According to our scenario, the elastic stress release by aseismic faulting and by the swarm is in any case of tectonic origin. So, a geothermal operation can contribute to releasing some of the elastic strain energy available to drive seismicity. Indeed, the seismicity rate is significantly reduced after the 2012 swarm, as similarly observed in the Ridgecrest aftershock suppression [Im et al., 2021]. Thus, in principle, it indicates the future seismic hazard is mitigated to some degree at the depth. However, our simulation also shows that the detailed strategy to avoid large seismic events during stress release should be carefully considered.
Our study lends support to the view that the reach of seismicity induced by a geothermal injection can be augmented substantially due to aseismic fault reactivation, as also suspected in other case examples [Lengliné et al., 2017]. Our study also points to the joint role of pressure and thermal contraction in driving aseismic slip. While it is well known that thermal effects can be a cause of anelastic deformation and permeability enhancement near boreholes [Sanyal et al., 1987;Ghassemi et al., 2008;Gan and Elsworth, 2014], the role of thermal stresses in driving inelastic deformation at the larger scale and in inducing seismicity is only being recognized [Im et al., 2021]. Table 1 CRediT authorship contribution statement Kyungjae Im: Conceptualization, Methodology, Software, Writingoriginal draft. Jean-Philippe Avouac: Conceptualization, Validation, Writingreview & editing.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  [Cooper and Simons, 1977] 2) Measured reservoir temperature: 149 -204 • C [Matlick and Jayne, 2008]