Rupture Directivity of the 2021 $ML$ 2.2 Gwangyang, Korea, Microearthquake: Toward Resolving High‐Resolution Rupture Process of a Small Earthquake

On 27 August 2021, $ML$ 2.2 earthquake occurred ∼7 km north–northeast of Gwangyang City on the southern tip of the Korean Peninsula (Fig. 1). Three small earthquakes with magnitudes of 1.1–1.7 followed within 26 hr, and two small earthquakes preceded the shock (Table 1). The largest shock of $ML$ 2.2 was felt by residents around the epicentral area with the maximum modified Mercalli intensity III, according to Korea Meteorological Administration (KMA). The sequence ceased on 28 August, and no additional shocks were detected. Although the sequence consists of small shocks, it poses a significant impact on earthquake hazards around the industrial town of Gwangyang, situated within 10 km from the epicentral area (Fig. 1). In the recent years, microearthquake swarms with M≤3 have been detected in southern Korea, for example, 2020 Haenam sequence (Son et al., 2021) and 2019–2020 Suncheon swarms (Kwak et al., 2022) (Fig. 1a). The frequency of reported microearthquake swarms has increased since 2016, perhaps due to the greatly expanded deployment of permanent seismographic stations that are over 300 by the end of 2021.

The epicentral area is a stable continental region with low seismicity, and the bedrock is Pre‐Cambrian gneiss with no known major faults or geological structure. In much of southern Korea, a compressional stress regime developed in the late Miocene (∼10 Ma), mainly due to plate‐driving forces acting on the downgoing slab along the Japan trench and inhibiting forces balancing them (Kim et al., 2022). Currently, this horizontal compression with an east–west trend (∼70°) prevails in southern Korea, including the epicentral area. Since the late 1990s, the southeastern Korean Peninsula has been seismically active with the 2016 $Mw$ 5.5 Gyeongju earthquake; however, during the twentieth century, significant earthquakes occurred in the southwestern region. On 3 July 1936, a moderate sized but damaging earthquake of magnitude 5.3 occurred ∼27 km north, followed by a magnitude 5.2 event ∼7 km east of the 2021 Gwangyang events (Fig. 1a). Hence, the 2021 Gwangyang microearthquakes may provide an opportunity to look at the long‐term seismicity rate in the region and to understand the seismogenesis of the past damaging earthquakes.

Our preliminary analysis of the seismic records from the 2021 Gwangyang sequence suggested that the microearthquakes were very well recorded at 55 stations within 100 km, and that the station azimuthal coverage and density are excellent (Fig. 1b). We observed that three‐component records showed varying high‐frequency contents depending upon the station azimuth (Fig. S1, available in the supplemental material to this article). This may suggest the effect of rupture propagation from the $ML$ 2.2 mainshock. Hence, we attempted to analyze the details of the 2021 Gwangyang microearthquake sequence, taking advantage of dense local records.

The previous studies of earthquake rupture processes commonly revealed as rupture directivity have been conducted for earthquakes in seismically active regions such as California and Japan due to dense earthquake monitoring (e.g., Tan and Helmberger, 2010; Wang and Rubin, 2011; Kane et al., 2013; Yoshida et al., 2019). In the recent years, detailed studies of rupture complexities and directivity of small‐to‐medium‐sized earthquakes that were induced by fluid injection have been reported largely due to the deployments of dense seismic stations and advanced analysis methods (e.g., Folesky et al., 2016; López‐Comino and Cesca, 2018; Király‐Proag et al., 2019; Wu et al., 2019). It should be useful to understand rupture directivity for small shocks so that the source‐scaling relations from small‐to‐large earthquakes can be more clearly examined.

In this study, we carry out a detailed analysis of rupture directivity for an $ML$ 2.2 microearthquake in a stable continental region to shed light on the rupture process of small‐sized earthquakes. We thus attempt to contribute to ongoing studies of the rupture process for small‐sized events that may suggest that the rupture directivity could be a size‐invariant characteristic of the earthquake rupture.

In the following sections, we report our waveform modeling for the $ML$ 2.2 mainshock with records from seismographic stations in high frequencies (0.3–2.5 Hz) and inversion for the focal mechanism. We located the sequence with small location uncertainties and well‐constrained focal depths using close‐in local stations. Precisely relocated microearthquakes align along east‐southeast–west‐northwest, suggesting the fault geometry that is consistent with the focal mechanism. Finally, we report our observation of the rupture directivity of the $ML$ 2.2 Gwangyang microearthquake. Our directivity analysis consists of (1) retrieving the apparent source time functions (ASTFs) of the $ML$ 2.2 microearthquake by employing the empirical Green’s function (EGF) technique, and (2) applying the waveform stretching method to determine the best‐fit parameters‐rupture propagation direction and rupture speed, assuming a unilateral rupture. The observed ASTF varies depending on the rupture direction and the rays leaving the source in the 3D focal sphere. The $ML$ 2.2 mainshock unilaterally ruptured in the slip‐parallel direction within the fault plane.

We retrieved waveforms from 55 stations within a 100 km radius from the $ML$ 2.2 mainshock (Fig. 1b). The sensors at those stations consist of 19 broadband, 1 short‐period (⁠$T0=2s$⁠; BAU), and 35 accelerometers, and the stations provided spatial coverage with an interstation distance of ∼15 km. The national seismic network, operated mainly by KMA, has been expanded rapidly since 2016 (primarily accelerometers as of 2021; see Data and Resources). All accelerometers are operated with a ±0.5g full‐scale setting, providing adequate‐quality data even for microearthquake studies.

We modeled $ML$ 2.2 mainshock waveforms recorded at six stations within ∼50 km to determine the focal mechanism and constrain its focal depth. We could model three‐component records in a relatively high‐frequency band of 0.4–2.5 Hz for five broadband and accelerographic stations, and 0.3–1.0 Hz for a short‐period station (BAU) with an adequate signal‐to‐noise ratio (SNR). Synthetic waveforms were calculated with the frequency–wavenumber integration method (Bouchon, 1981; Saikia, 1994) with a three‐layer crustal model of the southern Korean Peninsula (Kim et al., 2011). For station GMNA (15.4 km, azimuth = 95.6°) path, we used a two‐layer crustal model with a 2 km thick low‐velocity surface layer that fits the vertical and radial components that show S‐to‐P converted phase (Sp) arrival ∼0.4 s ahead of S (see Fig. 2 and Kim et al., 2022). We also tested global models, IASP91 and AK135 (Kennett et al., 1995), which include a 20 km‐thick upper crust and 15 km‐thick lower crust. These models produced synthetics with prominent reflected S at certain stations (e.g., MNDB) that were not observed from a likely source at ∼15 km depth; hence these models were not suitable for the study area.

We inverted the waveform data by allowing a small time shift between the observed and synthetics, and grid searched over double‐couple fault parameters, as given in Kim et al. (2010). We added P‐wave first motion data from eight stations to help constrain the focal mechanism. Figure 2 compares observed and synthetics calculated for the best‐fit focal depth of 16 km. The focal mechanism is strike‐slip faulting with a substantial oblique thrust component. We obtained a seismic moment of $M0=0.96±(0.40)×1013N·m$ (⁠$Mw$ 2.62) and the double‐couple fault parameters: strike = 114° ± 2°, dip = 70° ± 4°, and rake = 41° ± 4° (Fig. 2).

We located six earthquakes using the P and S arrival times with the regional model (Kim et al., 2011). Each event is located by using between 43 and 73 manually picked P and S arrivals, respectively, and four stations are situated at epicentral distances shorter or comparable to the focal depth. Hence, these microearthquakes are well located with small uncertainties and reasonably good depth control. Average root mean square residuals are 0.05–0.08 s; location uncertainties are 340–580 m horizontally and 670–980 m vertically at a 95% confidence level (Table 1). Focal depths range from 13.3 to 14.0 km. Events are generally located in the west‐northwest–east‐southeast direction, except for the smallest event 5 (Fig. 3a).

We performed double‐difference (DD) hypocenter relocation to delineate the causative fault by improving the relative location between events (Waldhauser and Ellsworth, 2000). We used 788 differential times (dt) between the event pairs derived from P‐ and S‐wave arrival times to solve for the relative event locations. The refined epicenters align along a nodal plane of the mainshock focal mechanism (Fig. 3b). The smallest event 5 seems to be slightly off the main cluster, possibly due to poor arrival picks on weak signals (Fig. 3b).

We hence performed waveform cross‐correlation (WCC) of P‐ and S‐wave arrival pairs to reduce the manual picking errors and improve the event location (Schaff and Waldhauser, 2005). We filtered waveforms in a 1–10 Hz band and used a 1 s‐long time window around both P and S arrivals. We retained 812 dt measurements with cross‐correlation coefficients (CC) exceeding 0.7. The DD relocation using both the catalog (phase) and WCC‐derived dt yields six events tightly aligned, and event 5 is now in the cluster with other events. When viewed in an across‐strike vertical section, hypocenters align along ∼70° SW dip direction (Fig. 3f), which is consistent with a mainshock nodal plane (Fig. 2). Good correspondence between a nodal plane striking 114° and event locations indicates that the nodal plane striking 114° and dipping 70° SW is the likely causative fault plane for the microearthquake sequence. Figure 3d shows the distribution of six hypocenters projected onto the inferred fault plane. Six events seem to migrate in the eastward and up–dip direction. The 95% confidence error bars from the bootstrap method (Efron, 1982) imply that the observed feature is robust with respect to the relative location uncertainty.

We first analyzed the mainshock rupture directivity by using the observed $τ(θi)$ at various azimuths. We assumed that the $τ(θi)$ is equal to the ASTF duration retrieved from observed records. We used the EGF method to retrieve the ASTF, because the EGF method minimizes site and path effects on observed source pulses (e.g., Hough, 1997).

We examined all the five small events as EGF for the target mainshock and found that event 1 (⁠$ML$ 1.3)—a foreshock, was the best EGF among the five candidate events. All the five EGF events have a common focal mechanism with the target event (Fig. S2). Event 5 is too weak with low SNR, and event 6 produced poorly defined ASTFs. Events 1, 2, and 4 yielded fairly good ASTF (Figs. S3 and S4) and are all good EGF events in terms of variance reductions (VRs). However, the shape and duration of the ASTFs from event 1 are better defined due to the adequate amplitude ratio between EGF and the target event (Table 1).

We used transverse‐component S in 2 s time windows, starting 0.5 s before the S arrival. We applied a 1 Hz high‐pass filter to suppress low‐frequency noise and performed multitaper deconvolution to obtain the ASTF of the target event (Prieto, 2022). We obtained well‐resolved ASTFs at 21 out of 55 stations examined (Fig. 4a). The ASTFs at azimuth range of 270°–310° show a long duration of ∼0.1 s, whereas ASTF durations at azimuth range of 80°–100° are the shortest (∼0.05 s). This azimuthal dependence of rupture duration indicates rupture propagation toward the east‐southeast.

The ASTFs at stations HAIA (42.9 km, azimuth = 97.6°), GACA (54.5 km, azimuth = 76.1°), and YGJA (69.1 km, azimuth = 86.9°), which are situated close to the presumed rupture direction (east–southeast), exhibit a small pulse immediately preceding the larger pulse (Fig. 4a). These may represent the existence of rupture initiation phase (e.g., Király‐Proag et al., 2019), but low SNRs of EGFs at these stations do not allow us to identify these pulses from the raw waveforms clearly.

We constrained that rupture propagates on the fault plane (strike 114°, dip 70°), which is well defined from both focal mechanism (Fig. 2) and event relocation (Fig. 3), and solved for $ϕ$⁠, $γ$⁠, and $Vr/Vc$⁠. Because the rupture vector is confined to the fault plane, $γ$ can be uniquely determined from $ϕ$⁠. Hence, the best‐fitting $ϕ$ and $Vr/Vc$ are found by grid search such that the L2 misfit between theoretical and observed $sij$ is minimized. Two ASTFs (HAIA and YGJA) with low VRs were excluded for the analysis (Fig. S3). We obtained $ϕ=98.5°±7.9°$⁠, $γ=36.3°±11.8°$⁠, and $Vr/Vc=0.72±0.06$ from the grid search (Fig. 4c,d), with bootstrap uncertainty (Fig. S5). Observed ASTFs aligned by angle $θ$ (angle between the rupture direction and ray takeoff vector) show a clear broadening of ASTFs with increasing angle, which is consistent with the assumed unilateral rupture model whose theoretical durations are plotted in Figure 4b.

In case when we grid‐searched rupture direction parameters without constraining the rupture on the fault plane, we obtained $ϕ=95.7°±6.7°$⁠, $γ=34.4°±17.2°$⁠, and $Vr/Vc=0.72±0.06$ (Figs. S5, S6). This is close to the constrained case, suggesting that we can identify the fault plane (as rupture plane) between two nodal planes of the focal mechanism by using the stretching method, even in the case with no knowledge of a well‐defined fault plane for a small earthquake such as $ML$ 2.2 Gwangyang microearthquake.

We calculated an average rupture length of 215 m (⁠$σ=38m$⁠) from $τ(θi)$ measurements at all stations, assuming $Vc=3.6km/s$⁠. The estimated rupture length roughly coincides with the spatial extent of relocated six events (∼150 m). Although we deal with very short ASTFs of duration less than 0.12 s, our directivity analysis is less affected by the limited temporal resolution, as most ASTFs (17 out of 19) are significantly broader than the minimum resolvable delta function (Fig. 4b).

Previous rupture directivity studies often assumed horizontal rupture. They analyzed the variation of ASTF duration with azimuth because of limited station coverage to resolve the vertical component of rupture and for simplicity (e.g., Kane et al., 2013; Folesky et al., 2016). In this study, we used stations with a range of azimuths and takeoff angles, and alleviated the assumption of horizontal rupture propagation. The rupture dip (⁠$γ$⁠) of 36.3° becomes 39.1° when it is projected on to the fault plane, and it is close to the rake angle of the focal mechanism (∼41°; Fig. 3d), implying that the rupture propagated in a nearly slip‐parallel direction (i.e., mode II rupture; Rubin, 2002).

The reported 2021 Gwangyang microearthquake sequence appears to define a line source (Fig. 3d) rather than a rupture area. Based on obtained rupture length, we attempt to estimate the possible rupture area of the mainshock. We obtained the average slip (u) of 1.53 cm by using an empirical relationship between rupture length and average slip for strike‐slip earthquakes in stable continental regions, that is $log10$$u=log10$L–4.149 (Leonard, 2014). Assuming a rectangular rupture, its width can be estimated from the seismic moment, $M0=μuA=μuLW$⁠, in which $μ$ is the rigidity (∼30 GPa), A is the rupture area (L × W), L is the rupture length (215 m), W is the rupture width, and u is the average slip (final offset, 1.53 cm). The width and the rupture area are ∼98 m and $0.021km2$⁠, respectively (Fig. 3d). It is likely that the rupture initiated at the lower‐left corner of the rupture area and propagated to the upper right side of the mainshock or the aftershocks may have occurred on a neighboring fault plane different from the mainshock rupture plane.

The recently expanded seismic network with good spatial coverage in the southern Korean Peninsula enabled a detailed investigation of the 2021 $ML$ 2.2 (⁠$MW$ 2.6) Gwangyang earthquake sequence. The mainshock focal mechanism determined by modeling high‐frequency waveforms at local distances indicates strike‐slip faulting with a substantial oblique thrust component. The best‐fit waveform data constrained the source at 16 ± 2 km depth. Precise DD relocation of six events in the sequence aligns the events along the west‐northwest–east‐southeast direction and suggests that the nodal plane striking 114° and dipping 70° to the southwest is the fault plane. Detailed ASTF analysis reveals that the mainshock unilaterally ruptured in the slip‐parallel direction within the causative fault plane. The rupture initiated at the lower left of the rupture area, and propagated eastward and up‐dip direction along the aftershock hypocenters. The estimated rupture length is ∼215 m based on a rupture speed of $0.72Vc$ (∼2.6 km/s), which agrees with the spatial extent of mapped hypocenters.

Waveform data of the Korea Meteorological Administration (KMA) were downloaded from the National Earthquake Comprehensive Information System website (https://necis.kma.go.kr) upon approval; data from stations GSU, HWSB, and JRB, operated by the Korea Institute of Geoscience and Mineral Resources (KIGAM) were acquired through written request. As of December 2021, KMA operates 282 stations, whereas KIGAM maintains 39 stations, which form the National Seismographic Network of Korea. We used Seismic Analysis Code (https://ds.iris.edu/ds/nodes/dmc/software/downloads/sac/), ObsPy (https://docs.obspy.org/) and Python library Multitaper (https://github.com/gaprieto/multitaper/) for deconvolution. All websites were last accessed in July 2022. The supplemental material for this article includes one table and six figures.

The authors declare that there are no conflicts of interest recorded.

This work was supported by Brain Pool Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Science and Information & Communication Technology (ICT) (2020H1D3A2A02047949). Y. Kim and M.‐S. Seo acknowledge support from the NRF grant funded by the Korean government (MSIT) (2022R1A2C1003006). The authors thank the reviewers for their valuable comments that helped us to improve the manuscript. The authors appreciate the staffs at the Korean Meteorological Administration, and the Korea Institute of Geoscience and Mineral Resources, who operated seismographic stations and provided data used in this study. The authors also thank colleagues at Seoul National University, J. Y. Park, S. Han, and Y. O. Son, who helped us analyze data.