Implementation of Seismic Ground Response Analysis in Estimating Liquefaction Potential in Northern Thailand

It has been known that northern Thailand is an active tectonic region in Southeast Asia. Some earthquakes with low to medium magnitudes had occurred in northern Thailand. The M w 6.1 Mae Lao Earthquake occurred on May 5 th , 2014 in Chiang Rai Province. The earthquake also resulted in the unique phenomenon of ground failure, which was known as liquefaction. Learning from the event, the liquefaction potential based on seismic ground response analysis was performed. Several site investigations including standard penetration test and seismic down-hole test in Chiang Rai Province were carried out. The next generation attenuation model was conducted to generate the ground motion for nonlinear seismic response analysis. The peak ground acceleration at the ground surface from seismic ground response analysis was used to analyze the empirical analysis of liquefaction potential. The results show that liquefaction could occur at the investigated locations during the earthquake. The results also confirm the liquefaction evidence found in Chiang Rai Province during the M w 6.1 Mae Lao Earthquake. This research can help the people to consider the earthquake impacts to northern Thailand.


Introduction
It has been known that northern Thailand is one of active tectonic regions in Southeast Asia (Mase et al., 2018a;Tanapalungkorn et al., 2020). This is due to the fact that several active faults exist in this region. The endogen energy could trigger the movement of active fault which can trigger earthquakes in northern Thailand (Mase et al., 2020a). Within the last two decades, this area has intensively undergone earthquake events. The recent strong earthquake occurring in Mae Lao, northern Thailand in May 5 th , 2014 ( Figure  1) is widely known as the Mae Lao Earthquake (Mase et al., 2020c). This earthquake has not only triggered the structural building collapse, but also triggered a unique phenomenon called liquefaction on the Mae Lao Basin area (Mase et al., 2020c). According to Soralump et al. (2014), the liquefaction during the Mae Lao Earthquake is known as the second liquefaction eyewitness during the modern era of Thailand. Learning from the Mae Lao Earthquake in 2014, intensive studies of liquefaction were started.
Several local researchers had reported and studied liquefaction impact during the Mae Lao Earthquake. Soralump et al. (2014) reported the massive liquefaction had been found on the basin area of Mae Lao. In this area, the sand boiled and massive cracks occurred. Ornthammarath and Warnitchai (2016) and Mase et al. (2020c) studied the interpretation of Mae Lao Earthquake and structural damage. Ornthammarath and Warnitchai (2016) reported that ground motion of Mae Lao Earthquake approached 0.3g at the epicentre. It has also exceeded the threshold of peak ground acceleration (PGA) which can trigger liquefaction, i.e. 0.1g (Kramer, 1996). Mase et al. (2020b) conducted a study of ground motion parameters and resonance effects that occurred during the earthquakes in northern Thailand. Tanapalungkorn and Teachavorasinskul (2015) performed the analysis of liquefaction potential in northern Thailand by using the multispring element model proposed by Iai et al. (1992). Tanapalungkorn and Teachavorasinskul (2015) also adopted the seismic ground response analysis to observe the soil behaviour, especially related to maximum excess pore water pressure ratio (r u max), which is also recommended by several researchers, such as Mase (2017a). Results show that in northern Thailand, liquefaction could happen. In general, previous studies concerned on the simulation of seismic ground response analysis to estimate the vulnerability of liquefaction. However, implementation of seismic response analysis, combined with empirical analysis based on site investigation data, is still rarely found.
This study was performed to investigate the liquefaction potential in northern Thailand based on one-dimensional seismic ground response analysis. The parameter of PGA at ground surface obtained from seismic ground response analysis was used as the parameter for empirical analysis of liquefaction. In this study, the empirical method to estimate liquefaction potential proposed by Idriss and Boulanger (2008) was implemented. Factor of safety (FS) is presented in this study. In general, the results of this study could give a better understanding of the implementation of seismic response analysis and empirical analysis of liquefaction. Furthermore, the results of this study could lead the local engineers to consider the liquefaction damage in northern Thailand.

Studied Area and Geological Condition
This study is focused on several sites spreading in Chiang Rai Province, northern Thailand. They are noted as BH-1 to BH-7 ( Figure 1). Those are capital cities of districts in Chiang Rai Province, northern Thailand. Those investigated points are surrounding the epicentre of Mae Lao Earthquake ( Figure 2). BH-1 is located at Mae Sai, whereas BH-2 is in Mae Chan. BH-3 and BH-4 are situated in Chiang Kong and Mueang, respectively. BH-5, BH-6, and BH-7 are located in Mae Lao, Phan, and Wiang Pa Pao, respectively. For those seven sites, the site investigation data including the standard penetration test (SPT) and seismic downhole data were collected. The example of site investigation data from the sites is presented in Figure 2. In this figure, the site investigation at BH-5, the closest investigated points to the epicentre, was selected. In general, the investigated sites are dominated by loose to dense sands. At the shallow depth, loose sandy soils classified as SM based on Unified Soil Classification System (USCS) were found up to 9 m depth. This layer has (N 1 ) 60 average of about 11 blows/ft and shear wave velocity (V s ) of about

Empirical Analysis of Liquefaction Potential
Liquefaction is a unique phenomenon resulting due to earthquake. According to Das and Luo (2016), liquefaction on sandy soils happened due to the excess pore water pressure (∆u) triggered by earthquake shaking. Excess pore water pressure significantly rises up which decreases the effective stress. Excess pore water pressure is also known as the main parameter of liquefaction (Mase, 2017a). During the liquefaction, sandy soils behave as a liquid material in which all structures standing on sandy soil layers could sink and tilt. Several researchers had proposed the method to estimate liquefaction potential. The common method to estimate liquefaction is the stress equilibrium method. This method was originally proposed by Seed and Idriss (1971).
The main concept of this method is to compare cyclic resistance ratio (CRR) and cyclic stress ratio (CSR). CRR is defined as the ratio against liquefaction, which is provided by soil resistance itself. Cyclic stress ratio (CSR) is defined as the stress ratio resulting from earthquake shaking. Idriss and Boulanger (2008) proposed the empirical method to estimate the liquefaction potential in sandy soils. These two researchers mentioned that liquefaction could occur when the factor of safety against liquefaction (FS) is less than 1. The formulation to derive FS is expressed in the following equation: Idriss and Boulanger (2008) also proposed the empirical formulation to determine CSR that was modified from the equation of Seed and Idriss (1971). The empirical formulation proposed by Idriss and Boulanger (2008) is expressed in Equations 2 -4b: where: CSR is cyclic stress ratio (no dimension), r d is depth reduction factor (no dimension), MSF is magnitude scaling factor (no dimension) (Idriss, 1999), K σ is overburden correction factor (no dimension) (Idriss and Boulanger, 2008), PGA max is maximum peak ground acceleration (m/s 2 ), g is gravity acceleration (m/s 2 ), σ v is total stress.

The depth reduction factor (r d ) in Equation 5a
is expressed in these following equations: where: z is the depth of the investigated point. (2008)  (N 1 ) 60cs is a corrected standard penetration value normalized by clean sand effect (in blow/ feet). It is calculated by these following equations: ER is the ratio of energy efficiency, N is measured SPT (in blow/feet), and FC is fine content (in percent).

One-dimensional Seismic Ground Response Analysis
Method of one-dimensional seismic ground response was developed based on seismic wave propagation through horizontally layered soils. The framework of seismic ground response analysis had been presented by several researchers, such as Mase et al. (2018b), Hashash et al. (2016), and Mase et al. (2017a and2018a). It was addressed to solve several cases on geotechnical earthquake engineering. In general, there are two main models implemented in seismic ground response analysis, i.e. equivalent linear and nonlinear models (Yoshida, 2015). Nonlinear model was proposed to solve the limitation of equivalent linear. Hashash et al. (2016) proposed the pressure-dependent hyperbolic model, which was intensively developed by Hashash and Park (2001). This model focuses on hysteresis loop during cyclic loading which has a backbone curve defined as a hyperbolic function. Nonlinear analysis was performed by defining the discrete time increments in time-domain on lumped mass system by Hashash et al. (2016). There are some improvements from the first-generation model, which is implemented in a pressure-dependent hyperbolic model, especially related to determination of appropriate model of nonlinear soil behaviour. Hashash and Park (2001) introduced the reference shear strain (γ r ), which was correlated to referenced confining pressure (σ vreff ) and fitting parameters from laboratory tests. In addition, very small damping ratio of material was considered. This parameter correlates to the dependency of strain equivalent, which has an important role in seismic ground response (Laird and Stokoe, 1993). The main results of one-dimensional seismic ground response analysis include the time-history of ground motions, the frequency content of spectral acceleration, and the hysteresis loop of shear stress-shear strain. In this study, the ground response parameter, especially PGA at ground surface, was used as the parameter for the analysis liquefaction potential, as elaborated in the previous section.

Analysis Framework
This study was initiated by collecting the site investigation data in the studied area, i.e. in Chiang Rai Province. The site investigation data collected in this study are standard penetration test (SPT), boring log, and seismic downhole test. The collected data were then studied to obtain the description of soil profiles in the studied area. From the preliminary study, the suspected layers to undergo liquefaction can roughly be estimated. After the data collection and preliminary knowledge were obtained, the ground motion analysis was implemented.
The first step is calculating spectral acceleration and peak ground acceleration on each site using the next generation attenuation (NGA) models proposed by Abrahamson et al. (2014), Boore et al. (2014), Campbell and Bozorgnia (2014), Chiou and Youngs (2014), and Idriss (2014). Several parameters, such as fault type, earthquake magnitude, and epicentre should be determined. For the Mae Lao Earthquake, the magnitude (M w ) is about 6.1 and the fault type is slip strike fault (Mase et al., 2020c). Furthermore, the spectral acceleration was calculated and the largest spectral acceleration from the models was selected as the target matching spectra.
To generate the ground motion on each site, the spectral matching method was implemented. It is because no ground motion records were available at the investigated sites. This method was proposed by Al Atik and Abrahamson (2010). The spectral matching method was performed to derive the ground motions, which were relevant to the local site conditions. The matched ground motion used to generate the artificial ground motion is the acceleration recorded at the closest station to the earthquake epicentre, i.e. Mae Chan Seismic Station (MEAJ) (Figure 1). This ground motion was obtained from Thai Meteorological Department or TMD (2019).
The generated ground motions were then used as the input motions to simulate seismic ground response analysis on each investigated site. In this study, the pressure-dependent hyperbolic model (Hashash et al., 2016) was employed to obtain the soil behaviour description during the Mae Lao Earthquake. Results, such as time-history of ground motion and spectral acceleration, are presented. Furthermore, PGA at ground surface from seismic wave propagation was used as the earthquake parameter in liquefaction analysis. To observe the liquefaction potential under conservative conditions, the ground water level is simply assumed to be located at the ground surface. The factor of safety against liquefaction (FS) was also studied in this research. Figure 3 presents spectral acceleration generated from NGA models analysis. As presented in Figure 3, five NGA models were implemented to determine the spectral acceleration on each investigated site. The first model is Abrahamson et al. (2014) or ASK14 model, and the second one is Boore et al. (2014) or BSSA14 model. The other models are Campbell-Bozorgnia Campbell and Bozorgnia (2014) or CB14 model and Chiou and Youngs (2014) or CY14 model, respectively. The last used NGA model is Idriss' (2014). Those attenuation models have considered the uncertainty in earthquake engineering problems, such as the magnitude of earthquake, the local site condition, the fault type, and the distance to rupture (Mase, 2018 andMase, 2017b). Generally, BSSA14 resulted in the largest value of spectral acceleration on each site. Therefore, BSSA14 model can be used as the target spectral acceleration to generate the artificial ground motion for the investigated sites.

Next Generation Attenuation Model Analysis and Spectral Matching Results
The results of spectral matching analysis from BSSA14 model are presented in Figure  4. In Figure 4, by using the spectra matching method, the generated spectra acceleration was derived. Generally, the tendency of artificial spectra acceleration on each investigated site is relatively consistent with the referenced spectral acceleration. The artificial spectral acceleration was then used as the input motion on the seismic ground response analysis. The input motion was applied at the bottom of the soil column. In other words, the bottom of soil layer can be assumed as the engineering bedrock for each investigated site. Another reason is because the bottom of investigated sites has V s ≥ 760 m/s (Mase et al, 2018a), which is also indicated as the engineering bedrock surface (NEHRP, 1998). During the seismic ground response analysis, several parameters such as time-history of ground motion and spectral acceleration at ground surface were observed. The detailed explanation of spectral  Figure 3. Results of NGA model analysis.
acceleration and time-history of ground motions is elaborated in the next section. Figure 5 presents the interpretation of onedimensional seismic ground response results during the simulation of Mae Lao Earthquake. In this figure, two main results including the spectral acceleration and time-history of ground motion at ground surface are presented. In general, the spectral acceleration from seismic ground response analysis tends to amplify at ground surface. Spectral acceleration on each site also presents the peak value at short-medium period (T < 0.5 s). It indicates that the ground motion tends to be more destructive for low to medium high-rise buildings. Mase et al. (2018aMase et al. ( , 2020b in the study of ground response analysis during the Tarlay Earthquake also reported that the general pattern of earthquake impact in northern Thailand tended to give more impacts to the low-medium high-rise building. The propagated ground motions on the investigated sites tend to enlarge at ground surface. Generally, the propagated ground motions could amplify up to 2.8 times.

One-dimensional Seismic Ground Response Results
The propagated seismic wave result also shows that PGA max at several sites, such as BH-4, BH-5, BH-6, and BH-7 could reach 0.3 to 0.4g. The results of this study confirm the previous I J O G Indonesian Journal on Geoscience, Vol. 8 No. 3 December 2021: 371-383 studies performed by Mase et al. (2020c) and Ornthammarath and Warnitchai (2016) who mentioned that bedrock PGA max of Mae Lao Earthquake recorded Mae Suai Dam was about 0.3g. This PGA max value could amplify at the ground surface up to 0.4g at ground surface (Mase et al., 2020c). Based on the results and Kramer (1996), it can be roughly estimated that several sites having PGA max > 0.1g could be possible to undergo liquefaction. The detailed explanation for the liquefaction potential on each investigated site is presented in the next section Liquefaction Potential in Chiang Rai Province Figure 6 presents the interpretation of liquefaction potential on each investigated site during Mae Lao Earthquake in northern Thailand. This figure presents the liquefaction potential corre- sponding to the depth and soil layer on each site investigation. For BH-1, the subsoils tend to be safe from liquefaction. This is because FS on each layer has exceeded the liquefaction threshold, i.e. FS equal to 1. For BH-2 and BH-3, the similar trend as BH-1 was also found. No liquefaction indication on this site could be caused by the relatively lower earthquake impact on this area since PGA at ground surface is generally lower. In addition, the soil resistance on those sites is relatively higher; therefore, the liquefaction potential can be reduced. For BH-4, the liquefaction indication was found at the first sand layer, i.e. SC-GC. This layer has FS less than one. For the other sand layers, such as the second SC and the third SC layers, FS is larger than 1. BH-4 is relatively close to the epicentre. Based on the seismic ground response analysis, PGA at ground surface is relatively higher because the distance to the epicentre is quite close. BH-5 is also indicated to undergo liquefaction, especially at the first and second layers. At those layers, FS is less than 1. The site is very close to the earthquake epicentre. During Mae Lao Earthquake, those sites are predicted to have PGA at ground surface about 0.3 to 0.4g. For the other sand layers, a higher soil resistance tends to play an important role in reducing the earth-quake energy to generate liquefaction. For BH-6, the liquefaction is indicated to happen on each investigated layer, because FS on each sand layer is less than one. In terms of the distance, BH-6 is not very close to the rupture. However, the soil resistance at this site is relatively lower than the other areas. Seismic ground response analysis noted that PGA at ground surface is predicted to be about 0.424g. It could be also the reason why the investigated layers of BH-6 are very vulnerable to undergo liquefaction. For BH-7, the first and second layers are very vulnerable to undergo liquefaction in the studied area. BH-1 is also predicted to have PGA of about 0.348g, which has already exceeded the minimum liquefaction threshold. Similar to BH-6, BH-7 is also not close to the earthquake epicentre, but the liquefaction could be possible. This may be caused by the soil resistance provided by the investigated sites, so liquefaction could happen in BH-7.

Conclusions
This paper presents the analysis of liquefaction potential in northern Thailand during the M w 6.1 Mae Lao Earthquake in 2014. NGA models were implemented to estimate the spectral acceleration at ground surface. The spectral matching method was implemented to generate ground motion at the investigated sites for seismic ground response analysis. Several key results, such as time-history ground motion and spectral acceleration at ground surface were used as the parameters to determine liquefaction potential in the studied area. Northern Thailand, especially Chiang Rai Province, which is dominated by sandy soils at shallow depth, could undergo the earthquake impact such as liquefaction. It is indicated by the variation of PGA max the at ground surface where the studied area is very possible to influence the earthquake damage. The most impacted sites are generally located close to the earthquake epicentre. This is due to the fact that a near location to the epicentre tends to undergo more shaking impact. Liquefaction is generally found at the first and second layers on each investigated site. However, at deeper layers, the other sand layers are also possible to undergo liquefaction. This may be caused by low resistance of soil layers at deeper depths. Another cause may be influenced by very large peak ground acceleration resulting during the earthquake. The results of this study would be recommended to the local engineers to consider liquefaction impact in northern Thailand that can be used as the reference of seismic hazard mitigation in northern Thailand.