Ground displacement analysis in the Sarulla geothermal field using Persistent Scatterer Interferometric Synthetic Aperture Radar (PS-InSAR) method

The Sarulla geothermal field is aggressively developed in Indonesia, reaching 330 MW in two consecutive years. Furthermore, the field is situated along the active Great Sumatra Fault. This condition has the potential to induce surface deformation phenomena attributed to exploitation activities, tectonic activities, or a combination of both. The Persistent Scatterer-Interferometric Synthetic Aperture Radar (PS-InSAR) method is a multi-temporal Synthetic Aperture Radar (SAR) image processing technique employed to detect surface deformations with millimeter-level accuracy. PS-InSAR has been applied in this paper over the Sarulla geothermal field, which comprises two exploited prospects (Silangkitang and Namora-I-Langit) and two unexploited prospects (Donotasik and Sibualbuali). The objective of this paper is to analyze the characteristics of ground displacement over the exploited and unexploited geothermal areas. Multi-temporal images of SENTINEL-1A with ascending and descending orbital modes from 2017 to 2021 were processed using PS-InSAR method. A series of interferograms from the selected primary image were paired with secondary images using SNAP program and Python-based program snap2stamps. The generated interferograms were used as input for the StaMPS program to obtain the persistent scatterers as the best coherent points. Ground displacement towards the satellite line of sight was observed within those two exploited geothermal areas adjacent to the production and injection zones. The descending images resulted in more significant PS points than ascending images. The exploited geothermal area shows a clustered pattern of PS points with a significant rate of ground displacement, whereas the unexploited area suggests a wide dispersion of PS points with a slight amount of displacement rate. Additionally, the detected ground displacement that occurred adjacent to the Great Sumatran Fault may indicate the active movement of the fault.


Introduction
The rapid growth of geothermal development implies rapid geothermal fluid exploitation that could potentially cause ground displacement due to changes in reservoir pressure, temperature, and underground water level [1][2] [3].The pressure change induces alteration in the bulk volume of the reservoir or subsurface space and affecting the reservoir displacement, that propagates to the surface, causing ground displacement.Numerous studies related to ground displacement due to geothermal fluid extraction have been conducted worldwide using leveling surveys such as those at Kamojang [4], Wairakei -Taupo [5], The Geysers [6], Cerro Prieto, Lardarello, and Travale [7].
The utilization of Synthetic Aperture Radar (SAR) has become widespread for monitoring ground displacement offering extensive spatial monitoring capabilities.Interferometric SAR is generated from the phase changes between two overlapping SAR images.The PS-InSAR method, applied in this paper, is one of many processing methods to detect the ground displacement over a wide area using a series of interferogram.This method enables us to address the decorrelation problem by identifying points or pixels with high-quality coherent targets that remain stable over long-time intervals, termed as persistent scatterers (PS) [8].
The objective of this paper is to analyse ground displacement using SAR data associated with exploitation and exploration purposes.Ground displacement induced by changes in subsurface properties due to fluid extraction and reinjection can provide an indication of spatially distributed active reservoir or injection area.Thus, the ground displacement monitored in this paper was used to identify the active production and reinjection zone.The characteristic includes the pattern and significance of ground displacement over the exploited area, serving as lessons to assess the exploration target over the unexploited area.
To achieve the objectives, the Sarulla geothermal field was selected as the study area due to its location adjacent to the major strand of the Great Sumatra Fault (GSF), the major right-lateral strike-slip fault that extends 1650 km along Sumatra Island (Figure 1a).The Sarulla geothermal field comprises four prospect areas within adjacent regions, with two of them already exploited.The two exploited areas also have a distinguished geothermal system, making this field more attractive for discussing the characteristic of ground displacement among these areas.In addition, we examine the ability of PS-InSAR method to obtain ground displacement map over the dense vegetation adjacent to the study area (Figure 1b).

Study Area
The Sarulla geothermal field is situated approximately 50 km south of Toba Lake, the largest volcanic crater lake in Sumatra (Figure 1).Based on exploration result conducted by Unocal from 1993 to 1998, there are three proven geothermal prospects area beneath the Sarulla contract area, i.e., Eastern Sibualbuali, Silangkitang (SIL), and Namora-I-Langit (NIL) with combined capacity of 330 MW.
Additionally, there is one undrilled prospect area at Donotasik [9].Currently, Sarulla Operations Limited operates this project under a joint operation contract with PT Pertamina Geothermal Energy This project significantly contributed to additional installed capacity around the 2017 -2018 period, reaching 330 MW [10].Commercial electricity generation commenced from Unit 1 (110 MW) in March 2017, located in the SIL area, followed by Unit 2 (110 MW) in October 2017 and Unit 3 (110 MW) in May 2018, both situated in the NIL area.To support the power plant capacity in the SIL area, four production wells are strategically positioned within the pulldown block along the Great Sumatra Fault (GSF), accompanied by nine injection wells located in the southeast and northeast of the production area.In the NIL area, the power plant is supported by 11 production wells concentrated in the northern part and nine injection wells situated in the eastern part, targeting the GSF [11] [12].
The NIL and SIL areas are situated between the major strands of the GSF and lie within a low gravity zone between these trands [9].These geothermal systems are associated with large volcanic debris basins that exhibit the structures of half graben in profiles normal to the two active faults strands of the GSF [13].Despite NIL and SIL having similar structural attributes, the geothermal system of these areas differ.NIL is associated with volcanic geothermal system, whereas SIL is a fault-controlled geothermal system associated to a deep heat source [11].
The NIL volcanic complex consists of Quaternary andesitic-rhyolitic lavas and tuffs, whereas SIL consists of Paleozoic metasediments considered as basement rocks, Tertiary deposits overlie the basement, followed by Mesozoic granites, and the upper part is recent alluvium (Figure 2).Donotasik is situated southeast of the SIL area and in the western part outside the Hopong Caldera.This prospect is the only area that has not been drilled to date.Several small domes of dacitic-rhyolitic were found adjacent to the eastern part of the shear zone towards the north [14].
The Sibualbuali area is located at the southeast end of the contract area, and its geothermal system is associated with the Sibualbuali volcano, bounded by the two major strands of the GSF.Four wells drilled targeting the eastern flank of Sibualbuali suggest that the geothermal system have strong lateral and vertical temperature gradients along the strands of the GSF [9].

SAR Dataset
A total of 49 C-band SAR images (at 5.3 GHz frequency corresponding to a 5.6 cm wavelength) from Sentinel-1A have been employed to investigate ground displacement over the study area, covering the period 2017-2021.The specification of the SAR images include Single Look Complex (SLC) with vertical-vertical (VV) polarization Interferometric Wide (IW) mode in both ascending and descending direction.Table 1 depicts the attribute data for each dataset.One single acquisition image was selected from each ascending and descending dataset as primary image to be considered as baselines in space and time.Subsequently, the primary image was coregistered with the remaining secondary images to form an interferograms.Figure 3 illustrates the chart corresponding to the baseline and temporal perpendicular for both ascending and descending dataset.

PS-InSAR Processing
The PS-InSAR has been applied to investigate the displacement, considering its ability to obtain a more accurate measurements compared to D-InSAR.On the other hand, when compared to the SBAS method [16], it incurs a high computational cost due to the large number of interferograms to be processed.Therefore, in this study, the PS-InSAR method is considered the most optimum for application.
The measured interferogram phase (ϕ) from two overlapping SAR images is given by equation [17]: where   is the phase attributed to the ground displacement () along the azimuth direction.λ is the SAR wavelength.  correspond to the linear displacement rate which derived from the shift in the orbital trajectory between the two images. ∥ is the distance between imaging positions parallel to the radar's look direction.  is the phase due to inaccuracy of external DEM,  is the height above the reference surface,  is the radar range, and θ is the incidence angle of the radar.  is the atmospheric phase delay.  is the noise due to temporal and geometrical decorrelation.
The PS-InSAR method in this paper was performed using three main open-source software, i.e., Sentinel Application Perform (SNAP) provided by ESA, the Stanford method for persistent scatterers (StaMPS) [18], and QGIS.SNAP program is used to generate interferograms.We used a set of Python scripts, namely snap2stamps [19] to call the SNAP routines to perform the TOPSAR interferometric processing automatically.Candidates of PS pixels were selected using mt_prep_snap as the UNIX command before processed in StaMPS.The selected PS candidates were ready to be processed using StaMPS to obtain the LOS displacement from each PS.
Figure 4 illustrates the coherence values, wrapped phase, and unwrapped phase of a sample interferogram from both ascending and descending directions.The coherence map depicts the pixel coherency value between two SAR images, with higher value indicating more coherence and a greater potential to be selected as PS candidates.The wrapped phase is obtained from the selected PS candidates, while the unwrapped phase is obtained after the PS weeding process, leading to a reduction in pixels.Positive value in both wrapped and unwrapped phases indicates a decrease in LOS range, corresponding to uplift deformation, whereas negative values indicate an increase in LOS range, corresponding to subsidence deformation.

PS-InSAR Deformation Results
In the PS-InSAR processing with a dispersion amplitude requirement of 0.41, the ascending and descending data stacks generated 62,315 and 82,802 PS, respectively.Figure 5 depicts the PS representing the displacement rate over the study area overlaid with the ESRI satellite imagery.The displacement rate relative to the Line-of-Sight (LOS) direction across the study area ranges from about -14.9 to 9.1 mm/year using ascending mode data and -26.4 to 23.1 mm/year using descending data.Positive and negative values follow similar rules as the phase changes.Although there is a slight difference in the LOS displacement rate between ascending and descending images, the spatial distribution appears similar in both directions.
Four areas were considered for analysis, covering NIL, SIL, Donotasik and Sibualbuali.The prospect boundaries from [9] for Donotasik and Sibualbuali area, and the updated prospect boundary pointed out by [11] for NIL and SIL areas were used.The PS clustered in exploited areas shows a more significant LOS displacement rate compared to the unexploited areas, which corresponds to geothermal exploitation activity.Additionally, the PS are distributed along the main GSF and its strands.Twelve points, labeled 6 A to L, associated with ground feature over the NIL and SIL area were denoted.Each point defines a region with a radius of 300 m, comprised of several PS to be analyzed.

Deformation Characterization at NIL Area
Four points associated with the production and reinjection area, as well as geological features, are depicted in Figure 6.Point A and C refer to the areas where the power plant and production well pad are located, respectively.Point B refers to the reinjection area, while point D covers an area over an outcropped alteration zone.The power plant and the vicinity of production wells area (point A and C) observed have a high number of PS and a negative LOS displacement rate, defining a subsidence phenomenon in the adjacent area (Table 2).This probably corresponds to the production activity in the vicinity of this area, which may be affected by fluid extraction from the reservoir.Another subsidence phenomenon is observed in the injection zone (point B) which is situated adjacent to the GSF.The negative correlation between fluid injection and subsidence phenomenon possibly occurred due to the fault activity playing a main role in ground displacement.According to the time-series LOS displacement (Figure 6), since the commencement of Unit III power plant, the ground displacement over these regions has constantly subsided until end of 2019 and became stable or slightly subsided after early 2020.
TUuplift displacement relative to the descending LOS is observed in the west part where the outcropped alteration zone is situated (point D).Meanwhile, from the ascending LOS, this region suggests subsiding at a near zero rate or remaining nearly stable.The uplift may probably correspond to the lateral and vertical movement of the GSF adjacent to the area.We observed three clustered PS over the SIL area denoted with points E, F, and G as depicted in Figure 7. Point E covers a region comprised of production wells and the power plant, while points F and G cover a region comprised of injectors.According to the PS statistical distribution (Table 3) from descending mode over the production area (point E) and the northwest injection area (point F), there is an uplift displacement.A slight subsidence displacement or nearly condition is observed over the southeast injection area (point G).In contrast, the ascending data depicts a subsidence trend suggesting that the direction of ground displacement is towards the satellite LOS.The existence of structure as the primary permeability in fault-controlled geothermal system plays a main role in the interplay between production and injection, Thus, production wells and injectors are targeted to intersect the GSF.The dominant uplift displacement over the production towards the northwest injection area probably indicates the existence of connectivity among these areas along the fault strand.

Deformation Characterization at the Unexploited Area of Donotasik and Sibualbuali
The PS are distributed widely over the entire area and show no significant displacement rate, as depicted in Figures 8 and Figure 9.The lack of PS is likely related to dense vegetation covering the area.The widely dispersion of PS suggests that ground displacement is not related to geothermal exploitation activity, as observed in NIL and SIL.The PS consistently exist in urban areas as the effect of stable pixels corresponds to man-made structures.At Donotasik, the PS were observed within the prospect boundary in a channeling pattern (Figure 9).This pattern coincides with several interpreted fault strands, which probably corresponds to the fault activity.The dispersed pattern of PS may occur due to the vegetated land cover, which produces distributed scatterers.Two points selected among the channeling pattern (H and I in Figure 9) depict a similar displacement trend both from ascending and descending directions.These selected points coincide around the fault strands, showing a slight uplift relative to the satellite sights (Table 4).
At Sibualbuali, the PS form a channeling pattern situated at the east part outside the Sibualbuali prospect boundary and east away from the major strand of the GSF (Figure 9).Within the prospect boundary, the PS disperse widely, mainly due to vegetated land cover, which produces distributed scatterers and depicts a uniform, slight displacement rate, suggesting no existence of geothermal exploitation over the entire area.Three points (J, K, and L in Figure 9) depict a relatively small significance of displacement except in the area of point L, which coincides with the fault at the southeast part of the prospect area (Table 5).

Conclusion
The PS-InSAR method was applied to the SENTINEL-1A dataset in two orbital modes over the Sarulla geothermal field.The method successfully generated 62,315 and 82,802 PS from the ascending and descending data stack, respectively, corresponding to the high-quality coherent pixels during the observation period.The descending orbital mode resulted in a more significant number of PS points rather than the ascending mode, mainly due to vegetation coverage in certain areas.The characteristic of PS distribution over the two exploited geothermal areas (NIL and SIL) revealed a clustered pattern of PS points.The processing of PS-InSAR using the descending dataset produced more significant ground displacement.This indicates that the direction of displacement towards the LOS of satellite in the descending orbital mode.Both NIL and SIL showed distinguished displacement trends.Predominantly, the ground displacement over the NIL area exhibited subsidence, whereas at SIL, uplift displacement seemed to dominate the entire exploited area, although not at a significant rate.
A negative correlation, where the injection area induces a subsidence, occurred in the injection area of adjacent to the strand of GSF.This subsidence is probably due to fault activity inducing topographical changes adjacent to the fault.In SIL, uplift displacement occurred adjacent to the production area, whereas the two injection areas indicated uplift in the northwest and a near stability in the southeast.The uplift displacement in SIL's production area indicated connectivity between the production area and northwest and southeast injectors through the fault, probably inducing pressure extension in the subsurface.
In unexploited areas, we observed unclear displacement patterns over some areas of Donotasik and Sibualbuali.The unclear PS pattern occurred due to the lack of C-band coherency in dense vegetation areas.In Donotasik, some of the PS formed a channeling pattern coinciding with the fault strand, indicating topographical changes corresponding to fault block movement.In Sibualbuali, the channeling pattern of PS formed at the east outside the prospect area and east away from the GSF main strand.However, at the southeast end of the prospect boundary, the channeling PS coincided with the fault.According to the displacement characteristics over the exploited areas, where the existence of the fault induced a more significant displacement rate, we concluded that this channeling pattern of PS points indicates the fault's existence as a high permeability zone over Donotasik and Sibualbuali prospects.

Figure 1 .
Figure 1.Location of the Sarulla Geothermal Field in Sumatra Island (a) and detailed location comprises four geothermal prospect areas superimposed on ESRI Satellite imagery depicting dense vegetation land coverage (b).

Figure 2 .Figure 3 .
Figure 2. Geological map of the Sarulla geothermal field with distribution of thermal manifestation denoted by red dot and stratigraphical units comprises Paleozoic to Quaternary lithologies compiled from [15] and [9].

Figure 4 .
Figure 4. Map of coherence, PS wrapped phase, and PS unwrapped phase (from left to right) of generated interferogram from paired images -primary and secondary images of ascending orbit (top) and descending orbit (bottom).

Figure 5 .
Figure 5. PS-InSAR LOS displacement from SENTINEL-1A images in ascending (a) and descending (e) orbits over the Sarulla Geothermal Field.Each PS presented in color dots is the LOS displacement rate in mm/year unit with the base map of ESRI image.The clustered PS consistently exist over the man-made structures, such as power plant facility (b, c, f, g) and urban area of Sipirok (d and h).

Figure 6 .
Figure 6.The PS-InSAR LOS displacement map in ascending (a) and descending (b) over the NIL area (top) and the time series of LOS displacement of each observed points (bottom).The yellow dashed line denotes the geothermal prospect based on resistivity surveys [11].

Figure 7 .
Figure 7.The PS-InSAR LOS displacement map (top) in ascending (a) and descending (b) over the SIL area and the time series of LOS displacement of each observed point (bottom).The yellow dash line denotes the geothermal prospect based on resistivity surveys [11].

Figure 8 .
Figure 8.The PS-InSAR LOS displacement map in ascending (a) and descending (b) over the unexploited area of Donotasik (top) and time series displacement (bottom).Prospect and faultscompiled from[9] and[20]

Figure 9 .
Figure 9.The PS-InSAR LOS displacement map in ascending (a) and descending (b) over the unexploited area of Sibualbuali (top) and time series displacement (bottom).Prospect and faultscompiled from[9] and[20].

Table 1 .
Specification of Sentinel-1A ascending and descending orbit data used in this study.

Table 2 .
Statistical distribution data of PS from ascending and descending mode within selected region adjacent the NIL exploited area.

Table 3 .
Statistical distribution data of PS from ascending and descending mode within selected region adjacent the SIL exploited area.

Table 4 .
Statistical distribution data of PS from ascending and descending mode within selected region adjacent the Donotasik prospect area.

Table 5 .
Statistical distribution data of PS from ascending and descending mode within selected region adjacent the Sibualbuali prospect area.