Ground motions induced by pore pressure changes at the Szentes geothermal area, SE Hungary

Excessive thermal water volumes have been extracted from porous sedimentary rocks in the Hungarian part of the Pannonian Basin. Thermal water production in Hungary increased significantly from the early 1970s. Regional-scale exploitation of geothermal reservoirs without re-injection resulted in basin-scale pressure drop in the Upper Pannonian (Upper Miocene) sediments, leading to compaction. This compaction resulted in ground subsidence primarily through poro-elastic coupling. We investigated surface deformation at the Szentes geothermal filed, SE Hungary, where the largest pressure decline occurred. Subsequently, hydraulic head recovery in the western part of the geothermal reservoir was initiated in the mid-1990s. We obtained data from the European Space Agency’s Envisat satellites to estimate the ground motions for the period of November 2002–December 2006. We applied inverse geomechanical modeling to estimate reservoir properties and processes. We constrained the model parameters using the Ensemble Smoother with Multiple Data Assimilation, which allowed us to incorporate large amounts of surface movement observations in a computationally efficient way. Ground movements together with the modeling results show that uplift of the Szentes geothermal field occurred during the observation period. Since no injection wells were operated at Szentes before 2018, and production temperatures remained relatively constant through the entire production period, we explain ground uplift with pore pressure increase due to natural recharge. The estimated decompaction coefficients of the reservoir system characterizing the elastic behavior of the Szentes geothermal reservoir varies between ~ 0.2 × 10–9 and 2 × 10–9 Pa−1. Compaction coefficients of the reservoir system corresponding to the earlier depressurization period, from ~ 1970 to the mid-1990s, may be significantly larger due to the potential inelastic behavior and permanent compaction of clay-rich aquitards. The improved parametrization enables better forecasting of the reservoir behavior and facilitates the assessment of future subsidence scenarios that are helpful for the establishment of a sustainable production scheme.


Introduction
Hungary is among the most favorable countries for geothermal development within Europe (e.g., Békési et al. 2018;Cloetingh et al. 2010;Horváth et al. 2015;Lenkey et al. 2021;Limberger et al. 2014Limberger et al. , 2018, with an average geothermal gradient of ~ 45 °C/ km, and a mean surface heat flow of 100 mW/m 2 (Lenkey et al. 2002). The highest values occur in the SE part of the country, corresponding to the thinnest parts of the crust and lithosphere (e.g., Horváth et al. 2006). The outstanding geothermal conditions of Hungary originate from Miocene extension of the Pannonian Basin, related to subduction and roll-back in the Eastern Carpathians, resulting in the thermal attenuation of the lithosphere (e.g., Horváth 1993;Horváth et al. 2006). The different timing of syn-rift sediment infill in the sub-basins indicates that extension in the Pannonian Basin migrated in space and time (e.g., Fig. 1b) (Balázs et al. 2016). Subsequent to the syn-rift phase, thermal subsidence of the Pannonian Basin initiated, accompanied by continuous post-rift sedimentation (e.g., Juhász 1991;Sztanó et al. 2013). From Late Miocene to recent times, basin inversion has occurred due to a push from the Adriatic microplate towards the Pannonian Basin (Bada et al. 2007;Horváth and Cloetingh 1996), forming its present-day basement geometry (Fig. 1a, b).  Haas et al. (2014) and modified after Békési et al. (2018). b Interpreted composite reflection seismic transect from the eastern part of Hungary, with the main tectonic and stratigraphic features of the area, after Balázs et al. (2016) Thermal water production in Hungary has already started in the nineteenth century for bathing purposes, and the country is still famous for its large amounts of hot springs and thermal spas. Geothermal wells most commonly target Upper Miocene and Quaternary porous sedimentary reservoirs (Fig. 1b). Fractured Mesozoic carbonate rocks also have significant geothermal potential, and are being utilized in several locations within Hungary (e.g., Goldscheider et al. 2010;Horváth et al. 2015;Mádl-Szőnyi et al. 2015;Szanyi and Kovács 2010). Additionally, fractured Mesozoic crystalline basement rock are also suitable targets for deep geothermal developments (e.g., Békési et al. 2018;Horváth et al. 2015;Vass et al. 2018), however, there has been no active exploitation of such resources yet. Until 2018, geothermal energy was exclusively used for direct heat purposes in Hungary. Since then, the first geothermal power plant has become operational in Tura, with an installed capacity of 3.35 MW e (Nádor et al. 2019). Compared to its vast potential, Hungary is still lagging behind in geothermal energy production. Further geothermal developments are crucial to increase the role of renewable energy resources within Hungary's (and Europe's) energy demand. However, the long-term sustainable production of geothermal energy requires cautious planning and regulation. Exploitation in excess of natural recharge can result in reservoir pressure decline, causing a decrease in production rates. Furthermore, such "overexploitation" of geothermal reservoirs may lead to compaction, land subsidence, or even induced seismicity (e.g., Allis 2000;Békési et al. 2019;Keiding et al. 2010;Maghsoudi et al. 2018;Trugman et al. 2014;Meer et al. 2014).
Previous studies have demonstrated the effectiveness of interferometric synthetic aperture radar (InSAR) based monitoring combined with (inverse) modeling of geothermal fields (Heimlich et al. 2015;Keiding et al. 2010;Parks et al. 2020;Trugman et al. 2014;Vasco et al. 2013Vasco et al. , 2002. In case of geothermal applications, inverse models constrained with ground motion data are commonly based on single or distributed volume/ pressure sources in the subsurface (e.g., Kiyoo 1958;Okada 1985;Segall 2010;Yang et al. 1988). Ground deformation due to seismic events occurring at geothermal areas has also be modeled using analytical methods in order to estimate source parameters (e.g., Békési et al. 2021). More complex numerical models have also been applied to predict ground motions at geothermal sites. For instance, Vasco et al. (2013) applied coupled numerical modeling to resolve the observed ground deformation at the Geysers Geothermal Field. However, coupled models with large number of model parameters are not always suitable for inversion exercises incorporating extensive calibration datasets such as detailed surface movements as they require abundant runtime. The advantages of probabilistic ensemble-based approaches for inversion of surface subsidence have already been shown in several studies outside the geothermal energy arena (Baù et al. 2015;Candela et al. 2017;Fokker et al. 2016). Such ensemble-based approaches have limitations in terms of model parameters, including only a limited number of reservoir properties, and cannot directly distinguish between different subsurface processes. Still, ensemble-based techniques are capable of dealing with large amounts of measurements and they can successfully estimate driving parameters of subsurface processes. A recently developed ensemble-based subsidence interpretation and prediction tool (ESIP) is capable to distinguish between different compaction behaviors, and allows to predict a wide range of subsurface parameters (Candela et al. 2022). However, ESIP requires detailed information on the local geology and pore pressure history of a study area for the whole operation period, that are not always available. In this study, we combine the Ensemble Smoother with Multiple Data Assimilation [ES-MDA, Emerick and Reynolds (2013)] with the analytical solution of Geertsma (1973), in order to estimate first-order properties of the Szentes geothermal area, without the requirement of high-resolution subsurface measurements, such as pressure time series. This study is the first application of mapping and modeling ground motions due to geothermal activities in Hungary, that aims to demonstrate the usefulness, requirements, and limitations of inverse geomechanical models for geothermal sites in Hungary and worldwide.

The Szentes geothermal field
The Szentes geothermal field is located in SE Hungary, in the northern part of the Makó trough ( Fig. 1a). Upper Miocene delta front and delta plain sediments (Upper Pannonian Újfalu Formation) and delta top and alluvial plain sediments (Upper Pannonian Zagyva Formation) host the geothermal reservoir, deposited in the confines of a progradational delta system (e.g., Fig. 1b) (Juhász 1991;Sztanó et al. 2013). The targeted Upper Pannonian formations are largely heterogeneous; built up by the alternation of sandstone, silty sand, silty clay and clay layers (e.g., Juhász 1991). Geothermal production primarily targets the Újfalu formation, with a mean thickness of approximately 600 m in the vicinity of the geothermal wells, and a mean top depth of ~ 1500 m inferred from well logs. The individual sandstone layers have small lateral extent, therefore, it is difficult to correlate the sand bodies drilled in different wells. Still, hydraulic connection between the individual sand bodies exist. Natural pore pressure in the Upper Pannonian formations of the Szentes area is slightly above hydrostatic with a gradient of 10.2 MPa/km, however, extreme overpressure in the underlying Lower Pannonian strata (71 MPa/km) is observed (Bálint and Szanyi 2015). Such high overpressure in the low-permeability Lower Pannonian sediments is present in several locations of the Great Hungarian Plain, and most likely explained by tectonic processes (Almasi 2002;Tóth and Almási 2001;Balen and Cloetingh 1995).
Szentes is the first area where geothermal energy has been produced for industrial application in Hungary, with the first well drilled in 1958 (Bálint and Szanyi 2015;Szanyi and Kovács 2010). There are 45 wells drilled in total (including the Szegvár and Fábiánsebestyén areas), and 32 of them were producing thermal water up to 90 °C in 2010. Starting from the early 1970s, thermal water production reached ~ 6.5 million m 3 / year, which has decreased to an average of 5.7 million m 3 /year from the early 1990s (Szanyi and Kovács 2010). Yearly extracted thermal water volumes remained approximately constant after the early 1990s, although the exact yearly amounts cannot be precisely estimated due to the different owners of the individual wells, where production data are not always documented or provided. Additionally, geothermal wells have not been continuously operated; production during the summer period has been terminated in several wells. Yearly thermal water withdrawal has also slightly varied depending on winter temperatures. The complete lack of reinjection wells resulted in significant reservoir pressure decline and hydraulic head decrease with up to 36 m (Figs. 2a, b and 3). Reservoir pressures have started to increase since the mid-1990s due to the decrease in production rates in the western part of the geothermal filed, accompanied by 4-8 m Regular monitoring of the Szentes geothermal wells from the beginning of thermal water withdrawal has not taken place, therefore, it is difficult to fully understand reservoir behavior through time. There are multiple factors responsible for the lack of regular measurements, for instance the different owners and operators of geothermal wells, inconsistencies in regulations of thermal water production in Hungary. We gathered downhole pressure measurements from the Szentes area ( Fig. 3), showing that in some cases 2-3 decades passed between measurements. During 2009-2010, an investigation of 20 geothermal wells took place in the framework of the Jedlik project (Bálint and Szanyi 2015), including reliable downhole pressure measurements, showing that pore pressures were still below natural conditions in 2009-2010 (Fig. 3). Pore fluid pressures in 2010 were up to ~ 0.6 MPa lower than in 1980 (Fig. 3, yellow circles connected with a solid line), and depletion has probably reached its maximum in the early 1990s. Bálint  Bálint and Szanyi (2015) and Szanyi and Kovács (2010). Please note that the original undisturbed water levels approximately agree with the estimated values for 1970. The largest decrease of hydraulic heads in 2000 occurred in the eastern part of the geothermal field, whereas hydraulic recovery of the western area has already started in the mid-1990s. The locations of well K-652 is labeled in a with white, where hydraulic head time series and InSAR-derived ground motions are plotted in Fig. 5 and Szanyi (2015) studied the possible effect of reinjection towards a more sustainable field operation, by investigating the most suitable location and depth for an injection well. Since then reinjection in a single well has initiated, while production using downhole pumps is still ongoing.

InSAR data
We used C-band radar acquisitions of the European Space Agency's (ESA) Envisat satellites in order to map the ground motions for the period of November 2002-December 2006, summarized in Table 1. We assessed the availability and quality of radar images acquired on both ascending and descending satellite orbits. We introduced a criterion on the perpendicular baselines, and we discarded scenes with perpendicular baselines over 600 m, temporal baselines above 1 year and interferograms with no visual interferometric coherence. This criterion has reduced the number of scenes in case of both ascending and descending geometries, and we had to exclude the ascending images from further analysis due to the low number of scenes that we found insufficient for time series analysis.
To create the individual interferograms on the descending geometries, we used the ESA SNAP software (S1TBX-ESA Sentinel-1 Toolbox, http:// step. esa. int). We applied orbital corrections on the individual scenes using DORIS Precise orbits. We selected a single master image using a criterion of minimizing the temporal and perpendicular baselines and we co-registered the images to the master geometry. We used the

1-arc-second (30 m) resolution Shuttle Radar Topography Mission Digital Elevation
Model (SRTM DEM) to compute the topographic phase that we subsequently removed from the interferograms. We obtained the time series of ground motions by Persistent Scatterer Interferometry (PSI), within the Matlab-based workflow of STAMPS (Hooper et al. 2007(Hooper et al. , 2012. In STAMPS, PS pixels are initially selected based on their amplitude dispersion, initially described by Ferretti et al. (2001). We selected a threshold of 0.4 for the amplitude dispersion index, allowing for a large number of pixels to be included in further analysis. Later, the number of PS was reduced during an iterative process based on phase stability. Interferograms were then corrected for spatially correlated and uncorrelated error terms as described by Hooper et al. (2012).
The PSI results indicate dominantly positive movements (~ uplift) at the Szentes geothermal field, with a deformation rate of ~ 5 mm/year (Fig. 4, outlined with a black rectangle) for the period of November 2002-December 2006. The largest positive LOS movements (~ uplift) is observed outside the geothermal field, in the north, with velocities up to 15 mm/year. Minor negative movements are observed south of the field, near the geothermal wells of Szegvár, and over the eastern and northeastern part of the Szentes field.
We plotted groundwater levels from well K-652 and LOS velocities of PS pixels within 500 m distance from the well (Fig. 5). Hydraulic head measurements from the Szentes area are rather limited; well K-652 was chosen because it has the most frequent observation history, and covers the duration of the InSAR monitoring. Hydraulic head data in the . Arrows indicate the flight direction of the satellite and the look direction with the corresponding incidence angles. The locations of the geothermal wells are plotted with black triangles. The modeling area that corresponds to the Szentes geothermal filed, excluding the Szegvár thermal wells that are located south of the field, is outlined with black rectangle K-652 well, located in the southern part of the geothermal field, show a decreasing trend until 1995, when the recovery of water levels started (Fig. 5a). The InSAR time series can be approximated with a linear trend (Fig. 5b), with 2 major negative and positive deviation in ~ January 2005 and August 2005. These anomalies might be explained by the fact that thermal water production was terminated for the summer period, allowing for the shortterm recovery of groundwater levels. In the vicinity of well K-652, the InSAR measurements indicate positive LOS motions (~ uplift) of the area, although groundwater levels do not show any significant variation (only a minor increase of 0.5 m) between 2002 and 2007 (Fig. 5a). It is important to note that InSAR observations are averaged within 500 m distance of the geothermal well and the open sections of the geothermal well only partly cover the vertical extent of the geothermal reservoir, since the aquifer units are separated with clay-rich aquitards with variable thickness and depth. Formation pressures might be different away from the boreholes, as the initially artesian wells required pumping after several years of production, when water levels declined and no longer reached the ground surface. If pressure changes are restricted to certain parts of the reservoir, the measured groundwater levels may not reflect these changes that might explain the discrepancies between the observed groundwater levels and LOS movements.

Inverse geomechanical modeling
We modeled ground movements using the analytical solution of Geertsma (1973), for a disk-shaped reservoir undergoing a uniform pore pressure increase/decline, embedded in a homogeneous elastic half-space, with a Poisson's ratio ν . In the Geertsma model, the vertical (u z ) and radial (u r ) components of movements at the surface due to pressure change, p ., attributed to a disk-shaped source with radius, R, and height, h, compaction coefficient or compressibility, c m , and burial depth, D, can be obtained as (Geertsma 1973): where e is the Euler's number, J 0 and J 1 are first-type Bessel functions of order 0 and 1, respectively, and r is the distance measured from the center of the source. The formula can be simplified by substituting c m h p = h , where h is the height change of the reservoir. This simple model cannot take into account complex geological settings and variations in subsurface geomechanical and petrophysical properties, but it can be efficient for a first approximation of ground deformation induced by injection/production activities. Due to its simplicity, the driving parameters of the Geertsma model can be estimated in a computationally inexpensive way by inversion techniques. Estimates of reservoir parameters and related uncertainties, for instance the compressibility, is highly relevant for more complex 3D reservoir models, to increase their reliability by choosing reasonable input parameters.
For the modeling of ground deformation, we selected an area of approximately 20 × 15 km outlined in Fig. 4, based on the observed surface movements and well locations at Szentes. We excluded the southernmost geothermal wells of Szegvár and the location of a major positive ground motion anomaly outside the geothermal field in the north (Fig. 4) from the modeling. It was necessary since the modeled ground motion pattern would have been too complex to capture the major ground motion anomalies that occur in the vicinity of the Szentes geothermal wells with a simple model. We primarily concentrated on the major positive LOS anomaly identified on descending satellite orbits for the period of November 2002 and December 2006 (Fig. 4). We first tested the model with one single Geertsma source with uniform pressure increase (source 1), selecting broad prior ranges of model parameters. We assumed a uniform pressure increase based on the InSAR observations and the hydraulic head profile (Fig. 2b), indicating increasing groundwater levels at the bottom of the geothermal reservoir (base of the Újfalu formation). Having only one source, an unmapped negative LOS anomaly remained ~ east of the main uplift signal (Fig. 6b, e). To account for this subsidence, we included an additional small source with negative pressure change (source 2) to the model (Fig. 6a, c) based on the InSAR measurements. (1) To estimate the model parameters and uncertainties, we applied ensemble-based probabilistic inversion, having the ascending and descending InSAR datasets as target observations. The Ensemble Smoother [ES, Emerick and Reynolds (2013)] estimates the model parameters by a global update, incorporating all data available. Therefore, inverse problems with large number of observations can be solved by the ES in a computationally efficient way. For non-linear forward models, the ES requires several iterations, where the predictions of the previous run is used as an input for the subsequent data assimilation step (ES-MDA, Emerick and Reynolds 2013). In case of the forward model based on the analytical solution of Geertsma (1973), only a weak nonlinearity between ground motion and certain model parameters (radius and depth of the disk-shaped sources) exist. Still, we choose ES-MDA as we expected a better performance than a single ES and to be able to evaluate the improvement of the solution with each DA steps.
The solution for a single data assimilation for the updated model ensemble is: In Eq. 3, M is the prior ensemble of model parameters, GM represents the result of the forward model working on all ensemble members, and GM′ is the difference between GM and its mean. N e is the number of ensembles, and D is an ensemble of data realizations, created by perturbing the measurements according to their covariance matrix ( C d ) . The mean of the ensemble is taken as the best estimate, which is used as input for the next update in case of ES-MDA. The number of data assimilation steps, N a must be selected a priori. The data covariances used for the update steps are increased by a multiplication factor, α i for i = 1, 2…, N a , and α i must be selected as N a i=1 1 α i = 1 (Emerick and Reynolds 2013). This is necessary to compensate for the effect of multiple application of an ES.
For the inverse modeling with ES-MDA, we fixed the depth of the sources to 1500 m, which was identified as the mean depth of the top of the reservoir from well logs and seismic data. We approximated the subsurface with an elastic half-space, with a Poisson's ratio ν = 0.28, that we selected based on average values for the area inferred from well logs. The model parameters include the location (X and Y), radius (R) and height change  ( h ) of the two disk-shaped sources, 8 parameters in total. We assigned prior values for the parameters assuming normal distribution, except for the height change of both sources, where a lognormal distribution was chosen (Table 2). We selected prior values for the height change to be positive in case of source 1, accounting for a positive pressure change at depth attributed to the major positive LOS displacement anomaly. The height change of source 2 was assumed to be negative, approximating pressure decline and corresponding subsidence. Measurement errors of the InSAR observations were set to 4 mm, after testing the influence of different values (2-8 mm) on the inversion procedure. We first included only one source in the parameter estimation procedure allowing for broad prior values, to get a first-pass model for the major positive pressure source. Using the mean posterior model parameters of the one source model, we estimated the parameters of the second source combined with the first source having fixed parameters. We iteratively fixed one source at a time, until no significant improvement compared to the previous model was achieved. Then, we assigned relatively narrow ranges to the prior parameters, with mean and standard deviation reported in Table 2. Having these prior statistics, we estimated all model parameters of the two sources through the same ES-MDA run, using all InSAR observations. The data assimilation procedure was performed with 100 ensemble members and 4 iterations or data assimilation steps.
Modeling results indicate the presence of an extensive positive pressure source, with a radius of 8755 ± 38 m (source 1, red dashed circle in Fig. 6a), and a smaller negative source with a radius of 1269 ± 34 m (source 2, blue dashed circle in Fig. 6a). Source 1 is located at the center of the Szentes geothermal field, covering the locations of most of the geothermal wells. The location of source 2 corresponds to the vicinity of production wells in the eastern part of the field. Source 2 is positioned almost entirely inside the large positive source (Fig. 6a), suggesting that source 2 is superimposed on top of the significantly larger source 1.
Radius of source 1 (m) Height change of source 1 (m) The model can adequately approximate the major structures of the ground motion observations (Fig. 6). Still, residual surface movements (Fig. 6c) are significant at certain areas, suggesting that the simple model could not entirely reproduce the InSAR time series displacements. The overall RMSE of the model is 7.5 mm. The average observed displacements in the model area is 20.3 mm with a maximum of 59 mm, and only 9% of the observations show movements > 30 mm. The average modeled displacement is 20.38, which is almost identical with the observed one, while the maximum modeled movement is 23.4 mm, which is significantly smaller than the observed maximum value. Local ground motion anomalies with large displacements (> 30 mm) can explain the relatively large RMSE, showing that the complexity of the subsurface cannot be captured by the model. Still, the motivation for such simple model was to have an estimate on the sources that are responsible for the major pattern of ground movements and to find reasonable intervals for compaction coefficients that are mapped within h 1 and h 2 and related uncertainties of the reservoir system. The resulting posterior model parameters and their standard deviations suggest that all model parameters are well constrained compared to their prior statistics (Table 2, Fig. 7). Figure 7 demonstrates the improvement of the uncertainties of model parameters through each data assimilation run. No clear correlation between model parameters have been identified except for the height change ( h 2 ) and the radius (R 2 ) of source 2 (Fig. 7d). The correlation between h 2 and R 2 is clearly shown by Fig. 8, identified as a linear relationship between the logarithm of the parameters with a slope close to − 0.5. Apparently, all posterior parameters can be constrained independently,  1 (a, b) and source 2 (c, d), including all 4 data assimilation steps. The final parameters are shown with red circles but the product h 2 R 2 can be estimated better than the two parameters separately. This is related to the fact that the subsidence of a relatively small source is in first order proportional to its compacting volume, for which π�h 2 R 2 is the measure.

Discussion
PS-InSAR time series of ground deformation together with the modeling results suggest that uplift of the Szentes geothermal field due to pore pressure recovery occurred between November 2002 and December 2006. Since no injection wells were operated at Szentes before 2018, and production temperatures remained relatively constant through the entire production period (Bálint and Szanyi 2015), we explain ground uplift with pore pressure increase due to natural recharge. The inferred pressure recovery of the Szentes geothermal area is in agreement with increasing water levels measured at the bottom of the Upper Miocene (Upper Pannonian) sediments between 2000 and 2010 by Bálint and Szanyi (2015) (Fig. 2b). Additionally, on top of the recharge of the whole area, subsidence of the northeastern part of the geothermal field is observed, although pressure and groundwater level measurements indicating pressure decline in the reservoir are not available.
With the Geertsma models for surface displacement due to two disk-shaped reservoirs embedded in an elastic half-space with uniform pressure increase and decline, we managed to reproduce the major anomalies observed on the ground movement maps (Fig. 6). However, significant misfits between observed and modeled surface deformation remained at certain locations (Fig. 6c). We attribute those to the limitations of the model, assuming constant reservoir properties and a perfectly circular shape. As a matter of fact, sedimentation in the Makó trough followed aggradational and progradational cycles, with various lithologies (e.g., Sztanó et al. 2013). The resulting local variations in sand and clay content of the reservoir cause heterogeneous properties, most importantly a highly variable compaction coefficient, causing variations in the magnitude of ground deformation mapped with InSAR.
The formula of Geertsma (1973) shows a linear relationship between ground motion and the product of reservoir compaction coefficient, reservoir height, and pressure change. Inversion techniques fail to constrain each of these parameters individually. Therefore, we constrained the product of the three parameters of both sources (Table 2 Fig. 7). Then, we plotted the decompaction coefficient associated with the source 1 as a function of reservoir pressure, with constant reservoir heights of 200, 400 and 600 m (Fig. 9), based on the mean posterior products of the three parameters ( h 1 ). We outlined intervals with reasonable values for pressure change based on water level variations for the section in Fig. 2b, groundwater level measurements in well K-652. The inferred decompaction coefficient of source 1 varies between ~ 0.2 and 2 × 10 -9 Pa −1 (Fig. 9), characterizing the elastic behavior of the Szentes geothermal reservoir.
The compaction coefficient attributed to the small source with pressure decline was not assessed due to the lack of pressure/groundwater level measurement supporting the pressure decline that could explain the observed subsidence, and the correlation between the radius and the height change of source 2, as shown in Figs. 7d, 8. We attempted to compare the resulting decompaction coefficient attributed to the Upper Pannonian sediments with a porosity of 21-28% at Szentes with various sandstones worldwide, to reveal if the resulting values are realistic. We only made comparisons with compaction coefficients of sandstones that were inferred from ground deformation observations. For instance, Thienen-Visser et al. (2015) estimated the compaction coefficient of the Rotliegend sandstones in the Groningen gas field with porosities of 12-22% to a maximum value of 2 × 10 -10 Pa −1 . This value agrees with the lower limit of the inferred decompaction coefficient estimated for the Upper Pannonian sandstones in Szentes. Since the lower limit of reservoir porosity in Szentes (21%) is almost identical with upper limit of the porosity of the Rotliegend sandstone in Groningen (22%), we can conclude that our estimates are realistic.
The decompaction coefficient of the Szentes reservoir system may be different from the compaction coefficient corresponding to the earlier depressurization periods, as the water-bearing sand bodies are often separated with aquitards composed of clay layers that often experience permanent, irreversible compaction (e.g., Pijnenburg et al. 2019). The individual compaction coefficient of clays are considered to be significantly higher than the sand bodies (e.g., Zimmerman 1990). As a result, they could contribute significantly to the overall compaction coefficient of the reservoir system. The compaction of aquitards may be responsible for the observed local subsidence in the northeastern part of the field, although this hypothesis is not supported by any pressure/groundwater level measurements. Compaction of highly compressible lithologies (mudstones) were identified as main sources of ground subsidence for instance at the Wairakei-Tauhara and Ohaaki geothermal areas in New Zealand, with up to 15 m (Hole et al. 2007). To properly assess the contributions of elastic and inelastic deformation at the study area and fully understand the compaction behavior of the reservoir, ground deformation observations and more frequent pressure measurements of the major depletion period are required.

Conclusions
We performed PS-InSAR time series analysis of Envisat SAR images acquired on descending satellite tracks, for the period of November 2002-December 2006, covering the Szentes geothermal area in SE Hungary. The ground movement observations were used in inverse models based on the analytical function of Geertsma (1973), for ground deformation induced by a constant pressure change in a cylindrical source at depth. We constrained the model parameters with an ensemble smoother with multiple data assimilation and we managed to reproduce the main features of the observed surface movements with the combination of a larger source with positive pressure change, and a smaller depleting source. The procedure yielded reasonable estimates of the model parameters and their uncertainties.
Ground movements together with the modeling results suggest that uplift of the Szentes geothermal field occurred between November 2002 and December 2006, which is explained by pore pressure recovery supported by natural recharge. Based on the InSAR data, inverse modeling, and groundwater level variations, we managed to provide reasonable estimates on the decompaction coefficient of the reservoir (~ 0.2 × 10-9-2.0 × 10 -9 Pa −1 ), characterizing the elastic behavior of the sandstone aquifers. Additionally, permanent compaction of aquitards represented by imbedded clay layers of the reservoir system may have occurred during the earlier depressurization period, from ~ 1970 to the mid-1990s. Compaction of aquitards may have continued in the northeastern part of the field, where a minor negative ground motion anomaly is observed.
The InSAR surface movement estimates, employed in an ES-MDA inverse modeling scheme, have shown to combine to an efficient tool to make a first-order estimate of the properties and processes at the Szentes geothermal area. The inferred properties and processes are highly relevant as an input of more complex reservoir models and for predictions for future production scenarios, targeting a sustainable production plan. To achieve more reliable estimates on the reservoir behavior, ground motion observation of the entire thermal water extraction period, and more frequent pressure/hydraulic head measurements would be required. Future projects combining InSAR with forward and inverse modeling with high-resolution subsurface datasets potentially have an important role in in the sustainable exploitation of the highly prospective geothermal resources in the Pannonian Basin.