Integration of Numerical Models and InSAR Techniques to Assess Land Subsidence Due to Excessive Groundwater Abstraction in the Coastal and Lowland Regions of Semarang City

Abstract

This study was carried out to assess land subsidence due to excessive groundwater abstraction in the northern region of Semarang City by integrating the application of both numerical models and geodetic measurements, particularly those based on the synthetic aperture radar interferometry (InSAR) technique. Since 1695, alluvial deposits caused by sedimentations have accumulated in the northern part of Semarang City, in turn resulting in changes in the coastline and land use up to the present. Commencing in 1900, excessive groundwater withdrawal from deep wells in the northern section of Semarang City has exacerbated natural compaction and aggravated the problem of land subsidence. In the current study, a groundwater model equivalent to the hydrogeological system in this area was developed using MODFLOW to simulate the hydromechanical coupling of groundwater flow and land subsidence. The numerical computation was performed starting with the steady-state flow model from the period of 1970 to 1990, followed by the model of transient flow and land subsidence from the period of 1990 to 2010. Our models were calibrated with deformation data from field measurements collected from various sources (e.g., leveling, GPS, and InSAR) for simulation of land subsidence, as well as with the hydraulic heads from observation wells for simulation of groundwater flow. Comparison of the results of our numerical calculations with recorded observations led to low RMSEs, yet high R² values, mathematically indicating that the simulation outcomes are in good agreement with monitoring data. The findings in the present study also revealed that land subsidence arising from groundwater pumping poses a serious threat to the northern part of Semarang City. Two groundwater management measures are proposed and the future development of land subsidence is accordingly projected until 2050. Our study shows quantitatively that the greatest land subsidence occurs in Genuk District, with a magnitude of 36.8 mm/year. However, if the suggested groundwater management can be implemented, the rate and affected area of land subsidence can be reduced by up to 59% and 76%, respectively.

1. Introduction

Land subsidence refers to a gradual or sudden vertical deformation of the ground surface, which is usually economically and socially detrimental since it often inevitably leads to structural damages to buildings and public infrastructure, as well as expands the flood inundation area [1,2,3]. Land subsidence may occur due to natural and anthropogenic (manmade) processes. The former (natural factors) includes fault compaction and tectonic movement, etc., whereas the latter (anthropogenic factors) includes excessive extraction of subsurface fluids (e.g., groundwater, oil, and gas), building loads on the ground surface, and so on [4,5,6,7,8,9].

Numerous studies have been extensively conducted to explore the physical mechanisms behind land subsidence arising from natural compaction. Gambolati and Teatini [10] employed a 1-D nonlinear finite element model to investigate soil compaction induced by groundwater flow through an isothermal sedimentary basin subjected to a continuous vertical sedimentation process to mimic the evolution of the accreting Quaternary column. Zoccarato et al. [8] used an adaptive finite-element mesh to analyze the development and evolution of the Mekong Delta in Vietnam and then described accretion and natural consolidation to characterize delta dynamics. Aside from natural compaction, another crucial natural factor that brings about land subsidence is tectonic movement. Tectonic subsidence is most common in subduction zones. This can be seen in the Pingtung Plain, Taiwan [11], which has a very active subduction zone and is prone to land subsidence. In addition, groundwater extraction frequently exacerbates land subsidence in this area. The impact of soil texture, cyclic loading, and gravitational body force on one-dimensional consolidation of unsaturated and saturated soils was quantitatively examined in a series of papers by Lo et al. [12,13,14]. Employing the theory of poroelasticity, Lo et al. [15,16] have recently formulated a boundary-value problem based on a set of coupled partial differential equations to numerically model the spatial and temporal distribution of excess pore water and air pressures in a two-layer soil system with an upper unsaturated zone and a lower saturated zone caused by external loads.

Land subsidence can be induced not only by natural sources, but also by human intervention. Since 1935, Beijing, China’s capital city, has been experiencing ground subsidence, with at least five significant sites. Yang and Ke [17] concluded that the rapid urban-area development in Beijing increased the density of multi-story buildings, which accelerated the occurrence of land subsidence, as evidenced by a time series of surface displacement data recorded from the Persistent Scatterer Interferometry Synthetic Aperture Radar (PSInSAR). Excessive groundwater pumping, on the other hand, is perhaps the most significant human intervention that gives rise to land subsidence. In fact, land subsidence resulting from groundwater removal has been a central problem in many regions over the world, including San Francisco Bay and the Florida Everglades in the United States [18], Romagna in Italy [5], central Taiwan [19], Esna City in Egypt [20], Bangkok in Thailand [21], central Saudi Arabia [22], Tehran in Iran [23], and Shanghai in China [24].

Land subsidence is an active study issue that involves various disciplines since it typically gives rise to large economic losses and disperses around the world. Several studies on land subsidence have been conducted using geodetic measurement techniques such as leveling, GPS, and satellite image (e.g., InSAR) to monitor surface deformations. Geodetic measurements are very useful in quantifying and interpreting the magnitude of land subsidence [25,26,27,28,29]. In addition to monitoring ground surface deformation, the technique of GPS and InSAR can also be used to assess the hydraulic head and aquifer system parameters for a regional area suffering from land subsidence [30,31,32,33,34].

Field experiments are another technique to gain physical insight into the behaviors of land sinking [35,36], but, they are often time-consuming and expensive. Therefore, numerical modeling that usually takes little time and is more cost-effective has become a preferred method for understanding and estimating the vertical deformation of the ground surface in local regions [37,38,39]. The integration of numerical simulation and field measurement is shown to be powerful because the implementation of both approaches can support the assessment process of land subsidence in a comprehensive framework that can systematically link its causes and effects.

The present study focuses on the northern part of Semarang City, which is located on the northern coast of Indonesia’s Java Island. Land subsidence is a common phenomenon there and has been monitored using different geodetic measurements, including leveling, GPS, and satellite images [40,41,42,43,44]. Natural compaction and excessive groundwater pumping are the main causes of land subsidence in Semarang City [43,45], which has been taking place for over a century [46]. Land subsidence due to natural compaction in the study area has been investigated by Sarah, et al. [47]. It is revealed that natural compaction is prevalent in the Demak region (east of the study area), but that this is not been the case in Semarang City. The rapid piezometric drop in Semarang city has drained the excess pore pressure developed during deposition so that the current land subsidence is entirely due to groundwater exploitation. Measurement of active tectonic role to land subsidence using GNSS data sets [48] revealed that during 2011–2017, the rate of tectonic subsidence in Semarang City was 3–5 mm/year. Considering that the subsidence rate exceeds the natural displacement rate by an order of magnitude, the calibration process has not been affected by natural displacements.

Several sub-districts in Semarang City have undergone over-exploitation of groundwater resources, which has evidently given rise to extreme declines in groundwater elevations and has thus been regarded as the main cause of land subsidence there [49]. Field observations have indicated that several housing estates in the northern part of Semarang City are found to be experiencing land subsidence. People in the area have no alternative but to move and leave their homes as a result of land subsidence. Figure 1 depicts several houses in Genuk District that have been affected by land subsidence. Figure 1a shows that in the areas where land subsidence is occurring, the ground elevation of the houses is lower than the road elevation, resulting in heavy inundation during the rainy season. The picture clearly points out that newer houses must raise their base elevation to be at or above the road level. Figure 1b illustrates the current condition of an older house whose window elevation is the same as the ground elevation due to land subsidence.

In Indonesia, the availability of groundwater observation data and land subsidence monitoring data is still very limited. The use of satellite imaging, leveling, and GPS in land subsidence monitoring and assessment is also still separate. Furthermore, most groundwater and land subsidence studies in Indonesia, particularly in Semarang City, have been limited to individual hydrogeological, geomatics, and soil compression models, with no study combining different models to test their relationship and predict changes in groundwater that cause variations in land subsidence. Furthermore, a well-structured regulatory mechanism that can effectively and significantly reduce the negative impacts of land subsidence is required. The current study was thus undertaken to start from establishing a physically-consistent framework for numerically modeling the phenomenon of land subsidence caused by groundwater abstraction that occurs in the northern part of Semarang City. Our model was verified by monitoring data of crustal deformation obtained from leveling, GPS, and Sentinel-1 InSAR in Semarang City. The InSAR analysis generates a time series of deformations and corresponding speeds in the study area, which are also used for calibration. Next, the prediction of future land subsidence with and without groundwater management measures was performed based on our well-validated numerical model. Lastly, the effectiveness of the proposed regulatory policies was quantitatively evaluated in terms of reductions in the affected area of land subsidence in spatial and temporal scales.

2. Study Area

2.1. Focus of Study Area

Semarang City is the capital of Central Java Province and also the fifth largest city in Indonesia. Semarang City is divided into 16 sub-districts: Mijen, Gunungpati, Banyumanik, Gajah Mungkur, Semarang Selatan, Candisari, Tembalang, Pedurungan, Genuk, Gayamsari, Semarang Timur, Semarang Utara, Semarang Tengah, Semarang Barat, Tugu, and Ngaliyan. Among these sub-districts, Tugu, Semarang Barat, Semarang Utara, Semarang Tengah, Semarang Timur, Gayamsari, Genuk, and the northern part of Pedurungan are located in the study area (i.e., the northern part of Semarang City). According to field observations, Genuk, the northern part of Semarang Utara, Gayamsari, and Semarang Timur are experiencing land subsidence. The study area, as shown in Figure 2, spans 23.5 km × 10.5 km and is situated between 6°50′–7°10′ S and 109°35′–110°50′ E, near to the Java Sea.

2.2. Geological Setting

Geologically, the northern part of Semarang city is situated in an alluvial plain that is bounded to the north by the Java Sea and to the south by the Semarang highland. The north plain in this area is composed of a thick sequence of Holocene alluvial deposits (Qa) and the south highland comprises Quartenary volcanic rocks (QTd) and Tertiary sedimentary formation (Tmpk) [50]. The geological formations in the northern area of Semarang City are strongly affected by the sedimentation process and shoreline changes. The heterogeneity of subsurface is influenced by the sea-level fluctuation from the last glacial maximum until the Holocene transgression, as explained in Sarah et al. [47]. The rapid advancement of the shoreline indicates the rapid deposition of the Semarang alluvial deposit, forming soft clayey sediment on the upper part of the subsurface. In the eighth century, Semarang City originated from a cluster of small islands off the coast of Pragota (now Bergota) [51]. The coastline changes result from the sedimentation process. Old topographical maps revealed that the Semarang shoreline advanced 8–12 m/year from 1695–1940. Figure 3a shows that an 884 m advancement took place over 1847–1991 [52]. It appears that accretion of 884 m occurred over 144 years, implying that sedimentation in the northern part of Semarang City was extensive during this time period.

Figure 3b provides an evolution of the coastal accretion and retreat processes for changes in shoreline from 1984 to 2016 based on satellite images. The coastal shoreline in 1984 was accreting towards the sea. This phenomenon was influenced by natural sedimentation and intensive reclamation for land development that began in the 1980s [53]. The shoreline position in 1994 was similar to that in 1984. Slight advancement was found near the west canal (Kanal Barat), possibly due to an increase of sedimentation discharge from the headwater. An examination of the shoreline in 2004 showed that land subsidence may have affected the coastal morphology. Shoreline retreat was observed at the western and eastern sections. The combined influence of sea-level rise and land subsidence caused the seawater to transgress the land; thus, the west and east coastal lands are susceptible to inundation. This was not seen in the center part, because some flood protection systems were built for the city center. The dyke, canals, and polder pumping system help prevent the city from the multiple disasters of subsidence, coastal flood (rob), and seasonal flood during the rainy season. In 2016, the shoreline retreated further in the western part, implying a more aggravating impact of land subsidence and sea-level rise. The 2016 shoreline position in the city center remained similar to 2014 and only some slight advancements were observed near the west canal, Semarang Port, and east Kaligawe industrial complex. The small shoreline advances were due to sediment discharge from the west and east canals, as well as reclamation work in the port, trade zone, housing, and industrial complex.

Figure 4 presents a geological map of Semarang City, which is one part of the Semarang-Demak groundwater basin [50]. The lithological point and axes are elaborated in Figure 5. The northern region of Semarang City is mainly composed of alluvial deposits (Qa) due to the existence of a sedimentation process that also causes changes in the coastline. The thickness of the alluvial deposit ranges from 20 to 100 m, becoming thicker to the north and east (Figure 5). The bulk of the sediments in the basin consists essentially of grey, very soft clay, and silt, with abundant calcareous and shell fragments. Tons of sands and gravels are found at various depths from 5 to 80 m (Figure 5a,b). At the southern highland, thin clay and silt overlay alluvial fans of sand and gravel (Figure 5b). Apart from Qa, Semarang City has Pleistocene sediments (in the form of the Damar Formation, Qtd), which were generated by sedimentations due to variations in sea level in at least the past 500 years or approximately in the middle of the 14th century [54], caused by transgression and regression processes that form deltas and tidal deposits [45].

The stratigraphy in Figure 5 is divided into three units, clay to silty clay (Unit 1), sand lenses (Unit 2), and volcanic sandstone (Unit 3). The clay unit has a soft to medium consistency with an N-SPT value of 1–9. The sand lenses are medium dense (n-spt value 17–22), and the volcanic sandstone is dense to very dense (n-spt value >30). The geotechnical properties are derived from borehole investigations in Terboyo and North Semarang, and laboratory analysis [45]. Compressibility of the aquitard units was derived from 1-D oedometer tests [55] and hydraulic conductivity was obtained from the falling head permeability test [56]. Properties of the coarse-grained sediments are characterized based on their grain-size composition [57]. The geotechnical properties for the subsurface units are presented in

Table 1. The geotechnical properties.

Stratigraphic Unit	Natural Unit Weight (γn) (kN/m3)	Initial Void Ratio	Compressibility Index	Coefficient of Recompression	Modulus of Elasticity (E) (kPa)	Hydraulic Conductivity (k) (m/s)
		(e0)	(cc)	(cr)
Clay to silty clay (Unit-1)	15–17	1.2–1.7	0.35–0.74	0.10–0.17	1468–2000	1.68 × 10−10–6.54 × 10−9
Sand lenses (Unit-2)	17–19	1.5–2.0	0.41–0.77	0.11–0.18	4000–5000	1.59 × 10−6–5.03 × 10−5
Volcanic sandstone (Unit-3)	21–25	1.56–2.16	0.45–0.8	0.13–0.2	6900–7000	1.20 × 10−6–2.20 × 10−5

.

2.3. Hydrogeological Setting

The groundwater basin of the northern region of Semarang City, a part of the Semarang-Demak groundwater basin, is composed of unconfined and confined aquifers. The former (unconfined aquifer) is near the ground surface, with the groundwater table in direct contact with the atmosphere. The latter (confined aquifer) is separated from the unconfined aquifer by an impermeable barrier and is in compressed or semi-stressed status, consisting of lenses of sand and gravel that are covered by a layer of clay or sandy loam [53].

A hydrogeological evaluation to quantify interstitial fluid movements in aquifer layers needs a numerical representation of the groundwater flow regime. The hydrogeological dataset being inputted into our numerical models was collected from a geoelectric field study conducted in February 2019, as well as borehole works from a variety of sources, including the Balai Besar Wilayah Sungai (BBWS) Pemali-Juana, and Marsudi [53]. The lithology of the Semarang groundwater basin is illustrated in Figure 5.

The confined aquifer in the northern part of Semarang City is made up of two formations, i.e., the sandstone and conglomerates of the Damar formation along with the alluvial fan, as well as lenses of sand and gravel of the Garang deltaic deposit. The Garang delta was developed from the former river channels during the deposition of the deltaic deposit. In the Damar delta, groundwater flows from the volcanic rocks in the southern hills to the sedimentary basin in the north Semarang. The Damar formation has long served as a reliable source of fresh groundwater, mostly exploited by industries using deep wells at depths of 60–180 m. The Garang delta is also exploited by wells, but to a lesser degree due to its limited lateral extent. The numbers of registered wells and their output capacity are presented in Figure 6.

In Semarang City, field investigation reveals that regional groundwater flows from the south-southwest to the north-northeast through the deposits of gravel lens and sand lens. As a result, the groundwater level tends to fall to the north-northeast, with the conical decline’s center pointing towards LIK Kaligawe. Pump wells have been considerably used to extract groundwater in Semarang City and have two types, i.e., dug wells (shallow wells) and boreholes (deep wells). Boreholes, with depths ranging from 60 to 180 m, are positioned in a confined aquifer and often utilized for industrial purposes, whereas dug wells are typically installed by locals for their daily needs. In 1900, the number of pump wells in Semarang City was first reported, but these data are only available until 2010. Figure 6 shows the production capacity and the total number of groundwater wells [53,58,59,60], indicating that both the number and capacity of wells are increasing every year. Previous groundwater monitoring showed that the groundwater condition in Semarang city is already stressed. Over-exploitation of groundwater resources caused the formation of a cone of depressions that was first observed in 1984, but widened southward in 2010. Groundwater over-withdrawal also leads to the lowering of the piezometric pressure, thus increasing the effective stress of the aquifer system. However, due to the low permeability of the aquitard layer, the dewatering process is delayed, thus giving rise to gradual land subsidence. When the pore pressure decrease is smaller than the preconsolidation stress of the aquitard, the phenomenon of small, elastic, reversible subsidence occurs [61]. On the contrary, when the pore pressure drop is larger than the preconsolidation stress, it induces large, inelastic, irreversible subsidence.

3. Ground Deformation

3.1. Leveling and GPS Survey for Land Subsidence Monitoring

Land subsidence measurements were taken in the northern part of Semarang City by several institutes. Land subsidence was measured at 137 locations intermittently from 1991 to 2019 [62]. However, because the measurements (leveling and GPS) were carried out by different agencies and not synergized with each other, the measurement data at each point were not recorded continuously, but instead only for very short periods of time. Accordingly, in the current study, the data from the monitoring points of land subsidence that are close to each other, which are measured at different years, are used for calibration of the numerical model. Figure 7 depicts the time scale and period of leveling and the GPS survey points for land subsidence monitoring.

3.2. Interferometric Synthetic Aperture Radar (InSAR) Data and Processing

The most widely applied technologies for monitoring ground subsidence may be leveling and GPS, although their usage is still restricted due to a lack of geographic samples and expensive cost [63], which accordingly makes it difficult to undertake long-term monitoring of land subsidence. InSAR technology has been developed in recent years to tackle this issue [25]. Although this satellite-derived radar imaging has the potential to become a formidable tool for providing low-cost, continuous ground movement data over a vast area [64], the method of assessing ground deformation using InSAR image suffers from low accuracy because of the spatial-temporal decorrelation of the distance between objects and satellite orbits, high land cover complexity, and signal interference caused by air conditions [63,64,65,66].

DInSAR (Differential InSAR) is a technique for generating a large-scale map of line-of-sight (LOS) components using highly precise displacements [67]. A multitemporal deformation map, as well as many differential interferograms of the same zone from separate tempo acquisitions, must be considered. The short baseline subset (SBAS) and the persistent scatterer interferometric synthetic aperture radar (PSInSAR) are two sophisticated DInSAR algorithms that can be used to monitor the deformation of the earth’s surface over time. SBAS is carried out by using a mixture of differential interferograms produced by data pairs with a modest orbital separation (baseline) [68]. PSInSAR is based on phase characteristics and detects low-amplitude pixels with phase stability that are not found by conventional amplitude-based methods. It also uses spatially-correlated phases rather than historical phases so that the variables can be examined across time [69]. Consequently, the SBAS and PSInSAR algorithms are applied to address the problem of inaccuracies that may occur while using InSAR to detect land subsidence.

The InSAR data used in this study was taken from Sentinel-1 between 2015 and 2020 (https://comet.nerc.ac.uk/comet-lics-portal/, accessed date: 25 June 2020), with frame ID 076D_09725_121107 covering Central Java Province. The boundary of the area covers 6°50′0″–7°10′0″ S and 110°12′0″–110°37′0″ E, and more than 100 InSAR images were used in the analysis. An open-source program called LiCSBAS was applied in the present study to execute the InSAR time series analysis using LiCSAR products [70]. Generic Atmospheric Correction Online Service (GACOS) was used to correct atmospheric inaccuracies appearing in SBAS-InSAR [71]. Due to the correctness and coverage of unwrapped data, as well as loop closure control, inaccurate unwrapped data must be recognized and removed in the time series analysis. GACOS atmospheric products employ the interpolation technique of iterative tropospheric decomposition (ITO) to remove elevation-related and turbulent signals from the zenith total delay (ZTO) and then provide high-resolution tropospheric delay maps for InSAR and other data [71].

SBAS-InSAR is an analytical method for the multi-image InSAR time series to generate an estimation of the deformation of the earth’s surface by combining interferometric pairs from small time-space baselines. Land subsidence detection and monitoring using SBAS-InSAR has now become predominant. The analysis starts with N + 1 SAR pictures taken at ordered times (t₀, …, t_N) and is based on the same area. Assuming that each acquisition can interfere with other images, each SB subset thus must have at least two acquisitions, giving rise to the following inequality for the number of potential differential interferograms M [68]:

$$\frac{N+1}{2}\le M\le N\left(\frac{N+1}{2}\right)$$

(1)

Using the estimated generic j-interferogram of the SAR acquisition at times $t_A$ and $t_B,$ the topographic phase component removal in the azimuth pixels and coordinate range ($x,r$) can be described as follows:

$${\Phi }_j\left(x,r\right)=\Phi \left(t_B,x,r\right)-\Phi \left(t_A,x,r\right)\approx \frac{4\pi }{\lambda }\left[d\left(t_B,x,r\right)-d\left(t_A,x,r\right)\right]$$

(2)

where $d\left(t_B,x,r\right)$ and $\left(t_A,x,r\right)$ are the cumulative line-of-sight (LOS) deformations at time t_A and t_B with respect to the reference instant t₀ and λ is the transmitted signal’s center wavelength. Consequently, one can obtain $d\left(t_0,x,r\right)\equiv 0$, and it is fair to identify $d\left(t_i,x,r\right)$ as the required deformation time series, with i = 1, …, N. Assume $\Phi \left(t_i,x,r\right)$ to be the corresponding phase component; therefore, we have $\Phi \left(t_i,x,r\right)\approx \frac{4\pi d\left(t_B,x,r\right)}{\lambda }$.

Based on a series of displacement results, an SB inversion was performed on the interferogram network to estimate the velocity of a surface pixel over time. It is assumed that an M-unwrapped interferogram stack $d={\left[d_1,\dots , d_M\right]}^T$ was generated from N images acquired at $\left(t_0,\dots ,t_{N-1}\right)$ incremental displacement vector $m={\left[m_1,\dots , m_{N-1}\right]}^T$ (i.e., m_i is the incremental displacement between time t_i−1 and t_i), and can be extracted by solving Equation (3):

$$d=Gm$$

(3)

where G is an M × (N − 1) zero architecture matrix representing the interferogram network relationship with incremental displacements, and the unwrapped interferogram (difference between two acquisitions) is the sum of the corresponding incremental displacements [72]. Cumulative displacements (i.e., displacement time series) are calculated by adding the incremental displacements for each acquisition. The mean displacement velocity is then computed based on at least the quadrature of the cumulative displacements.

The NSBAS approach [73] was used, which imposes a temporary limitation to obtain a more practical time series of the displacement even with a disconnected network.

$$\left[\begin{array}{c}d\\ 0\end{array}\right]=\left[\begin{array}{c}\left[\begin{array}{ccccccc} & & \hspace{1em}G& \hspace{1em}& \hspace{1em} & 0 \hspace{1em}& 0\end{array}\right]\\ \left[\begin{array}{ccccccc}1& 0& \cdots & \cdots & 0& -t_1& -1\\ ⋮& \ddots & \ddots & & ⋮& -t_2& ⋮\\ 1& \cdots & 1& \ddots & ⋮& ⋮& ⋮\\ ⋮& & ⋮& \ddots & 0& ⋮& ⋮\\ 1& \cdots & 1& \cdots & 1& -t_{N-1}& -1\end{array}\right]\end{array}\right]\left[\begin{array}{c}m\\ v\\ c\end{array}\right]$$

(4)

A linear displacement (d = vt + c) is assumed if “γ“ is the scaling (weighting) element in the temporal constraints. The low time limit has little effect on the solution within the network’s linked components (e.g., 0.0001). As a result, the time restriction component only affects the connection via network gaps. Equation (4) can be used for pixels with fully connected networks as well as pixels with gaps.

Due to the limitation of interferometric phases, phase unwrapping was applied to quantify the absolute value of the phase with respect to a reference point inside the interferogram [74]. In the loop closure and mask time series phase, the mask was created utilizing multiple noise indices acquired at earlier stages for the time series and speed of displacement. If the value of any noise indicator for that pixel was greater than or less than the threshold set, that pixel was masked. The generated time series also includes several noise-related conditions, such as residual tropospheric noises, ionospheric noises, and orbital errors. A space-time filter can be utilized to isolate these components from the displacement time series (i.e., high-pass time and low-pass space) [75].

3.3. Land Deformation Mapping

High-precision subsidence mapping can be obtained using a satellite-based technique based on SAR using the DInSAR method. The DInSAR method is based on analyzing SAR images to identify surface changes to sub-centimeters along the sensor’s line of sight to the target, or Line of Sight (LOS), in order to calculate the LOS displacement value. LOS is the surface displacement between the satellite and the ground pixel in the azimuth direction along the flight path. The method is commonly used to estimate vertical displacement is to divide the LOS displacement obtained from DInSAR by the cosine of the incidence angle, assuming no horizontal movement occurs [76]. Furthermore, by using the least-squares algorithm, which is a matrix approach that aims to estimate the unknown parameters, during the inversion process, the displacement time series (in millimeters) can be obtained. After obtaining the displacement time series, a regression is performed to determine the annual average rate of land subsidence.

Before mapping land deformations, InSAR vertical displacement must be validated with GPS station measurements. There are several GPS points in the InSAR image frame, but many of them do not have enough temporal overlap with the InSAR 2015–2020 time series. Only the GPS points of K371, KOP 8, SMK 3, N259, BM11, 259, SFCP, and 1114 have temporal overlap with InSAR time series. Figure 8 depicts an InSAR validation for 2015–2020 with displacements monitored using GPS for 2015–2019. When compared to GPS monitoring results, it appears that in Figure 8, the InSAR displacement is in good agreement with an RMSE of 0.810 cm and R² of 0.949.

Figure 9 gives a map of the deformation velocity generated by the InSAR analysis in the northern section of Semarang City from 2015 to 2020. Inspecting Figure 9, it can be noted that the most significant deformation occurs in the north-east part of Semarang City, with a rate above 80 mm/year, whereas the west and southern parts of Semarang City are dominated by an uplift rate up to 20 mm/year. Figure 10 illustrates a time series of the deformation with and without a spatial time filter at the black dot in Figure 9 in more detail. The time series of the black dot in Figure 9 reveals that the magnitude of land subsidence was roughly 300 mm from 2015 to 2020.

A plot was constructed from the base map to determine the position of land deformation and uplift in northern Semarang City, creating a deformation map as shown in Figure 11. By inspecting Figure 11, one can observe that the highest deformation occurs in Genuk District, followed by Semarang Utara District and parts of Semarang Timur and Gayamsari Districts. Meanwhile, Semarang Timur and Gayamsari, as well as practically all of Tugu, Semarang Barat, Semarang Tengah, Semarang Selatan, Pedurungan, Gajahmungkur, Candisari, Tembalang, and Ngaliyan’s sub-districts experience the uplift status. The vertical uplift in the southwest part of Semarang City may be related to thrust faults in Semarang City, which is the part of the Baribis–Kendeng active fault [77].

4. Groundwater and Land Subsidence Numerical Model

4.1. Groundwater and Geotechnical Subsidence Equation

Groundwater and land subsidence behaviors have been evaluated extensively and quantitatively using numerical modeling. In the current study, the groundwater flow analysis was performed with the flow package in MODFLOW, while the land subsidence analysis was carried out with the SUB (Subsidence and Aquifer-System Compaction) package. The SUB package can be applied to simulate both elastic (recoverable) and inelastic (non-recoverable) interbed compactions with either no or delayed drainage.

The partial differential equation that combines Darcy’s Law with mass balance for describing a three-dimensional groundwater movement in MODFLOW takes the form:

$$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)-W=S_s\frac{\partial h}{\partial t}$$

(5)

where x, y, and z express Cartesian coordinates; K_xx, K_yy, and K_zz designate the tensor components of hydraulic conductivity in the x, y, and z axes; W represents the volumetric flux of water sources and (or) sinks per unit volume; S_s denotes specific storage; and h signifies hydraulic head.

The mechanism of vertical soil deformation induced by groundwater pumping is similar to that of soil consolidation acted upon by vertical load compression, in which the cause of consolidation is an increase in vertical effective stress and is then followed by a reduction in pore volume. However, in the former, excessive groundwater extraction leads to a decrease in pore water pressure and in turn raises effective stress, thus bringing about a decrease in pore volume [33,78].

According to Terzaghi’s effective stress concept, this can be mathematically written as:

$${\sigma }_{ij}^{\prime }={\sigma }_{ij}-{\delta }_{ij}u$$

(6)

where ${\sigma }_{ij}^{\prime }$ is the effective stress tensor component; ${\sigma }_{ij}$ is the total geostatic stress tensor component; ${\delta }_{ij}$ is the Kronecker delta function where ${\delta }_{ij}=1 if i=j or {\delta }_{ij}=0 if i\ne j$; $u$ is the fluid pore pressure; i and j represent the Cartesian coordinates x, y, and z. For 1D vertical stress problems, the equation reduces to:

$${\sigma }_{zz}^{\prime }={\sigma }_{zz}-u$$

(7)

The general equation that evaluates the thickness of compaction or expansion, Δb, between intervals t_n−1 and t_n can be written as

$$\begin{gathered} \Delta b=\frac{0.434b_0}{\left(1+e_0\right)\sigma \prime }\left[C_n\left({\sigma }_n^{;}-{\sigma }_{c,n-1}^{;}\right)+C_r\left({\sigma }_{c,n-1}^{;}-{\sigma }_{n-1}^{;}\right)\right] \\ {C}_n=\{\begin{array}{c}{C}_c,{\sigma }_n^{;}>{\sigma }_{c,n-1}^{;}\\ C_r,{\sigma }_n^{;}<{\sigma }_{c,n-1}^{;}\end{array} \end{gathered}$$

(8)

where e₀ is an initial void ratio, b₀ is initial thickness, σ′ is effective stress, ${\sigma }_{n-1}^{;}$ and ${\sigma }_n^{;}$ are the effective stresses at t_n−1 and t_n, respectively; ${\sigma }_{c,n-1}^{;}$ is the preconsolidation stress at t_n−1. The relationship of ${\sigma }_n^{;}$ to ${\sigma }_{c,n-1}^{;}$ is used to decide whether the value of $C_n$ is $C_c$ or $C_r$. The equation provides a quantitative estimate of settlements in over-consolidated sediments, those in normally-consolidated sediments, and those in sediments under transition from the over-consolidated state to the normally-consolidated state.

4.2. Groundwater and Geotechnical Subsidence Equation

The land subsidence process involves the hydro-mechanical coupling between the flow of groundwater and the deformation of the solid matrix in the aquifer system. Our models are rigorously established in a physically-consistent manner with the land subsidence monitoring results from GPS and InSAR. A predictive model for future development of land subsidence is then constructed based on groundwater usage regulations adopted by Lo et al. [49]. In contrast to Lo et al. [49], who used a grid measuring 250 m × 250 m to perform a numerical simulation of the groundwater model and its management, the simulation in this study used a grid measuring 100 m × 100 m. Furthermore, more aquifer parameter data were collected and input for modeling land subsidence. Since the northern part of Semarang City is part of the Semarang groundwater basin, a conceptual model for the groundwater and land subsidence model was made for the entire Semarang groundwater basin. Figure 12 is a schematic of the computational cell of a groundwater model for the Semarang City groundwater basin. The height of the cell is equivalent to the vertical depth of the soil layer at each site. Following the geological investigation presented in Figure 5, the model’s z-axis complies with the stratigraphic condition that is made up of three layers, consisting of (1) upper aquitard intercalated with (2) sand and gravel lenses, and (3) the Damar Formation at the base.

The boundary conditions for the numerical simulation of the Semarang groundwater basin are prescribed as (1) a no-flow boundary in some southern parts that reflect the existence of impermeable rocks and reverse faulting, (2) a groundwater divide boundary in the southwest part (extending from the south to the north), and (3) a constant head boundary in the north, west, and east parts due to sea and rivers.

The topographic data used in the numerical simulation are cited from DEM produced by DEMNAS Indonesia with a resolution of 0.27 arc-sec ≈ 8.325 m, while the data for the coordinates of the well are taken from the Mining and Energy Agency of Central Java Province, and the tidal data is the average tide of 1.07 m from the Meteorology and Geophysics Agency.

Semarang, similar to other Indonesian cities, has two seasons, i.e., dry and wet (rainy). The dry season lasts from April to September, while the wet season begins in October and ends in March. The average annual temperature in Semarang City is 28.08 °C, with a humidity of 76.61%. In Semarang City, annual rainfall ranges from 1500 to 3000 mm. Rainfall during the rainy season recharges groundwater and also causes river overflows and floods [79]. However, we did not include the effects of river overflow and flooding as specified head conditions due to insufficient data. Rainfall data recorded by the Indonesian Water Resource Agency, Public Work Department, are used for recharge input data from 1970 to 2010. Mangkang Barat/West Mangkang River, Mangkang Timur/East Mangkang River, Garang River, and Canal Timur/East Canal River are all included in the model and are prescribed as the boundary conditions (IBOUND). The piezometric level of the Damar aquifer in Semarang City drops rapidly and has surpassed the aquitard past maximum effective stress, and the demand for groundwater from the Damar aquifer was assumed constant during the dry and wet season.

The aquifer of Semarang City is separated into 20 hydrogeological polygons to input the hydrogeological parameters into the numerical model, as shown in Figure 12. The hydrogeological polygons separating sections were initially obtained by borehole and laboratory data in Table 1. The values of the elastic skeletal storage coefficient (Sfe) and inelastic skeletal storage coefficient (Sfv) are used for the no-delay interbed layer, whereas vertical K (hydraulic conductivity), elastic particular storage, and inelastic specific storage are integrated into the delay interbed layer. For modeling purposes, the initial properties were adjusted to satisfy the model calibration, as summarized in

Table 2. Hydrogeological model parameter.

Stratigraphic Unit	Hydraulic Conductivity (k) (m/s)	Sfe	Sfv
Clay to silty clay (Unit-1)	1.68 × 10−10–6.54 × 10−9	5.56 ×10−3	1.31 ×10−2
Sand lenses (Unit-2)	1.59 × 10−6–5.03 × 10−5	2.45 ×10−3	2.45 ×10−3
Volcanic sandstone (Unit-3)	1.20 × 10−6–2.20 × 10−5	1.42 ×10−3	1.42 ×10−3

. In this modeling, a delayed bed was assigned for Unit-1 as justified by its behavior as an aquitard. Unit-2 and 3 were assigned non-delay beds and Sfe and Sfv were set equal to establish elastic behavior.

4.3. Groundwater Numerical Model for the Steady-State Flow Model

The steady-state flow model was utilized to ensure that the hydrogeological input data are rational before they are next applied as the initial condition for the transient flow model. In the steady-state flow model, the values of appropriate input parameters for all grids, such as recharge and hydraulic conductivity, are determined to guarantee that the computed head in the model is able to match the observed head accurately. The steady-state flow model was performed based on the period of 1970 to 1990. Since the steady-state conditions do not take into account the temporal changes, the calibration of the steady-state process used the average observed data from 1970–1990 obtained from observation wells to assess whether the developed model is in agreement with field measurements. There are 54 observation wells in Semarang City. Unfortunately, much recorded data at observation wells have been lost. To this end, the model in this work was calibrated and verified using data from six observation wells, i.e., Prawiro Jaya Baru (O1), PRPP (O2), SMKN 10 (O3), Wot Gandul (O4), Santika Hotel (O5), and LIK. Kaligawe (O6).

To reduce the difference between the calculated and measured heads, a number of trial-and-error attempts were made to calibrate the hydraulic conductivity and recharge within the range of their theoretical values. The calibration results between the computed and observed values are depicted in Figure 13. It can be seen that values are close to the 45° line on the graph, mathematically pointing out a perfect correspondence between the computed head and the observed head. The colored bars represent the magnitude of the error at each observation well, with green color indicating that the error value is less than 1 m, yellow color indicating that the error value is between 1 m and 2 m, and red color indicating that the error value is greater than 2 m.

Figure 13 shows that the correlation value is very close to the 45° line at all observation wells. The value of Root Mean Square Error (RMSE) computed for all observation wells is 1.270 m. The individual discrepancy in water head for Prawiro Jaya Baru (O1), PRPP (O2), SMKN 10 (O3), Wot Gandul (O4), Hotel Santika (O5), and LIK Kaligawe (O6) are 0.618 m, 1.805 m, 1.812 m, 0.995 m, 0.980 m, and 0.597 m, respectively. Among them, the values at Wells O2 and O3 are relatively larger. Considering that the distance of the O2 and O3 observation wells to the boundary conditions is 0.29 km and 1.56 km, respectively, it is possible that the specified head boundary has a strong influence on the observation wells.

4.4. Groundwater and Land Subsidence in the Past

After the steady-state flow model was precisely calibrated, simulation and calibration were then undertaken for the model of transient flow and land subsidence to analyze the effect of groundwater flow due to pumping on land subsidence. The recharge rate utilized in the transient flow model was based on precipitation data from 1990 to 2010. Calibration of the transient flow model was conducted based on the period from 1990 to 2010, whereas the period from 2010 to 2020 was subsequently applied for validation. As achieved in the steady-state flow model, the computed groundwater level was thus compared with the observed level in the historical record of six observation wells from 1990 to 2010. The calibration of land subsidence magnitude was performed with the recorded data at the GPS measurement stations around the four observation wells (O2, O3, O4, and O6) since the other two observation wells are far from the GPS measurement points. Land subsidence data from 1990 to 2014 were utilized for land subsidence calibration. In addition to the GPS measurement, InSAR data were also used for calibration and validation. The use of InSAR can bridge the gaps between the leveling and GPS methods, and can provide the latest deformation because the satellite image is available until 2020. Figure 14 shows the subsidence contour map at the end of 2010 based on the simulation result, as well as the comparison of the computed and observed groundwater heads from 1990 to 2010. Figure 15 presents a graphic representation of the simulation result of land subsidence at O2, O3, O4, and O6, as well as the calibration with various geodetic measurements.

Figure 14 reveals that the most significant drop in land subsidence and groundwater level in 2010 occurred in Semarang Utara and Genuk sub-districts, whereas there was a vertical uplift in the southwest part of Semarang City. The values of RMSE for the calibration of groundwater level calculated in the transient flow model at O1 to O6 are 0.42 m, 1.91 m, 1.82 m, 1.76 m, 0.86 m, and 1.1 m, respectively, while the coefficient of determination (R²) takes the value of 0.6, 0.84, 0.91, 0.8, 0.75, and 0.98 at O1 to O6, correspondingly. One can also find from Figure 15 that the highest land subsidence in 2020 still took place in the northeast region of Semarang City. The RMSE values for calibrating the simulation of land subsidence at O2, O3, O4, and O6 are 5.88 cm, 7.34 cm, 7.32 cm, and 2.17 cm, respectively, while the R² values at these four observation wells are 0.91, 0.97, 0.99, and 0.99, correspondingly. This implies in the mathematical sense that the land subsidence model is in good agreement with field measurements. InSAR was also utilized to calibrate the land subsidence model when no land subsidence monitoring data from leveling or GPS was available, as done in Figure 15.

Figure 15 also reveals the fact that O2 (PRPP) at the northwest region of Semarang City and O6 (Kaligawe) at the northeast suffered from the greatest land subsidence. It can also be observed that the linear portion of the subsidence curves in Figure 15 represents the fast subsidence rate during 1990–2000. The subsidence patterns follow the groundwater withdrawal trend as observed in Wotgandul and Kaligawe. Different patterns are shown in the northwest area (PRPP and SMKN); the fast subsidence rate occurred from 1990 to 2000, but subsidence appears to decelerate afterward while groundwater levels were still falling. These differences are possibly due to the heterogeneity of the subsurface stratigraphy and hydraulic conductivity. In addition, a physical explanation behind this phenomenon can be explored by revisiting Figure 5, which shows that the sand lenses within the alluvial deposit are more intensive in the northwestern part. These sand layers act as drainage paths that accelerate pore pressure dewatering of the clayey soil, causing fast subsidence during the early period of a piezometric drop. Further decrease in piezometric levels thus creates smaller subsidence due to the dissipation depletion of excess pore water pressure. Therefore, land subsidence alleviated after 2010, partially because further lowering of piezometric levels did not surpass the soil preconsolidation pressure, in turn leading to a smaller subsidence rate. Preconsolidation pressures of the subsurface clay were measured by oedometer tests from undisturbed samples retrieved from borehole SMG-01 (northeast Semarang) and SMG-02 (northwest Semarang). The clay layers in SMG-01 are underconsolidated to normally consolidated, while the clay layers in SMG-02 are overconsolidated up to 40 m depth and normally consolidated at the depth of 50–60 m. The critical groundwater level required by the clay layer to undergo inelastic compaction was calculated from its preconsolidation pressure. It was shown that the piezometric groundwater level in the northeast has surpassed the critical groundwater level; therefore, inelastic compaction occurs. Meanwhile, in the northwest, the piezometric groundwater level is still above the critical groundwater level; hence, the subsidence is small.

Figure 15 also depicts that in Wot Gandul and Kaligawe, land subsidence was found to accelerate after 2010. This can be physically illustrated by the fact that groundwater piezometric levels have dropped beyond the recent preconsolidation pressure, thus causing fast, irreversible subsidence. A detailed examination of Figure 14 indicates that every 1 m drop of a piezometric level can result in an approximately −0.1 m vertical settlement.

To validate the simulation results, simulated subsidence was compared to subsidence based on InSAR analysis at points other than observation points. Figure 16 depicts sample points for detecting land subsidence using InSAR and simulated subsidence.

According to Figure 16, the simulated subsidence versus InSAR in the northern part of Semarang is close to the 45° line and has a small RMSE value (below 2 cm), whereas the RMSE values in the east-south area (Points 4 and 6) are quite large (above 2 cm) and quite far from the 45° line. The east-south compaction rate has a high RMSE value, which could be attributed to the onset of the subsidence process and the difference in time spans used in the two types of data sets in the area. The land subsidence simulation was carried out by taking a period of 30 years starting in 1990 and ending in 2020, whereas InSAR only examined the period 2015–2020. Based on the land subsidence simulation, the piezometric pressure in the east-south part of the simulation was drastically reduced at the start of the simulation. This is consistent with the piezometric contours of the Damar Formation aquifer [80], as well as the very large decrease in piezometric pressure of the Damar Formation aquifer until 1997, followed by a slower decrease [45]. However, when compared to the GPS measurement results for subsidence, the simulation results in the area have similar values.

4.5. Future Projection of Land Subsidence

This section concerns whether the problem of land subsidence will worsen if groundwater over-extraction continues in the northern part of Semarang City. A further assessment was thus conducted to predict the future development of land subsidence, after which the countermeasures for managing the problem of drops in regional groundwater level and resulting land subsidence are recommended. Lo et al. [49] have recently conducted a numerical study to provide future predictions of groundwater levels under different scenarios based on the implementation of various management measures. Their results show that among three proposed management measures (i.e., TS2, TS3, and TS4), a regulation strategy to reduce by 10% both the number and production capacity of deep wells from 2035 to 2050 has the most significant effect on groundwater restoration.

The same two scenarios (TS3, and TS4) for controlling groundwater pumping used by Lo et al. [49] were adopted and then performed in the land subsidence simulation in this study. The TS3 and TS4 scenarios were compared with the TS1 scenario without any management measure to quantify the effectiveness of these two management strategies on land subsidence. The TS1 scenario is a direct projection of land subsidence from the trend line (calculated in Lo et al. [49]) of the output capacity of deep wells from 2010 to 2050. The TS3 and TS4 scenarios are the representations as a consequence of a reduction in the number of deep wells and their production capacity by 5% and 10%, respectively, annually from 2025 to 2034, and then maintaining their number and capacity from 2035 to 2050. Figure 16 depicts the simulation results of land subsidence from 2010 to 2050 at O2, O3, O4, and O6. The contour maps of land subsidence in 2050 are presented in Figure 16 for TS1, TS3, and TS4, respectively.

A comparison of Figure 15 and Figure 17 reveals that the subsidence trend at PRPP (O2) and SMKN 10 (O3) subject to the TS1 scenario until 2050 takes a similar path to Figure 15, with rates of 10.8 mm/year and 11.98 mm/year, respectively. Higher rates of land subsidence were found at Wot Gandul (O4) (25.83 mm/year) and Kaligawe (O6) (36.8 mm/year). When the control measures (i.e., TS3 and TS4) are put in place, the subsidence degree at PRPP (O2) shows a rebound of around 10 to 27 cm at the end of 2050 as compared to 2020. A slowing down of subsidence rates from 2020 to 2050 can be also observed at SMKN 10 (O3) (26.03% and 51.42% under the TS3 and TS4 scenarios, respectively), Wot Gandul (O4) (27.27% and 53.03%, respectively), and LIK Kaligawe (O6) (27.87% and 59.08%, respectively). However, at PRPP (O2), we found that variations in effective stress due to changes in groundwater level do not surpass its preconsolidation pressure; therefore, the reduction in subsidence rates when the groundwater level is raised under the TS3 and TS4 scenarios would not be significant, while in other places, the rebound of the groundwater level is capable of reducing the subsidence rate.

Particular attention should be directed again to Figure 17 and Figure 18, where two groundwater management strategies are demonstrated to indeed mitigate the further development of land subsidence. A quantitative analysis based on these figures is presented in

Table 3. Reduction in the affected area of land subsidence.

Region	The Affected Area of Land Subsidence (km2)			Reduction in the Affected Area TS3–TS1 (%)	Reduction in the Affected Area TS4–TS1 (%)	Land Use
	TS1	TS3	TS4
1	0.47	0.12	0.12	−74.85%	−75.60%	Government offices, city-scale trade zone, and low-density housing
2	10.91	10.55	9.53	−3.30%	−12.64%	National-scale trade zone, city-scale trade zone, neighborhood-scale trade zone, industry, government offices, high-density housing
3	1.06	0.43	0.29	−59.55%	−72.38%	Sub-city-scale trade zone, neighborhood-scale trade zone
4	6.32	6.16	5.60	−2.47%	−11.28%	Industry, higher education, medium-density housing

, summarizing the results of the specified analysis for each scenario in Figure 15 and Figure 16. The affected area of land subsidence is divided into four regions: Region 1 is around Prawiro Jaya Baru (O1), Region 2 is around PRPP (O2), SMKN 10 (O3), and Wot Gandul (O4), Region 3 is around Hotel Santika (O5), and Region 4 is around LIK Kaligawe (O6). Among these regions, Region 2 has the greatest subsidence area. The percentage of decrease in the affected area (an area that has subsidence of more than 1 cm) of land subsidence is 2–75% if the TS3 scenario is taken, while the TS4 scenario is more effective with the percentage being 11–76%. Irrespective of the TS3 and TS4 scenarios, the greatest decrease in the affected area of land subsidence among all regions occurs in Region 1, with values up to 74.85% and 75.60% under the TS3 and TS4 scenarios, respectively. However, reducing land subsidence in Region 2 is critical and can yield the greatest economic benefits, since many commercial centers and government agencies, such as national-scale, city-scale, and neighborhood-scale trade zones, industry, and government offices, as well as high-density housing, are located in this region.

5. Conclusions

The northern part of Semarang City that is formed by alluvial deposits is prone to land subsidence, which has been exacerbated by over-drafting groundwater from deep wells in recent decades. To quantitatively analyze and model this problem, a systematic investigation of land subsidence caused by excessive groundwater abstraction was carried out in the current study using MODFLOW’s flow and subsidence packages integrated with the recorded data obtained from groundwater level observation, leveling, GPS, and InSAR. The precise calibration and validation of our model were achieved excellent agreement with field measurements, and it was then applied to project future developments of subsidence and groundwater level in 2050.

To mitigate the negative consequences of land subsidence in the future, two groundwater regulatory strategies are proposed. Our results indicate that the subsidence rate can fall from 26% to 59% in different regions with the introduction of these two measures. It is also shown that although Genuk sub-district (LIK Kaligawe) has a land subsidence rate of 36.8 mm/year, these measures can reduce the subsidence rate by up to 59%. The maximum reduction in the affected area of land subsidence was demonstrated to occur in Region 1 (around Prawiro Jaya Baru), where the affected area can be decreased by up to 75%. Since the northern part of Semarang City is experiencing land subsidence and is located in the coastal area, we strongly suggest that evaluation and prediction of the risk of coastal flooding due to land subsidence and global sea-level rise should be at the forefront of future research.