Hazard Risk Geospatial-Based Assessment of Groundwater Depletion and Land Subsidence for Water Conservation

1. Introduction

In recent decades, groundwater has become one of the most important sources of water among the available resources throughout the world [1,2]. The importance of groundwater as a water resource is increasing as a result of the continuous growth of the global population, which is predicted to exceed 11 billion in 2100 [2]. Nevertheless, many cities are currently facing severe water resource instability [3], with developing and poor urban areas being the most vulnerable to such instability. A total of 21 of the world’s 33 megacities are located on the coast, and in Asia, for example, more than one billion people live within 100 km of the sea [4]. The impact of climate change on water resources and their management has already been seen in various regions around the globe [5]. In 2010, one-fifth of the world’s population, roughly 1.2 billion people, were living in areas with water scarcity. Following projected urban and industrial water use, this number will double to 2.4 billion by 2050 [4]. From these figures, it is clear that sustainable water resource management, particularly of groundwater, is critical to maintaining water demand. In many areas, a lack of management measures for groundwater utilization has led to the extensive over-exploitation of this resource, inducing groundwater depletion and thereby becoming a threat to the sustainable development of socioeconomic and ecological health [1,3].

Concerning groundwater resource management, the development and integration of pilot case studies on hazard–damage relationships at local or regional scales are essential to provide information for relevant agency managers and policy makers [6,7,8]. Despite its criticality, the relationship between hazards and damage is poorly studied [6]. Regarding the inter-correlation of water-related hazards with damage to land, such as land subsidence, it is considered advantageous as a parameter in assessing the phenomenon of groundwater depletion. Land subsidence can be defined as a geological phenomenon that causes a slow lowering of ground surface elevation [9]. One of the factors responsible for land subsidence is groundwater extraction [10,11]. Land subsidence also impacts groundwater, in that it reduces the capacity of aquifers to store water [12].

Land deformation measurements extracted from interferometric synthetic-aperture radar (InSAR) data processing [13,14] provide a method to estimate groundwater level, particularly in confined aquifers below the land surface [15,16]. Land deformation is the process of soil compression as a result of the increase in effective stress [17]. Groundwater extraction removes the water pressure that increases the effective stress developed by the soil skeleton, resulting in the compaction of aquifers and aquitards [17,18]. The elevation change of the groundwater surface is comparable to the total compaction throughout the confined aquifer system [9]; put simply, groundwater draw-down produces a time constant that delays compaction, which is expressed as subsidence of the ground surface above the aquifer. Therefore, land deformation measurements above a confined aquifer over time enable the estimation of the change in water level over that specific period [18]. Due to the lack of conventional groundwater level measurement owing to inherent shortcomings, InSAR-derived land deformation data present a novel technique that is efficient for producing information necessary for groundwater-related estimation.

The utilization of the InSAR technique, which facilitates measurements with centimeter-scale accuracy, makes it possible to assess deformation due to groundwater pumping activities. For instance, elastic storage coefficients have been effectively estimated using collocated InSAR and well measurements [19,20]. Furthermore, it was possible to perform a long-term analysis of aquifer system responses to pumping based on InSAR time-series analysis using the Persistent Scatterer (PS) method [21,22] and Small Baseline Subset (SBAS) method [18,23,24,25]. These methods highlight the advantages of the advanced InSAR technique for obtaining groundwater levels based on time-series land deformation measurements. Moreover, regarding the evaluation of InSAR results, centimeter-scale elevation accuracy can be achieved through measurement of continuous satellite navigation technology such as a Global Navigation Satellite System (GNSS) [26,27,28].

Groundwater is primarily supplied for domestic usage in urban areas and for irrigation in agricultural areas [2]. The application of InSAR can be challenging in different areas, such as urban and agricultural regions. Conventional InSAR processing encompasses decorrelation of the interferometric signal. Moreover, the characteristic electrical changes of scatterers lead to the appearance of temporal and geometric decorrelations. Time-series SAR interferometry represents an advanced technique of differential InSAR. It involves processing multiple SAR images acquired at different times over the same area through an effective data processing and analysis procedure to extract information on phase displacement from the other phase components [29]. The analysis enables a decrease in atmospheric and topographic errors to reduce signal decorrelation and residuals. The PS [30,31,32] and SBAS [33,34] methods of time-series InSAR have been mainly applied to the study of the temporal evolution of land deformation. These methods offer a temporal solution based on the principle of dominant reflectivity of centers of permanent scatterers to dilute the uncorrelated phenomena [31]. The PS method mainly records Persistent Scatterer (PSs) reflectors, such as those found in urban areas, whereas the SBAS method responds to Distributed Scatterers (DSs), which are typically found in agricultural regions [29]. Thus, the extracted measurement of subsidence distribution can be employed to evaluate groundwater level information in areas with wide area coverage.

In this study, the utilization of time-series InSAR to retrieve land displacement is proposed to assess the change in groundwater level in urban and agricultural regions on the island of Bali, Indonesia. Despite its importance, groundwater has become an issue on the island due to risks related to large-scale groundwater extraction [35,36,37,38,39]. Appropriate water management is crucial for sustainable groundwater development. Information regarding the groundwater condition, such as the water level, is fundamental for resource development strategies in the area [39]. In this regard, in this study, the distribution of land displacement generated from time-series InSAR processing using the PS and SBAS methods is integrated with the distribution of groundwater level interpolated from monitoring well site data in order to evaluate the correlation of land subsidence with groundwater depletion in the area. Geospatial information on groundwater depletion throughout Bali Province, Indonesia, is presented based on the analysis, with this information being essential for the development of groundwater utilization strategies for sustainable water management.

2. Materials and Methods

2.1. Study Area

The study area of Bali (Figure 1) extends from 8°03′140” to 8°50′48” S and 114°25′53” to 115°20′2” E. Administratively, as a province, Bali covers an area of 5632 km², including surrounding smaller islands, i.e., Nusa Penida, Nusa Lembongan, Nusa Ceningan, Serangan, and Menjangan. The area consists of Denpasar City as the provincial capital and eight regencies. The land use comprises 61.3% dryland agriculture, 23.4% forest, and 14.3% wetland rice fields, with alluvial, Regosol, Andosol, Latosol, and Mediterranean soil types. The alluvial soil is distributed on the western side of the island, for example, in the Jembrana Regency. Latosol is generally found in the southern regions, such as the Badung Regency and Denpasar City. Regosol is present in the northern region of the Karangasem Regency, in addition to the Singaraja Regency [40].

2.2. Data Sources and Instruments

Well data and satellite data were used as data sources for analysis. The well data consists of groundwater level measurements at monitoring wells. A report on groundwater level measurement from 2008 to 2013 at 18 monitoring wells in the study area was provided by the Ministry of Mining and Natural Resources (ESDM) of Indonesia in Bali Province (see Table 1). The monitoring wells (see Figure 1) are distributed within the groundwater basin boundaries. A groundwater level measurement was conducted in 2017 at the monitoring wells through field observation by utilizing a water-level measurement tool. Additionally, a handheld GPS was used to record and check the position information of the wells.

The primary data used for generating land subsidence information was satellite data, i.e., SAR data of the PALSAR from the Advanced Land Observing Satellite (ALOS) satellite, Digital Elevation Model (DEM) data from the Shuttle Radar Topography Mission (SRTM), and GNSS-CORS data. Open-source computer software was employed for the processing of satellite data. InSAR processing was conducted using GMT5SAR software (Scripps Institution of Oceanography, University of California, San Diego, CA, USA); this software runs on the UNIX operating system using C codes [41,42]. RTKLIB software (Tokyo University of Marine Science and Technology, Tokyo, Japan) was specifically utilized for GNSS positioning based on the CORS data [44,45]. Moreover, the QGIS image processing software (Spatial Planning and Information, Department of Land Affairs, Eastern Cape, South Africa) was used to view, edit, and analyze geospatial information from the processed satellite data [45,46].

Seven ALOS PALSAR track acquisitions cover the entire study area (see Figure 1), with a total of 63 scenes from 2007 to 2010. The tracks consist of three paths, i.e., path numbers 422, 423, and 424, with two, three, and two rows, respectively (see Table 2). The data characteristics used for the analysis were fine mode, single polarization (HH), processing level of 1.0, 236 mm wavelength, 46-day orbital revisit period, ascending, 34.3° incidence angle, and 13.1 km critical baseline [47,48,49,50]. Based on the look-direction characteristic of PALSAR, the azimuth (compass) direction of the radar illumination is perpendicular to and orientated toward the right of the satellite track [49]. Put simply, based on Figure 1, the illumination azimuth direction is East North East (ENE). Moreover, one scene of the SRTM DEM data was utilized for topographic phase removal of InSAR processing. The data was SRTM 1 Arc-Second Global elevation data that extends to latitudes of +/−60° and has a resolution of 1 arc-second, or approximately 30 m [51,52]. In order to evaluate the InSAR-derived land displacement, data recorded at three CORS stations, i.e., CSRJ, CDNP, and CPBI (see Table 3), that provide GNSS data consisting of carrier phase and code range measurements was used to generate ground surface elevation in subsided areas. GNSS enables autonomous global geospatial positioning [27,53]. CORS networks have been rapidly installed in Indonesia to provide practical measurement capabilities, in addition to high-accuracy positioning [54].

2.3. Research Scheme

A research scheme (see Figure 2) was developed for the hazard risk assessment of groundwater depletion due to land subsidence. The method was applied to the study area of Bali, Indonesia, which represents a small-island developing site with limited water resources. Information on InSAR-derived subsidence distribution was integrated with groundwater depletion information calculated from monitoring well data. Specifically, the InSAR result was compared with GNSS-CORS data measurement. Lastly, a map was produced to present the conservation zones that were found to be vulnerable to land subsidence as a result of groundwater utilization, based on spatial information of groundwater depletion. Such measures are essential for sustainable groundwater management.

2.4. Groundwater Depletion Monitoring

Groundwater depletion was calculated using the collected secondary data of groundwater level from 2008 to 2013 and field observation data from 2017. The percentage of groundwater depletion (Δh) was determined using groundwater level in the initial year (h₁) and groundwater level in the final year (h₂) of the data in 2008 and 2017, as follows:

$$∆\mathrm{h}=\left({\displaystyle \frac{{\mathrm{h}}_1-{\mathrm{h}}_2}{{\mathrm{h}}_1}}\right)\times 100\mathrm{\%}$$

(1)

The total groundwater depletion was determined by subtracting the final groundwater level from the initial level in the period 2008–2017 for each monitoring well. The values of total groundwater depletion (in percent) were input into the monitoring well points as vector data. The information was then interpolated to create a map of the distribution of groundwater depletion throughout the study area.

Spatial interpolation involves the use of vector points with known values to estimate values at locations where the value is unknown to create a raster surface covering an entire area. The interpolation result is typically a raster layer. It is important to find a suitable interpolation method to optimally estimate values for unknown locations [46]. For the distribution of groundwater depletion, the interpolation was performed in the QGIS application using the Inverse Distance Weighting (IDW) method. In this method, the sample points are weighted during interpolation such that the influence of one point relative to others declines with distance from the unknown point that will be created. Weighting is assigned to sample points through the use of a weighting coefficient that controls how the weighting influence will decline as the distance from the new point increases. The greater the weighting coefficient, the less influence distant points from the unknown point will have during the interpolation process. As the coefficient increases, the value of the unknown point approaches the value of the nearest observational point [55]. In general, values greater or equal to 1 are used for the coefficient. The weighting coefficient was set to 2, which represents the default value known as the inverse distance squared weighted interpolation. The coefficient was adjusted to the number suitable for analysis. The set value resulted in a few different surfaces at a larger distance from nearby points.

2.5. Land Displacement Measurement

Measurements of land subsidence were mainly derived from InSAR analysis. Fundamentally, the InSAR technique combines two SAR scenes acquired at different times in the same area, which are then used to generate an interferometric pair [56]. For the measurement of the phase difference caused by ground movement, the differential interferometric (DInSAR) phase for a pixel is composed of phase components as follows [56,57]:

ϕDInSAR = ϕdefo + ϕres_topo + ϕatmo + ϕnoise

(2)

More recently, advanced time-series InSAR has been developed to overcome some of the limitations of conventional InSAR. The technique involves the use of multi-temporal SAR scenes to reduce atmospheric effects by stacking multiple interferograms. The most common time-series InSAR techniques are PS and SBAS, which are able to provide millimeter-scale precision in deformation mapping [30]. PS is appropriate for retrieving PSs with a single master for baseline configuration and pixel selection criteria based on amplitude dispersion [29,31,32,58]. Conversely, SBAS is suitable for deriving DSs by executing a small baseline configuration and a coherence-based pixel selection criterion [33,34]. In this study, the PS and SBAS methods were combined (see Figure 3) in order to enable displacement to be observed both in urban and agricultural regions. The data processing step conducted using GMT5SAR consists of preprocessing, focus and alignment, pair selection, topographic phase removal, phase filtering and unwrapping, and layer stacking. The time-series technique in GMT5SAR is generated by performing a layer-stacking process whereby the input interferograms are processed using batch processing. The characteristics of the data type affect the parameter values for processing [48].

The first step of the InSAR process was the preprocessing of raw SAR data and orbital information to extract orbit and baseline estimation, i.e., orbital position, velocity, near range, and Doppler centroid. Secondly, the images were focused to create two single look complex (SLC) images, and the repeat images were aligned with the reference image. The third step was pair selection of the interferograms based on the PS and SBAS methods. One scene was selected as a super master (reference scene) that was common to the other scenes employed as slaves (repeat scenes) following the PS method. Each surrogate master needs to be aligned with the super master before it can provide a surrogate. Conversely, certain maximum baselines are required for pair selection in the SBAS method. The baselines for SBAS pair selection were determined at 2500 days for the temporal baseline and 1000 m for the perpendicular baseline. Considering the long wavelength of ALOS PALSAR data [47,48,49,50], the large baselines are acceptable for the process.

The subsequent step after pair selection was DEM generation for topography removal and to transform the DEM from longitude, latitude, and topography into range, azimuth, and topography. The next step involved handling the reference and repeat SLCs using the DEM generated from the topographic phase removal, followed by filtering, unwrapping, and geocoding. The processes of interfering, filtering, and unwrapping computed amplitude, phase, and coherence images of the interferograms without the effect of topographic phase. Specifically, the geocoding process transformed all of the products from the range/azimuth coordinate system to longitude and latitude. The parameters of the Gaussian filter wavelength, perform land mask, and unwrap and geocode threshold that needed to be set in the configuration file were Gauss 300, 1 (perform), and 0.45, respectively, for both the PS and SBAS interferograms. Subsequently, layer stacking was conducted to stack the generated interferograms. For SBAS layer-stacking, the other parameters for the smoothing factors, namely, wavelength, incidence angle, range, and atmospheric correction, were 0.2362 m, 34.3°, 850014–848965, and 1 (perform), respectively.

A comparison of the InSAR result with RTK-GNSS measurements as a reference was carried out for evaluation. Fundamentally, RTK is able to perform measurements with up to centimeter-level accuracy by employing phase measurements of the signals’ carrier wave and relying on a single reference station or interpolated virtual station to perform real-time corrections [57,59]. The time required for a signal to travel from the satellite to the receiver is used to calculate the distance between the receiver and the satellite [26]. The data format of GNSS-CORS data downloaded from a receiver was converted into ASCII format of Receiver Independent Exchange Format (RINEX).

In regard to the data in the study area (see Table 3), the observation type of L1, L2, C1, P2, S1, and S2 recorded every 30 min can be processed using RTKLIB software. The input for the data consisted of rover, base station, navigation, precise clock, and Earth Orientation Parameters (EOPs). Parameters for running the calculation through the RTK-post tool of RTKLIB, i.e., positioning mode, frequencies/filter type, elevation mask/SNR mask, ionosphere correction, troposphere correction, and satellite ephemeris/clock, were set into static, L1+2 (combined), 15, estimate Slant Total Electron Content (STEC), estimate Standard Deviation (STD), and precise, respectively [43,44].

3. Results

Hazard risk spatial information of groundwater depletion and land displacement was extracted from the data analysis. The distribution of land displacement was constructed based on a mean velocity map of InSAR time-series displacement for approximately three years (February/October 2007 to February/December 2010). The result was evaluated with the trend of GNSS-CORS land surface elevation measurements at stations located in subsided areas. Furthermore, a map of groundwater depletion percentage was established to represent the changes in groundwater level based on the observed land subsidence. In order to present essential information for water conservation measures, vulnerability levels of high, moderate, and low were determined, corresponding to the hazard risk assessment.

3.1. Spatial Information of Groundwater Depletion

The percentage of total groundwater depletion from 2008 to 2017 at 18 monitoring wells was interpolated and mapped (see Figure 4) to present spatial information of the distribution of groundwater depletion throughout the study area. The total percentage of groundwater depletion ranged from −66.24 to 25.24% (positive values indicate an increase in groundwater level and negative values indicate a decrease).

3.2. Distribution of Land Subsidence

Subsidence was defined based on the InSAR time-series displacement. The InSAR products were estimated using multi-temporal interferograms generated through time-series analysis using the PS and SBAS methods. In regard to the PS pairs, eight interferograms for each track (providing a total of 58 interferograms for the seven tracks) were generated, with scenes 05750, 09343, and 05565 (see Table 1) as the super masters for paths 422, 423, and 424, respectively. The scenes were chosen with the criteria of being able to cover the area of interest and minimize perpendicular baselines between interferometric pairs [49,51]. In correspondence to the scatter chart of the baseline distribution (see Figure 5), the super master scene was located in the middle with a value close to 0.

In the case of SBAS pairs, a total of 84 interferograms were extracted within 1000 m of the perpendicular baseline and within 950 days of the temporal baseline. A different number of interferogram pairs was generated for each track, i.e., 31, 31, 25, 27, 28, 28, and 29 interferogram pairs for path 422 row 7010, path 422 row 7020, path 423 row 7000, path 423 row 7010, path 423 row 7020, path 424 row 7010, and path 424 row 7020, respectively. The time-series of cumulative displacement and mean velocity of the InSAR results indicated upward (+) and downward (−) displacement from 2007 to 2010 in the entire study area (see Figure 6).

The calculated mean displacement velocity (see Figure 7), which ranged from –8.4 to 8.4 mm/yr, was used to estimate subsided and non-subsided sites. Areas with upward (+) displacements (i.e., with no subsidence) were observed in the northern and eastern regions of Bali. Downward (−) displacements (i.e., subsidence) were distributed in the southern coastal areas, including an area of Denpasar City and the Badung Regency. These areas have been the center of development activity on the island, i.e., settlements and tourism accommodation. Groundwater is mainly utilized to meet domestic water needs [60]. Likewise, subsidence was also measured on the western side of the island. This area contains irrigated paddy fields that use groundwater as their main irrigation source due to limited water resources [40].

Furthermore, the InSAR measurements were evaluated with GNSS-CORS measurements as reference data. Referring to the line-of-sight (LOS) displacement equation [50], the vertical displacement was calculated by dividing the LOS value of InSAR-derived time-series cumulative displacement by the cosine of the incidence angle. Considering that the LOS displacement is proportional to the vertical displacement, regardless of the effect of horizontal displacement, the downward displacement was considered sufficient to represent subsidence in the study area. Therefore, the calculated vertical displacement was employed for further evaluation with GNSS measurements, in addition to analysis of groundwater depletion.

Regarding the GNSS-CORS measurement, the height that represents the land surface elevation in the subsided area was estimated through the RTK-post tool of RTKLIB processing program using available and processable data. The position of a rover station in the subsided area was measured based on the position of a base station located in the non-subsided area. The distribution of the GNSS-CORS stations is shown in the mean displacement velocity map (see Figure 7), which signifies stations in areas with downward (subsided) and upward (non-subsided) displacements. As shown on the map, there are two stations in the subsided area, i.e., CDNP and CPBI, and one station in the non-subsided area, i.e., CSRJ. Based on the distribution of the GNSS-CORS stations, stations CDNP and CPBI were designated as rover stations, with station CSRJ as the base station for both rover stations. The reference or base position of the selected base station was set by referring to the latest recorded information of latitude, longitude, and height of the station (see Table 3). The base station also meets the distance criteria of a base station to rover station distance of less than 80 km [50].

The significant decreasing trend of height changes at stations distributed in the subsided area, i.e., CDNP and CPBI (see Figure 8), corresponded with the subsidence trend of the cumulative displacement time series. In regard to the evaluation result that highlighted subsidence at stations on the southern side of Bali, the InSAR-derived displacement was confirmed to represent the subsidence distribution.

3.3. Hazard Risk Geospatial-Based Assessment of Groundwater Depletion and Land Subsidence

The values of the average percentage of groundwater level change integrated with the values of downward displacement at each monitoring well were plotted on a graph (see Figure 9) to present a quantitative estimation of the hazard risk correlation between groundwater depletion and land subsidence. Monitoring wells with a decreasing trend of groundwater level changes exhibited downward displacement, denoting subsidence. Based on monitoring well MW04-BD04, which is located in a subsided area, the average groundwater depletion between 2008 and 2017 was 8.03%; in comparison, the average subsidence rate between 2007 and 2010 was 5.58 mm/yr. Conversely, in non-subsided (uplifted) areas, monitoring wells were found to exhibit no trend of decreasing groundwater level. For instance, monitoring well MW12-BL01, which had a negative average groundwater depletion of –0.60%, exhibited a subsidence rate (i.e., uplift) of 3.3 mm/yr.

Furthermore, hazard risk assessment for the zonation of groundwater conservation was performed based on the values of groundwater depletion and land subsidence. The changes in groundwater level related to groundwater depletion and downward displacement, which were considered as subsidence, were classified into four conditions, i.e., safe, vulnerable, critical, and damage zones. Concerning a publication by the local agency of Public Works (DPU) related to the zonation of groundwater conservation [56], the criteria for the zonation (see Figure 9) were determined by considering the significance levels of the decreasing trends of groundwater level or the values of the average percentage of groundwater depletion and the presence/absence of land subsidence.

Based on the groundwater conservation zones derived from the hazard risk for land subsidence caused by groundwater depletion (see Figure 10), areas with no groundwater depletion or subsidence, including non-groundwater basin areas, were classified into safe zones. A settlement area with low population density dominated by vegetation such as forests in the northern region of Bali was assigned to this category, for instance, the areas of Tejakula Sub-District, Tejakula District, Buleleng Regency. Regarding vulnerable zones, areas with subsidence and low (<10%) groundwater depletion in the western region of Bali were classified as vulnerable areas, for instance, a paddy field area in Manis Tutu Sub-District, Melaya District, Jembrana Regency. In 2017, the total area of paddy fields reached 6,757 ha in the Jembrana Regency. Groundwater is mainly utilized in the area for irrigation systems, with a total requirement of 10,135.50 L [40].

Subsided areas with moderate (10–25%) and high (>25%) groundwater depletion were categorized as critical and damage zones, respectively. Such areas were observed in the southern and eastern regions of Bali, such as Kesiman Sub-District, East Denpasar District, for the critical zone and Sanur Sub-District, South Denpasar District, for the damage zone that was located in Denpasar City. Besides being developed as a settlement center, the area has also been extended for tourism accommodation, e.g., hotels and villas. The average water requirement for a hotel is 1500 L/room/day, and that for a villa is 2500 L/room/day [28].

4. Discussion

Water resources in the province of Bali, Indonesia, are quite abundant, consisting of groundwater, river water, lakes, reservoirs, and springs. Groundwater in groundwater basins represents the second-largest water resource in Bali, supplying roughly 1598 million m³/year of water [61]. This factor signifies the importance of groundwater as an alternative water resource in the area. The usage of groundwater for water supply has increased following an increase in population and economic activities in Bali. The population of Bali exhibited an annual growth of 1.23% between 2010 and 2015. Agriculture, industries, and services, including tourism, are the main economic activities in the area [47]. Groundwater has largely been used to meet domestic needs, in addition to irrigation for agricultural areas with limited water resources [40].

Several problems have arisen due to the over-extraction of groundwater in the region of Bali, i.e., the risk of land subsidence and saline water intrusion [27,28,38,40,60,62]. Groundwater depletion was observed in the southern region of Bali from groundwater level monitoring from 2001 to 2015 at wells located in four different area classifications, i.e., tourism, river basin, settlement area, and wetland [38]. Seawater intrusion has also been shown to increase groundwater salinity in the Candidasa region in the eastern part of Bali [39]. In addition to climate change, groundwater depletion and land subsidence are considered as causes of seawater intrusion [60,63,64,65,66].

Based on the results of the InSAR analysis, the shortcomings of the technique can be identified from several inconsistencies in the displacement on the InSAR images. InSAR essentially measures object directions and vertical movements with an angle based on the distances to different objects on the Earth’s surface [63]. The inconsistencies indicate that the technique may be inadequate for indicating exact amounts of subsidence, as it also estimates the object movements [30]. The accuracy of the result can be improved by employing more datasets in both ascending and descending orbits. Research has shown the effectiveness of the utilization of multi-track data for extracting a substantial value of vertical displacement [64,65,66]. Moreover, the trend of height changes at GNSS-CORS stations, i.e., CDNP and CPBI, demonstrated good agreement with the subsidence areas identified by means of InSAR. In another study, continuous GPS measurement at stations CDNP and CPBI showed negative deformation in the northern and southern regions of Bali; in comparison, a positive trend was shown in the central region. Station CSRJ was also used as a reference, considering that the area is free from earthquake effects and bias [66].

Furthermore, spatial information on the distribution of groundwater depletion can be derived through data interpolation and mapping that cover the entire study area. The groundwater level measurement from institutional data was provided for a five-year period from 2008 to 2013; however, there was a lack of data covering the subsequent period. Regarding this matter, continuous measurement is necessary to monitor the changes in groundwater level. Concerning the lack of conventional groundwater level measurement due to its shortcomings, InSAR-derived land displacement could be employed to estimate groundwater level with higher temporal resolution [13,14,15,16,17,18].

The application of geospatial analysis for hazard risk assessment offers rapid and low-cost processing [67,68,69]. Moreover, the information can also be integrated with other parameters as input for groundwater utilization plans. As detailed in the groundwater utilization plan provided by the Ministry of Public Works in Bali Province in 2013, designated protection and utilization zones for groundwater have been developed. The protected zone covers groundwater recharge areas, whereas the utilization zone consists of safe, vulnerable, critical, and damage zones. The amount of groundwater discharge is limited to 150 m³/day in the safety zone and 72 m³/day in the vulnerable zone, and groundwater use is restricted in the critical and damage zones [61,62]. Based on the results of this study, groundwater conservation zoning can be used to guide the planning of utilization zones.

5. Conclusions

In this study, a hazard risk assessment was performed using spatial information of land subsidence distribution and groundwater depletion. The risk of subsidence due to groundwater depletion can be assessed and confirmed according to the conformity of the groundwater level changes based on monitoring well data and subsidence sites identified through InSAR analysis. This type of information is important in supporting the planning of sustainable groundwater development in the area, as it considers the condition of the groundwater and the risk related to it. Concerning its advantages of quick and low-cost processing, InSAR is considered appropriate for establishing continuous spatial information of groundwater conservation zoning based on the risk of land subsidence caused by groundwater depletion. In supporting sustainable groundwater development, the information can be considered appropriate as input for related water management planning, such as groundwater conservation zonation.