The Using of GIS to Delineate the Liquefaction Susceptibility Zones at Yogyakarta International Airport

: Spatial analysis is performed to delineate liquefaction susceptibility zones at Yogyakarta International Airport(YIA). The low to medium cohesionless soil consistency is predominantly observed on the upper subsoil. A shallow groundwater level and low fines content have also enlarged the likelihood of earthquake-induced liquefaction. An SPT based liquefaction triggering procedure is adopted in this study to indicate the Factor of Safety ( FoS ) whereas the Liquefaction Severity Index ( LSI ) is used to measure the severity of liquefaction by presuming its manifestation. Inverse Distance Weighted (IDW) interpolation in QGIS is chosen to produce a map with 50 m × 50 m grid size. The analysis results show the YIA’s area is prone to undergo liquefaction at various depths. However, thin liquefied layers may not generate sufficient artesian flow pressure to eject water or sand. The LSI analysis concludes that YIA area is categorized as a non-liquefied to moderate severity where the West side is the governing area.


Introduction
During the Yogyakarta International Airport (YIA) planning stages, the assessment of the liquefactionprone area was a paramount task in geotechnical earthquake engineering. The YIA is located near the confluence of Eurasian and Indo-Australian tectonic plates thus it can be subjected to earthquake-induced liquefaction. Geological Agency of the Ministry of Energy and Mineral Resources [1] examined that the YIA soil deposits were vulnerable to be liquefied with severe manifestations. Liquefaction events generally occur in loose sand with low fines content below groundwater level [2]. Otherwise, the existence of an overlying low-kv layer typically with soil behavior type (Ic > 2.6) impedes the upward seepage and prevents its manifestations [3]. Moreover, a previous study [4] suggested that the existence of thick non-liquefied soil near the ground surface reduces the damaging effect if the underlying liquefied layer is thin. Subsequently, the Rockworks software will be used to generate soil lithology.
The analysis of liquefaction was obtained from the SPT based liquefaction triggering procedure by Boulanger and Idriss [5] as an update from former study of simplified procedure by Seed and Idriss [6].
The method identifies the liquefiable depth and the Factor of Safety (FoS) at each depth, thus it helps determine the suitable soil improvement method. Other findings by Sonmez [2] concluded that a factor of safety is feasible to indicate the condition of soil layer during and after the ground shaking whether it is liquefied or not. However, it is not a practical index to prepare susceptibility maps, thus Liquefaction Severity Index (LSI) was adopted to measure its severity. This index estimates the liquefaction severity by taking into account the probability of liquefaction as a function of depth. Subsequently, the variability of each soil layer and its vulnerability can therefore be described as a single index.
A spatial variation is usually discovered in the liquefaction-prone area, thus the Geographic Information System (GIS) becomes a very useful tool [7]. It helps interpolate the LSI attribute into areas or rasters to a desirable extent. Open sources QGIS software was used for this paper. Prior to this, several studies on various locations [8][9][10] conclude that the availability of a liquefaction susceptibility map help third parties to arrange the soil improvement plan. However, additional data collection for instance the Probabilistic Seismic Hazard Analysis (PSHA) and Site Specific Response Analysis (SSRA) should be conducted to enhance the results of liquefaction analysis as stated in Hartono [11].

Research Significance
This paper aims to perceive geological hazards associated with earthquakes. The YIA project is situated in an unfavorable location since the subduction and shallow crustal earthquakes have been recorded. The initial presumption indicates the existence of loose sand with low fines content that potentially generate an upward seepage and hence exacerbates the forthcoming hazards. Consequently, the predefined LSI method was performed to measure the severity of liquefaction and presume its manifestation. However, the results are less likely practical unless depicted in a form of map.
The open sources QGIS software will delineate the liquefaction susceptibility zones by utilizing the LSI attribute. As the aforementioned map for the assessment of ongoing project or soil improvement plan is unavailable, this paper is necessary.

Methodology
The YIA project is located at the Southern side of Kulon Progo Regency as shown in the aerial photographs of Figure 1. It is situated at 396280 m E and 9126640 m N thus categorized as zone 49S based on UTM WGS-84. Region of interest for soil improvement is composed of approximately 90 ha airside and landside zone [12].

Site and Subsurface Condition
The geological setting of the YIA was obtained from a map published by BAPPEDA DIY [13]. It is underlain predominantly by alluvial sands comprised of loose gravel, sand, silt, and clay that have been deposited by waterways in landforms.

Data Collection and Analysis
A geotechnical investigation was conducted in 2017-2018. It was manifested in 84 boreholes data and several laboratory tests as reported in the soil investigation reports [14,15]. The borehole drilling was mostly conducted up to 14 m below the ground surface due to the existence of very dense layers on the last three samples of SPT. The boreholes data did not contain any information regarding the surface elevation and reference level. As a result, this paper retrieved DEMNAS data for elevation information which refers to EGM2008 datum. The soil consistency of the YIA based on Terzaghi et al. [16] classified primarily as loose, medium, dense until very dense sand and the average groundwater level was noticed at 3 m below the surface level. In order to depict the soil lithology of YIA, Rockworks17 software was properly used.
The soil data consist of a large dataset of soil parameters thus the programming becomes a useful tool. Subsequently, the liquefaction analysis was calculated using an open-source scientific environment written in Python named Spyder. The use of Python [17] as a data handler was used in this work to enhance the analysis processes.
The spatial analysis was conducted to delineate the spatial information regarding which area is prone to liquefaction. The geoprocessing activities through Inverse Distance Weighted (IDW) interpolation method was chosen for this paper. The project boundary will be inputted as a vector layer hence LSI value will be set as an interpolation attribute. Moreover, the distance coefficient was applied as its default value whereas the size of raster is restricted to become a 50 m × 50 m square grid. Consequently, LSI map will become comprehensible for engineering practices and soil improvement planning.

Design Criteria
The design criteria for this project refers to the Indonesian National Standards and the supporting documents [18][19][20]. Futhermore, the site class of YIA was modeled as stiff soil (SD) in accordance with SNI 1726:2019 [19] by using borehole data provided in the soil investigation report [14]. The site class depicted the soil responses against earthquake ground motion, and hence determined the amplification factor of ground acceleration. However, the commonly used parameter for categorizing site class based on predefined code [19] solely incorporates the upper 30 m S-wave velocity (Vs30) and N-SPT data. As a consequence, when the distance of surface layer to the engineering bedrock does not extend up to 30 m, it may lead to an incorrect conclusion. Nevertheless, since engineering bedrock was not captured in the soil investigation reports [14,15] thus the SD class was then used in this paper.
A former study [21] defined the liquefaction zone based on the grain size distribution and uniformity coefficient (Uc). Cyclic liquefaction is associated with sediment ejecta which occurs in poorly graded and untreated cohesionless soil when subjected to seismic loading. Subsequently, a total of six boreholes samples were previously examined and shows a uniformity coefficient (Uc) of 3.9.
Index properties of subsoil throughout the depth were observed to be typical [11] and therefore the grain size distribution of soil samples at 6.0 -6.5 m depth from ground surface was plotted in Figure 2. The result indicates that the YIA subsoil is widely prone to be liquefied. Liquefaction analysis requires the number of fines content as the increment of ( 1 ) 60 that affecting 7,5 value. According to grain size distribution [15] the average fines content (FC) is approximately 3.0% thus it is categorized as clean sand.
Earthquake loading for YIA was set for 1000 years return period, and hence the infrastructure complies with 75 years design lifetime as well as a 7% probability of exceedance. Since neither the Probabilistic Seismic Hazard Analysis nor Site Specific Response Analysis was available, the Peak Ground Acceleration (PGA) and earthquake magnitude were obtained from reference [20]. Therefore the design criteria were modeled to indicate project susceptibility against liquefaction with PGA bedrock of 0.35g, amplification factor (SD) of 1.25, and Mw of 7.2.

Liquefaction Analysis
The SPT-based liquefaction triggering procedures have been developed over the years. It began in Japan and achieved the landmark work of Seed et al [6] which has been the standard in engineering practice for over two decades. Afterward, an update of CPT and SPT based liquefaction triggering procedures was developed by Boulanger and Idriss [5]. Several factors are speculated to govern the liquefaction process and its manifestation since the YIA soil deposits consist of medium sand and low fines content. A shallow groundwater table and the absence of a cohesive layer also exacerbated its processes. Day on [22] found that the hazard associated with soil liquefaction during earthquakes was also encountered in typical soil deposits. Consequently, when liquefaction occurs, the effective stress of the soil decrease whilst the pore water pressure increase abruptly and causes the soil to lose its ability to withstand any forces.

SPT Based Liquefaction Triggering Procedure
The SPT based liquefaction triggering procedure [5] was selected to determine the forthcoming liquefaction hazard as the recent development of former study [6]. The Dimensionless Cyclic Stress Ratio (CSR) values computed by assuming that a soil column of width and length was moving horizontally as a rigid body in response to the maximum acceleration exerted by the earthquake [22]. However, the soil column did not behave as a rigid body and most likely deformable. Consequently, the dimensionless values of depth reduction factor ( ) were introduced to comply with the aforementioned soil behavior [23]. The SPT-based case histories were reprocessed using the revised MSF relationship and showed a good agreement with the previous study.
Where = maximum horizontal acceleration at ground surface induced by the earthquake (m/s 2 ), = acceleration of gravity (9.81 m/s 2 ), 0 = total vertical stress at a particular depth where the liquefaction analysis was performed (kPa), ′ 0 = vertical effective stress at the same depth in soil deposit where 0 was calculated (kPa), = depth in meter below the ground surface where liquefaction analysis was performed, and = 1 − 0.012 z.
Soil resistance is derived from Cyclic Resistance Ratio ( ) values. Skempton [24] observed the different field test procedures affect N-SPT value used in the design. Furthermore, the consideration of confining pressure for soil layers beneath was also necessary as given in the formulas below.
for anticipated (M = 7.5) earthquake namely 7,5 was determined by Boulanger and Idriss (2014) to indicate soil resistance against ground shaking, the 7,5 was retrieved from the deterministic chart provided in reference [5]. Accordingly, the duration effects of loading cycles and amplitudes were applied by using the revised dimensionless Magnitude Scale Factor (MSF) obtained by combining (1) laboratory-based relationships and (2) correlations of equivalent uniform loading cycles with earthquake magnitude.
Where ( 1 ) 60 = dimensionless corrected N-SPT values in accordance with field test and fines content, = earthquake magnitude on the area of interest, FC = fines content of soil samples. Liquefaction will be observed if the cyclic stress induced by the earthquake exceeds the capacity of soil to resist. The Factor of Safety (FoS) is defined as a preliminary indication of liquefaction on the design level.

Liquefaction Severity Index (LSI)
Several limitations due to the use of simplified procedures were addressed by Sonmez [2]. To overcome some limitations of FoS, liquefaction potential index (LPI) and its severity categories were initially proposed by Iwasaki in 1982. However, nonliquefied areas could not be divided based on LPI. Furthermore, the moderate category of LPI consists of a broad range of values and thus, is uncertain although the low and high categories of LPI are properly defined. Sonmez [2] modified the ( ) term in the LPI equation by considering the threshold value of 1.2 between the non-liquefiable and marginally liquefied categories.
The different degree of susceptibility classes and nonliquefied area was introduced to express a severity map by using historical data as a starting point. The FoS is associated with the probability of liquefaction ( ). Overburdened pressure of soil deposit is also considered as a depth factor ( ) = 10 − 0.5 z for < 20m, where = depth in meter below the ground surface where liquefaction analysis was performed.
The severity classification terms given in Figure 3 were verified by observing the Chi-Chi earthquake in 1999 and revealed to have a good performance. Furthermore, the threshold of a depth of 20 m was used to take into account the liquefaction phenomenon where the previous studies indicated that liquefaction was less likely to occur beneath.

Soil Lithology
The long profile of YIA subsurface was created using Rockworks17 as given in Figure 4. Existing geological setting of Yogyakarta shows that the project location is composed of alluvial sand deposits. The observed borehole data also shows a typical cohesionless layer throughout the depth thus the soil lithology seems uniform.
Representative boreholes represent typical soil lithology layers that consist of very loose to loose sand at surface level up to the 6 m depth, medium sand that varies from 4 m to 8 m depth, dense sand from 6 m to 10 m depth, and underlaid by very dense sand that typically observed below the 8 m depth. The depth refers to a vertical distance from the ground surface to a certain elevation. The figure also shows that the West and East side consists of a relatively thick very loose to medium sand layer that is prone to liquefaction, whereas the middle side consists of a thin medium sand layer overlying the very dense sand. However, the deep groundwater elevation was observed predominantly in the east side, thus it creates a relatively thick unsaturated soil condition at the upper subsoil. This unsaturated layer is not affected by the increment of artesian water pressure thus the soil could maintain its minimum effective stress to cope with seismic loading and hence be categorized as a non-liquefied layer. The SPT captures soil as interval data thus the identified soil layers only represent certain depths. Consequently, the upper subsoil was not properly identified in this paper thus its resistance to deal with artesian flowinduced sediment ejecta can not be assessed.

Factor of Safety Analysis
The result shows a high likelihood of earthquakeinduced liquefaction on the YIA. Liquefaction potential analysis using SPT based liquefaction triggering procedure indicates that 64 out of 84 boreholes locations were prone to be liquefied at various elevations. It was determined by using the FoS = 1 declared in the previous section as a threshold. Analysis of two representative boreholes DB-31 and DB-38 are given in Table 1. The groundwater level was situated at 2.25 m and 2.06 m depth for DB-31 and DB-38 respectively. Unsaturated soils above the groundwater are categorized as a non-liquefied layer. Moreover, it shows that medium sand gives insufficient strength to resist the seismic load. On the other hand, dense and very dense sand provide sufficient cyclic resistance.
In order to portray resistance of each soil consistency against liquefaction, the SPT samples were summarized in Table 2. Loose soil was prone to undergo liquefaction due to its low resistance. However, it consists of both saturated and unsaturated soil since the average groundwater level was noticed at 3 m depth from the ground surface. The existence of unsaturated loose soil made several samples of loose sand categorized as non-liquefied samples. Moreover, the medium sand is foreseen to be the transition zone between liquefied and non-liquefied layers since it consists of both layers regardless of its saturated or unsaturated conditions. Factor that governs liquefaction on the medium sand is the overburden pressure. The existence of thick overlying soil gives sufficient pressure to the medium sand beneath that increase ( 1 ) 60 .

Liquefaction Severity Index Analysis
Liquefaction is associated with the increment of artesian water pressure and may eject the water or sand out of the soil layer. However, the thickness of liquefied layer governs the sediment ejecta occurrence. The amount of ejecta observed at the ground surface after the earthquake indicates the severity of liquefaction so-called LSI. Although 64 out of 84 boreholes locations were prone to be liquefied, only several boreholes with relatively thick saturated soil generate severe manifestations. Therefore, it can be seen in Table 1 that although two meters thick liquefied layer was observed in the 4 m depth, the LSI shows a very low status. LSI results show that there are 14 locations classified as non-liquefied zone dispersed mostly in the runway area. A total of 34 locations have a very low LSI status that dominated the most of YIA area. These locations may undergo liquefaction but presumably have a low manifestation at the ground surface based on historical data. Low LSI status was observed at the 29 locations both on the west and east side of YIA. The other 7 locations have relatively thicker loose sand and shallow groundwater level thus classified as moderate status. The average LSI value is 14.8 and elucidates that in general liquefaction will occur most likely with a low manifestation based on the historical data stated in [2]. However, the historical data is limited to the literature and it will be necessary to compare this work with observed liquefaction manifestation on the typical site. Figure 5 shows the compilation of N-SPT, FoS, and LSI results as a function of depth for all boreholes. The surface elevation of YIA project was obtained from DEMNAS data thus the depth of each borehole was substituted by the elevation towards geoid/ EGM2008. The YIA's ground surface is relatively flat with an approximately 1:500 slope. The results indicate that the liquefaction predominantly occurs at 1 m -4 m elevation whilst rarely occurring above 4 m and below 1 m elevations. It is also concluded that the top layers below the average 7 m ground surface elevation were dominated by unsaturated soil and hence considered to be non-liquefied. Several improvement methods outlined in reference [18], for instance the application of stone columns take into account the liquefiable depth thus the embedded depth of stone column could be defined appropriately using the results of this study. Improvement methods such as dynamic compaction are also suitable for  [12] d Total of SPT samples that are being analyzed to be liquefied or non-liquefied [23] cohesionless soil projects. The depth of improvement for the dynamic compaction depends on observed liquefiable depth. Required energy can be fulfilled by designing the drop height as well as the weight of pounder.

Spatial Analysis
The groundwater level map ( Figure 6) was used to depict the distribution of saturated soil thickness within the project area. It is found that there is a comparatively deep groundwater level observed in the east side of the YIA whereas other sides indicated shallow levels. Although Figure 6 gives several insights regarding the possible thickness of saturated soil, it should be compared with LSI map to conclude the conformity between them.
The distribution of LSI was illustrated in Figure 7. A relatively thick saturated very loose to loose soil observed at the west side of YIA shows the low to moderate LSI status thus it becomes the governing area. It gives a reasonable agreement with the existence of shallow groundwater level. Whereas the middle and east side of YIA give mainly very low LSI status that correspond to the thin dense sand overlying very dense sand at middle side and relatively deep groundwater level as shown in Figure 4.
A comprehensive map for engineering practices is necessary to enhance improvement works. As a result, the interpolated attributes was depicted in 50 m × 50 m raster size. Improvement works both using stone column or dynamic compaction as declared in analysis and discussion usually require grid area in order to monitor its progress or geotechnical instrument. The approval of improvement work can be taken based on each grid compliance towards certain threshold.   This paper does not take into account other factors that may govern liquefaction, for instance changes to soil properties during earthquakes. A clean sand layer with high fines content will be naturally compacted when subjected to a horizontal force as the result of a shaking motion. Moreover, earthquake-induced artesian water pressure also causes a sand boil that is observed at the relatively flat ground surface. However, it requires a dynamic finite element and effective stress analysis using robust constitutive soil models to execute. It can precisely back-analyze wellinvestigated liquefied sites and provide realistic models in conformity with the site condition. Therefore, for the next research, it is suggested to model the earthquake-induced artesian water pressure. Furthermore, liquefaction analysis using CPT-based method is desirable and has a number of advantages compared with SPT-based method used in this paper.

Conclusions
Spatial analysis is performed to delineate the liquefaction susceptibility zones at Yogyakarta International Airport which is located near the confluence of the tectonic plates. The project area underlaid alluvial sand deposits that have predominantly low to medium consistency in the upper 8 m. The grain size distribution shows that the cohesionless soil samples of the YIA project have a large uniformity coefficient (Uc) of 3.9 and situated in the liquefaction prone zone. Moreover, the existence of shallow groundwater level of approximately 3 m depth and low fines content governs earthquake-induced liquefaction. The liquefaction analysis results show that 64 out of 84 locations of YIA are prone to undergo liquefaction at various depths based on FoS values and since it did not presume the liquefaction manifestation, the LSI method was performed. YIA area is categorized as a non-liquefied to moderate severity status based on LSI. West side of YIA is the most governing area due to the thicker saturated very loose to loose sand observed beneath the surface. The LSI map depicts the vulnerability of YIA existing condition towards liquefaction and hence provides guidance for engineering practices. Due to several limitations of this paper, effective stress analysis using robust constitutive soil models is suggested to be performed for executing a realistic model in conformity with the site condition.