Complex rupture process of the 2014 Thailand MW 6.2 earthquake on the conjugate fault system of the Phayao fault zone

Tira Tadapansawut, Graduate School of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan tira.tadapansawut@gmail.com Yuji Yagi, Faculty of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan yagi-y@geol.tsukuba.ac.jp Ryo Okuwaki, Faculty of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan Mountain Science Center, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan COMET, School of Earth and Environment, University of Leeds LS2 9JT, UK rokuwaki@geol.tsukuba.ac.jp Shinji Yamashita, Graduate School of Science and Technology, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan luxemburg1386@gmail.com Kousuke Shimizu, Graduate School of Science and Technology, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan smzksk@geol.tsukuba.ac.jp


Introduction
The seismicity of Thailand is relatively low: less than 10 earthquakes with a magnitude greater than 5 have been registered since the 1970s (ISC, 2021) (Fig. 1). Although situated in a low seismicity zone, Thailand is surrounded by major active faults, such as the Sagaing Fault in Myanmar and the major Aliao Shan-Red River fault north of Thailand (Morley et al., 2011) ( Fig. 1). These faults are subject to a progressive clockwise strain rotation caused by the motions induced by the escape tectonics from the Tibetan Plateau to SE Asia and the Sumatra-Andaman subduction zone (e.g., Huchon et al., 1994;Leloup et al., 2001;Morley, 2002;Morley, 2007;Morley et al., 2011). Thailand has complex geological structures that include multiple active fault zones (e.g., Morley, 2002;Morley et al., 2011) (Fig. 1). According to Charusiri and Pum-Im (2009) and Morley et al. (2011), many active fault zones in Thailand are part of a strike-slip fault system trending northeast-southwest and northwest-southeast.
These trends are a result of the development of the major Cenozoic rift basin that is subject to a north-south compression and east-west extension. Geological records suggest that there is historical seismicity since the Late Quaternary at the northern part of Thailand associated with the active fault zones (e.g., Fenton et al., 2003;Pailoplee and Charusiri, 2016). One of the largest historical earthquakes in 1545 collapsed an immense pagoda in Wat Chedi Luang temple in Chiang Mai province (Kázmér et al., 2011) (Fig. 1). Even though the seismicity of Thailand is relatively low compared to that of the surrounding region, seismic hazard study (e.g., Fenton et al., 2003) shows the likelihood of earthquakes recurring on the order of 50 years with anticipated intensity in the area of active fault zones of II to VII on the Modified Mercalli scale (Pailoplee et al., 2009;Pailoplee and Charusiri, 2016).  (Styron and Pagani, 2020). (b) The wider view of the study region. The rectangle is the map region of (a). The bathymetry/topography is from GEBCO (2020).

Figure 2:
The study area of the 2014 Thailand earthquake. The yellow star shows the epicenter. The red beachball shows the GCMT solution of the 2014 Thailand earthquake (Dziewonski et al., 1981;Ekström et al., 2012). The black beach balls show the focal mechanism of the relocated aftershocks with the moment magnitude larger than 4.2 (Pananont et al., 2017). The dashed lines highlight the aftershock lineations. The solid lines are the active faults (DMR, 2016, Kanthiya et al.,2019 A relocated hypocenter of the 2014 Thailand earthquake was obtained by Pananont et al. (2017) using the double-difference algorithm HypoDD (Waldhauser, 2001): it is located at 19.733°N and 99.689°E with a 5 km hypocentral depth that is between the middle of the MLF and the top of the PF (Fig. 2). The centroid moment tensor solution shows the nodal planes oriented NNW-SSE and ENE-WSW with MW 6.2 (GCMT, 2014). The relocated aftershocks during the first week (Pananont et al., 2017)  the N-S trending aftershocks (NSTA) and the ENE-WSW trending aftershocks (EWTA). The regional moment tensor solutions of the aftershocks (Pananont et al., 2017) are located along the NSTA and EWTA, with their strike directions aligned with the trends of the NSTA and EWTA (Fig. 2).
The regional moment tensor solutions of the mainshock and the aftershocks determined by Noissagool et al. (2016) show that the principal compressive stress orientation is NNE-SSW (N18E) that is consistent with the regional stress orientation in northern Thailand (Heidbach et al., 2016;Tingay et al., 2010;Simons et al., 2007). The high shear stress zone is related to the strike orientation of the active MLF that is close to the EWTA: N30E-N50E. This high shear stress zone contributes to the initiation of slip based on Mohr-Coulomb failure criteria (Noisagool et al., 2016). Pananont et al. (2017) studied the aftershocks sequence occurring within hours by analyzing the changes in the stress field due to the rupture, for which they computed the Coulomb stress changes: they suggested that the mainshock occurred on the right-lateral faulting along the NSTA. They argued that the complex rupture process that has produced the complicated pattern of the aftershock distribution has a more elaborate geological and geomechanical origin. However, the source mechanism of the 2014 Thailand earthquake has not been clearly understood; whether the rupture evolves along the apparent conjugate fault system inferred from the aftershock distribution. The detailed imaging of the source process of the 2014 Thailand earthquake should be a critical basis to illuminate the causative relationship between the rupture evolution and the geometric complexity in the fault system for the smallerscale, M6-class earthquake, which has been difficult to investigate in a means of finite-fault inversion.
A possibility of the complex fault geometry can be expected from a simple observation of the teleseismic waveforms. In general, the waveform shape at the stations at close azimuth and epicentral distance become similar if an earthquake occurs along a single, simple fault plane, which shares the similar radiation pattern among those stations. However, as shown in Figure 3, the waveforms of the TIP and ARU stations at the nearby similar azimuth and epicentral distance show the different waveform shape and amplitude, which is unexpected if the earthquake rupture propagates along a single flat plane with a constant slip vector (Tadapansawut et al., 2021). This may imply that the mainshock mechanism may involve geometric complexity. In addition, the aftershock distribution with two major trends of the NSTA and EWTA (Fig. 2) may suggest the complexity of fault geometry of the mainshock.
To resolve the possible complex fault geometry, we apply a new framework of the flexible finite-fault inversion algorithm for teleseismic body waveforms . We introduce a relative weight smoothness constraint that is proportional to the components of each basis moment tensor (Kikuchi and Kanamori, 1991) into the potency density tensor inversion of Shimizu et al. (2020) which can mitigate the effect of the modelling errors originating from the uncertainty of Green's function (Yagi and Fukahata, 2011) as well as the uncertainty of fault geometry (e.g., Ragon et al., 2018) (see details in the Method section). This method can simultaneously estimate the distribution of the focal mechanism and the slip along the assumed model plane; it enables the reconstruction of complex rupture processes, including those occurring along faults containing fault bends and those consisting of multiple subevents, without a priori assumption of the fault geometry (Tadapansawut et al., 2021;Yamashsita et al., 2021). The improved flexible finite-fault inversion framework has been applied to large earthquakes such as the 2020 MW 7.7 Caribbean earthquake (Tadapansawut et al., 2021) and the 2018 MW 7.9 Gulf of Alaska earthquake , but it has never been applied to smaller-scale M6-class earthquakes like the 2014 Thailand earthquake. In the following sections, we first perform a numerical test to evaluate the resolvability of the rupture process for a moderate-sized event. Then we apply the algorithm to the teleseismic body waves of the 2014 Thailand earthquake. We estimate the spatiotemporal distribution of both the slip and the fault geometry and discuss the source process of the 2014 Thailand Earthquake, which is heavily controlled by the geometric complexity of the fault system and the associated local stress field.

Method: Relative weight potency-density inversion
To construct a rupture model of the 2014 Thailand earthquake, we apply the flexible teleseismic finite-fault inversion of Shimizu et al. (2020). The method can resolve the fault geometry and the slip along the flexible assumed model plane without a priori fault geometry assumption; it represents the shear-slip vectors with five basis double-couple moment tensor components (Kikuchi and Kanamori, 1991). The method is based on the novel finite-fault inversion of Yagi and Fukahata (2011), which can mitigate the modelling error originating from the uncertainty of Green's function. The observation equation is defined as where ! is the teleseismic waveform at station . "! is Green's function for the th component of the basis double-couple moment tensor at station , and "! is the error of Green's function.
̇" is the potency-rate density function for the th component of the basis double-couple moment tensor at the source location of the assumed model plane ( ). #! is the background and instrument noises.
Although method of Shimizu et al. (2020Shimizu et al. ( , 2021 can resolve both fault geometry and slip, the source focal mechanism change within subevent is still difficult to reveal because of the spatiotemporal smoothing constraint which was introduced by the Gaussian with a same covariance into the potency-rate density function without distinguish for all five basis doublecouple components. This may introduce bias because the potency-rate density of the dominant basis component becomes smoother than those of the minor components. To mitigate the bias due to the smoothing constraints, we applied a new framework of the relative weight smoothness constraint developed by Yamashita et al. (2021): it adds an inverse relative weight parameter (1/ " ) to the standard deviation of each basis component that is proportional to the double-couple component of the GCMT moment tensor solution (Fig.   S3). To avoid the instability of the solution due to the extremely small relative standard deviation, we set the minimum weight smoothness constraint to 5% of the maximum relative standard deviation. The sensitivity analysis of the minimum percentage of the relative weight smoothness constraint showed that when it is set to 5%, 10% and 20%, the solutions do not visibly differ (Figs. S1 and S2). This new framework has been proven efficient for the analyses of the source process of the 2018 Gulf of Alaska earthquake  and the 2020 Caribbean earthquake (Tadapansawut et al., 2021). There it solves the problem of oversmoothing the slip distribution of the major components, and allows more clear capture of the rupture propagation.
The GCMT solution (Dziewonski et al., 1981;Ekström et al., 2012) for the 2014 Thailand earthquake, shows dominant strike-slip faulting associated with the pure strike-slip M1 and M2 moment tensor components (Kikuchi and Kanamori, 1991) (Fig. S3). Therefore, the relative weight smoothness constraint on the strike-slip components should become dominant over other components like the vertical slip components M3, M4 and M5 (Fig. S3).
A comparison of the results between the non-weighted smoothing constraint  and the relative weighted smoothing constraint  methods is presented in Fig. S4. With the non-weighted smoothing constraint method, the moment-rate function is a smooth function with two peaks and the reverse fault component is dominant; in contrast, with the relative weighted smoothing constraint method, the moment-rate function has several distinct peaks and the strike-slip component is dominant in all regions.

Numerical test
In this study, the new finite-fault inversion method is applied to a small-scale M6-class earthquake for the first time. To evaluate the resolvability of this approach, we perform a numerical test using synthetic waveforms. The synthetic rupture process is assumed to take place along the dipping planes of two faults: Fault 1 and Fault 2 (Fig. 4). We generate the synthetic waveforms at 25 teleseismic stations (Fig. 3), which are the same stations that we use for the analyses of the 2014 Thailand earthquake (see the Data and Model setting section).
Each fault model consists of a set of 1×1 km 2 sub-fault grids situated along the strike and the dip. The total fault length for Fault 1 is 9 km and for Fault 2 is 19 km; the width of both faults is 7 km. The input source model is composed of a spatially variable slip distribution (Fig. 5a).
The input total seismic moment is 0.58 × 10 19 Nm (MW 6.4). We assume that the rupture initiates at the hypocenter on Fault 1. Then the rupture-front propagates bilaterally along Fault 1 at a constant speed of 3.6 km/s (Fig. 5e). The rupture on Fault 2 initiates at its eastern edge at 1.4 s hypocentral time and propagates unilaterally towards the southwest (Fig. 5e) at a constant speed of 3.6 km/s. We generate the synthetic waveform by using Green's function incorporating random Gaussian noise and then adding the background noise (Fig. S10). Then, we invert the synthetic waveforms using the new finite-fault inversion method adopted in this study. The model plane is a non-rectangular horizontal model plane covering the input sources.
The sub-fault has a dimension of 2 × 2 km 2 . The moment-rate function for each sub-fault is represented as a linear B-spline function with a duration of 7.2 s. The total rupture duration is set at 9.0 s. The maximum rupture velocity assumption is set at 3.6 km/s. The output total seismic moment is 0.57 × 10 19 Nm (MW 6.4). The moment-rate function (Fig. 5d) is consistent with that of the input (Fig. 5b) and shows minor moment release at ~3 s and main moment release at 4-6 s. The resolved spatiotemporal rupture evolution shows that the rupture originates at the hypocenter and propagates southwards along Fault 1 during the first 2 s. Then a rupture occurs on the east end of Fault 2 and propagates in the southwest direction (Fig. 5f).
The strike orientation of the spatiotemporal distribution of the potency-rate density tensor at each time window is consistent with the strike orientations of the input faults. Overall, the spatiotemporal distribution of the potency-rate density agrees with that of the assigned input model. The results of this numerical test confirm that our new framework of finite-fault inversion can resolve complex ruptures even for cases of smaller-scale events consisting of multiple rupture segments and with geometric changes in the fault system.

Data and model setting
For the analysis of the 2014 Thailand earthquake rupture, we use 25 vertical components of the globally observed teleseismic waveforms (Fig. 3b) obtained from the Global Seismographic Network and the Federation of Digital Seismograph Network provided by IRIS-DMC. The data are selected based on signal-to-noise ratio high-enough to distinguish the Pwave arrival and to ensure azimuth coverage (Fig. 3) Although the earthquake magnitude is small, the waveform at the first 10 s can be distinguished from the noise level (Fig. 6).
We manually pick the first arrival of P-wave, and convert it into a velocity waveform to remove the instrument response. Then we resample it at 0.2 s. Following Kikuchi and Kanamori (1991), we calculate Green's function at 0.1 s sampling rate for the components of each basis moment tensor. The finer sampling of Green's function with respect to the observed waveform sampling ensures sufficient resolution for the time shift relating to location of each subfault to the hypocenter. After this, we resample Green's function at 0.2 s, which is the waveform sampling rate. The simplified 1-D near-source structural velocity model from Wongwai et al. (2013) is applied to calculate the Haskell propagation matrix for Green's function (Table S1). The sensitivity of the structural velocity model is tested against the global model of the Earth's crust CRUST1.0 (Laske et al., 2013; Table S2) and the regional structural model developed by Wongwai et al. (2013) (Figs. S6 and S7). The results show that the slip pattern is robustly resolved for the different velocity models. However, the amplitudes of the moment-rate functions by applying both of these structural velocity models are slightly different, probably due to the different velocity and thickness of the first layer (Tables S1 and   S2). We do not apply a low-pass filter to the waveforms, as we want to preserve all complexities of the rupture process contained in the waveform data.
The assumed model plane is confined by the relocated aftershock distribution and covers the NSTA and EWTA that are the expected rupture fault planes. The method we use allows the assumption of a horizontal model plane that is independent of the actual fault plane(s); however, such a supposition can produce a very smooth solution that will impair the interpretation of the rupturing path or the fault geometry. This problem is distinct in the conjugate strike-slip fault earthquakes with multiple fault planes because if the model space is wide and covers unnecessary space where the slip is unlikely to occur, then the unnecessary slip is squeezed out from the actual slip due to the smoothing effects. To mitigate this issue, Yamashita et al. (2021) restricted the horizontal model plane only to the aftershock region, and obtained a non-rectangular plane; in that way, the rupture propagation is captured in detail and the solution is more stable . We assume the model plane to have a strike of 60° and a dip of 0°. The model plane is a non-rectangular horizontal model plane with a maximum total length of 25 km along the EWTA and 12 km along the NSTA (Fig. 7). The sub-fault has a dimension of 2 × 2 km 2 and lies along the strike and the dip. The moment-rate function for each sub-fault is represented as a linear B-spline function with a duration of 7.2 s.
The total rupture duration is set at 9.0 s. The maximum rupture velocity is set at 3.6 km/s and is approximated from the preliminary rupture duration around 7 s at a distance of 25 km from the assumed model fault plane. The approximated value is also equal to the first layer shear wave velocity (V s ) of the simplified structural velocity model (Wongwai et al., 2013). We tested the sensitivity of the model to maximum rupture velocity by setting it to 3.2, 3.6, and 4.0 km/s (Figs. S8 and S9). The analysis shows that if we assume the rupture-front speed to be 3.2 km/s, then the early rupture behaviors at <3 s are masked and not well resolved (Fig. S9); in contrast, the model becomes consistent even if we assume the rupture-front speed to be 3.6 or 4.0 km/s (Fig. S9). For the maximum rupture-front speed, we choose to use the value of 3.6 km/s. As the initial rupture point, we use the relocated hypocenter with coordinates 19.733°N, 99.689°E at 5 km depth (Pananont et al. 2017).

Results
The total moment tensor solution, calculated by the integration of all potency-rate density tensors, exhibits strike-slip faulting with the two nodal planes striking at 249° (ENE-WSW) and 339° (NNW-SSE) (Fig. 7a). The moment-rate function shows at least two rupture episodes. One is between 0 and 1.5 s with a low moment-rate and the other is between 1.5 and 4.5 s with a high moment-rate. The highest moment-rate occurs at around 3.5 s (Fig. 7b). The total seismic moment is 0.36 × 10 19 Nm (MW 6.3), which is slightly larger than the GCMT solution (MW 6.2) and the USGS W-phase moment tensor (MW 6.1). The larger seismic moment in our work is probably due to our model covering a wider area that includes the aftershock distribution along the NSTA and EWTA.
The static distribution of the potency density reveals two large potency zones located in the middle of the EWTA and one in the middle of the NSTA. The larger potency density in the EWTA is around 1.8 m and the potency density in the NSTA is around 1.3 m (Fig. 7c). The nodal plane distribution of the potency density shows that the strike orientation rotates clockwise along the EWTA from 240° to 265° (Fig. 7c) and along the NSTA from 7° to 24° (Fig. 7c).
The spatiotemporal distribution of the potency-rate density exhibits two rupture episodes, one along the NSTA and the other along the EWTA (Fig. 8a). The initial rupture of the mainshock originates at the hypocenter in the first 1.5 s and propagates south along the NSTA at a rupture speed of ~3.0 km/s. The second rupture occurs at the eastern edge of the EWTA between 1.0 s and 1.5 s and propagates southwest along the EWTA at a rupture speed of ~3.5 km/s (Fig. 8a). The second rupture has the highest potency-rate in the middle of the EWTA between 2.0 and 3.0 s and terminates at the west end of the EWTA at 4.5 s. These two rupture episodes coincide with the dominant peaks seen in the moment-rate function (Fig. 7).
The spatiotemporal distribution of the moment tensor solution shows two dominant patterns of strike-slip faulting, both with smaller-scale fluctuation of the fault geometry in each lineation (Fig. 8b). One with a strike at NNE-SSW near the epicenter occurred between 0.5 and 1.5 s and the other with a strike at ENE-WSW northwest from the epicenter occurred between 1.5 and 4.0 s. The nodal plane distribution extracted from the resultant spatiotemporal potencyrate density tensor (Fig. 8b) exhibits clockwise strike rotation from ~18° to ~33°. For the rupture propagating from the epicenter towards the south along the NSTA during the first 1.5 s, this rotation coincides with the timing when the large potency-rate density is observed. The clockwise rotation of the strike also occurs from the middle to the west end of the EWTA from ~218° to ~250°; it is associated with the second rupture arising at the eastern edge of the EWTA, propagating west during the period between the 2.0 and 3.5 s and having the highest potency-rate of around 1.2 m/s (Fig. 8b). These spatiotemporal changes of fault geometry are robust for different assumptions of the regional structural velocity model and of the maximum rupture velocities (Figs. S6 to S8).

Discussion
Our finite-fault source model of the 2014 Thailand earthquake distinguished two rupture episodes that show a dominant strike-slip faulting consisting of different rupture lineations along the NSTA and EWTA (Fig. 8). The nodal plane distribution along the NSTA shows nodal strikes in the NNE-SSW direction and the auxiliary plane in the ESE-WNW direction (Fig. 7). The nodal plane distribution along the EWTA shows nodal strikes in the ENE-WSW direction and the auxiliary plane in the NNW-SSE direction (Fig. 7). The consistency between the nodal plane distribution (Fig. 7c)  spatiotemporal potency-rate density distribution (Fig. 8) facilitates identification of the possible fault geometry. The striking plane along the NSTA is determined to be in the NNE-SSW direction and is associated with the rupture propagating towards the south. The striking plane along the EWTA is determined to be in the ENE-WSW direction and is associated with the rupture propagating towards the southwest. The obtained two dominant fault planes along the NSTA and EWTA are consistent with the two distinct trends of the relocated aftershock distribution (Pananont et al., 2017). The first is the N-S trend (~180° from north) along the NSTA located near the epicenter and the second is the ENE-WSW trend (~60° from north) along the EWTA located northwest from the epicenter. Although the geometry of our model, designed to cover the aftershock distribution area, is non-rectangular, the potency density and the potency-rate density of each sub-fault are estimated independently from the assumed model geometry.
The strike orientation of the potency density tensor shows geometric bends; since they are changes in the strike direction of the rupture-hosting fault, they play an important role in the earthquake rupture process (e.g., Kase & Day, 2006;Ulrich et al., 2019). We observe distinct fault geometries along the NSTA and the EWTA: one with clockwise strike orientation rotation from the north to the south along the NSTA (7° to 24° of the reference points in Fig.   7) and the other toward the west along the EWTA (240° to 265° of the reference points in Fig.   7). These lineations of the strike directions in the NSTA and EWTA coincide with the spatial pattern of the aftershock distribution. In addition to these general lineations of the fault geometry, the strike orientation of the spatiotemporal potency-rate density tensor distribution ( Fig. 8b) also shows dynamic changes of the fault geometry. During the first 1.5 s, the strike orientation rotates clockwise as it propagates from the northern to the southern edges of the NSTA (Fig. 8b). Then, at around 1.5 s as the rupture migrates from the NSTA to the EWTA, the fault strike direction changes from NNE-SSW at the northern edge of the NSTA to ENE-WSW at the eastern edge of the EWTA, which implies a conjugate fault with fault planes lying in the NSTA and EWTA. Next, between 2.0 and 3.5 s the second rupture propagates along EWTA from its eastern edge towards the south-west and terminates at around 4.5 s at its western edge. In this process, the time at which the strike orientation rotates clockwise corresponds to the time of the largest potency-rate density. According to the surface fault lines (DMR, 2016;Noisagool et al., 2016), the orientation of the known active conjugated strikeslip faults of the PF and MLF shows striking at N5E-N13E and N30E-N50E; this is consistent with our findings that at the northern edge of the NSTA, the striking is in the NNE-SSW direction and at the eastern edge of the EWTA, in the ENE-WSW direction. The multiple subevents at the conjugated strike-slip fault system are possibly due to the complex rupture evolution (e.g., Yamashita et al., 2021;Meng et al., 2012). Therefore, we conclude that the rupture evolution of the 2014 Thailand earthquake is characterized by multiple sub-events in the conjugated strike-slip fault system of the PF and MLF. Furthermore, studies using the flexible teleseismic finite-fault inversion have shown that complexities in the faulting system, like geometric bends, can cause the non-smooth rupture propagation of the mainshock (e.g., Okuwaki et al., 2020;Tadapansawut et al., 2021;Yamashita et al., 2021). The dynamic rupture simulation demonstrates that rupture perturbation could have occurred from the bend along a strike-slip faulting (Duan & Oglesby, 2005;Kase & Day, 2006). Thus, we suggest that the rupture along the NSTA and EWTA exhibits the complexity of the fault geometry that includes a bend. This complicated fault system is the reason for the fluctuation of the rupture front. It also can act as a barrier for the termination of the rupture propagation at the southern edge of the NSTA and the western edge of the EWTA.
The rupture evolution of the 2014 Thailand earthquake displays two distinct rupture episodes with rupture directions along the NSTA and EWTA. These rupture episodes reveal two perpendicular planes that coincide with the aftershock distribution pattern along the NSTA and EWTA (Fig. 8). The aftershock distribution shows a spatial gap of around 5 km located between the northern edge of the NSTA and the eastern edge of the EWTA (Figs. 7 and 8). In this gap, our source model shows the lowest potency-rate density for the entire rupture duration (Fig. 8). The agreement between the low potency-rate density area and the spatial gap in the aftershock distribution suggests that the two conjugate faults are not connected. The discontinuity of the co-seismic slip across the gap in the conjugated fault system suggests that the second rupture episode initiated at the eastern edge of the EWTA may have been triggered by the initial rupture episode along the NSTA.
The spatial distribution of the P-axis (or the maximum compressive stress axis) azimuth, extracted from the potency density tensor for each sub-fault ( Fig. 9) exhibits clockwise rotation from the northern edge to the southern edge of the NSTA, from 52° to 69° azimuths, and from the eastern edge to the western edge of the EWTA, from 15° to 39° azimuth.
This clockwise rotation of the P-axis azimuth is in accordance with the aftershock lineation along the NSTA and EWTA (Fig. 9a). The histogram of the P-axis azimuth distribution displays two peaks, one at 10°-30° and the other at 50°-70° (Fig. 9b). Since most of the aftershocks of magnitude 4 and above, for which a focal mechanism solution is estimated by Noisagool et al. (2016), occur in the EWTA (e.g. Pananont et al., 2017), the principal compressive stress axis orientation of 18° obtained from the mainshock and the aftershocks of the 2014 Thailand earthquake (Noisagool et al., 2016) is most likely reflecting the one in the EWTA domain, which coincides with our estimates of the P-axis azimuth distribution along the EWTA (Fig. 9). Whilst, the direction of the P-axis azimuth along the NSTA obtained in this study (~50°, Fig. 9) is rotated clockwise by about 30° from the principal compressive stress axis orientation obtained by Noisagool et al., (2016). As a result, the strike-slip direction with strike of 24° at the southern part of the NSTA obtained from this study (Fig. 7), is opposite to that expected from the principal compressive stress axis orientation of 18° obtained by Noisagool et al., (2016). On the other hand, if Coulomb's friction factor is a typical value of 0.6, the two peaks of our P-axis histogram (Fig. 9) can be naturally explained as a shift of the P-axis of the conjugate fault plane (e.g. Iio, 1997), which leads to ~35° principal stress axis; about 15° clockwise from the principal stress axis orientation obtained in Noisagool et al. (2016). We should mention, however, the focal mechanism solutions obtained in this study are affected by dynamic changes in the stress field due to seismic waves, and estimation of the principal stress axis is beyond the scope of this study. Our results suggest that further investigation of the stress field in this region is needed, taking into account the spatial bias of aftershock distribution, which affects the estimates of the principal stresses for the conjugate fault earthquake.

Conclusion
We construct a source model for the 2014 Thailand MW 6.2 earthquake that occurred within the Phayao fault zone in northern Thailand, by applying a new framework of the flexible teleseismic P waveform finite-fault inversion and resolved both the fault geometry and the slip.
Our source model exhibits complex rupture evolution consisting of two rupture episodes along a conjugated strike-slip fault system that comprises two distinct fault planes. These planes coincide with the relocated aftershock distribution. The initial rupture originates at the hypocenter and propagates southward along the north-south oriented fault plane near the epicenter. Then the second rupture episode is triggered north of the epicenter at the eastern edge of the conjugated east-west oriented fault plane and propagates southwestward until the rupture terminates. Our source model shows not only the conjugate fault geometry but also the fault bends that are related to the smaller-scale features of the aftershock lineation in each rupture episode. Our model also suggests that the conjugate fault system of the 2014 Thailand earthquake is not connected at the junction; the observed spatial gap (~5 km) may account for the triggering of the second rupture episode. The spatial variation of the principal stress axis inferred from our finite-fault model suggests an in-situ stress state of the Phayao fault zone, which is responsible for the complex rupture evolution of the 2014 Thailand earthquake.

Availability of data and materials
Teleseismic waveforms were obtained from the following networks: the Global Seismograph