Design of Stone Column to Mitigate Soil Liquefaction: Cases Study of Yogyakarta International Airport

A low-to-medium cohesionless soil with low fines content was predominantly observed at the surfaces of Yogyakarta International Airport (YIA). The condition exposed subsoil of YIA to Liquefaction in addition to its location on a high seismic zone which has increased the likelihood of massive ground shaking. This means it is necessary to improve soil condition and vibro-replacement using stone column was selected as the appropriate method due to its recent popularity for the enhancement of sandy ground. Stone column has the ability to reduce the Cyclic Stress Ratio (CSR) of liquefiable soil and can be reliably evaluated. Therefore, this study was conducted to evaluate the risk of Liquefaction at YIA by adopting the SPT-based Liquefaction triggering procedure and presuming its manifestation using Liquefaction Severity Index (LSI). It is pertinent to state that the theoretical approach introduced by Priebe was used to design the geometry and center-to-center distance of stone column. The results were presented in the form of maps with a 50 m × 50 m grid size which include the cut and fill, LSI before and after improvement, stone column spacing, as well as stone column depth. It was discovered that the triangular spacing required for stone column ranged from 1.25 m to 2.5 m while the maximum depth was found to be 6 m. More-over, stone column inclusion efficiently reduced the severity of Liquefaction from medium to very low for the areas studied. However, stone column has several limitations and this means a combination of soil improvement methods needs to be applied to areas with moderate LSI.


INTRODUCTION
Yogyakarta International Airport (YIA) was found to be vulnerable to Liquefaction at different depths (Hartono, 2021).It was also discovered from the geology setting that the project site was underlain by alluvial sand deposits, has loose to medium consistency, low fines content with 3%, and a relatively shallow groundwater level, thereby, increasing its proneness to Liquefaction.Moreover, the existence of active faults in the vicinity increased its vulnerability to earthquakeinduced Liquefaction.This led to the application of the SPT-based Liquefaction triggering procedure to determine Factor of Safety (FoS) of native subsoil towards Liquefaction (Hartono and Fathani, 2022).Liquefaction Severity Index (LSI) was also adopted to measure Liquefaction's severity in terms of its manifestation.Furthermore, a previous study showed 64 out of the 84 boreholes data contributed to the risks that led to the mod-erate LSI on the west side as shown in Figure 1 (Hartono and Fathani, 2022).It was concluded that YIA generally requires soil improvement to mitigate Liquefaction risk and this was done in late 2018 using dynamic compaction (Setiaji et al., 2018).
The aforementioned study was initiated based on the native soil conditions in 2017 but the determination of design level is also mandatory before the construction phase and this is the focus of this study.The achievement of design level required a cut and fill process to induce the re-arrangement of YIA subsoil followed by another Liquefaction analysis and then an alternative design of stone column to reduce the risk of liquefaction.
Stone column has become the most popular technique to improve the sandy ground.The process Figure 1 The modification of the 2017 LSI map of existing conditions from Hartono and Fathani (2022) involves filling up vertical boreholes with gravel and compacting them using a vibrator (Castro, 2017).Its main function is to act as inclusions with higher stiffness, shear strength, and permeability than native soil.This means it has the ability to improve the bearing capacity, reduce the settlement, and dissipate the water during ground shaking.Practical cases have also shown the impact of stone column superposed with the densifying effect due to the probe vibration.The determination of the value of stone column added and not covered by dynamic compaction makes it important to design an alternate solution to improve the native soil using this method with the support of engineering tools such as GIS and Python.It is recommended that further studies focus on comparing the cost and performance of the two solutions in order to gain tremendous insights.

Data collection
YIA project is located on the Southern side of Kulon Progo Regency and situated at 396280 m E and 9126640 m N, thereby, leading to its categorization as zone 49 S based on UTM WGS-84 coordinate system.A total of 84 boreholes survey was conducted without any information regarding the surface elevation (PT.Nur Straits Engineering (NSE), 2017) and the drilling mostly ended up at 14 m below the ground surface due to the existence of a very dense layer on the last three samples of SPT.Moreover, the laboratory test report was also used to indicate the percentage of fines content for the According to the aerial photo presented in Figure 2, the existing condition of YIA was observed to be aquaculture during the borehole drilling in the middle of 2017.This shows that the existing surface had several elevations and this led to the retrieval of surface elevation data from DEMNAS in line with the EGM2008 datum due to the absence of topography data during the period of conducting this study.

Liquefaction analysis
Simplified procedure is the method widely recognized for liquefaction analysis (Seed and Idriss, 1971) and has also been applied to new cases histories (Boulanger and Idriss, 2014).Therefore, the SPT-based liquefaction triggering procedure was selected to analyze liquefaction hazard.It is important to note that Cyclic Stress Ratio (CSR) assumed soil column to be a rigid body moving horizontally in response to maximum acceleration by earthquake (Day, 2012).Meanwhile, soil column did not behave as a rigid body and was also found to be deformable, thereby, leading to the introduction of the depth reduction factor r d as follows where a max = maximum horizontal acceleration at the ground surface (m s -2 ), a max = acceleration of gravity (9.81 m s -2 ), α v0 = total vertical stress at (kPa), α ' v0 = effective vertical stress (kPa), z = depth in meters below the ground surface where liquefaction analysis was performed, and r d =1 -0.012 z.
The resistance of native soil was considered using the Cyclic Resistance Ratio (CRR).It was also important to adjust the observed N-SPT value by applying a reduction factor (Skempton, 1986).Moreover, the existence of fines content within the sandy soil had the ability to increase soil density and this led to the adoption of adjustments for fines content ∆(N 1 ) 60 .The resistance of the native soil toward the anticipated M W 7.5 earthquake (CRR 7.5 ) was retrieved using a deterministic chart or Equation 2. Accordingly, Magnitude Scale Factor (MSF) was used to account for the duration effects of loading cycles and amplitudes. (2) where (N 1 ) 60CS = corrected N-SPT values in accordance with field test and fines content and M = earthquake magnitude.Liquefaction generally occurs when the cyclic stress exceeds the cyclic resistance of native soil and this is the reason it was determined by measuring FoS as follows.
There are several limitations to the application of the SPT-based liquefaction triggering procedure.For example, FoS is not a practical parameter to prepare liquefaction severity maps (Sonmez and Gokceoglu, 2005) and this led to the introduction of liquefaction Potential Index (LPI) to measure the severity of liquefaction (Iwasaki et al., 1982).Meanwhile, it was impossible to divide nonliquefied areas based on LPI.It is also important to note that the moderate category of LPI consists of a broad range of values.Therefore, Sonmez and Gokceoglu (2005) introduced LSI with 6 degrees of susceptibility.FoS was associated with the probability of liquefaction P L (z) and the overburdened pressure of soil deposit was considered as a depth factor w(z) = 10 -0.5 z for z < 20 m, where z = depth in meters below the ground surface analyzed for liquefaction.A threshold depth of 20 m was used for liquefaction because previous studies indicated that liquefaction was less likely to occur beneath the value.

Stone column
The ground improvement method required a considerable number of column, thereby, indicating a complex modeling process of real geometry and this led to the application of simplified geometrical models (Castro, 2017).This study adopted the concept of a unit cell in axial symmetry due to its suitability and correspondence to the functional use of airport.It involved the uniform distribution of a great number of columns in a wide area under a uniform load.
Stone column are usually distributed in triangular or square grids with their influence, known as the tributary area, ideally shaped as a hexagon or a square as shown in Figure 3.The tributary area was transformed into a circle which was equal to d e = 1.05 -1.13 s for triangular and square grids respectively in order to allow axial symmetry conditions.The s is the center-to-center distance between column.
An adaptable method was developed to measure the performance of stone column based on a theoretical basis (Priebe, 1995).The limitations of this method include (1) the assumption of column to be a rigid layer and uncompressible and (2) neglect of  (Castro, 2017) the bulk density of column and soil, thereby, indicating the ability of column not to fail in the end bearing.There was also an initial pressure difference between column and native soil which led to bulging.This led to the introduction of the basic improvement factor (n 0 ) to evaluate column material towards shear force while surrounding soil reacted elastically using the following equations.
A c = 0.25πD 2 (10) where K ac is the active earth pressure (kPa), Φ c is the internal friction of stone column material (degree), D is the diameter of stone column (m) while A C and A is the area of a single column and the area of a unit cell respectively (m 2 ).In most cases, the application of the Poisson ratio (ms = 1/3) is adequate for the state of the final settlement.Moreover, the chart provided in Figure 4 can also be used to determine the n 0 .
Stone column material was practically compressible and this induced the increment of area ratio.Therefore, reduced improvement factor (n 1 ) was proposed by Priebe (1995) to address the compressibility of column material as indicated in the following equations. ) Frictional resistance of column was able to carry the external load and provided stabilizing effects on the nearby soil.Moreover, stone column had the capacity to reduce CSR in a seismic event and this led to the introduction of the reduction factor, α, as follows.

Calculation
A large dataset of soil parameters was used in this study and a back-calculation was required to derive the minimum spacing for stone column.This means it would be tiresome to use manual calculation, thereby, Python in geotechnical engineering practices was applied as the calculator, data handler, and visualizer (Yogatama and Adi Tirta, 2021).Subsequently, an open-source scientific environment written in Python known as Spyder was used to address several complex calculations.

Geographic Information System (GIS)
The spatial analysis using GIS assisted the engineers to delineate liquefaction susceptibility maps (Hartono and Fathani, 2022).This method has also been confirmed in previous studies to be a very useful tool especially when facing spatial variability.It is also important to note that this study focused on analyzing the cut and fill, planning of stone column's spacing and depth, and visualizing liquefaction analysis as a map.This led to the selection of open-source software known as QGIS to produce maps.Moreover, the delimited text data containing information were interpolated using Inverse Distance Weighted (IDW) interpolation as the geoprocessing tool.This was followed by setting the raster's size at 50 m × 50 m grid to ensure the practicality for construction work.
2.6 Design criteria

Design elevation
It was noted in Section 2.1 that the existing condition of YIA had several elevations ranging from +5 m MSL to +8 m MSL.For operational purposes, the filling process was required up to a certain elevation that can endure throughout design lifetime with due consideration for the settlement, land subsidence, and creep.This condition is commonly recognized as construction or fill level while the surface level after ground-lowering activities is known as design elevation.It is important to note that the settlement analysis was not part of this study, hence, a design elevation set at +7.4 m MSL with due consideration for the fluctuation of sea water level, surge, and safety from seawater overflow was used.

Design water elevation
Groundwater level position usually influences the effective vertical stress.This was observed from the theoretical belief that a shallow groundwater level reduces the effective vertical stress on each soil layer because of the water pressure generation.The condition usually leads to an increment in CSR value and a reduction in FoS.Therefore, the fluctuation of groundwater level was neglected and assumed to be situated at a certain elevation for design phase.The groundwater was set to be Hartono (2021) previously identified the subsurface condition of YIA in several cross-sections.Hartono and Fathani (2022) also stated that the low-to-medium sandy soil consistency was predominantly observed on the upper subsoil.Moreover, the layers underlain by very dense sandy soil are presented in Figure 5.The absence of index properties led to the definition of several parameters by correlation.It was discovered that the correlation by Bowles (1996) was adopted to estimate saturated unit weight of soil as indicated in Table 1.
Several sieve tests were conducted in 2018 and a quick review showed that many boreholes had typical grain size distributions and the average fines content (FC) was estimated at 3.0% and categorized as clean sand.

Fill material
A cut-and-fill process was conducted before soil improvement phase using fill material found within the project boundary.The saturated unit weight and fines content was based on the information presented in Section 2.6.3.Moreover, the filling material was compacted to satisfy the standard requirement and mitigate liquefaction in the middle.It is also important to note that the target N-SPT usually varies depending on the thickness but the value was set to 20 for the filling material used in this study.This study aimed to determine the minimum spacing between stone column to provide sufficient improvement factors to the subsoil.A back calculation was conducted in order to define several items of stone column before the analysis.The diameter was set to be 0.8 m as a common practice globally and as a correspondence to the minimum spacing between stone column which was set at 1.25 m with due consideration for the area replacement ratio.Moreover, the spacing was proposed in several categories with an increment of 0.25 m.

Seismic load
The return period of an earthquake is usually derived from SNI 8460:2017 (Badan Standardisasi Nasional, 2017) and depends on the type of infrastructure but the value for airport is not clearly discussed in the code.The use of a high return period is not feasible in terms of cost and improvement effort because YIA possessed massive risks of ground shaking.Therefore, a return period of 1000 years was used in this study and this corresponds to a 7% probability of exceedance for the 75 years design lifetime.
Probabilistic Seismic Hazard Analysis (PSHA) or Site-Specific Response Analysis (SSRA) was not available for liquefaction analysis.This led to the derivation of Peak Ground Acceleration (PGA) at the bedrock of 0.35 g from earthquake hazard map by National Center for Earthquake Studies (2017).Meanwhile, the boreholes data showed that the upper 30 m of YIA subsurface was classified as medium soil (SD) and the amplification factor was set at 1.25 based on SNI 8460:2017.The magnitude of controlling earthquake is usually derived from earthquake de-aggregation but there was none in this study, hence, historical data was used to determine the magnitude.The National Center for Earthquake Studies (2017) collected data on earthquake experienced near YIA and this was subsequently used to calculate the magnitude.The values recorded varied but the highest was found to have led to a lower MSF with subsequent reduction in CRR and low SF.Therefore, megathrust earthquake with an Mw of 7.2 were considered in developing the conservative stone column design.

Cut and fill analysis
Cut and fill analysis was conducted to estimate the volume of fill material and prepare the landside for soil improvement efforts.The DEMNAS data also showed the maximum cut and fill thickness to achieve design elevation of +7.4 m MSL were 0.6 m and 2.4 m respectively.The plan for the cut and fill works was presented as a map with a 50 m × 50 m grid area as indicated in Figure 6.
The layers of subsoil changed after reaching design elevation and this means additional soil investigations needed to be conducted ideally.However, it was assumed that the middle of the fill layer had the ability to contribute to Liquefaction analysis due to the absence of borehole data after the cut and fill process.The layer created additional vertical stress on soil beneath but the excavation area was favorable due to the reduction in the vertical stress.LSI of the existing condition was determined using the actual groundwater level to represent the bumpy surfaces as indicated in Figure 1.It was reasonable to explain the risks in the existing condition, especially in the year 2017, but this was not desirable to plan for soil improvement activities.The new LSI map before improvement was re-produced based on design criteria stated in Section 2.6.
The representative boreholes used for Liquefaction analysis were DB-11 and DB-13 located on the west and east side respectively.The cut and fill analysis showed that the DB-11 would be filled with sandy material while DB-13 was to experience the opposite.Table 2 presents Liquefaction analysis for representative boreholes.The fill material overlying the loose sand on the DB-11 was observed to be subjected to overburden pressure, thereby, increasing CSR value of the loose sand.FoS was also relatively low at this elevation and this led to a moderate LSI value.Moreover, the groundwater level situated at +5.40 m MSL made some parts of the fill materials behave unsaturated and were removed from Liquefaction analysis.Conversely, some parts of medium sand at DB-13 were excavated during the earthworks and this affected FoS of subsoil towards Liquefaction.This was due to the fact that there was a reduction in the height of the medium sand layer during the earthworks.LSI value only reached very low categories with no indication of liquefied layer and this was used as the starting point in designing stone column.It is also pertinent to note that soil A total of 84 boreholes were examined and only 9 were discovered to have sufficient FoS to resist upward seepage induced by sediment ejecta, hence, had an LSI value of zero.This was most likely noticed in the excavation area.Another 17 boreholes were also observed to have a very low LSI value and this means some of their layers had a FoS slightly below the requirement.Moreover, 42 and 16 boreholes had low and moderate LSI values respectively.The map provided in Figure 7 summarized LSI value before improvement of YIA and the west side was found to be more vulnerable to Liquefaction in line with the findings of previous studies in Figure 1 even though different criteria were applied.

Stone column design
The ordinary approach to achieve minimum FoS after stone column inclusion was through the determination of the spacing input by trial and error.However, this study applied back calculation to the minimum reduced improvement factor (n 1 ) to achieve an FoS of 1.1 It was possible to set FoS at different values but this was restricted by the spacing distance.The friction angle of stone column material (φ c ) was determined to be 45°while the other parameters were based on the information presented in Section 2.6.5.Moreover, stone  (Terzaghi et al., 1996) b Corrected N-SPT value c Unsaturated soil layers were removed in FoS analysis; following equation ( 5) column were planned to be distributed in triangular grids and this led to the definition of the tributary area as d e = 1.05 s.It is pertinent to note that the back calculation was executed in Spyder using scipy interpolate enhanced by enormous correlation data between s and n 1 as the library.Stone column was used for improvement up to the bottom of the liquefied layer (FoS < 1) and was automatically determined for each borehole by the Spyder.
LSI calculation results for the representative boreholes are tabulated in Table 3 and it was discovered that DB-11 required a 1.25 m spacing stone column from the surface elevation of +7.4 m MSL to +1.58 m MSL.Meanwhile, a reduction factor of 0.35 was applied on the loose and medium sand layer which led to an increment in FoS.The findings also showed that the loose sand did not reach an FoS of 1 because the spacing had a minimum distance of 1.25 m but LSI value was successfully reduced from moderate to low.This analysis was already conservative because the ability of stone column to dissipate the excess pore water pressure was not considered.In common cases, after a ground-shaking occurred, the fines particle within the native soil could be trapped inside stone column and this has the ability to interrupt the function of stone column as a drainage.This is the reason this function was neglected in this study.Meanwhile, the DB-13 was found to be safe from Liquefaction even without the inclusion of a stone column.Some other boreholes also did not require stone column and this allowed the engineer to avoid the higher cost of design.
The maps were produced in 50 m × 50 m grid size to portray the necessity of stone column inclusion in each area.Figure 8 shows that stone column spacing ranged from 1.25 m to 2.50 m with an interval of 0.25 m.Moreover, the analysis showed that several boreholes, including DB-13, did not need stone column and this led to the omission of several grids from the interpolation and map layer.(Terzaghi, Peck, and Mesri, 1996) b Corrected N-SPT value c Unsaturated soil layers were removed in FoS analysis; following equation ( 17) Figure 9 also indicates the depth of stone column to be embedded beneath the ground surface was mostly 4 m.

Performance of stone column
The inclusion of a stone column theoretically reduced CSR value.It was also discovered that the native soil was densified and stone column provided a reinforcing effect to the surrounding.Meanwhile, the secondary benefit that was not considered in this study was the drainage of excess pore water pressure during ground shaking.The results presented in Figure 10 showed that the application of stone column to improve soil shifted FoS of upper soil layers from the liquefied to the not-liquefied zone even though several soil layers were still liquefied.This is one of the limitations of using stone column due to its minimum spacing.It is important to note that several soil layers used in this study required < 1.25 m spacing to comply with the threshold of FoS ≥ 1 but this was not feasible.

Alternatives for soil improvement
Section 4.1 was used to discuss the limitations of applying stone column to mitigate soil Liquefaction.This was further confirmed by Figure 10 that  and Bray, 2019).Moreover, the existence of a thick non-liquefied soil layer near the ground surface has the ability to reduce the damaging effect of Liquefaction (Ishihara, 1985).This means the existence of liquefied layer at low to moderate LSI area after stone column inclusion does not imply the possible manifestation of Liquefaction at the ground surface.This insight served as the basis for geotechnical engineers to design the map provided in Figure 11.There is a maximum LSI determined based on engineering judgment usually allowed in a common project and later provided to the owner.
A situation mandating additional improvement requires using Figure 12 to examine the eligible improvement to be applied to the study area before or after stone column inclusion.The sieve tests conducted on the representative boreholes showed that YIA subsoil was generally classified as a sand material and this means dynamic compaction or vibratory probe is suitable to compact the loose sand layer before stone column installation.

Python for geotechnical calculation
The efforts to optimize soil improvement are often overlooked in general practices.This is due to the fact that contractors are primarily concerned with minimizing the total time required during soil improvement design.Therefore, Python was utilized to speed up calculations and minimize human errors in this study.This is necessary because the incorporation of optimization during design phase can lead to several benefits such as cost and material savings.However, some limitations were encountered when using Python to analyze Liquefac- This means it was merely an empirical approach based on historical data.Therefore, dynamic finite elements and effective stress analysis methods are necessary to provide a more realistic model.They can be used to produce a robust constitutive soil model in order to represent soil behavior at the site more accurately.Meanwhile, it is important to note that the application of these methods means Python is no longer useful for running the model.

CONCLUSION
This study showed that the cut-and-fill process affected the resistance of subsoil to endure earthquake-induced Liquefaction.Most of the excavated areas showed a very low to low LSI status while the areas with fill materials were at moderate levels as shown on the maps.It is perti-nent to note that the diameter of stone column was 0.8 m diameter, the center-to-center distances ranged from 1.25 m to 2.50 m with an interval of 0.25 m, and the maximum depth was 6 m.The installation of stone column reduced CSR value of the DB-11 representative borehole by up to 35%.However, stone column had several limitations associated with the inability to improve the minimum center-to-center distances that make FoS any further.This means the moderate LSI areas needed to combine other soil improvement methods.Therefore, the dynamic compaction or vibratory probe methods were proposed to be conducted before installing stone column in order to compact the upper loose sand layer.Dynamic finite elements and effective stress analysis were also found to be mandatory for future studies in the process of modeling the realistic soil behavior when the ground shakes and this means there will be no need for Python to run the model again.

DISCLAIMER
The authors declare no conflict of interest.

Figure 2
Figure 2 Aerial photo of YIA during soil investigation

Figure 3
Figure3Simplification of the unit cell to axial symmetry condition(Castro, 2017)

Figure 6
Figure 6 Cut-and-fill map

Figure 7
Figure 7 LSI map before improvement

Figure 8
Figure 8 Stone column depth map

Figure 9
Figure 9 Stone column length map

Figure 10 N
Figure 10 N-SPT, FoS, and LSI graph in a function of depth

Table 2 .
LSI calculation before improvement for representative boreholes a Cohesionless soil consistency was sorted based on

Table 3 .
LSI calculation after improvement for representative boreholes a Cohesionless soil consistency was sorted based on