Source parameters of the 8 February 2016, Mw=4.2 Los Humeros earthquake by the inversion of InSAR-based ground deformation

We investigate ground movements induced by the 8 February 2016, Mw=4.2 earthquake at the Los Humeros Geothermal Field (Mexico) using Sentinel-1 radar interferometry. Previous estimated focal mechanism solution based on seismic data with a hypocentral depth of 1900 m could not resolve the measured coseismic surface deformation pattern. In this study, we applied inverse elastic dislocation models to estimate the source parameters of the seismic event. Our models suggest the reverse reactivation of the Los Humeros normal fault at a shallower depth (


Introduction
Superhot Geothermal Systems (SGHS) with their high enthalpies are important geothermal reservoirs. These systems, with temperatures often exceeding 350 • C, are frequently marked by recent or even active tectonic or volcanotectonic deformations. This results in complex reservoirs in which fault structures may act as prolific permeability pathways for geothermal fluids. Conversely, exploitation of such fields may yield considerable risks for induced seismicity. Indeed, seismic events can cause significant damage in the vicinity of the geothermal sites and may even lead to the termination of production (e.g. Deichmann and Giardini, 2009;Kim et al., 2018). The understanding of fault geometry and behavior is therefore crucial for a safe and sustainable exploitation of these important resources.
Synthetic Aperture Radar Interferometry (InSAR) can be used to estimate ground movements and the increased temporal resolution of recent satellite missions bolsters its deployment as a cost-effective method for long-term monitoring of SHGS. Previous studies based on InSAR time series analysis, for instance at the Reykjanes Geothermal Field, Iceland Receveur et al., 2018) and the Los Humeros Geothermal Field, Mexico (LHGF, Békési et al., 2019a) demonstrate the added value of surface deformation studies to understand reservoir dynamics and induced seismicity. The interpretation of ground deformation together with other geological and geophysical measurements can contribute to understand the complex behavior of fault structures controlling a geothermal field and provide valuable information for geothermal exploitation.
We have employed radar interferometry using Sentinel-1 images to constrain fault geometry and kinematics at the Los Humeros Geothermal Field in Mexico (LHGF, Figure 1a). We have studied the coseismic deformation due to the Mw=4.2 earthquake on February 8, 2016. The epicenter of this event had been located along the trace of the Los Humeros fault (Figures 1a-b and 2) (Lermo Samaniego et al., 2016), which had previously been interpreted as a normal fault, as indicated by the topography and displacements of geological units (Figure 1c). However, surface motions mapping the coseimic deformation and the focal mechanism solution based on seismological data indicated reverse movement with a left lateral component (Figure 1b). Initial results showed that the proposed mechanism of the moment tensor solution does not resolve the measured coseismic surface deformation pattern (Békési et al., 2019b). In this study, we have revised our earlier elastic dislocation models by inverting the InSAR-derived surface movements with constraints from existing geological and geophysical models. Our findings for the source parameter of the earthquake have been corroborated by an InSAR time-series analysis, indicating prolonged surface motions after the event, and an investigation of the underlying processes.

The Los Humeros Geothermal Field
The LHGF has been exploited since 1990 by the national Mexican electrical company (Comisión Federal de Electricidad (CFE)). 25 wells are presently producing a total of 6 Mt of steam yearly, with an installed capacity of ca. 95 MW of electric power (Gutiérrez-Negrín, 2019) ( Figure 1a). Further exploitation of the geothermal field is being planned, in particular in the area northwest of the current production zone, where temperatures locally reach 400 • C at 2.5 km depth (Gutiérrez--Negrín, 2019).
The geothermal field is situated within the Los Humeros Volcanic Complex (LHVC) in the Eastern sector of the Trans-Mexican Volcanic Belt, Mexico. It consists of two nested calderas (Figure 1). The Los Humeros Caldera was formed ~160 ka ago, with the younger (70 ka) Los Potreros Caldera nested in the western part, hosting the geothermal field  Norini et al. (2019). The Los Potreros caldera rim (blue line), fault structures (orange and red lines), the epicenters of the most significant (Mw > 3.0) earthquakes (red and yellow stars, (Lermo Samaniego et al., 2016) and the production and injection wells of the Los Humeros Geothermal Field are presented. (b) The focal mechanism solution of the 08 February 2016 earthquake (from Lermo Samaniego et al., 2016). (c) W-E cross-section along the location of the 8 February 2016 earthquake showing the geometry and kinematics of the fault structures based on Norini et al. (2019) and this study. The trace of the section is shown in (a). . The geothermal reservoir is built up by 10.5-1.5 Ma years old pre-caldera fractured andesites with an average depth of 1.5 km (Old Volcanic Succession unit, Figure 1c). The cap rock of the hydrothermal system is mostly composed of welded ignimbrites and post-calderas volcanic deposits (rhyodacite, andesite, and basaltic andesite) of the LHVC unit ( Figure 1c, Norini et al., 2019;Norini et al., 2015). The LHGF is controlled by the recent and still active fault system in the central part of the Los Potreros caldera (Figure 1, Norini et al., 2015) by providing permeability pathways for geothermal fluids. These faults are considered to be resurgent structures that might be related to shallow magmatic intrusions (Norini et al., 2019).
The topographic expression of the resurgence fault system is represented by rectilinear and curvilinear prominent fault scarps displacing by tens of meters Upper-Pleistocene-Holocene volcanic deposits in the center of the Los Potreros caldera (Figure 1a-c). The main fault swarm, represented by the Maxtaloya fault, the Los Humeros fault and some subparallel fault strands, is slightly curvilinear and runs NNW-SSE for ~8 km (Figure 1). Multiple N-S, NE-SW and E-W curvilinear splays depart from the main NNW-SSE fault swarm, depicting a complex deformation pattern (e.g. Arroyo Grande, Las Viboras and Las Papas faults, Figure 1a, Calcagno et al., 2018;Norini et al., 2019). All these faults invariably displace Holocene pyroclastic deposits, with vertical displacement rates of up to 10 mm/yr (Norini et al., 2019;Norini et al., 2015).
Reverse and normal displacements along these faults have been interpreted as the expression of changing stress field in the subsurface, as also documented by well logs and geophysical data (Carrasco-Núñez et al., 2017;Corbo-Camargo et al., 2020;Lermo et al., 2008;Lermo Samaniego et al., 2016;Lorenzo-Pulido, 2008;Norini et al., 2019;Norini et al., 2015;Rocha et al., 2010). The most probable origin of the local stress field has been suggested to be the pressurization of the shallow magmatic and/or hydrothermal system of the caldera complex (Montanari et al., 2017;Norini et al., 2019;Norini et al., 2015).
The Los Humeros fault exhibits a prominent fault scarp with a maximum elevation of 70-80 m above the caldera floor, showing a westward dip and dominant dip-slip normal displacement of a hydrothermally altered pumice fall deposit and underlying lava flows with a minor component of left lateral motion (Norini et al., 2019;Norini et al., 2015). Partly buried, eastward dipping, NNW-SSE reverse faults have also been identified 0.5-1 km to the west of the Los Humeros fault scarp by means of fieldwork, well logs, magnetotelluric sounding and seismic data (Corbo-Camargo et al., 2020;Lermo Samaniego et al., 2016;Norini et al., 2019) (Figure 1a (Figure 1c), although the subsurface geometry of these structures is rather uncertain, due to the limitations of the geological and geophysical datasets.
Changes in pressure of the shallow magmatic/hydrothermal system may have induced cyclic inversion of the dip-slip resurgence faults (cycling shifting between uplift and subsidence) (Norini et al., 2019). Also, the operation of the LHGF, with extraction and injection of fluids in the crust, may have induced ground motions superimposed on the long-term geological deformation signal (e.g. Békési et al., 2019a). Indeed, a shift between normal and reverse movements along the resurgence faults may have been induced in the shallow subsurface (down to about 1500 m of depth, i.e. the mean production depth, Gutiérrez-Negrín, 2019) ( Figure 1c) by changing pressure in the geothermal reservoir due to the natural evolution the hydrothermal/magmatic system and the extraction and injection of fluids, as discussed in the next sections.

InSAR processing
We used SAR images from ESA's Sentinel-1A and Sentinel-1B satellites. To map the coseismic deformation due to the February 8, 2016 earthquake, we used one ascending and one descending image pair (with acquisition dates of 29 January 2016 -10 February 2016 and 7 February 2016 -19 February 2016, respectively). We processed the individual interferograms using the GAMMA software (Werner et al., 2000). We applied a 30-m resolution Shuttle Radar Topography Mission Digital Elevation Model (USGS 1 ARC-second SRTM DEM, https://doi. org/10.5066/F7DF6PQS) to remove the topographic phase and obtain the differential interferograms. The interferogram pairs spanning the coseismic deformation showed unwrapping errors along the surface rupture due to the lack of coherence; therefore, we masked the surface rupture area of the interferograms for further modeling and interpretation (white area in the center of the fringe pattern Figure 3 a,d). Other InSAR data with larger wavelength (such as ALOS-2 observations) were also considered to overcome the decorrelation of Sentinel-1 interferograms, but there were no acquisitions mapping the studied coseismic deformation.
We further performed a time-series analysis of 20 Sentinel-1 images of descending track after the coseismic deformation, between 19 February 2016 and 16 May 2019. For the time series analysis we were only interested in the long-term deformation signal, therefore we picked Sentinel-1 images acquired on descending orbit every two months, making sure to avoid the co-seismic displacement due to the 8 February 2016 seismic event. We selected a single master image acquired on 13 February 2017 based on the minimization of the temporal and perpendicular baselines. We followed the Persistent Scatterer (PS) method using the MATLAB-based workflow of STAMPS (Hooper et al., 2012). PS candidates were selected based on their phase characteristics, and the number of PS were reduced in an iterative process through the estimation of phase noise. The wrapped phase of the final PS selection was corrected for spatially uncorrelated DEM error. After phase unwrapping, the spatially correlated error terms were estimated and removed from the interferograms (Hooper et al., 2007). The time series displacements were finally obtained from the corrected interferograms.

Methodology
The ground deformation due to the 8 February, 2016 earthquake suggests movements along a single rupture (Figure 3a, d). The ascending interferogram shows movements towards the satellite line of sight (LOS) at both sides of the fault, up to 42 mm and 87 mm in the eastern and western side of the rupture, respectively ( Figure 3a). The descending interferogram shows movements away from the satellite in the eastern block, with a maximum displacement of -105 mm, and movements of the western block towards the satellite with up to 40 mm ( Figure 3d). The focal mechanism solution based on seismological data by Lermo Samaniego et al. (2016) indicated reverse movement with a left lateral component at 1900 m depth, but the proposed mechanism does not resolve the measured coseismic surface deformation pattern (Békési et al., 2019b). We performed a Bayesian inversion with a significantly shallower single fault model and uniform slip (Model 1) on a rectangle (Okada, 1985) to reveal the source parameters of the seismic event.
Since the volcanic structures of the area commonly have a curved geometry , we also tested fault models with varying dip. We approximated the curved geometry with a segmented fault based on two rectangles with uniform slip in an elastic half-space (Model 2, Model 3).
To find the source parameters and uncertainties of the fault, we inverted the Sentinel-1 ascending and descending interferograms simultaneously using the Bayesian approach implemented in the Geodetic Bayesian Inversion Software (GBIS, Bagnardi and Hooper (2018)). The algorithm adopts the Markov Chain Monte Carlo method with automatic step size selection using the Metropolis-Hastings algorithm, to find the posterior probability density functions (PDFs) of the model parameters. Prior to the inversion, the InSAR data were subsampled using an adaptive quadtree sampling algorithm (Bagnardi and Hooper, 2018;Decriem et al., 2010), aiming to reduce the computational costs. Data errors required for the Bayesian inversion were simulated using an exponential function with nugget, fitted to the experimental semi-variograms calculated from the data (Bagnardi and Hooper, 2018). In each inversion procedure we performed, the PDFs were sampled through 10 6 iterations.
In case of the single-fault inversion, we used the already existing rectangular dislocation model with uniform slip (Okada, 1985) in GBIS. The 9 parameters in case of the single fault model include the location (X and Y coordinates and depth of the midpoint of the top edge), length, width, strike, dip, strike-slip, and dip-slip of the dislocation (Table 1). Subsequently, to approximate a curved fault geometry, we added a new forward model in GBIS. The model for the segmented fault consists of two rectangular dislocations connected to each other at depth, with the Table 1 Parameters, uncertainties (2.5% -97.5% confidence intervals), and RMS misfits of the single fault model (Model 1), the segmented fault models with identical slip (Model 2), and independent slip (Model 3), and the seismological solution of Lermo Samaniego et al. (2016) for the 8 February, 2016, Mw=4.2 Los Humeros earthquake. The uncertainties of the hypocenter locations by Lermo Samaniego et al. (2016) are reported to be <300 m in all 3 directions, however the exact values are not known. The coordinates of the faults are also reported in UTM zone 14N, and the rakes are calculated for comparison with the focal mechanism solution.   *In case of this study, it marks the coordinates of the center of the top edge of the rectangular plane. **Negative dip angles are due to convention signs in GBIS (dislocations dipping to the west have negative sign, Bagnardi and Hooper (2018)), and strike values are therefore oriented northwards. In the main text we converted negative dip values and corresponding strike values to match the focal mechanism solution for easier comparison. ***In case of this study, these are the UTM14N coordinates of the center of the top edge of the (upper) rectangular plane. ****Mw values of this study are calculated with the optimal fault parameters and a shear modulus of 8 GPa.
upper edge of the lower dislocation connected to the lower edge of the upper dislocation. In order to ensure the connection between the dislocations, we assumed that the two active fault segments have identical length and strike. We inverted the model with identical and different slip values for the upper and lower dislocations (Model 2 and Model 3, Table 1). Under these assumptions we had 11 model parameters with the identical and 13 model parameters with different dip-slip and strike-slip at the two segments. The model parameters include the location (X and Y coordinates and depth of the midpoint of the top edge), length, width, strike, dip, strike-slip, and dip-slip of the upper dislocation and the width, dip, strike-slip, and dip-slip of the lower dislocation. For all models we approximated the subsurface with an elastic halfspace with a Poisson's ratio ν = 0.25. This value is considered as a lower bound for the Poisson's ratio of the geological units of the area, calculated based on the mechanical rock property database of Los Humeros (Weydt et al., 2018;Weydt et al., 2021). We approximated the Poisson's ratio of the LHVC units (subdivided into a post-caldera unit, caldera unit, and pre-caldera unit) between 0.265-0.27, with and uncertainty of ±0.02. The Poisson's ratio of the underlying limestone basement was estimated to be 0.28±0.02. In order to account for the uncertainties of the Poisson's ratio and its effect on the resulting model parameters and their uncertainties, we tested the single fault model (Model 1) with Poisson's ratios of µ=0.25 and µ=0.3. We found that there is negligible influence of the variation of Poisson's ratio on the resulting model parameters, the uncertainties of the resulting model parameters largely overlap ( Figure S1a). Additionally, we found an elastic half-space a reasonable approximation, since the geological units do not have significantly different Poisson's ratios. Vp/Vs ratios based on seismic tomography  near the activated fault segment also do not show any significant variations within the LHVC sediments, supporting the applicability of an elastic half-space for this study.
We set up constraints based on the 3D geological model of Los Humeros as constructed through the use of geological and geophysical observations, and partly on the focal mechanism solution from seismological data (Lermo Samaniego et al., 2016) to find reasonable model parameters (Table 1). We selected lower and upper bounds of the model parameters in order to account for a wide range of possible solutions (Table 1) (Table 1). In case of the single fault model, we inverted for all the 9 model parameters simultaneously. In case of the segmented fault models, we followed the methodology of Funning et al. (2005) for two-fault models. We fixed one fault segment and inverted for the other one iteratively. We performed a parameter search for the location, strike and dip of the faults based on prior information, keeping the rest of the parameters of both faults fixed to obtain a first model reproducing the geometry of the observed ground deformation pattern. Then, we updated the model bounds and inverted for all the parameters.

Model 1 -single fault
The upper and lower bounds and inversion results of parameters of the single fault model (Model 1) are listed in Table 1. Contrary to the focal depth of 1900 m derived from seismological data (Lermo Samaniego et al., 2016), the observed surface movements could only be reproduced by a fault activation at a shallower depth, extending from the surface down to ~1000 m depth (Figure 4a). The inferred fault has reverse kinematics with a minor left lateral component (strike=161 • ). A good fit with both the ascending and descending InSAR data was achieved; the model mismatch is quantified in a root-mean-square (RMS) misfit of the order of 3 mm (Table 1).

Model 2 -segmented fault, identical slip
The segmented fault model approximating the curved geometry of the inferred fault structure was first tested with identical slip values for the upper and lower dislocation (Model 2). Model 2 shows a slightly improved fit with the InSAR data compared to Model 1 ( Table 1). The fault kinematics in case of Model 2 agree with the single fault solution. The width of the upper dislocation is significantly smaller than the lower one (208 m and 852 m, respectively, Figure 4b, Table 1), and the active fault segments extend from the surface down to ~1000 m depth. The dip angle increases with depth, with 57 • in case of the upper dislocation, and 67 • for the lower dislocation.

Model 3 -segmented fault, independent slip
We also performed the inversion leaving the dip-slip and strike-slip component of the lower fault segment independent from the values of the upper segment. This configuration allows us to mimic the change in the slip components with varying dip angles and depth. We set up the bounds of the strike-slip components of the two dislocations in order to allow only left lateral movements. It was important to avoid inconsistent lateral movements along the two connected dislocations. The parameters and uncertainties of the best-fitting model (referred to as Model 3) are presented in Figure 3 and Table 1. Model 3 has mostly similar parameters as Model 2, although only the lower fault is having left lateral kinematics ( Figure 4c, Table 1). Another difference with the results from Model 2 is that the two fault segments have similar widths (480 m and 442 m, Figure 4c, Table 1). The misfit of Model 3 has slightly improved compared to the previous models ( ). Due to the similar rake values to the seismological solution and the lowest overall misfit, Model 3 represents the most realistic source parameters solution of the 8 February 2016 earthquake.
The strike-slip components of the dislocations are marked by relatively large uncertainties (up to a variation of 30% -37% within 2.5% -97.5% confidence intervals) compared to the dip-slip components (up to a variation of 4% -16% ) in case of all models (Table 1, Figure S1). This can be explained by the fact that the N-S oriented movements, and therefore the strike-slip component, are less constrained by the InSAR data due to the satellite view geometries.

Deformation between February 2016-May 2019
The velocity map of the study area covering the LHVC with average yearly LOS movement rates after the coseismic deformation, between 19 February 2016 and 26 May 2019 is presented in Figure 5. Movements away from the satellite (~subsidence) are concentrated within the Los Potreros caldera with up to 8 mm/year near the inner-caldera structures (Figure 5 a,b). The 8 mm/year rate is in agreement with previous studies based on Envisat data between 2003 and 2007 (Békési et al., 2019a;Cigna et al., 2019). The Los Humeros fault in the north, and the Maxtaloya fault in the south appear to constrain the deforming area in the west. Subsidence is in general consistent with the location of the active production wells (Figure 5a).
An inconsistency between deformation and geothermal production is observed near the Las Papas fault. Subsidence of this area with comparable rates was not observed previously by Békési et al. (2019a), but it is possible that they mapped the signal within the topographic error term and removed it from the interferograms. Results of this study show that the center of the Los Potreros caldera floor experiences deformation up to 7 mm/year (Figure 5b Figure 5a). Furthermore, the observed subsidence is consistent with the previously studied coseismic motion of the eastern fault block (Figure 3d). Afterslip on the seismogenic fault and the viscoelastic relaxation after the main seismic event must contribute to the deformation signal. However, additional long-term volcano-tectonic activity may be partly responsible for the observed subsidence.

Discussion
We modeled the coseismic deformation due to the 8 February 2016, Mw=4.2 earthquake at Los Humeros with the effect of a single dislocation (Model 1), and also with two connected rectangular dislocations representing two fault segments to approximate a curved fault geometry (Model 2, Model 3). We jointly inverted the ascending and descending interferograms using a Bayesian approach in order to constrain the fault parameters and associated uncertainties. We based our models on prior information on the geometry of the mapped and inferred fault structures after Calcagno et al. (2018) and Norini et al. (2019). Our models are able to resolve the observed ground deformation pattern and suggest fault reactivation at very shallow depth (<1 km). Additionally, our current models are in agreement with the surface trace of the mapped fault structures (Figure 3).
The observed surface movements together with the modeling results indicate a reverse kinematics with a left lateral component (Table 1), interpreted as the reactivation of the Los Humeros fault (Figure 3, 4). The reverse movement and the geometry of the Los Humeros fault are in good agreement with the focal mechanism solutions of Lermo Samaniego et al. (2016) Table 1).    is shown with blue (the shading indicates the uncertainty of reservoir extent). The trace of the section is shown in Figure 3b.
However, Lermo Samaniego et al. (2016) observed a significant left lateral component. We tested the segmented fault models with identical slip of the upper and lower dislocation (Model 2) and we also performed the inversion with different slip values (Model 3). A significant left lateral movement (rake=58 • , Table 1) could only be reproduced by Model 3, and is only attributed to the lower fault segment. The inconsistency between the strike-slip components of the geodetic models and the seismological solution can partly originate from the uncertainty of the N-S oriented movements due to the satellite orbits. To better constrain the strike-slip kinematics through the InSAR inversion, an independent observation with a different view geometry would be required. Still, the strike-slip components necessary to reproduce the rake values of the moment tensor solution (42 • ) fell outside the confidence intervals of the fault models (Table 1, Figure S1). The mismatch in the magnitude of strike-slip component could further be explained by the uncertainty assigned to the focal mechanism inversion of Lermo Samaniego et al. (2016), limited to the polarities of the P-wave arrivals at the vertical components of five seismic stations. Additionally, errors in the moment tensor solution may occur due to uncertainties in the velocity model adapted by Lermo Samaniego et al. (2016) (e.g. Vasyura-Bathke et al., 2021. The inferred independent slip values of the upper and lower dislocation (Model 3) could also suggest that the kinematics of the Los Humeros fault varies with depth and dip angle.
There is no specific indication of the hypocenter location error of the 8 February 2016 earthquake by Lermo Samaniego et al. (2016), however they report a maximum error of 300 m in all directions (longitude, latitude and depth) for a subset of seismic events. The error of up to 300 m can be considered relatively low, but still the uncertainties of the horizontal coordinates inferred from the seismological model largely overlap with the results of this study (Table 1).
The InSAR data can be fitted with the reactivation of the Los Humeros fault extending from the surface down to ~800-1000 m depth, but no displacements can be inferred at larger depths. One explanation of the depth extent of the reactivation is that the Los Humeros fault is connected to a reverse fault at this depth, indicated with the black dashed line in Figure 4. While this fault was not identified by several studies (Carrasco-Núñez et al., 2017;Ferriz, 1982;Norini et al., 2015) the existence of the buried reverse fault structure is suggested by e.g. Lermo et al. (2008) and Norini et al. (2019). The occurrence of the earthquake at a maximum depth of 1 km would require the presence of stress-resistant material in shallow depth. Seismic tomography results show a significant positive Vp anomaly near the seismic event, delineating a relatively rigid body between ~500-1500 m depth , that could build up stresses for a Mw=4.2-4.35 rupture. Toledo et al. (2020) estimated the top of the pre-caldera group (hosting the geothermal reservoir) to ~600 m depth in the vicinity of the earthquake, which is significantly shallower than the average reservoir depth of approximately 1000 m. This would suggest that the reactivated fault segment inferred from our study extends also within the upper part of the reservoir (Figure 4). The very shallow fault reactivation inferred from the InSAR inversion could suggest that the earthquake hypocenter is located near the bottom of the rupture that is often observed (e.g. Gaudreau et al., 2019;Karasözen et al., 2016). Compared to our models, the earthquake was originated significantly deeper (1900 m) according to the seismological solution (Table 1, Figure 4). The uncertainty of the hypocentral depth based on the seismological model is estimated to be as low as 300 m (Lermo et al., 2008), and the depth of the fault plane is estimated to be up to 1100 m (within 97.5 % confidence interval in case of Model 2, Table 1) according to the geodetic models, showing that the uncertainties of the two models do not overlap. Such a large difference in source depth of ~900 m cannot be explained by the limitations of our models, for instance the assumption of an elastic half-space. Lermo Samaniego et al. (2016) used a 1D P-wave velocity model (Lermo et al., 2008) and first arrivals from five seismic stations. The velocity model of Lermo et al. (2008) adapts significantly different values from from the results of more recent studies, especially at shallow depth (Löer et al., 2020;Toledo et al., 2020). Also, due to the largely heterogenous subsurface of the study area, the hypocentral depth might be biased as a result of differences in seismic velocities from the 1D model. All these limitations of the seismological model may explain the differences in estimated focal depths.
The calculated moment magnitude from all three slip models based on the InSAR data is approximately 4.35, using a shear modulus of 8 GPa (Table 1). We estimated the shear modulus of the volcano-sedimentary units of the LHVC based on the mechanical rock property database of Los Humeros (Weydt et al., 2018;Weydt et al., 2021). The estimated moment magnitude of 4.35 is larger than the Mw=4.2 value reported by Lermo Samaniego et al. (2016), corresponding to an approximately 1.7 times larger seismic moment (3.8 x 10 15 Nm in case of model 3, compared to 2.2 x 10 15 Nm of the seismological solution). Using a shear modulus of 5 GPa would result in a modeled moment magnitude of 4.2, although such value for the shear modulus would underestimate the actual properties. Still, the uncertainties of the rock mechanical properties of different rock types can be as large as 4 GPa (Weydt et al., 2018), leading to relatively high uncertainties (±0.15) of the estimated moment magnitudes. Weston et al. (2012) found that there is a slight trend for an overestimation of the moment magnitude for thrust events studied using InSAR. This discrepancy can be explained by modeling errors due to the absence of InSAR observations in the near-fault zone due to the lack of coherence and/or the limitations of the seismological solution, such as the influence of a curved fault geometry on the estimation of moment magnitude (e.g. Contreras-Arratia and Neuberg, 2020). Another factor is the mapping of potential deformation from aftershocks, afterslip and viscoelastic relaxation in the interferograms. Additionally, we approximate the subsurface with an elastic half-space instead of a layered media. This can also contribute to the overestimation of the moment magnitude, especially due to a shallower source depth compared to the seismological solution (Weston et al., 2012).
Our model results show remaining misfits with both the ascending and descending data (Figure 3 c, f). These can be explained by the model assumptions that we applied. The fault surfaces have significant irregularities, and the interferograms suggest movements along several fault structures not present in our models (Figure 3 a, d). For instance, both the ascending and descending data show misfits with the modeling results northeast to the surface trace of the faults. This is explained by the unmapped dislocation along the NE-SW striking Arrayo Grande fault. Also, the termination of the Los Humeros fault and its potential connection to the sub-parallel Maxtolya fault in the south might be responsible for a more complex slip distribution. The dip angle of the faults in the area is also likely to change with depth, as suggested by both geological and geophysical data (e.g. Carrasco-Núñez et al., 2018;Norini and Groppelli, 2020;Urbani et al., 2020).
Over geological time scales, the Los Humeros fault with normal kinematics has induced the lowering of the western fault block, creating the present-day topography. However, active ground movements based on the InSAR observations indicate the reverse reactivation of the Los Humeros structure. A possible source of subsidence of the eastern fault block observed on the coseismic interferogram and prolonged surface motions (Figure 3d, Figure 5a, b) is the depressurization of the hydrothermal system. This pressure source could also explain the local radial stress field observed by Norini et al. (2019). Connection between the shallow fault reactivation and the deep pressure source is ensured by the buried reverse structure inferred from Norini et al. (2019) and this study (Figure 4). Norini et al. (2019) and Urbani et al. (2020) suggested the resurgence of the Los Potreros caldera floor east of the Los Humeros fault, that is supported by the presence of intrusive-like bodies on seismic tomography images . Although the current subsidence of the caldera floor indicate depressurization at depth (Figure 5a, b), resurgence of the caldera floor can still exist in a larger time scale. In such a case, recent subsidence might be related to an episodic pause in the overall resurgence process.
The correlation between the 08 February 2016 earthquake and the increase of the injection rate in the H-29 well is obvious ( Figure 2). Since the earthquake occurred seven days after the onset of injection, it is likely that poro-elastic stress changes due to fluid injection are partly responsible for the triggering of the event. Additionally, the depressurization of the hydrothermal system due to production from several other geothermal wells may also contributed to the occurrence of the earthquake. However, it is difficult to distinguish between the different factors that are responsible for the observed ground motions, also due to the lack of interferometric coherence observed in the near-fault zone.
Coseismic deformation mapped with radar satellite data provides independent observations to constrain the location, geometry, and kinematics of the seismic event. Our study demonstrates that ground movements mapped with InSAR and guided by other geological and geophysical observations can significantly contribute to our understanding of the fault structures and processes that control the geothermal reservoir.

Conclusions
We reproduced the coseismic surface deformation of the 8 February 2016, Mw=4.2 earthquake at Los Humeros, with a two-component rectangular dislocation model approximating the curved geometry of the fault structures of the LHVC. We constrained the model parameters and the underlying uncertainties with a Bayesian inversion. The results imply fault reactivation at very shallow depth, down to 1 km.
Modeling results suggest that fault kinematics changes with depth, with a purely reverse kinematics of the upper dislocation, and significant left lateral component of the lower dislocation (rake=58 • ). Models reproduce the observed ground movements with high precision, and fits with the surface trace of the Los Humeros structure, that was previously interpreted as a normal fault. The reactivated Los Humeros fault in the north and the Maxtaloya fault in the south and their potential continuation to a buried fault structure in larger depth act as major boundaries of the subsiding area observed after the coseismic deformation, for the period of February 2016 and May 2019. A potential explanation of the reverse reactivation of the Los Humeros fault and corresponding downward movement of the eastern fault block can be the depressurization of the whole hydrothermal system. Such depressurization might occur due to the exploitation of the geothermal field and/or due to natural pressure/temperature changes related to magmatic activity.
Our study demonstrates the added value of employing InSAR monitoring and inverse modeling under geological and geophysical constraints. It can help improve the understanding of complex processes of geothermal sites. Such information is highly valuable for future operations of the field and can contribute to our understanding of SHGS in general.

Data Availability
Sentinel-1 data are accessible through The Copernicus Open Access Hub (https://scihub.copernicus.eu) after registration. The open-source Geodetic Bayesian Inversion Software (GBIS) used to construct the models for the co-seismic deformation is available at https://comet. nerc.ac.uk/gbis/.

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.