Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6, 3816-3831 2024 Publisher: Learning Gate DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate © 2024 by the authors; licensee Learning Gate * Correspondence: y5708211@ms18.hinet.net The source characteristics and ground motion parameters measure with structural damages of Mw6.4, 2018 Hualien earthquake Taiwan from GPS data Boi-Yee Liao1*, Yue-Tan Zhang2, Sen Xie3, Wen-Hua Chi4 1Institute of Engineering Technology Management, International College, Krirk University, 3 Ram Inthra Rd, Anusawari, Bang Khen, Bangkok 10220, Thailand; y5708211@ms18.hinet.net (B.Y.L.). 2Department of intelligent Engineering, Guangzhou Nanyang Technological Vocational College, 1123 Huanshi East Road, Conghua District, Guangzhou City, Guangdong Province, China; sciencechn@foxmail.com (Y.T.Z.). 3Department of Business Administration, International College, Krirk University, 3 Ram Inthra Rd, Anusawari, Bang Khen, Bangkok 10220, Thailand; xie.sen@email.krirk.ac.th (S.X.). 4Department of Early Childhood Education, International College, Krirk University, 3 Ram Inthra Rd, Anusawari, Bang Khen, Bangkok 10220, Thailand; m0423@ms12.hinet.net (W.H.C.). Abstract: On February 6, 2018, a moderate earthquake with a magnitude of ML6.2, known as the 0206 earthquake, struck Hualien in eastern Taiwan, resulting in significant destruction, 17 fatalities, and over 300 injuries. This research investigates the kinematic source model of the earthquake by inverting coseismic Global Positioning System (GPS) displacement recordings from the area surrounding the epicenter. The inverted source model indicates that most asperities above the hypocenter are minor, with large slips occurring at depths shallower than 10 km, which contributed to the severe damage in Hualien City. Additionally, we calculated the Coulomb stress changes at depths ranging from 6 km to 14 km, revealing that most aftershocks within a month occurred within a Coulomb stress range of >0.1 bar, extending from shallow to deep and from southwest to northeast of the epicenter. Notably, despite high Coulomb stress, few aftershocks were generated in the western part of the epicenter, suggesting a seismic gap that may lead to a more significant earthquake in the future. Two years later, two moderate earthquakes (ML=5) occurred at depths of 17.53 and 12.1 km in the region where Coulomb stress had increased, validating our predictions. Furthermore, by applying the derived seismic source model and a stochastic semi-empirical technique to assess five parameters related to building damage, we found that all damaged buildings were within the danger alert level, with some located along the Meilun Fault, indicating that the destruction was due to the combined effects of the 0206 earthquake and the fault's induced movement. Keywords: Coulomb stress change, GPS, Semi-empirical stochastic method, Source model, Structural damage parameter. 1. Introduction Taiwan's eastern region is known to be prone to destructive earthquakes due to its location at the convergence of the Philippine Sea Plate and the Eurasian Plate [1]. Unfortunately, this has resulted in severe losses in both lives and property. The high-risk seismic zone in eastern Taiwan is a cause for concern. Earthquake activities focus on plate boundaries, which only add to local society and infrastructure challenges. In early February 2018, a series of earthquakes occurred off the coast of Hualien. At least 16 earthquakes with a magnitude greater than 4.5, including one magnitude of 6.1, distributed in the depth range of 3 to 15 kilometers [2]. The major earthquake with a magnitude of Mw 6.4 struck the northern coastal area of Hualien City on February 6, 2018, resulting in 17 fatalities, over 300 injuries, and the collapse of four more significant buildings [3]. From this point forward in the paper, the main seismic event will be the "0206 earthquake." In recent years, many researchers have conducted in-depth investigations into various observed data to understand this earthquake https://orcid.org/0009-0002-3393-1701 3817 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate comprehensively. For example, according to [4], the 0206 earthquake was likely caused by three faults, with the most significant impact on the west-dipping fault. This event eventually triggered the shallower Meinong Fault, leading to surface rupture. Wen et al. (2019)[5] combined the teleseismic body waves with GPS observations to investigate the kinematic source of the 0206 earthquake. After extensive computational analysis and fitting with actual observational data, they believed that two complex fault models caused the 0206 earthquake. They suggested that the mechanism of moderate- sized events that occurred in eastern Taiwan is on the west-dipping fault. Based on on-site geological and topographical surveys [6] and apparent source time function [7], it was concluded that the rupture of the 0206 earthquake propagated from the north to the southern Milun Fault. This result triggered the movement along the Milun fault and resulted in structural deformation, causing severe damage to buildings. In addition, the aftershock distributions offered some insights related to the structural characteristics. As [1] indicated, the aftershocks affected an area that extends from a depth of 5 to 15 kilometers towards the southwest of the epicenter and has reached the Taiwan Longitudinal Valley. Meanwhile, it is interesting to note that the 0206 earthquake occurred at the precise location with a significant change in the seismic velocity and the Poisson's ratio in the horizontal direction. Further analysis revealed that this location corresponds to the Taiwan Longitudinal Valley fault and its branches [8]. By employing various actual data and methodologies, novel questions regarding the mechanism of the 0206 earthquake have gradually come to light, contributing significantly to our understanding of the seismic structures and mechanisms. However, these findings typically take some time to fully manifest after an earthquake event. As Hakami et al. (2013)[9] mentioned, inadequate disaster preparedness measures can result in a second wave of casualties and property losses when a disastrous earthquake occurs. In earthquakes, according to the findings outlined in Chiu et al. (2020)[10], the survival rate of victims is not only inversely proportional to time but also rapidly diminishes. When an individual is trapped for five days following an earthquake, the reported survival rate can be as low as 6%. Therefore, this study explores the possibility of using GPS data to rapidly calculate corresponding main parameters shortly after an earthquake, providing valuable insights for disaster prevention and mitigation efforts. There have been some excellent studies recently, for example, Alif et al. (2021)[11], Duan et al. (2022)[12], and Peidou et al. (2024)[13], that have utilized GPS data for seismic source inversion. These studies have led to many significant achievements and have helped advance our understanding of earthquakes and their effects. Furthermore, extracting meaningful parameters from seismic records to gain a deeper understanding of seismic characteristics and subsequently assess the damage caused by earthquakes to buildings has always been a direction actively pursued by engineering seismologists. Liao et al. (2020, 2022) [14, 15] used seismic spectral intensity parameters (SIs) to assess earthquakes with magnitudes greater than 6.0 occurring at various locations in Taiwan. They found that SIs could effectively evaluate the damage to buildings of different heights caused by earthquakes. According to Campbell and Bozorgnia (2010)[16], using the PEER-NGA database to formulate Ground Motion Prediction Equations (GMPEs) for Cumulative Absolute Velocity (CAV) is an effective method for assessing potential structural damage following an earthquake. CAV serves as a robust parameter for swift damage evaluation, facilitating prompt responses to likely hazards to building safety. Based on the maximum deformation energy stored within seismic waves, Danciu and Tselentis (2007)[17] found that the characteristic intensity (Ic) is a valuable indicator associated with seismic damage. Massumi and Gholami (2016)[18] have suggested that the role of Ic in the inter-story drift of buildings is significant, which makes it applicable for seismic damage assessment of buildings with different periods. It is interesting to note that Zhai et al. (2007)[19] conducted a study on the impact of Maximum Incremental Velocity (MIV) on inelastic displacement ratio spectra (IDRS). According to their findings, the influence of MIV tends to decrease as the period ranges of structures increase. This fact highlights the significance of MIV as a critical parameter that needs to be considered when designing buildings and structures to withstand seismic activity. Recently, Liao et al. (2023)[20] integrated Seismic Intensity (SIs), Characteristic Intensity (Ic), Cumulative Absolute Velocity (CAV), and Maximum Incremental Velocity (MIV) to assess the seismic damage caused by the 2016 Meinong earthquake (ML=6.4) in 3818 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate Taiwan. They discovered correlations between damaged buildings, the critical thresholds of these parameters, and between the parameters themselves through regression analysis. This study is structured into three key steps. The first step involves the use of GPS displacement data to invert the seismic rupture process. Subsequently, the seismic source parameters, the corresponding parameters of strong ground motion in the respective area, and the Coulomb stress changes are calculated. Finally, these results are integrated with building damage assessments to evaluate the accuracy of the calculations. This structured approach ensures a comprehensive and accurate evaluation of the seismic damage caused by the 0206 earthquake. 2. Methodology 2.1. Source Inversion One surface observation point in a homogeneous medium environment will experience displacement due to an underground fault slip. Assuming the fault plane is divided into M*N sub-faults, the surface displacement Ui at one observation point contributed by all the sub-faults can be represented as [21] Ui = 1/F ∬ Δuj[δjkΣ ∂ui n ∂εn + μ( ∂ui j ∂εk + ∂ui k ∂εj )]νkdΣ. (1) where ui(x1, x2, x3) is the displacement vector at observation point, Δuj(ε1, ε2, ε3), is the dislocation vector across a fault surface Σ, ν𝕜(ε1, ε2, ε3) is the normal vector to the surface Σ, ui j (x1, x2, x3; ε1, ε2, ε3) is the i-th component of the displacement at the point (x1, x2, x3) due to a force F in the j-th direction placed at the point (ε1, ε2, ε3) on the surface. μ represents an elastic constant of the medium. If an observation point k is located on a surface with its GPS coseismic vertical displacement recording Uk, then according to (1) and the numerical method described by Okada (1992)[22], given a suitable seismic source model Δuij (i=1~N, j=1~M) distributed on the fault plane, the corresponding theoretical vertical ground displacement, represented as Uk, can be calculated well. In this study, we utilize the Genetic Algorithm (GA) algorithm to derive an appropriate seismic source model with the object function to minimize the differences between observed data and theoretical values of n locations. That is written as follows min (∑ |Uk n k=0 − Uk|). (2) The GA is an optimization algorithm based on natural genetics and biology concepts. It simulates the process of biological evolution by using a fitness function to measure the fitness of individuals and employing genetic operations (such as crossover, mutation, and selection) to generate new individuals, gradually optimizing the solution. Liao et al. (2008, 2013)[23, 24] applied GA to solve seismic inversion and derived the main characteristics of the source models compared with other studies. Therefore, we still use GA in this research to detect the source model of 0206 earthquakes. 2.2. Coulomb’s Stress Changes The Coulomb’s stress change acting on the target failure plane is denoted as ∆CFS and can be represented by Ganas et al. (2008)[25] .'CFS n += (3) where ∆𝜏 is the shear stress change on the receiver plane, ∆σn is the change in normal stress, and μ is the effective coefficient of friction. The parameters ∆σn and ∆𝜏 can be calculated as ∆𝜎𝑛 = ∑𝑖𝑗𝑛𝑖𝒔𝒊𝒋𝑛𝑗. (4) ∆𝜏 = ∑𝑖𝑗𝑙𝑖𝒔𝒊𝒋𝑛𝑗. (5) where ni and li are the normal and shear slip vectors, respectively. In this research, the shear modulus was assumed as 4*1010 Pa, the Poisson ratio was 0.25, and the μ value was 0.4. The stress tensor sij is represented as sij = 2μeij + λδij∑kekk. (6) where δij are the components of the identity matrix, eij are strain tensors, and μ and λ are Lamé 3819 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate parameters. According to Harris et al. (1998)[26] that the Coulomb’s stress change greater than 0.1 bar may even trigger the occurrence of other earthquakes when the background stresses are intense, meaning if the Coulomb’s stress change is above zero, it suggests that there is additional stress being added to a fault, which can push it towards brittle failure. On the other hand, if the change decreases, it can help prevent a rupture from happening. 2.3 Semi-Empirical Stochastic Method If the fault plane is divided into N*M sub-faults, the seismic wave Acc k(t) to an observation station is represented as 𝐀𝐜𝐜𝐤(𝐭) = ∑ ∑ 𝐀𝐢𝐣 𝐤 (𝐭 − 𝐭𝐢𝐣).𝐌 𝐣=𝟏 𝐍 𝐢=𝟏 (7) where tij is the delayed time, aij(t) is the waveform from ij-th subfault, k=1~3 mean north-south and east-west and vertical directions, respectively. The sub-fault nij with seismic moment M0 to a station at a distance r can create the corresponding acceleration spectrum aij(f) at one direction of a seismic wave (Boore , 1983)[27] r effPffSCMf Qfr c )/( max0ij ),(),()(a − = (8) 34 PRFR C  = (9) where R is the radiation pattern that is affected by rupture parameters, F is the free surface amplification effects, PR is a constant,  and  are the density and the shear wave velocity. ),( c ffS , represented as a source spectrum, which is calculated as a 2− model. S(𝑓, 𝑓𝑐) = ξ 1+( 𝑓 𝑓𝑐 )2 . (10) where fc is the corner frequency and  is a constant to estimate for the free surface effect. The acceleration time domain aij(t) can be obtained by performing a Fourier inverse transform of eq. (8). To reduce the errors from the underestimations of low-frequency portions (f≦1Hz) of aij(t), Joshi et al. (2012)[28] employed the correction function F(t) convolving with the aij(t) to obtain the seismic wave F(t) from the ij-th sub-fault. The F(t) is represented as F(t) = δ(t) + [ N−1 r(1−exp(−1))) ] ‧exp (− t r ). (11) Aij(t)=aij(t)*F(t). (12) Finally, the seismic waveform in a specific station can be computed by summing up the contributions of seismic waves from all sub-faults as described in Equation (7). 2.4. Parameters of Intensity Given the inverse seismic source model and in combination with the seismic wave simulation method described in section 2.3, seismic waves can be computed for different locations. Consequently, important parameters can be extracted from these seismic waves to assess their correlation with building damage. In this study, five parameters were computed: seismic spectral intensities acceleration (SIa) and velocity (SIv), standardized version of cumulative absolute velocity (CAVstd), characteristic intensity (Ic), and maximum incremental velocity (MIV). The SIs can evaluate seismic damages of the buildings with different heights (Liao et al., 2020)[29] and are defined as follows: SIa(ε) = 1 T2−T1 ∫ Sa(ε, T)dT T2 T1 . (13) SIv(ε) = 1 T4−T3 ∫ Sv(ε, T)dT T4 T3 . (14) where Sa(ε, T) and Sv(ε, T) are spectral acceleration and velocity, and T1~T4 in this study are 0.1, 0.6, 0.6, 1.6 corresponding to the buildings with 1~7 and 7~21 floors in Taiwan. The CAVstd means filtering out some low-amplitude, non-damaging ground motions (under 0.025g) to be a potential 3820 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate intensity measure related to damages by [16]. CAVstd can represented as CAVstd = ∑ (H(PGAi − 0.025) ∫ |Acc(t)|dt i i−1 N i=1 ) (15) where N is the number of non-overlapping one-second time intervals, PGAi is the peak ground acceleration (g) in time interval i, and H(x) is the Heaviside step function. The Ic is expressed as Ic = arms 1.5 ts 0.5. (16) arms = √ 1 tr ∫ [a(t)]2dt tr 0 . (17) where a(t) is an acceleration time history, tr means the total duration of the ground motion, and ts represents significant duration time which is the interval between the times at which 5% and 95% of the Arias intensity[17]. The MIV is defined as the maximum area under the acceleration time history between two successive zero-crossings, which is suitable to access the seismic damages due to the long- duration pulse by source directivity effect in near-fault regions [24]. MIV is represented as MIV = max (∫ |a(t)|dt) ti+1 ti . (18) 3. Results and Discussions 3.1. Source Model of 0206 Earthquake This study employs the inverse analysis of the coseismic vertical displacement recorded at GPS stations to infer the rupture process of the 0206 earthquake. The details are presented in Table 1. Using the coordinate difference between the three-day average before the earthquake and the one-and-a-half- day average after the earthquake, the coseismic displacement from continuous GPS stations can be accurately estimated [5]. The comparisons between the GPS records and the simulated results are depicted in black and red arrows, respectively, in Figure 1. Meanwhile, Figure 2 presents the distribution of the inverted seismic source. Table 1. Coseismic offsets recorded at continuous GPS stations. Station Lon.(deg) Lat.(deg) U(mm) CHUN 121.3863 23.4542 6.0000 DNFU DSIN 121.4755 23.6864 23.6325 9.0000 121.3912 99.9000 FLNM 121.5126 23.5998 69.1000 FONB 121.4466 23.7476 -2.4000 HNSN 121.5141 23.5995 8.7000 HUAL 121.3013 24.3390 13.7000 HUAN 121.6067 23.9767 70.4000 HUAP 121.2658 24.1448 4.9000 JSUI 121.7427 24.3103 29.7000 JYAN 121.4171 23.4933 3.5000 KNKO 121.2195 24.2438 -0.1000 NAAO 121.4989 23.4735 10.9000 NDHU 121.8034 24.4506 16.9000 NSAN 121.5440 23.8985 36.4000 PEPU 121.3760 24.4295 2.7000 SCHN 121.6035 24.0192 98.1000 SHUL 121.6448 24.1291 -7.4000 SICH 121.5559 23.7889 10.2000 SLIN 121.6476 24.1270 -7.6000 SOFN 121.4346 23.8132 5.6000 WULI 121.5914 23.8716 14.4000 YENL 121.3017 24.3536 18.4000 3821 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate E415 G970 G972 G979 121.5950 23.9048 -34.4000 121.5662 24.0058 50.9000 121.5972 23.9267 -52.0000 121.5887 23.9265 163.6000 GA37 121.5285 23.9463 7.4000 GA40 121.6013 23.9849 -77.6000 GE46 121.5749 24.0049 40.8000 GE63 121.5964 23.9880 -73.9000 GE64 121.6120 24.0082 120.9000 GA41 121.6021 24.0157 40.5000 E549 121.5853 23.8455 -40.6000 GA43 121.5338 23.8880 162.2000 GE56 121.5129 23.8856 -11.4000 GE57 GE08 121.4916 23.7280 -41.6000 121.4608 23.7557 -27.0000 GE55 121.4380 23.7313 -37.5000 Figure 1. 3822 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate Comparisons of the recorded vertical GPS displacements (black arrow) with the simulated vertical surface displacements (red arrow) based on the inverted source model. The red star is epicenter of the 0206 earthquake. Figure 2. The source model of 0206 earthquake inverted by combining GA method with GPS displacements listed in Table 1. Based on the results shown in Figure 1, it can be observed that there are some larger ground displacement values near the epicenter, but the ground displacements gradually decrease with increasing distance from the epicenter. This illustrates the geometric attenuation phenomenon of seismic energy transmission which is relevant to distance. Moreover, the theoretically computed ground displacements from the inverted seismic source fit better near the epicenter, but the relative errors increase as the distance from the epicenter increases. This phenomenon suggests that the calculation method proposed by Okada (1992)[22] can provide better estimates of ground displacements near the epicenter where the influence of the seismic source is greater, under the assumption of homogeneous geological structures. However, as the distance from the epicenter increases, the effects of geological and site-specific conditions become more significant, making simple structures insufficient to evaluate the surface displacements. According to the calculations made by USGS, the Hualien 0206 earthquake had a nodal plane 1 with a strike, dip, rake of (209.o, 73.o, 22.o), while the other nodal plane had a strike, dip, rake of (112.o, 69.o, 121.o). Moreover, a study by Lee et al. (2018)[30] suggested that the earthquake may have been caused by a west-dipping fault movement, implying that the Hualien fault plane is part of the Central Range fault. This study initially adopts a westward dipping fault plane to simulate a seismic source model. 3823 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate Utilizing the method described in section 2.1., theoretical ground displacements are computed and compared with actual observed GPS recordings, and the seismic source model was iteratively adjusted using the genetic algorithm until the error between theoretical and observed displacements converged. From the rupture process depicted in Figure 2, several characteristics can be discerned. Firstly, the maximum slip displacement is approximately 1.08 m, with an average of about 0.41m. Secondly, larger asperities are distributed above the hypocenter and within a range closer to the surface (<10km) which is agreeable to the results in [4], explaining why the moderate-sized 0206 earthquake can induce significant seismic activity and surface displacement. When there is a significant slip on a fault plane closer to the surface, the energy transmitted through seismic waves experience less dissipation due to the shorter propagation distance. Consequently, it leads to more powerful ground shaking and larger surface displacement. Finally, it is noteworthy that there is a large asperity distributed along a south- southwest direction, consistent with the study of [2], which may serve as a pressure source for the aftershock sequences occurring in the Taiwan Longitudinal Valley and trigger the movement of Meilun Fault. 3.2. Coulomb’s Stress Changes at Different Depths of 0206 Earthquake We have successfully computed the Coulomb stress changes for depths ranging from 6 km to 12 km using the inverted seismic source model and present the distribution of aftershocks (3650 gal) are mainly distributed in the southwest region extending from the epicenter to the southeast, covering an area of approximately 370 km2. This phenomenon may be attributed to the directivity effect of the seismic source. Our speculation is agreeable with the study of [2]. The buildings (white squares in Figures) damaged by the 0206 earthquake are all located within this area where SIa exceeds 650 gal, falling within the critical threshold, and are distributed along the Meilun Fault. It reveals that when assessing the seismic damage of buildings ranging from 1 to 7 stories, it is appropriate to use a critical threshold of 400 gal for the SIa parameter. This is because the risk of damage to buildings increases when the SIa exceeds this threshold, especially for structures of these heights. Having a critical threshold of 400 gal provides a reference for evaluating the extent of damage that earthquakes may cause to these buildings, and can assist in formulating appropriate disaster 3827 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate prevention and response measures. It is worth noting that from the distribution of SIv in Figure 6, values exceeding 30 cm/sec are also observed, mainly located in the southern region of the epicenter. However, only above the Meilun Fault line on the right side, SIv exceeds the critical threshold of 30 cm/sec. Most of the damaged buildings accumulated around the south of Meilun Fault line, including some hotels exceeding 7 stories in height, are situated at locations where SIv is approximately 28 cm/sec. There could be two explanations for this phenomenon. The first explanation is that adopting a critical threshold of 30 cm/sec is based on historical earthquake records and may require an adjustment to account for variations around this threshold. The second possibility is that the destructive force of the 0206 earthquake has pushed these buildings, ranging from 7 to 21 stories in height, to the edge of their seismic resilience. Additionally, the increase in Coulomb stress generated by the 0206 earthquake's seismic source induced the movement along the Meilun Fault. This combination of seismic force and fault movement may have resulted in the damages being particularly concentrated along the fault line of the Meilun Fault. The distributions of Ic and CAVstd are depicted in Figures 7 and 8, respectively. Figure 7. Ic distribution calculated using the inverted source model. The dashed line indicates areas with an Ic value of 316cm1.5/sec2.5; damaged buildings were located within the contour line. The red star represents the hypocenter of the 0206 earthquake, and the white squares represent the buildings damaged by this earthquake. The red line means Meilun Fault line. 3828 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate Figure 8. CAVstd distribution calculated using the inverted source model. The dashed line indicates areas with an CAVstd value of 418cm/sec; damaged buildings were located within the contour line. The red sta represents the hypocenter of the 0206 earthquake, and the white squares are the buildings damaged by this earthquake. The red line means Meilun Fault line. From these two figures, it can be observed that damaged buildings are located at positions where Ic is greater than 316 cm1.5/s2.5 and CAVstd greater than 418 cm/s. According to previous studies of [17- 19], from a physical perspective, the deformation and structural damage caused by earthquakes are expected to be reflected through Ic and CAVstd. As the values of both parameters increase, it indicates that the buildings experience deformation beyond their capacity to withstand due to the destructive force of the earthquake, resulting in severe internal structural damage. This ultimately leads to the collapse or severe damage of the buildings. Figure 9 shows the distribution of MIV caused by 0206 earthquake. 3829 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate Figure 9. MIV distribution calculated using the inverted source model. The dashed line indicates areas with an MIV value of 30 cm/sec; damaged buildings were located within the contour line. The red star represents the hypocenter of the 0206 earthquake, and the white squares are the buildings damaged by this earthquake. The red line means Meilun Fault line. From Figure 9, it can be observed that the destroyed buildings are almost located within or near the range of MIV=30.0 cm/sec. This analysis result is consistent with [20] findings, indicating that setting the boundary value of MIV to 30.0 cm/sec is an effective indicator. Furthermore, according to the research of [33], some nearby strong-motion stations recorded pulse-like velocity ground motion waveforms on both sides of the Meinong Fault line. The velocity pulse waveforms, whether caused by source directivity or site effects, are a primary cause of structural damage to buildings. Additionally, these waveforms are also responsible for the higher MIV near the Meinong Fault. 4. Conclusion In the event of an earthquake, immediate recognition of the distribution patterns of aftershocks and areas where buildings may be damaged is absolutely critical. This enables the swift allocation of disaster relief resources and prevents secondary disasters caused by aftershocks. As such, it is paramount for the effective implementation of earthquake disaster prevention and mitigation engineering. This study conducted research on the Hualien area of Taiwan, focusing on the ML 6.4 earthquakes that occurred on February 6, 2018. The researchers used GPS displacement data and the genetic algorithm (GA) method to inverse the source rupture process and calculate Coulomb stress changes. They then observed the distribution of aftershocks at different depths within one month. Seismic waves by employing the semi- empirical stochastic method were also calculated at different locations to extract five important parameters, which were then correlated with the damaged buildings from the earthquake to gain a better understanding of the relationship between parameter thresholds and their application in 3830 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate earthquake damage assessment. The study found that the source directivity is toward to south-west direction and that larger slip occurred in the shallower portions of the fault plane of the seismic source model, leading to larger surface displacements near the epicenter. After analyzing the changes in Coulomb stress, it was observed that aftershocks happening within a month after the main shock were all within the range of Coulomb stress increase of 0.1 bar. Additionally, as the depth increased, the distribution of aftershock sources shifted from the west-southwest direction of the main shock epicenter to the north-northeast. Although numerous blocks in the western part of the epicentral area showed an increase in Coulomb stress, only a few aftershocks occurred. However, these areas are likely to become seismicity-prone regions in the future, as evidenced by two moderate earthquakes (ML=5.0) occurring there two years later. Finally, the study found that the damaged buildings were located within or in very close proximity to the critical range of various parameters. This clearly indicates that the critical values of different parameters can be used as strong indicators of the distribution areas of damaged buildings. It is noteworthy that, based on the results of SIv, this study suggests that the higher buildings (7-22 floors) that were damaged by the 0206 earthquake may have been a result of a combination of seismic forces and induced movement of the Milun Fault. Acknowledgments: The authors would like to express their gratitude to the anonymous reviewers for their diligent work. Your valuable feedback has contributed significantly to the meaningfulness of this paper, enriching its content and depth. Additionally, we extend our heartfelt thanks to Krirk University in Thailand for providing research funding, which facilitated the completion of this study. Copyright: © 2024 by the authors. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). References [1] K. C. Hao, Z. K. Guan, W. F. Sun, P. Y. Jhong and D. Brown, "Aftershock Sequence of the 2018 Mw 6.4 Hualien Earthquake in Eastern Taiwan from a Dense Seismic Array Data Set," Seismological Research Letters, vol. 90, pp. 60– 67, 2018. DOI: https://doi.org/10.1785/0220180233 [2] P. R. Jian, S. H. Hung and L. Meng, "Rupture behavior and interaction of the 2018 Hualien earthquake sequence and its tectonic implication," Seismological Research Letters, vol. 90, pp. 68–77, 2018. DOI: https://doi.org/10.1785/0220180241 [3] J. H. Nieh, T. H. Hsu, H. C. Cheng, K. C. Chong and P. F. Lai, "2018 Taiwan Hualien earthquake-disaster lessons we learned in the emergency department of a tertiary hospital," J Acute Med., vol. 10, pp. 149-155, 2020. DOI: 10.6705/j.jacme.202012_10(4).0003. [4] M. H. Huang and H. H. Huang, "The complexity of the 2018 Mw 6.4 Hualien earthquake in east Taiwan," Geophysical Research Letters, vol. 45, pp. 13249–13257, 2018. DOI: https://doi.org/10.1 029/2018GL080821. [5] Y. Y. Wen, S. Wen, Y. H. Lee and K. E. Ching, "The kinematic source analysis for 2018 Mw 6.4 Hualien, Taiwan earthquake," Terr. Atmos. Ocean. Sci., vol. 30, pp. 1-11, 2019, DOI: 10.3319/TAO.2018.11.15.03 [6] S. Y. Huang, J. Y. Yen, B. L. Wu, I. C. Yen and R. Y. Chuang, "Investigating the Milun Fault: The coseismic surface rupture zone of the 2018/02/06 ML 6.2 Hualien earthquake, Taiwan," Terr. Atmos. Ocean. Sci., vol. 30, pp. 311-335, 2019. DOI: 10.3319/TAO.2018.12.09.03 [7] R. D. Hwang, C. Y. Lin, W. Y. Chang, T. W. Lin, Y. L. Huang and J. P. Chang, "Multiple-event analysis of the 2018 ML 6.2 Hualien earthquake using source time functions," Terr. Atmos. Ocean. Sci., vol. 30, pp. 367-376, 2019. DOI: 10.3319/ TAO.2018.11.15.01 [8] G. Toyokuni, D. Zhao and K. H. Chen, "Structural control on the 2018 and 2019 Hualien earthquakes in Taiwan," Physics of the Earth and Planetary Interiors, vol. 312, pp. 106673, 2021. DOI https://doi.org/10.1016/j.pepi.2021.106673. [9] A. Hakami, A. Kumar, S. J. Shim and Y. A. Nahleh, "Application of soft systems methodology in solving disaster emergency logistics problems," International Journal of Industrial Science and Engineering, vol. 7, pp. 783–790, 2013. DOI: https://doi.org/10.5281/zenodo.1089433. [10] Y.Y. Chiu, H. Omura, H. E. Chen and S. C. Chen, "Indicators for post-disaster search and rescue efficiency developed using progressive death tolls," Sustainability, vol. 12, pp. 8262, 2020. DOI:https://doi.org/10.3390/su12198262. [11] S.M. Alif, E.I. Fattah, M. Kholil and O. Anggara, "Source of the 2019 Mw6.9 Banten intraslab earthquake modelled with GPS data inversion," Geodesy and Geodynamics, vol. 4, pp.308–314, 2021. DOI: https://doi.org/10.1016/j.geog.2021.06.001. https://creativecommons.org/licenses/by/4.0/ 3831 Edelweiss Applied Science and Technology ISSN: 2576-8484 Vol. 8, No. 6: 3816-3831, 2024 DOI: 10.55214/25768484.v8i6.2826 © 2024 by the authors; licensee Learning Gate [12] H. Duan, J. Chen, S. Zhang, X. Wu and Z. Chu, "Coseismic fault slip inversion of the 2013 Lushan Ms 7.0 earthquake based on the triangular dislocation model, " Scientific Report, vol. 12, pp. 3514, 2022. DOI: 10.1038/s41598–022– 07458–z. [13] A. Peidou, D. F. Argus, F. W. Landerer, D. N. Wiese and M. Ellmer, "GPS displacement dataset for the study of elastic surface mass variations," Earth System Science Data, vol. 16, pp. 1317–1332, 2024. DOI: https://doi.org/10.5194/essd–16–1317–2024. [14] B. Y. Liao, C. H. Liu, T. W. Sheu and Y. T. Yeh, "Earthquake hazard in central Taiwan evaluated using the potentially successive 2013 Nantou Taiwan earthquake sequences," Geomatics, Natural Hazards and Risk, vol. 11, pp. 678–697, 2020. DOI: https://doi.org/10.1080/19475705.2020.1745901. [15] B. Y. Liao, H. C. Huang and S. Xie, "The source characteristics of the Mw 6.4, 2016 Meinong Taiwan earthquake from teleseismic data using the hybrid homomorphic deconvolution method," Applied Sciences, vol. 12, pp. 494, 2022. DOI: https://doi.org/10.3390/app12010494. [16] K. W. Campbell and Y. A. Bozorgnia, "Ground motion prediction equation for the horizontal component of cumulative absolute velocity (CAV) using the PEER-NGA database," Earthquake Spectra, vol. 26, pp. 635–650, 2010. DOI: https://doi.org/10.1193/1.4000012. [17] L. Danciu and G. A. Tselentis, "Engineering ground-motion parameters attenuation relationships for Greece," Bulletin of the Seismological Society of America, vol. 97, pp. 162–183, 2007. DOI: 10.1785/0120040087 [18] A. Massumi and F. Gholami, "The influence of seismic intensity parameters on structural damage of RC buildings using principal components analysis," Applied Mathematical Modelling, vol. 40, pp. 2161–2176, 2016. DOI: https://doi.org/10.1016/j.apm.2015.09.043. [19] C. Zhai, S. Li and L. Xie, "Study on inelastic displacement ratio spectra for near-fault pulse-type ground motions," Earthquake Engineering and Engineering Vibration, vol. 6, pp. 351–355, 2007. DOI:https://doi.org/10.1007/s11803– 007–0755–x. [20] B. Y. Liao, S. Xie and T.S. Hsieh, "Correlations between ground motion parameters measures and structural damages of the Mw6.4, 2016 Meinong Taiwan earthquake using hybrid simulation method," Civil Engineering and Architecture, vol. 11, pp. 1372–1382, 2023. DOI: 10.13189/cea.2023.110321. [21] J. A. Steketee, "On volterra's dislocation in a semi-infinite elastic medium," Canadian Journal of Physics, vol. 36, pp. 192–205, 1958. [22] Y. Okada, "Internal deformation due to shear and tensile faults in a half-space," Bulletin of the Seismological Society of America, vol. 82, pp. 1018–1040, 1992. [23] B. Y. Liao and H. C. Huang, "Rupture process of the 2002 Mw 7.9 Denali earthquake, Alaska, using a newly devised hybrid blind deconvolution method," Bulletin of the Seismological Society of America, vol. 98, pp. 162–179, 2008. [24] B. Y. Liao, Y. T. Yeh, T. W. Sheu, H. C. Huang and L. S. Yang, "A rupture model for the 1999 Chi-Chi earthquake from inversion of teleseismic data using the hybrid homomorphic deconvolution method," Pure and Applied Geophysics, vol. 170, pp. 391–407, 2013. DOI:10.1007/s00024–012–0498–0 [25] A. Ganas, A. Gosar and G. Drakatos, "Static stress changes due to the 1998 and 2004 Krn mountain (Slovenia) earthquakes and implications for future seismicity," Natural Hazard Earth System Sciences, vol. 8, pp. 59–66, 2008. [26] R. A. Harris and R. W. Simpson, "Suppression of large earthquakes by stress shadows: a comparison of Coulomb and rate-and-state failure," Journal of Geophysical Research, vol. 103, pp. 24439–24451, 1998. [27] D. M. Boore, "Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra," Bulletin of the Seismological of Society American, vol. 73, pp. 1865–1894, 1983. [28] A. Joshi, P. Kumari, S. Singh and M. L. Sharma, "Near–field and far–field simulation of accelerograms of Sikkim earthquake of September 18, 2011 using modified semi-empirical approach," Natural Hazards, vol. 64, pp. 1029–1054, 2012. DOI: 10.1007/s11069–012–0281–7 [29] B. Y. Liao, C. H. Liu, T. W. Sheu and Y. T. Yeh, "Earthquake hazard in central Taiwan evaluated using the potentially successive 2013 Nantou Taiwan earthquake sequences, Geomatics," Natural Hazards and Risk, vol. 11, pp. 678–697, 2020. DOI: https://doi.org/10.1080/19475705.2020.1745901. [30] S. J. Lee, T. C. Lin, T. Y. Liu and T. P. Wong, "Fault–to–fault jumping rupture of the 2018 Mw 6.4 Hualien earthquake in eastern Taiwan," Seismological Research Letters, vol. 90, pp. 30–39, 2018. DOI: 10.1785/0220180182.j [31] B. Y. Liao and H. C, Huang, "Coulomb stress changes and seismicity in central Taiwan due to the Nantou blind- thrust earthquakes in 2013," Journal of Asian Earth Sciences, vol. 124, pp. 169–180, 2016. DOI:10.1016/j. jseaes. 2016.05.001 [32] B. Y. Liao, S. Xie and T. S. Hsieh, "Using Coulomb stress changes and seismic spectrum intensities evaluation of the seismic hazard potential of Meishan earthquake in central Taiwan," Civil Engineering and Architecture, vol. 9, pp. 2196–2204, 2021. DOI: 10.13189/cea.2021.090709. [33] Y. Y. Lin, H. Kanamori, Z. Zhan, K. F. Ma and T. Y. Yeh, "Modelling of pulse-like velocity ground motion during the 2018 Mw 6.3 Hualien earthquake, Taiwan," Geophysical Journal International, vol. 223, pp. 348–365, 2020. DOI: https://doi.org/10.1093/gji/ggaa306