Introduction

Northern Thailand and adjacent regions are characterized by widespread intraplate seismicity (Fig. 1) associated with strike-slip and extensional faulting developed within a complex continental interior setting (Morley and Wang 2023). Although earthquake magnitudes are generally moderate, seismic activity is distributed across several major fault systems, making this region one of the most tectonically active areas in Thailand. Prominent structures include the Mae Chan Fault Zone (MCFZ), Mae Lao Fault Zone (MLFZ), and Phayao Fault Zone (PYFZ), which together accommodate ongoing deformation within the northern Indochina block (Fig. 2b). Both historical and instrumental records indicate that these faults remain active and capable of generating damaging earthquakes (Pananont et al. 2017; Tadapansawut et al. 2022; Sethanant et al. 2025).

Fig. 1
Fig. 1
Full size image

Major active faults in Thailand and adjacent regions (modified from DMR (2023) and Wang et al. (2014)) and epicentral distribution from the years 2010 to 2024 (TMD, 2024). Red beach balls are earthquakes of No. 1: the 2011 Mw6.7 of Tarlay earthquake (Shi et al. 2018), No. 2: the 2019 of Mw6.4 Saiyabuli earthquake (U.S. Geological Survey, Earthquake Hazards Program (U.S. 2019)), No. 3: the 2014 Mw6.3 of Mae Lao earthquake (Pananont et al. 2017), No. 4 and No. 5 are the 2022 Mw6.4 and the 2023 Mw6.4 of Keng Tung earthquake (U.S. Geological Survey, Earthquake Hazards Program (U.S. 2023)), respectively. CM Chiang Mai, LP Lampang, PY Phayao, CR Chiang Rai

Fig. 2
Fig. 2
Full size image

Tectonic setting, seismicity, and geological cross-section in northern Thailand. a Tectonic subdivision of Thailand showing major crustal blocks, including the South China Block, Indochina Block, Nan Back-arc Basin, Sukhothai Zone, Inthanon Zone, and Sibumasu Block. b Geology and active fault zones in the northern region encompassing Chiang Rai (CR), Phayao (PY), and Wang Nuea. Red lines indicate active faults from DMR (2023), while black lines represent the Latest Miocene—Recent faults within the Wiang Pa Pao Basin (Mankhemthong et al. 2020b). Earthquake epicenters (2010–2024, TMD data) are marked by circles scaled and colored by magnitude, including significant events such as the 2014 Mw 6.3 Mae Lao and 2019 Mw 4.9 Wang Nuea earthquakes. Background topography highlights structural and fault-related lineaments. c Simplified geological cross-section showing stratigraphic units and fault geometries across the region, including bedding orientations and fault types. The cross-section emphasizes the relationship between surface geology, fault activity, and subsurface structure (Mankhemthong et al. 2020b). MTFZ Mae Tha Fault Zone, PYFZ Phayao Fault Zone, MCFZ Mae Chan Fault Zone, MLFZ Mae Lao Fault Zone

Over the past decade, seismic activity along the MLFZ and PYFZ has attracted increasing attention following the occurrence of several moderate to large earthquakes. Among these, the 2014 Mw 6.3 Mae Lao earthquake represents the largest instrumentally recorded event in northern Thailand. This earthquake nucleated at a depth of approximately 7 km, produced some of the strongest ground shaking documented in the country, and caused widespread damage to more than 15,000 buildings, accompanied by surface fissuring (Noisagool et al. 2016; Titaya 2016). The spatial distribution of aftershocks and isoseismal patterns showed a clear NE–SW alignment, consistent with the regional strike of the fault system extending from northern Thailand into Laos (Chang et al. 2022). As the epicenter was located near the northern termination of the PYFZ, the Mae Lao earthquake provided important constraints on fault interactions and regional stress transfer, and its source characteristics have since been incorporated into national seismic hazard models (Tadapansawut et al. 2022).

On 20 February 2019, the Mw 4.9 Wang Nuea earthquake occurred within the central segment of the PYFZ, approximately 40 km south of the Mae Lao rupture. According to the Thai Meteorological Department, the mainshock was shallow, with a reported focal depth of about 5 km, and was followed by more than 70 aftershocks (TMD, 2019). Despite its proximity to the Mae Lao earthquake, the Wang Nuea event displayed distinct characteristics, including the absence of clear surface rupture and a diverse set of aftershock focal mechanisms. These differences suggest that fault segments within the same regional stress field may respond in markedly different ways, depending on fault orientation, inherited structural fabrics, and depth-dependent mechanical properties of the crust.

A particularly striking feature of the Wang Nuea earthquake sequence is the coexistence of reverse, strike-slip, and normal faulting within a relatively small spatial domain. Such variability in faulting style raises important questions regarding the vertical and lateral heterogeneity of the stress field, the role of inherited structures, and the mechanical stratification of the crust beneath northern Thailand. Previous studies in the region have suggested that intraplate seismicity may exhibit depth-dependent changes in faulting behavior, reflecting transitions from transpressional deformation in the upper crust to more extensional regimes at greater depth (Noisagool et al. 2016; Pananont et al. 2017). However, the extent to which individual earthquake sequences record systematic depth-dependent seismotectonic behavior remains insufficiently constrained.

Independent geophysical constraints on crustal structure in this region come from magnetotelluric (MT) surveys. Boonchaisuk et al. (2017) used 3-D MT imaging to resolve the resistivity structure of the Phayao Fault Zone, identifying a conductive zone interpreted as saline fluid in the source region of the 2014 Mw 6.3 Chiang Rai (Mae Lao) earthquake. More recently, Amatyakul et al. (2025) extended MT-based resistivity imaging across northern Thailand (1994–2024), linking fluid-related conductive anomalies to the spatial distribution of intraplate seismicity in the Wang Nuea region. These studies provide an independent, non-seismological line of evidence for the presence of fluids or mechanically weak material within the PYFZ, a factor relevant to the depth-stratified faulting behavior discussed later in this study.

In this study, we investigate the 2019 Mw 4.9 Wang Nuea earthquake and its associated aftershock sequence to better understand fault geometry, rupture processes, and depth-dependent seismotectonics along the PYFZ. By integrating full-waveform inversion, aftershock focal mechanism analysis, and morphotectonic observations, we aim to clarify how inherited fault systems and crustal heterogeneity control intraplate earthquake behavior in northern Thailand. The results provide new insights into vertically stratified deformation within an obliquely convergent continental interior and have broader implications for understanding intraplate seismic hazards in Southeast Asia.

Regional tectonic and seismological framework of northern Indochina

Far-field effects of the India–Eurasia collision

The present-day tectonic regime of northern Thailand is primarily influenced by the ongoing convergence between the Indian and Eurasian plates. Since the Cenozoic, this collision has resulted in large-scale crustal thickening in Tibet and lateral extrusion of continental blocks toward Southeast Asia, accommodated by major strike-slip fault systems and distributed intraplate deformation (Tapponnier et al. 1982; Morley 2007). Finite element modeling and stress inversion studies indicate that compressional stresses associated with the India–Eurasia collision are transmitted hundreds to thousands of kilometers into the continental interior (Müller et al. 2015), producing a dominant NE–SW to NNE–SSW compressional stress field across much of Indochina (Meyers et al. 2018).

Deformation in mainland Southeast Asia is therefore not localized along a single plate boundary but distributed across a broad and mechanically heterogeneous region. Within this setting, inherited crustal structures exert a first-order control on fault kinematics, allowing strike-slip, reverse, and normal faulting to coexist within relatively small areas. Variations in fault orientation, depth, and local rheological conditions strongly influence the style of deformation and earthquake rupture behavior (Morley and Wang 2023). This far-field stress regime provides the regional boundary conditions for fault reactivation along major intraplate fault systems in northern Thailand, including the Phayao Fault Zone examined in this study.

Tectonic setting of northern Thailand

The present-day deformation of northern Thailand is strongly controlled by inherited crustal structures associated with the Indosinian Orogeny. The PYFZ is located within the north–south-trending Chiang Mai–Chiang Rai Suture Zone (CMCRSZ), which marks the tectonic boundary between the eastern margin of the Sibumasu Terrane, referred to as the Inthanon Zone, and the Sukhothai Zone to the east (Metcalfe 2017; Sone and Metcalfe 2008; Fig. 2a). This suture records the closure of the Palaeo-Tethys Ocean during the mid-Triassic Indosinian Orogeny, as evidenced by the widespread deposition of post-orogenic non-marine sediments across both terranes (Hada 1997; Sone and Metcalfe 2008). The CMCRSZ includes a well-developed accretionary complex composed of oceanic crustal materials, seamount carbonates, volcanic rocks, and mélanges, which were emplaced during subduction and collision processes between the Sibumasu and Sukhothai terranes (Hara et al. 2017; Mankhemthong et al. 2020a). Present-day fault orientations indicate that the suture closed in a west–east direction, resulting in north–south structural fabrics in northern Thailand (Sone and Metcalfe 2008). However, structural studies within the CMCRSZ show that thrust and shear zones have predominantly north–south trends but include significant variations, possibly reflecting a curvature in the Sibumasu margin during collision (Hara et al. 2017; Mankhemthong et al. 2020a).

The Wang Nuea area in Lampang Province contains a wide variety of rock units distributed across distinct geographic zones (Fig. 3). Silurian–Devonian metasedimentary rocks are extensively exposed in the westernmost part, particularly along the highlands surrounding the Wiang Pa Pao area (Wongwanich and Boucot 2011). These represent the oldest rock units in the area and form part of the regional basement. To the east of this zone, forming a central-western belt, Carboniferous sandstones and Carboniferous–Permian sedimentary sequences are present, extending southwards across Highway No. 120, which connects Wiang Pa Pao to the southern boundary of the study area. In the southeastern part, Permian carbonate rocks dominate, forming a continuous north–south-trending unit adjacent to Triassic sedimentary rocks, where the well-developed Wang Nuea segment of the PYFZ forms a contact. In the Wang Nuea and Wiang Pa Pao areas, both Triassic sedimentary rocks and Triassic granite intrusions are reported, with the igneous plutonic bodies often surrounded by their associated contact metamorphic zones (Baum and Hahn 1977; DMR 2006). The western highland zones, particularly west of Wiang Pa Pao, are underlain by Cretaceous granite, forming elevated terrain that contrasts with the surrounding sedimentary sequences. Tertiary sediments and Quaternary alluvial deposits are concentrated in low-lying areas, especially near major streams and rivers as well as around urban centers (DMR 2006).

Fig. 3
Fig. 3
Full size image

Geologic and structural maps of the study area. a Simplified geologic map showing the distribution of sedimentary and igneous rocks, delineated by stratigraphic age. b Structural map showing fault and bedding orientations across the region. Red lines represent active fault zones identified by DMR (2023), while black lines show faults of the latest Miocene to recent age as mapped by Mankhemthong et al. (2020b). Bedding information is derived from DMR datasets (black from Baum and Hahn 1977, blue from DMR 2006). MJ Mae Jaem, WP Wiang Pa Pao, WN Wang Nuea

Geomorphologically, the PYFZ generally trends north–south and bends toward the south-southeast, exhibiting predominantly dextral strike-slip motion along a contractional bend. Morphotectonic landforms along the Wang Nuea segment, such as fault scarps, triangular facets, shutter ridges, offset streams, and linear valleys, further support the interpretation of dextral strike-slip deformation along the PYFZ. The distribution of bedding orientations, fault traces, and associated geomorphic features reflects the active tectonics shaping this structurally complex and geologically varied region. However, the surface expression of the PYFZ within the Tertiary cover in the Wang Nuea area is poorly exposed, which limits structural interpretation of hidden fault segments.

In the Wiang Pa Pao area, the Wiang Pa Pao Basin (WPB) marks a significant lithological and structural boundary between Triassic granitoids of the Sukhothai Arc Terrane and Palaeozoic metasedimentary rocks of the Sibumasu Terrane (Mankhemthong et al. 2020a). The basin developed in association with Cenozoic multi-phase rifting. During the Late Oligocene–Miocene, regional east–west extension produced major rift basins along pre-existing crustal anisotropies, including NW–SE, NE–SW, and WSW–ENE structural trends. In the latest Miocene to Recent, dextral strike-slip deformation reactivated these structures, forming a north–south-trending asymmetric half-graben in the WPB (Mankhemthong et al. 2020b).

Subsurface interpretation using gravity data, combined with surface geological mapping, indicate that the WPB is underlain by an inherited thrust fault related to the Indosinian Orogeny, which juxtaposes the Sibumasu and Sukhothai terranes (Fig. 2c). This concealed thrust is inferred from the sharp gravity gradient and the abrupt basement lithology contrast across the basin, indicating tectonic emplacement of metasedimentary rocks against granitic bodies along a crustal-scale compressional structure (Mankhemthong et al. 2020a). The orientation and reactivation of this pre-Cenozoic thrust system have exerted a fundamental control on subsequent fault segmentation and basin evolution. According to Mankhemthong et al. (2020a), whose study proposed a structural model of the Mae On Depression and the Mae Tha Fault Zone (MTFZ), and Mankhemthong et al. (2020b), the intersection of the PYFZ with these inherited structures is particularly relevant for interpreting the complexity of faulting observed during the Wang Nuea earthquake sequence.

Present-day stress and seismicity in northern Thailand are governed by a regional regime of NE–SW compression (Fig. 1), resulting from the ongoing collision between the Indian and Eurasian plates (Morley 2007). This stress orientation influences the reactivation of pre-existing faults with varying geometries, leading to diverse fault kinematics across the region. Active faults exhibit a range of orientations, from NE–SW striking left-lateral faults such as the MLFZ to N–S or NNW–SSE trending normal and strike-slip faults in WPB and Chiang Mai Basin. Earthquakes in the region typically occur at shallow depths, often less than 10 km, though some events have been recorded at depths exceeding 15 km, suggesting activity not only within the upper crust but possibly at the brittle–ductile transition (Morley 2015). The variation in fault orientation and depth highlights the complex interaction between regional stress and local structural inheritance, making northern Thailand a key area for understanding intraplate seismic hazard.

Regional seismicity and depth-dependent faulting

Instrumental seismicity across northern Thailand, eastern Myanmar, Laos, and southwestern China is dominated by shallow earthquakes, with most hypocenters occurring above 15 km depth, although events extending to depths of approximately 25 km have been reported (Noisagool et al. 2014; Sethanant et al. 2025). Focal mechanism solutions reveal substantial variability in faulting style, including strike-slip faulting along major fault zones, reverse faulting in transpressional settings, and normal faulting at greater depth or within extensional basins.

Several studies have documented systematic depth-dependent changes in faulting style within earthquake sequences in northern Thailand. For example, the 2014 Mw 6.3 Mae Lao earthquake involved predominantly strike-slip rupture at intermediate depths, whereas deeper aftershocks exhibited normal faulting consistent with a rotated stress field and reduced differential stress (Noisagool et al. 2016; Pananont et al. 2017). Regional seismic tomography further identifies low shear-wave velocity zones within the mid-crust, interpreted as rheologically weak layers that may limit brittle failure and promote stress decoupling between the upper and lower crust (Wu et al. 2022).

These observations indicate that intraplate earthquakes in northern Indochina commonly occur within a vertically stratified crust. In this framework, the upper crust accommodates transpressional and strike-slip deformation, while deeper levels transition toward extensional or transtensional regimes. This depth-dependent seismotectonic behavior provides a critical context for interpreting the mixed focal mechanisms and rupture characteristics of the 2019 Wang Nuea earthquake sequence analyzed in this study.

Data and methods

This study integrates full-waveform seismic inversion with morphotectonic analysis to characterize the source parameters and faulting behavior of the Mw 4.9 Wang Nuea earthquake, which occurred in northern Thailand in 2019. Seismic data were obtained from broadband stations operated by the Thai Meteorological Department (TMD) and span the period from 20 February to May 2019, covering the mainshock and its significant aftershocks.

Waveform data were first processed using the Seismic Analysis Code (SAC; Goldstein et al. 2003; Goldstein and Snoke 2005) to identify and verify P- and S-wave arrival times. These arrivals served as input for subsequent hypocentral parameter refinement and waveform modeling. All seismograms were deconvolved to remove instrument responses, detrended, demeaned, and filtered using a four-pole Butterworth bandpass filter of 0.18–0.9 Hz for both Pnl and surface waves. This relatively high-frequency band was chosen because of the small-to-moderate magnitudes of the events (Mw 3.15–4.91), their shallow depths, and the regional station distances (Supplementary Table S1).

To determine the earthquake’s moment magnitude, focal depth, and focal mechanism, we employed the Cut-and-Paste (CAP) inversion method (Zhu and Helmberger 1996), which partitions the waveform into five segments: two Pnl segments and three surface wave segments. These segments are inverted independently to optimize sensitivity to depth-dependent and structural variations within the crust. Synthetic waveforms were computed using the frequency–wavenumber (F-K) method (Zhu and Rivera 2002), incorporating a hybrid velocity model composed of CRUST 1.0 (Laske et al. 2013) and a modified version of the global CRUST 1-D model (Dziewonski and Anderson 1981). The one-dimensional velocity model adopted for waveform inversion in this study is based on the CRUST 1.0 global model (Laske et al. 2013; black lines in Fig. 4). To evaluate the sensitivity of focal mechanism solutions to velocity model choice, CAP inversions were performed using four models: (1) CRUST 1.0 (Laske et al. 2013), (2) the local minimum 1-D model of Saetang and Duerrast (2023), (3) the receiver-function-derived model of Pananont et al. (2017), and (4) the average Vp model derived from stations CRAI and CMMT of Buphu et al. (2025), combined with an assumed Vp/Vs ratio of 1.75. As shown in Fig. 4, the Buphu et al. models exhibit notably higher Vp (~ 6.5–7.0 km/s) in the 3–8 km depth range compared to CRUST 1.0 (~ 6.1 km/s), a difference most consequential for the shallow events (Mainshock, LP1, LP2, LP5, and LP6). The Moho depth of the four models is about 37 km for CRUST 1.0, more than 40 km for Saetang and Duerrast (2023), about 33–35 km for Pananont et al. (2017), and about 32–36 km for Buphu et al. (2025) (Fig. 4). The model of Saetang and Duerrast (2023) uses 5-km layers and hypocenters that are almost all shallower than 20 km, so its Moho is not directly resolved. Independent receiver-function estimates near the study area give a crustal thickness of 33.3 km at station LAMP and 32.3 km at station CMMT (Noisagool et al. 2014), and joint inversion of surface-wave dispersion and receiver functions gives 36.5 km at CMMT (Tadapansawut 2012). The CRUST 1.0 Moho therefore lies at the upper end of the regional range of about 32–37 km. All focal depths obtained in this study are well above the Moho in all four models. Comparative focal mechanism solutions across all four models are presented in Table 1, using variance reduction as the quantitative basis for final model selection. All four models give the same faulting style for every event, so the depth-dependent pattern of reverse, strike-slip and normal faulting is robust to the choice of velocity model. Fault-plane orientation is less well constrained, particularly for the smallest events (Table 1). The CRUST 1.0 model provides a stable reference framework for generating F-K Green's functions and was adopted as the primary model for the final solutions reported here, while acknowledging that local models (Saetang and Duerrast 2023; Buphu et al. 2025) provide better-constrained upper-crustal velocity structure for this region.

Fig. 4
Fig. 4
Full size image

The 1D velocity models of the crust used for Green’s function calculations. Shear-wave velocity (Vs) and compressional-wave velocity (Vp) profiles are plotted against depth for various 1D models. The global reference model CRUST 1.0 in the study area (black) is compared with locally derived models for the Chiang Mai–Chiang Rai region: Saetang and Duerrast (2023) (green), Pananont et al. (2017) (blue), and Buphu et al. (2025) (orange for CRAI station, red for CMMT station). Vp and Vs are represented by solid and dashed lines, respectively. Depth ranges of the Conrad and Moho discontinuities are from Saetang and Srisawat (2025). The Moho depth of each model is about 37 km (CRUST 1.0), more than 40 km (Saetang and Duerrast 2023), about 33–35 km (Pananont et al. 2017) and about 32–36 km (Buphu et al. 2025), and is marked on the figure

A grid search inversion approach was applied to minimize the waveform misfit and identify the best-fitting double-couple solution. The azimuthal distribution of the stations used for the mainshock, including station codes, geographic coordinates, epicentral distances, azimuths and takeoff angles, is provided in Supplementary Table S1 and illustrated in Supplementary Fig. S8. The azimuthal gap is explicitly reported and discussed for each event in Supplementary Information. The CAP grid-search inversion approach was specifically adopted because it is robust under sparse or asymmetric station distributions (Herrmann et al. 2011). Additionally, a sensitivity analysis was conducted to assess the impact of small variations in hypocentral coordinates on the inversion output. This step confirmed that the derived focal mechanism and depth estimates were stable for minor shifts in the assumed epicenter location. The final inversion result is illustrated in Fig. 5, showing the best-fitting waveform solution (example from the mainshock) across selected stations along with the depth-misfit curve, which identifies the optimal focal depth. Equivalent plots for the mainshock and for each of the six Mw ≥ 3.0 aftershocks are given in Supplementary Figs. S1–S7.

Fig. 5
Fig. 5
Full size image

Best solution and representative waveform fit for the mainshock. a CAP inversion on local data (black lines are observed, red lines are synthetic) and b Depth misfit curve for CAP inversion

Aftershock relocation was performed using the double-difference earthquake location algorithm HypoDD (Waldhauser and Ellsworth 2000). A total of 107 aftershocks from the Thai Meteorological Department (TMD) catalog were selected for relocation. The relocation used 6 stations (TM.CMPR, TM.SUAB, TM.LAMP, TM.CRMM, TM.NAN, TM.NANS), a subset of the stations listed in Supplementary Table S1, at distances of 45–135 km from the epicenter with a maximum azimuthal gap of 107°. Phase arrival times of P- and S-waves were manually picked from broadband waveforms and used as cross-correlation differential travel times between closely spaced event pairs. The relocation procedure utilized the CRUST 1.0 velocity model (Laske et al. 2013) as the reference 1-D velocity structure, consistent with the model used for waveform inversion. The hypoDD relocation used only the CRUST 1.0 model, and the effect of the choice of velocity model on the relocated depths was not tested, which is a limitation of this study. Events were linked using a conjugate gradient (LSQR) algorithm, with a minimum of four common stations per event pair required for inclusion in the relocation network. Sixty-nine of the 107 relocated events share epicentral coordinates with at least one other event, in 17 groups; for 20 of these events no cross-correlation data were available, and for the remainder, the weighting given to catalog relative to cross-correlation differential times in the final iterations, together with the damping parameter used, was insufficient to separate events that shared identical input catalog coordinates. This is noted as a limitation of the relocation, and the fine-scale relative positions of these events are interpreted with caution and are not used to infer fault orientation. The relocated hypocenters constrain the spatial distribution of seismicity associated with the 2019 Wang Nuea earthquake sequence, and the results are presented alongside the CAP focal mechanisms in Fig. 6. The relocated position and the residuals of every event are listed in Supplementary Table S2.

Fig. 6
Fig. 6
Full size image

Relocated seismicity of the 2019 Mw 4.9 Wang Nuea earthquake sequence. a Map view showing epicenters of relocated aftershocks (circles scaled and colored by local magnitude) relative to the Wang Nuea segment of the Phayao Fault Zone (PYFZ; red lines from DMR, 2023). The mainshock epicenter and its focal mechanism are shown, and the Mw ≥ 3.0 aftershocks (LP1–LP6) are marked by ringed circles. Beach balls for these events are plotted in b and c, labelled by moment magnitude. b E-W depth cross-section showing the vertical distribution of events; LP3 and LP4 are clearly separated at ~22 km from the shallow cluster at 1–3 km and the mainshock at ~7 km. c N-S depth cross-section confirming the same three-group depth separation. d Magnitude-time plot of the relocated sequence; LP3 (Mw 4.56) on 13 March 2019 (21 days after the mainshock) falls within the active aftershock period. Relocations were performed using the HypoDD double-difference algorithm (Waldhauser and Ellsworth 2000) with the CRUST 1.0 reference velocity model

Principal stress orientations for each depth group (Table 2) were determined by stress inversion of the focal mechanism solutions using the method of Michael (1984), which minimizes the misfit between observed slip directions and the predicted shear stress resolved on each fault plane. We acknowledge that stress inversions based on only two focal mechanisms per depth group carry inherently limited statistical reliability; the results should therefore be treated as indicative rather than definitive, and alternative stress configurations within the uncertainty bounds are discussed in the interpretation section.

Table 1 Focal mechanism solutions (strike, dip, rake, depth) derived from CAP waveform inversion using four 1-D velocity models: CRUST 1.0, Saetang and Duerrast (2023), Pananont et al. (2017), and Buphu et al. (2025). Solutions are shown for the mainshock (MS) and six Mw ≥ 3.0 aftershocks (LP1–LP6).
Full size table
Table 2 Orientation of principal stresses with depth
Full size table

A morphotectonic study was conducted using a high-resolution ALOS PALSAR digital elevation model (12.5 m) and satellite imagery from Google Earth to delineate lineaments and fault traces. Ground deformation was performed through field investigations conducted by the Department of Mineral Resources (DMR) between February and May 2020. These field surveys aimed to characterize surface expressions of active deformation, lithologic boundaries, and fracture patterns associated with the seismogenic structures. Structural data, including joint orientations and brittle deformation features, are used to characterize the geometry of pre-existing fractures and active faults in the area.

Seismic waveform inversion and locations of the aftershocks

According to our full-waveform inversion results, the rupture process of the 2019 Mw 4.9 Wang Nuea earthquake can be interpreted as a shallow strike-slip event on a steeply dipping fault. The mainshock nucleated at a depth of approximately 7 km, which is slightly deeper than the previously reported 5 km by TMD (2019), with a strike of 83°, dip of 73°, and rake of +7° (Aki–Richards convention), indicating predominantly strike-slip motion. Our CAP solution is consistent with the USGS W-phase moment tensor within reported uncertainties; discrepancies with the TMD solution are attributed to differences in velocity model and station coverage. The mainshock mechanism has two nodal planes, an E–W plane (strike 83°, dip 73°, rake +7°) and a nearly N–S auxiliary plane (approximately 351°/83°/163°), and the mechanism alone cannot distinguish between them. Both are plausible structures within a dextral transpressional principal displacement zone. The near N–S plane, with dextral slip, would represent a strand sub-parallel to the PYFZ itself, whereas the E–W plane, with sinistral slip, would represent an antithetic cross-structure linking adjacent strands. The aftershocks are concentrated 5–10 km west of the mapped Wang Nuea segment, indicating that the sequence occurred on an unmapped structure within the PYFZ system and not on the mapped N–S trace. We therefore do not identify a preferred fault plane.

The comparison shows that focal mechanism solutions are broadly consistent across all four models. For the mainshock, strike varies from 55°–115°, dip from 69°–77°, and rake from 5°–15°, with all solutions indicating predominantly strike-slip motion. For LP3 and LP4, all four models consistently yield negative rake values (–20° to –60° and –40° to –60°, respectively) and depths of 15–22 km, confirming the normal-faulting character of these events regardless of velocity model choice. LP1 and LP2 consistently show reverse to oblique-slip mechanisms (rake +25° to +80°) at shallow depths (1–3 km) across all models. These results confirm that the main conclusions regarding depth-dependent faulting style are robust to velocity model choice (Table 1).

The aftershock distribution, based on relocated earthquakes from Buphu et al. (2025) and this study, confirms distinct spatial and depth-related variations in seismicity, with aftershocks concentrated west of the Wang Nuea segment of the PYFZ. HypoDD relocation of 107 aftershocks (ML 1.2–4.2) reveals a compact cluster concentrated approximately 5–10 km west of the Wang Nuea segment of the PYFZ (Fig. 6), spanning latitudes 19.23–19.30°N and longitudes 99.51–99.71°E. The relocated seismicity is predominantly shallow, with most events occurring at depths of 1–6 km; a secondary cluster occurs near 10 km depth. The spatial pattern defines a roughly NNW–SSE trending zone, broadly parallel to the regional strike of the PYFZ, and supports the inference that seismicity was nucleated on a network of sub-parallel fault segments rather than on a single continuous rupture plane. Aftershocks with Mw ≥ 3.0 can be grouped into three depth groups based on their location. The first depth group of aftershocks (LP1 and LP2), occurring at shallow depths (~ 1 km), are characterized by reverse to oblique-slip faulting, with steeply dipping planes and rake values of approximately +80° and +45°, respectively. In contrast, LP3 and LP4, recorded 21 and 22 days after the mainshock respectively, exhibit normal faulting mechanisms with negative rakes (–48° and –60°) at ~22 km depth (a depth determination that is consistent across all four tested velocity models), suggesting a different stress regime at greater depth and a likely triggered origin, though background seismicity cannot be entirely excluded. The final aftershock depth group (LP5 and LP6), located at 2–3 km depth above the mainshock, exhibit strike-slip mechanisms with steeply dipping planes and rake values of –6° and +21°, respectively. This distinct pattern of aftershock distribution and depth-dependent faulting styles suggests that the Wang Nuea earthquake sequence was governed by a vertically and laterally segmented fault system, likely influenced by interactions with a pre-existing fault network.

Source characteristics and their implications for seismotectonics

All moderate-sized events in the 2019 Mw 4.9 Wang Nuea earthquake sequence (Mw ≥ 3.0) occurred within the continental crust, above the Moho discontinuity, which lies at about 32–37 km depth in the Wang Nuea area (Noisagool et al. 2014; Tadapansawut 2012; Wongwai 2011). These events were relocated and exhibited well-constrained focal mechanisms and focal depths (Table 1). The aftershocks are distributed 5–10 km west of the Wang Nuea segment of the PYFZ, where active fault-related geomorphology is not clearly expressed in satellite imagery or DEM data (Fig. 7). This distribution confirms that seismic activity was clustered around the mainshock but with variations in depth and fault orientation.

Fig. 7
Fig. 7
Full size image

West-to-east cross-section and vertical distribution of earthquakes along the line marked in Fig. 3b, crossing the Wiang Ka Long, Wang Nuea and Pa Faek fault segments. Rock units follow the key at left, red lines are faults, and the topographic profile is drawn with a vertical exaggeration of two. Grey circles are the relocated earthquakes of this study projected onto the section, and the north-looking beach balls are colored by focal depth. The histogram at right gives the focal depth distribution. The black histogram is this study and illustrates the seismogenic thickness, while the dark grey and light grey histograms are the 2014 Mw 6.3 Mae Lao earthquake from Sethanant et al. (2025) and Pananont et al. (2017), respectively

A comparison with the 2014 Mw 6.3 Mae Lao earthquake sequence (Pananont et al. 2017; Noisagool et al. 2016) reveals a similar pattern of depth-dependent faulting and stress rotation in northern Thailand. While strike-slip motion is prevalent at intermediate depths, such as the Wang Nuea mainshock at ~7 km, the post-mainshock events LP3 and LP4 at ~22 km depth exhibit normal faulting with moderately dipping planes and a near-vertical σ₁ axis. These events are interpreted as potentially triggered by the Wang Nuea mainshock, based on their temporal proximity (21–22 days post-mainshock) and spatial co-location within the sequence footprint, while acknowledging that certainty of association is limited by the moderate mainshock magnitude. We therefore describe these events as deep post-mainshock events potentially triggered by the Wang Nuea earthquake, while noting the possibility that they represent background mid-crustal seismicity. Statistically, both earthquake sequences show a drop in earthquake occurrence at ~8 km and ~16 km depths (Fig. 7), coinciding with low shear-wave velocity (Vs) anomalies interpreted as weak, potentially ductile zones that may suppress brittle rupture (Wu et al. 2022). These findings highlight the complexity of crustal stress fields and support a model of depth-stratified seismotectonics, in which the upper crust accommodates transpressional or strike-slip faulting, while the lower crust transitions into an extensional regime. Regional tectonic factors, including inherited thrust faults and oblique basement structures, likely influence this variation in faulting styles. These interpretations further refine the conceptual frameworks proposed by Morley and Wang (2023) and Mankhemthong et al. (2020b), which describe multi-phase intraplate deformation shaped by both strike-slip and extensional processes.

Field investigation

The dominant north–south lineaments of the Wang Nuea and Wiang Pa Pao segments are concentrated along the margins of the Wang Nuea area and mainly follow rock contacts, while the NE–SW-trending lineaments are primarily located within zones of exposed basement rocks to the west (Fig. 3). However, the epicenters of the mainshock and aftershocks of the 2019 Mw 4.9 Wang Nuea earthquake are predominantly located away from the major north–south fault segments, suggesting that deformation may also be accommodated by distributed fault networks within the interior of the Wang Nuea area.

The exposure of Triassic sandstone along the Wang Nuea segment, at Wang Thong Waterfall in Wang Kaew Subdistrict, Wang Nuea District, Lampang Province (19° 16′ 16.5″ N, 99° 39′ 50.4″ E), was investigated (Fig. 8). This locality indicates variations in basement fabric that deviate from the dominant regional north–south structural trend. The natural outcrops of sandstone exhibit evidence of NNE- to NE-trending folding (Fig. 8a), which may have formed in association with compressional stress during the mid-Triassic Indosinian Orogeny. Such structural variations, including the strike and dip of foliations in Carboniferous and Permian rocks, have also been documented in the geological survey report of Wang Nuea and Wiang Pa Pao by the Department of Mineral Resources (1989). Therefore, it can be inferred that older rocks farther west may also have experienced deformation under a similar compressional regime and may host pre-existing fault networks. However, according to the morphological interpretation, the boundaries of exposed rocks still correspond predominantly to the active strike-slip fault segments in this region. The rocks at the waterfall site are cut by multiple unfilled fractures that intersect at various angles, suggesting a more recent phase of brittle deformation in this area (Fig. 8b). These fractures include steeply dipping planes with measured orientations of 195/55°, 140/64°, 150/60°, and 259/62°.

Fig. 8
Fig. 8
Full size image

Triassic sandstone outcrops at Wang Thong Waterfall (19° 16′ 16.5″ N 99° 39′ 50.4″ E), illustrating structural deformation features. a Outcrop showing folding within the sandstone. b Close-up of fractures and joints cutting across the sandstone, with a stereonet indicating the orientation of the measured fractures (see location in Fig. 3a)

Further south along the Wang Nuea segment, two Quaternary quarries were investigated for evidence of recent deformation. The quarry near Wat Ban Mae Ma (19°11′ 05.4″ N, 99° 40′ 00.9″ E) consists of well-stratified, semi-consolidated sediments of alluvial origin. The northwest-facing cliff illustrates faulted beds. Fault orientations are difficult to determine in this eroded, strike-parallel section. Fault dip angles vary (Fig. 9a). The second alluvial quarry is located at a roadcut in Ban Mae Yen (19° 07′ 53.3″ N, 99° 40′ 40.2″ E). Overall, the sediments consist of reddish-brown sandy to gravelly deposits with large clasts, including boulders. A discontinuity is observed where the pebbly sand layers are subvertically truncated by poorly sorted boulder conglomerates (Fig. 9b). There is no evidence of erosion or drag near the contact, which may have formed by a strike-slip fault.

Fig. 9
Fig. 9
Full size image

Semi-consolidated sediment outcrops showing deformation related to tectonic activity in the area. a Sediment exposure at a quarry near Wat Ban Mae Ma (19° 11′ 05.4″ N 99° 40′ 00.9″ E), displaying well-defined sedimentary bedding with a clear fining-upward sequence. The red-dashed lines highlight faults, indicating reverse faulting associated with tectonic activities in the Wang Nuea area. b Road-cut outcrop at Ban Mae Yen (19° 07′ 53.3″ N 99° 40′ 40.2″ E), where reddish-brown fine sand and coarse gravel (including boulder-sized clasts) were observed. The red-dashed line marks an abrupt change in sediment layers, with a visible discontinuity between the two sides of the outcrop and gravel oriented vertically along the fault. This feature represents a strike-slip fault that cuts through and offsets the sedimentary layers (see location in Fig. 3a)

Together, these structural and geomorphic observations highlight the overprinted deformation in the Wang Nuea area. The presence of both ductile and brittle features in basement rocks, as well as evidence of recent deformation in Quaternary sediments, indicates a prolonged, multi-phase tectonic history that continues to influence present-day seismicity. These findings emphasize that seismic hazards may arise from subtle or unmapped fault networks.

Depth-dependent seismotectonics in northern Thailand

Following the conceptual framework outlined above, the spatial and depth distribution of the Wang Nuea earthquake sequence provides clear evidence for vertically stratified deformation within the continental crust of northern Thailand. The observed variation in faulting style with depth reflects differences in stress orientation, mechanical strength, and the availability of pre-existing structures capable of reactivation at distinct crustal levels.

The 2019 Mw 4.9 Wang Nuea earthquake and its associated Mw ≥ 3.0 aftershocks (LP1–LP6) reveal a systematic spatial and depth-dependent variation in faulting style and stress orientation within the upper to mid-crust of northern Thailand (Fig. 10). Shallow aftershocks (LP1 and LP2), occurring at approximately 1 km depth in Carboniferous rocks, display predominantly reverse faulting. This corresponds to a compressional to transpressional stress regime along the western margin, characterized by a near-vertical σ₃ and a horizontally oriented σ₁ trending roughly northwest (Table 2). The reverse-faulting style of LP1 and LP2 may reflect localized reactivation of an inherited thrust system previously identified in the Wiang Pa Pao Basin region (Mankhemthong et al. 2020a). Although originally formed during the Indosinian Orogeny, this crustal-scale compressional structure appears to remain mechanically significant under the present-day transpressional stress regime. The orientation and location of LP1 and LP2 suggest that these events may have exploited pre-existing weaknesses along this structure, supporting the interpretation that reactivation of inherited compressional fabrics continues to influence fault segmentation and shallow crustal deformation in the Wang Nuea area.

Fig. 10
Fig. 10
Full size image

Conceptual lithospheric cross-section showing depth-dependent seismotectonics of intraplate earthquakes in northern Thailand (modified from Morley 2015), extending across the western to eastern parts of northern Thailand, including the Wang Nuea area. The section illustrates crustal deformation, major faults, and sedimentary basins. The profile extends from Doi Inthanon in the west, through the Chiang Mai Basin, Mae Tha Fault, and Wiang Pa Pao Basin, to the Phayao Fault Zone (PYFZ) in the east. The Wiang Pa Pao Basin and the earthquake locations are projected onto the section, as their actual positions lie north of the line of section. Subsurface structures include a low shear-wave velocity zone from Wu et al. (2022) and the Moho, with depth constraints derived from receiver function studies (Noisagool et al. 2014; Saetang and Srisawat 2025). The lower detachment surface and the estimated present-day brittle–ductile transition zone are also illustrated. NF Normal Fault, SS Strike-slip Fault, RF Reverse Fault (see location in Fig. 1)

Beneath the Wang Nuea area, deeper aftershocks (LP5, LP6) at 2–3 km depth are characterized by strike-slip faulting, possibly within Permian rocks. Their orientations are consistent with regional NNE–SSW compression and reflect a stress field influenced by the long-term far-field effects of the India–Eurasia collision (Meyers et al. 2018). The mainshock, at approximately 7 km depth, shows near-pure strike-slip motion on a steeply dipping plane, indicating reactivation of a steep fault segment within the PYFZ. The observed spatial variability in faulting styles suggests that pre-existing faults with varied orientations are responding to a transitional stress field. This highlights the complex nature of fault reactivation in an intraplate strike-slip system, where multiple fault segments are activated rather than rupture occurring along a single continuous fault plane.

In addition, LP3 and LP4, occurring at approximately 22 km depth, exhibit normal faulting mechanisms characterized by moderately to steeply dipping fault planes. Stress inversion results suggest a dominantly extensional regime, with a subvertical σ₁ and a horizontally oriented σ₃ trending east–west (Table 2), though these results should be treated as indicative given the small sample size (n = 2 per depth group). The transition from shallow strike-slip to deeper apparent normal faulting is consistent with a progressive decrease in differential stress with depth, combined with increasing lithostatic pressure and thermal weakening of crustal materials. At mid-crustal levels, such conditions may promote gravitational collapse and flexural downwarping, particularly in structurally inherited zones of weakness (Morley and Wang 2023). The occurrence of brittle failure at ~22 km depth suggests that, despite thermal weakening, the crust may remain within or near the brittle regime at this depth. We note, however, that LP3 and LP4 may alternatively represent background mid-crustal seismicity coincidentally recorded during the aftershock monitoring window; the available data do not allow this possibility to be conclusively excluded.

Comparable depth-dependent stress rotation has also been observed in the 2014 Mw 6.3 Mae Lao earthquake sequence, where strike-slip faulting predominated at depths of 12–14 km, whereas events deeper than 22 km exhibited normal faulting with variably oriented stress axes (Noisagool et al. 2016; Pananont et al. 2017). These findings reinforce a depth-stratified model of crustal deformation in northern Thailand, wherein the upper crust accommodates NNE–SSW horizontal compression, while the mid- to lower crust transitions into an extensional or transcurrent regime with rotated principal stress axes. Earthquake hypocenters in both the Mae Lao and Wang Nuea areas tend to decrease in number within low shear-wave velocity (Vs) anomaly zones interpreted from regional tomography (Wu et al. 2022). Independent constraints from magnetotelluric (MT) imaging further support the depth-stratified model. Boonchaisuk et al. (2017) and Amatyakul et al. (2025) imaged the resistivity structure beneath the PYFZ and Wang Nuea region, respectively. A conductive zone at ~8–10 km depth identified by Amatyakul et al. (2025) spatially coincides with the low-Vs anomaly from Wu et al. (2022) at the same depth, jointly suggesting the presence of fluids or mechanically weak material that may suppress brittle failure and contribute to the observed seismicity gap between the shallow cluster (LP1, LP2, LP5, LP6 at 1–3 km) and the mainshock (~ 7 km). At greater depths, the MT resistivity structure is consistent with a transition to a more ductile or fluid-weakened mid-crust, providing a physical mechanism for the inferred extensional faulting behavior of LP3 and LP4 at ~22 km. Integrating these geophysical constraints, we suggest that the depth-stratified faulting pattern reflects both the mechanical stratification of the crust and the distribution of pore fluids, rather than stress rotation alone. This spatial pattern reflects the interplay of multiple tectonic processes, including the long-term far-field effects of the India–Eurasia collision (Meyers et al. 2018), lower crustal flow facilitated by elevated temperatures and rheologically weak layers (Morley 2018; Wu et al. 2022), and inherited anisotropic lithospheric fabrics (Mankhemthong et al. 2020a; Morley et al. 2013).

Conclusions

This study provides new insights into the source characteristics, faulting styles, and seismotectonic framework of the 2019 Mw 4.9 Wang Nuea earthquake in northern Thailand. Full-waveform inversion and aftershock analysis reveal that the mainshock occurred as a shallow event on a steep, unmapped fault segment of the Phayao Fault Zone (PYFZ), with subsequent aftershocks distributed both vertically and laterally across multiple crustal levels. The contrasting focal mechanisms of aftershock depth groups which range from reverse faulting at shallow levels to normal faulting at depths of approximately 22 km, highlight a depth-stratified deformation regime shaped by inherited structural fabrics and variations in crustal rheology.

The coexistence of brittle and ductile deformation features, observed both in the field and in focal mechanism solutions, underscores the role of pre-existing thrust faults, lithological contrasts, and low shear-wave velocity zones in controlling fault segmentation and rupture behavior. The observed seismicity patterns and stress rotations are consistent with a broader regional model of transpressional deformation in the upper crust, which transitions into an extensional regime at depth. This model is further supported by comparative analysis with the 2014 Mw 6.3 Mae Lao earthquake sequence and reinforces the concept of a vertically heterogeneous crustal stress field in northern Thailand.