Natural Compaction of Semarang-Demak Alluvial Plain and Its Relationship to the Present Land Subsidence

Land subsidence is the lowering of ground surface due to natural and/or anthropogenic processes. Land subsidence in the Semarang-Demak plain has been going on for more than thirty years, however the contribution of natural and anthropogenic causes is relatively unknown. The Semarang-Demak plain has only been formed recently, as a result of rapid sedimentation during the Holocene. The basin mainly consists of underconsolidated thick clay, vulnerable to excessive settlement due to its own weight and additional pressures. The rate of natural subsidence is quantified by modelling the delayed dissipation of measured overpressure and the resulting vertical deformation, resulting in natural compaction rate of less than 0.8 cm/year in Semarang City and more than 0.8 cm/year in Demak Regency. The subsidence computed for parts of the Semarang-Demak plain were compared to the measured geodetic rate, and the relative contributions of natural and anthropogenic causes are derived. Modelling results show that natural subsidence is more significant at the eastern part of the plain (Demak region) with compaction rate reaching 0.9 - 2.2 cm/year that counts for 48 - 92% of the total land subsidence.


Introduction
Land subsidence experienced by major citiesacross the world has caused physical and material losses, and threatened to impede the sustainable regional development (Gumilar, 2013;Erkens et al., 2015). Land subsidence is defined as the lowering of ground surface due to subsurface movement. The hydrogeology of alluvial deposit or basin filled multi-aquifer-aquitard layers is prone to land subsidence, particularly when the aquitard is mainly soft clay (Bakr, 2015). Land subsidence can occur due to natural causes (e.g. natural compaction, tectonics), anthropogenic causes (e.g. overexploitation of groundwater, building loads, reclamation loads), and the combination of both (Gambolati et al., 1999;Hutasoit and Pindratno, 2004;Teatini et al., 2011;Zeitoun and Wakshal, 2013). Land subsidence due to anthropogenic factors has been widely studied, for example by Marsudi, 2001;Cui and Tang, 2010;Erban et al., 2014;Xu et al., 2016, and more. Meanwhile, land subsidence due to natural cause, such as natural compaction is seldom studied due to the difficulties in direct monitoring, and the lack of information on sediment properties and depositional history (Zoccarato et al., 2018).
Previous studies indicated that natural compaction of Holocene deposit is the main driver of land subsidence in the coastal areas (Meckel et al., 2006;Törnqvist et al., 2008;Stanley and Corwin, 2013;Zoccarato et al., 2018). Natural compaction is the process of reducing sediment volume due to the load of overlying sediment (gravitational compaction). Gravitational load of the overburden causes the underlying sediment to compact due to the expulsion of excess pore pressure and rearrangement of grain-to-grain contact (Skempton, 1970). When clay is deposited rapidly or reaches great thickness due to continuous sedimentation, pore pressure in excess of hydrostatic (overpressure) could be developed (Skempton, 1970 andGoulty, 2010). Overpressure condition has transient nature which would dissipate towards the steady state condition (hydrostatic). The delayed dissipation of overpressure can result in ongoing compaction long after sedimentation has stopped. The occurrence of overpressure is widely recognized in deep sediments (Ramdhan and Goulty, 2011;Satti et al., 2015, and more), while for shallow recent sediment some observations have been identified (Skempton, 1970;Lafuerza et al., 2008;Liu et al., 2014;Zoccarato et al., 2018). Knowledge of overpressure profile enables us to model the rate of natural subsidence due to delayed dissipation.
Uncompacted Quaternary sediments underlying municipal areas are prone to land subsidence, for example the alluvial deposits in the coastal cities of North Java and Bandung City underlain by lake sediments (Taufiq, 2010 andAbidin et al., 2012). Monitoring studies showed that the cities across North Java coastal plain are experiencing land subsidence at various rates, from 1 to more than 10 cm/year, such as in Jakarta (Abidin et al., 2011), Indramayu (Andreas et al., 2017), Semarang-Demak (Chaussard et al., 2013), and Surabaya (Kurniawan, 2011). Recent sedimentation formed the plain of Semarang-Demak at least 500 years ago (Bemmelen, 1941;in Marsudi, 2001), as reflected from the progradation of the shoreline towards the Java Sea (Tobing et al., 2000;Marfai et al., 2008). The rapid and recent sedimentation of the Semarang-Demak plain suggests that natural compaction of the subsurface sediment is still in progress. It is thought that the ongoing land subsidence results from the combined effect of natural and anthropogenic drivers. The land subsidence in Semarang-Demak has caused many adverse impacts such as a wider expansion of coastal flooding areas, and cracking of buildings and infrastructure . Understanding the relative contributions of natural and anthropogenic factors to land subsidence rate is important to assign the appropriate mitigation measure. This paper is aimed to quantify the proportion of the natural subsidence rates by numerical modelling, and the anthropogenic subsidence rate by comparing to the measured geodetic rate.

Geology and Hydrogeology
The studied area is located between coordinates of 110.36 o -110.64 o E and 6.92 o -07.01 o S, consisting of coastal lowland of the north Java alluvial plain and hilly area in the south. The studied area administratively belongs to the Semarang City and Demak Regency. The Semarang-Demak coastal plain consists of alluvial deposit (Qa) that lies unconformably on the top of Damar Formation (QTd), and the southern hilly area comprises Quaternary volcanic rock (QTd) and Tertiary rock (Tmpk) (Figure 1). The geological structures present in the studied area were formed during the Pliocene to Pleistocene, composed of thrust fault in the south and a strike slip fault parallel to Garang River (Thaden et al., 1996;Poedjoprajitno et al., 2008).
The Semarang -Demak alluvial plain was formed by the recent sedimentation as influenced by the sea level changes since the Late Pleistocene up to now (van Bemmelen, 1949;Schaeck Mathon, 1975). Transgression and regression of the sea formed the deltaic and tidal deposits of the Semarang-Demak alluvium. The lower deltaic deposit is made up of lenses of gravel and coarse sand deposited within the former channels of Garang River, and clay-silt and lenses of medium to fine sand deposited on the top of it. At the upper part, the deltaic deposit interfingers with the tidal deposit (Figure 2a and 2b). The tidal deposit consists of thick layers of calcareous and shell-bearing clay and silt, intercalated with lenses of fine sand. The thickness of the highly compressible silty clay to clay layer at the upper part varies from 25 to 60  Suwarti andWikarno, 1992, andThaden et al., 1996).

Explanation:
Tufaceous sandstone, volcanic breccia, conglomeratesk (Qtd) Electric conductivity (EC) (mS/cm) 600 Data Source: Geological Agency and LIPI reports Data Source: Geological Agency and LIPI reports m, becoming thicker towards the east (Demak). The upper clay layer plays a major part to the subsidence process in this area. Van Bemmelen (1941;in Marsudi, 2001) stated that muddy sedimentation occurred in the Semarang-Demak plain at least 500 years ago, as reflected from the advancement of Semarang coast line of 8 -12 m/year during 1695-1940. Investigation by Tobing et al. (2000 showed that the shore advancement still occurred during the years of 1847 -1991 as much as 884 m. These facts are in line with the historical proof that the Semarang port 300 years ago was further south in Bergoto and the width of the alluvial plain of about 3.5 km is in agreement with the shoreline advancement (Marsudi, 2001) (Figure 1). The sea level change in the Sunda Shelf (Saito et al., 1989;Yulianto, 2014) is referred to explain the sea level fluctuation in the studied area. At the last glacial age about 18,000 years ago, the sea level dropped to 60 -100 m, drying parts of the Java Sea. During this time, the gravels and coarse sands were deposited filling the Garang River channels (Figures 2a and 2b). For the period of 18,000 -17,000 years ago, the sea level rose gradually up to the current zero sea level at 6,500 year. During these times, transgression-regression occurred depositing alternating fine-grained and coarse-grained sediments. From 6,500 -4,500 years ago the sea level rose to the 3.75 -5.0 m above the present zero level. Marsudi (2001) estimated that the thick clay in the middle of Semarang City was formed during this period. The sea level gradually fell back to the present zero sea level during 1,500 -750 years ago, which is thought that the medium to fine sand lenses on the top of the clay sediment were deposited during this time. The sea level continued to drop until 0.5 m below the present level, where the beach ridges were deposited near shore. From 500 years ago up to now the sea level tends to increase reaching the present sea level, and during these times the upper clay filling the Muria valley to the Semarang plain was formed.
The recent and rapid sedimentation of the Semarang-Demak plain causes the development of pressures in excess of hydrostatic, particularly within thick clay dominated sequence, where the porewater pressure has not got sufficient time to attain equilibrium during rapid sedimentation. The occurrence of shallow overpressure prompts delay dissipation to happen, causing natural compaction of the sediment. The natural cause of land subsidence in Semarang-Demak plain is thought to be natural compaction. The contribution of tectonics to land subsidence is omitted, due to the unnoticeable role of the subsurface tectonic to the downward movement. A gravity study by Wardhana et al. (2014) indicated that the tectonic movement in the studied area has slowed down or even halted. It is therefore considered that natural driver of land subsidence in this area is attributable to the natural compaction of the recent alluvial deposit.
The hydrogeology of Semarang-Demak consists of unconfined and confined aquifer systems. The unconfined aquifer is formed near the surface, consisting of sandy silt, sand, and gravels with water table fluctuating to precipitation. The annual rainfall in Semarang-Demak varies from 1,500 -3,000 mm with wet season from November to April, and dry season from May to October. The water level of unconfined aquifer varies from 0.3 -3.0 m below the ground surface, and the dug well depths are from 6 -12 m (Marsudi, 2001;Sudaryanto and Wibawa, 2013). The confined aquifer in Semarang-Demak comprises Quaternary deposit aquifers (Garang delta aquifer and Quaternary marine aquifer) (Taufiq, 2009).
Transgression and regression due to the sea level fluctuation influenced the groundwater condition (Post, 2005), i.e. the variety of groundwater condition from fresh to brackish/saline as the result of depositional process. Measurement electrical conductivity in the aquifers and aquitards of the Semarang-Demak alluvial plain (Taufiq, 2010;Sarah et al., 2018) showed the variety of groundwater quality as shown in the subsurface profiles in Figure 2a and 2b. Both figures show that the electrical conductivity of Damar Formation is low (660 μS/cm), and so is the Garang delta on top of the Damar Formation (EC of 650 -756 μS/cm).
The aquifer and aquitards at the depth of 10 -35 m below the ground surface has high electrical conductivity (EC 3,950 -17,867 μS/cm) indicating the brackish-saline condition that was formed during the deposition (the Quaternary marine aquifer system). At the depths of less than 10 m, the electric conductivity in the aquifer and aquitard show low EC value at the east (Sriwulan and Gemulak), and a quite high EC value at the west (Bandarharjo) (see Figure 2a). The low EC value in unconfined aquifer in the east is possibly due to leaching during tidal deposition, while the high EC in the west is likely due to the contamination of seawater during tidal flooding. Figure 2b shows that at A2 site to Miroto, former beach ridge at the surface was formed that has fresh water quality (EC 1050 μS/cm). The former beach ridge is located 3.5 km away from the present shoreline, and was likely to be the former port of Semarang about 300 years ago. Van Bemmelen (1941;in Marsudi, 2001) stated that advancement of the shoreline towards the north during the year 1695 -1940 was 8 -12 m/year, which is in agreement with the distance of the former shore line (Miroto) to the present shoreline. Figure 2b shows that within the Quaternary marine aquifer system, from Poncol to Bandarharjo, the groundwater quality changes to brackish condition (EC 2200 -3950 μS/cm); the interfingering between the Garang delta aquifer and Quaternary marine aquifer is seen.
Exploitation of groundwater mostly involved the Garang delta aquifer and the Damar Formation aquifer of fresh water quality. The most widely exploited aquifer comes from the Damar Formation aquifer, where excessive groundwater utilization in industrial estates and densely populated areas shows the decrease of piezometric levels and is reflected in the depression of piezometric groundwater contours (Figures 3 and 4). Figure 3 shows the changes of piezometric levels in Semarang  Figure 4. It shows that the piezometric levels in 1984 were mostly positive, with some negative pressures started to form near the north shore, small cone of depression was also noticed in the southwest of the alluvial plain. Large cone of depression was seen in 2010, covering almost all of the studied area. A higher piezometric head is observed in the east part (Demak) and concentrated cone of depression is also observed in the south at the Mranggen industrial complex. The lowering of piezometric pressure adds extra pressure to the overlying aquitards, accelerating the naturally present compaction.

Methods
Methods employed in this study include acquisition of pore pressure data and modeling of pore pressure dissipation to compute the natural compaction rate. The measurement of pore pressure, governing equations, model setup, and material properties are explained as the following.

Pore Pressure Data Acquisition
The fast rate of sedimentation in the Semarang-Demak plain resulted in the occurrence of shallow overpressure (i.e. pore pressure in excess of hydrostatic condition). The identification of shallow overpressure was made directly by CPTu dissipation tests following procedures from Lafuerza et al. (2008) and Liu et al. (2014), and indirectly from laboratory consolidation tests on undisturbed samples retrieved from boreholes (Lafuerza et al., 2008;Long et al., 2011). Thirty six CPTu measurements and four boreholes (SMG01, SMG02, JTK, GEM) distributed across the Semarang-Demak plain were used in this study (Figure 1). Direct measurement of pore pressure by dissipation tests was made by halting the CPTu measurement at any depths and monitored the pore pressure variation with time until equilibrium pore pressure was achieved. Partial dissipation method was employed where the equilibrium pore pressure was obtained using 1/√t method (Flemings et al., 2008;Liu et al., 2014). The equilibrium pore pressure in excess of hydrostatic signifies the occurrence of overpressure, while the equilibrium pore pressure lower than and in line with hydrostatic represents the underpressure and hydrostatic state respectively. Since dissipation tests can not be carried out for all depths, indirect approach was made by using laboratory consolidation test data. Preconsolidation stress (σ′ p ) obtained from consolidation test was used to predict the overpressure condition using the equation (Lafuerza et al., 2008;Long et al., 2011): Equation (1) shows that the overpressure can be estimated as the difference between the present effective overburden pressure with the maximum past effective pressure (preconsolidation pressure). Preconsolidation pressure (σ′ p ) can also be estimated from CPTu data (Demers and Leroueil, 2002) as follows: u 0 : hydrostatic pressure (kPa) σ vo : overburden stress Studies by Lafuerza et al. (2008) and Liu et al. (2014) showed that the direct measurement of pore pressure by dissipation test and the indirect method by preconsolidation stress data resulted in similar results.

Governing Equations to Compute Natural Compaction Rate
Computation of the natural compaction rate was made by solving groundwater flow equations from transient condition (overpressure) to steady state condition (hydrostatic), and the resulting vertical deformation by one-dimensional consolidation theory. Numerical model of uncoupled consolidation was applied using finite element packages of Seep w/ (Geostudio, 2012) and Sigma w/ (Geostudio, 2013 where u a : atmospheric pressure atmosferik u w : pore pressure γ xy : shear strain ε x,y,z : strain at the x,y,z-direction τ : shear stress σ x,y,z : stress at the x,y,z-direction E : modulus of elasticity H : ν : Poisson ratio Analytical computation was carried out to validate the numerical modelling, using the formula for one-dimensional consolidation for underconsolidated clay (Day, 2010) : s c : consolidation settlement due to loading (m) c c : compression index e 0 : initial void ratio H 0 : initial clay thickness (m) σ′ vo : initial vertical effective stress (kPa) Δσ v : effective stress increase due to one-dimensional loading Δσ′ v : effective stress increase due to self-weight consolidation Time for consolidation was computed by (Day, 2009) :

Model Set-up
The numerical model was set up to match the overpressure profiles obtained from CPTu measurements at the maximum depth of 30 m. The subsurface stratification up to the depth of 30 m is dominated by clay intercalated with thin layers of sand ( Figure 5). The seepage model comprised 2501 nodes and 2400 quadrilateral mesh elements of 0.5 m x0.5 m size, and no flow boundary was set along the base and sides of the model ( Figure  5). Pressure heads were assigned at the centre of the model (Figure 5), the value of each head node represented the measured pore pressure (t 0 ) ( Figure 6). Groundwater modeling was carried out to dissipate the transient overpressure at t 0 to achieve the hydrostatic condition (t 1 ). The same geometry was used for stress-strain modeling, and 2542 nodes and 2440 quadrilateral mesh elements of 0.5 m x0.5 m size were used in this model (Figure 7). The displacements were set to zero in both the -x and -y directions along the bottom of geometric mesh, and along the vertical boundary of the geometric mesh (i.e. at both the left and right sides). The soil cannot move in the -x direction (U x =0), but is free to move in the -y direction. Along the exposed ground surface, the soil was free to move in both the -x and -y directions.
The engineering properties of the Semarang and Demak clays are influenced by the groundwater salinity, giving different compressive properties and hydraulic conductivity for freshwater and saline clays. The fresh to saline groundwater varies over depth, the material properties were applied according to its groundwater quality as had been outlined by Sarah et al. (2018). The material properties used in the numerical and analytical modeling were determined from laboratory test results as summarized in Table 1.

Pore Pressure Profile
Two typical pore pressure profiles have been identified in the Semarang-Demak plain, representing the thick and thinly overpressured layers. Figure 8 shows the pore pressure profile Dry unit weight, γ dry (kN/m 3 ) 9.0 -11.0 9.0 -11.0 Hydraulic conductivity, k (m/s) 3,35x10 -5 -1,09x10 -6 5,03x10 -5 -1,59 x10 -6 Modulus of elastisticity, E (kPa) 4,000 -5,000 4,000 -5,000 piezometric pressure of the Damar Formation aquifer below the alluvial deposit. The exploitation of Damar Formation aquifer results in the drop of piezometric pressure, as shown by the cone of depression in the north of Semarang City (Figure 4). The piezometric head drop of the Damar aquifer forces the overlying aquifer and aquitard to be in equilibrium with the aquifer pressure, inducing stress transfer from the pore water to grain-to-grain contact (i.e. subsidence) at the upper layers. This profile (Figure 8) is the typical pore pressure found in the north part of Semarang City, showing thinly overpressured layer. It is therefore evident that natural compac-tion due to the delayed dissipation of overpressure is hardly present in Semarang City, and the natural subsidence rate should be small. Figure 8 shows that the pore pressure profile up to the depth of 10 m varies from underpressure to slightly overpressure. The underpressure condition corresponds to the shallow aquifer pressure found at the upper part. Clay at the depth of 10 to 35 m is generally overpressured, some underpressures observed belong to the thin intercalating sands. The pore pressure profile from CPTu 32 site in Gemulak, Demak Regency (Figure 9) represents the characteristic of the pore pressure profiles at the east of the plain (Demak Regency). The thicker

I J O G
Natural Compaction of Semarang-Demak Alluvial Plain and Its Relationship to the Present Land Subsidence (D. Sarah et al.) 283 overpressure found in this region indicates that little anthropogenic stress is affecting the subsurface, and natural compaction is likely to occur.

Rate of Natural Compaction
The results of the seepage analysis are presented in the form of porewater pressure contours representing initial condition (overpressure) being dissipated into hydrostatic condition ( Figure 10). The seepage analysis results were imported to the stress-deformation model to compute the vertical displacement rate induced by pore pressure change. The finite element modeling was further validated by analytical model using Terzaghi 1-dimensional equation. The results of the natural compaction modeling presented in Figure 11, show compa-  rable results for finite and analytical models. The subsidence rate from the finite element modeling is slightly higher than the analytical method, most possibly due to the exclusion of the sand layer settlement in the analytical calculation. Figure 11 shows that the natural subsidence rate for CPTu13 site located in Demak Regency is about 1.74 cm/year, in which 90% of compaction would be achieved in 307 years. The compilation of natural compaction rates for Semarang-Demak plain is presented in Figure 12, showing various rates from less than 0.3 to 2.2 cm/year. It can be seen that the naturally driven subsidence rates are very small in Semarang City (< 0.3 cm/year), as the pore pressures are being depleted by rapid groundwater exploitation. Some sites in Semarang City show slightly higher rates of natural compaction, due to the less groundwater exploitation at the southwestern part and thicker clay sequence in the eastern part. As expected, the natural compaction rate is quite high in Demak region, varying from 1.1 -2.2 cm/year. The rate of natural compaction in Demak region is comparable to the natural compaction rate in other Holocene deposits such as in Mekong delta (2.0 cm/year) (Zoccarato et al., 2018) and Po delta (1,5 cm, 90% compaction achieved in 300 years) (Teatini et al., 2011).

Comparison of Natural and Anthropogenic Subsidence
Comparison of natural and anthropogenic subsidence was made by comparing the natural subsidence rate to the measured geodetic rate. The interferometric synthetic aperture radar (InSAR) results from Chaussard et al. (2013) was used to derive the major driver of land subsidence in all sites ( Figure 13).
Locations A and B in the north coast ( Figure  13) correspond to the Semarang Port and industrial complex of Terboyo (SMG01) that indicate the land subsidence rates of 7.8 and 6.0 cm/year. The natural compaction rates in areas A and B are 0.07 and 0.19 cm/year respectively, that indicate the natural compaction only contributes 1 -3% of the land subsidence rates. Meanwhile, location C in Figure 12 is located in the Mranggen industrial complex that corresponds to the high concentrated piezometric drop in the south (Figure 4). Therefore, it is presumed that the land subsidence in area C is entirely caused by anthropogenic forces.   Figure   The detailed comparison between natural compaction rate to the observed land subsidence rate is given in Table 2. Table 2 shows that the driving factor of land subsidence in Semarang City is mainly due to anthropogenic causes (91 -100%), while natural compaction contributes merely a minor fraction (1 -9%). Towards the east of the plain, the contribution of natural compaction to land subsidence increases from 28 -37% at the border of Semarang City, and much higher to 48 -92% in Demak Regency. The strong signal of natural subsidence at the east part of the Semarang-Demak plain suggests that further development of this area could result in even higher subsidence rate. The rates of land subsidence in Semarang-Demak are comparable to those taking place in Bandung and Jakarta. In North Jakarta and Bandung City, the average rate of land subsidence is about 7.2 cm/ year (Chaussard et al., 2013), which comes from mostly from anthropogenic causes. Taufiq (2010) calculated that anthropogenic factors contributed to 30 -70% of the total land subsidence rate in Bandung, while Pranantya et al. (2017) concluded that the anthropogenic factors could account for 30% of the land subsidence rate in North Jakarta. The contribution of natural factors to land subsidence in Bandung and Jakarta are not yet fully understood, hence it should be a point of interest for further research.

Conclusions
The ongoing land subsidence in the Semarang-Demak plain is caused by natural and anthropogenic factors. The rate of natural compaction was attained by modeling the 1-dimensional compaction due to delayed dissipation of overpressure. Quantification of the main drivers of land subsidence in Semarang-Demak plain has been achieved by comparing the geodetically measured land subsidence rate to the computed natural compaction. Spatial variability of land subsidence drivers is found across the Semarang-Demak plain where land subsidence in Semarang City is mainly governed by anthropogenic factor (91 -100%), and the natural compaction rate governs the land subsidence in the eastern part (Demak Regency) (48 -92%). The analysis presented in this study shows that the natural compaction of the Recent sediment has not been completed, therefore naturally driven subsidence will continue to occur for several decades to come. These results have significant implications for the regional development of this area, such that increasing anthropogenic drivers would speed up the already high subsidence rate. This scientific information is important for policy makers, as it would serve as a guidance for regional spatial planning and regulations, and should be used to prevent further environmental damages due to the impacts of land subsidence. Further research on the natural compaction of Recent sediments shall be extended to other areas in Indonesia, such as other cities on the north coast of Java, East Sumatra, and Kalimantan to understand the mechanism and the rate of natural compaction.