Radiated seismic energy from the 2021 ML 5.8 and ML 6.2 Shoufeng (Hualien), Taiwan, earthquakes and their aftershocks

We calculated radiated seismic energy (ES), seismic moment (M0), and moment magnitude (MW) and then determined the ES–ML, ES–M0, M0–ML, and MW–ML relations for the 2021 Shoufeng earthquake sequence (2.5 < ML < 6.3), where ML is the local magnitude. Notably, a crossover magnitude was detected as ML = 4.0 for the earthquake sequence. For ML < 4.0, we obtained logM0 ∝ ML, MW ∝ 0.67ML, and a low ES/M0, indicating a low average stress drop; for ML > 4.0, logM0 ∝ 0.67ML, MW ∝ ML, and a high ES/M0 were present, and then there was a high average stress drop. These derived relations implied that source duration (T) is independent of M0 for ML < 4.0. Moreover, the M0 ∝ T3 relation seemed able to interpret those relations for ML > 4.0. Nevertheless, the ES–ML relation remains logES ∝ 2.0ML for 2.5 < ML < 6.3. From this study, the derived relations could not predict the source parameters for ML > 6.3 events. This might indicate that ML saturates beyond ML 6.3. Through such analyses, we not only established the relations among source parameters but also elucidated the basic physics of the earthquake sequence. logEs is proportional to 2.0ML for the Shoufeng earthquake sequence in Taiwan. A crossover magnitude detected at ML = 4.0 divides the MW–ML relation into two parts. For ML > 4.0, logM0 ∝ 0.67ML and MW ∝ ML; for ML < 4.0, logM0 ∝ ML and MW ∝ 0.67ML. logEs is proportional to 2.0ML for the Shoufeng earthquake sequence in Taiwan. A crossover magnitude detected at ML = 4.0 divides the MW–ML relation into two parts. For ML > 4.0, logM0 ∝ 0.67ML and MW ∝ ML; for ML < 4.0, logM0 ∝ ML and MW ∝ 0.67ML.


Introduction
Since the 2018 MW 6.4 (M L 6.2) Hualien earthquake occurred (cf. Rau and Tseng 2019), high seismicity in the Hualien region has persisted for a while (also refer to the Central Weather Bureau (CWB); https:// scweb. cwb. gov. tw/ en-US). More than a thousand aftershocks (M L ≥ 2) followed the 2018 event within 3 months. Close behind the 2018 event, the Xiulin earthquake with M W 6.2 (M L 6.3) occurred approximately 18 km southeast of the 2018 event, on April 18, 2019 (cf. Lee et al. 2020 Fig. 1). Approximately 1 year later, on April 18, 2021, an M L 5.8 earthquake occurred again in the Shoufeng region, Hualien; 3 min later, an M L 6.2 earthquake (refined to be 6.26) occurred almost at the same location as the M L 5.8 event (Fig. 1).
Subsequently, approximately 500 aftershocks occurred. At first, the aftershocks were distributed around the M L 5.8 and M L 6.2 events; later, the aftershocks were distributed eastwards and became shallower (≤ 10-km depth). The two moderate-sized earthquakes seemed to feature a similar plane along a NE-SW strike based on the focal mechanisms (Fig. 1). The cross-section of the earthquake sequence along line BB′ in Fig. 1 indicated a west-dipping plane, which might be the fault plane for the two 2021 events. In this study, we did not analyze the rupture features of the M L 5.8 and M L 6.2 Shoufeng earthquakes but investigated the source parameters of the earthquake sequence, including the seismic moment (M 0 ), radiated seismic energy (E S ), and moment magnitude (M W ). Examining the M L 5.8 and M L 6.2 Shoufeng earthquakes' rupture feature would be a separate issue. For earthquakes, E S is one of the key source parameters needed to understand dynamic rupture because E S is strongly related to static stress drop (Δσ), dynamic stress drop (Δσ d ), and fracture energy (E g ) (cf. Kanamori and Heaton 2000;Kanamori 2001). Here, Δσ and Δσ d are defined as σ 0 -σ 1 and σ 0 -σ f , where σ 0 , σ 1 , and σ f are the initial stress, final stress, and frictional stress, respectively; Eg is the erergy that extends the fault during the earthquake. A derivative parameter called scaled energy, which is defined as the ratio of E S to M 0 , is often used to indicate the friction drops during the earthquake rupture (Kanamori and Heaton 2000) and also represents the average stress drop on the fault plane (Houston 1990;Kanamori et al. 1993). For large earthquakes, E S / M 0 is almost equal to a constant, 5 × 10 −5 (Kanamori 1977); however, Ide and Beroza (2001) reported a E S /M 0 of 3 × 10 −5 for a wide M 0 range. For small earthquakes, E S /M 0 is approximately 10-100 times smaller than large earthquakes (Abercrombie 1995; Mayeda and Walter 1996;Kanamori and Heaton 2000). Generally, E S can be calculated through the integral of velocity-squared seismograms after certain appropriate corrections are made (cf. Boatwright and Choy 1986;Newman and Okal 1998;Kanamori et al. 1993;Venkataraman et al. 2006). Although E S only represents a portion of strain energy (W) released during an earthquake and is approximately 10% less than W (McGarr 1999), it is a key source parameter for understanding earthquakes' dynamic ruptures. Gutenberg and Richter (1956) constructed a relationship between surface magnitude (M S ) and logE S as logE S = 1.5M S + 4.8 (E S in Nm). This indicates that logE S is proportional to 1.5M S . Such a relation is capable of being theoretically verified (Kanamori and Anderson 1975). However, most studies have revealed that logE S is proportional to 2.0M L (Thatcher and Hanks 1973;Seidl Fig. 1 (Left) Map showing the tectonic setting of the Taiwan region. Several moderate and large earthquakes compared with the Shoufeng earthquake sequence are also plotted. (Middle) Seismicity in the Shoufeng region, Hualien County, Taiwan. The red-tone circles denote aftershocks occurring before June 8, 2021; the blue-tone circles show aftershocks occurring after June 8, 2021. The stars indicate the two 2021 Shoufeng earthquakes and several moderate-sized earthquakes occurring around the source area of the 2021 events. Also included are the focal mechanisms reported from the BATS (IES 1996). (Right) AA′ and BB′ profiles reveal that these events are distributed at different depths (< 25 km). The cross-section of these events along line BB′ shows a west-dipping plane, which might be the fault plane of the two 2021 events and Berckhemer 1982; Kanamori et al. 1993;Dineva and Mereu 2009). Deichmann (2018a) further confirmed that logE S scales as 2.0M L by using numerical simulations. In Taiwan, Huang (2003) first used regional broadband seismic data to establish the relationship between logE S and M L , indicating that logE S ∝ 2.0M L . Chan et al. (2020) also obtained logE S ∝ 2.0M L for the 2019 Xiulin earthquake sequence. In addition, the observed scale between M W and M L (M W :M L ) ranges from 1:0.67 to 1:1.5 (Hanks and Kanamori 1979;Margaris and Papazachos 1999;Wu 2000;Wu et al. 2005;Chen et al. 2009;Grünthal et al. 2009;Sargeant and Ottemöller 2009;Bethmann et al. 2011;Zollo et al. 2014;Deichmann 2017;Munafò et al. 2016;Malagnini and Munafò 2018;Mereu 2020). For Taiwan's earthquakes, Wu et al. (2005)  Regarding the Shoufeng earthquake sequence, what is the relation between M W and M L ? Hence, in this study, we not only calculated E S for the 2021 Shoufeng earthquake sequence from near-field broadband seismograms but also constructed the logE S -M L relationship to infer whether logE S ∝ 2.0M L does or does not hold. Furthermore, several relations between source parameters were also discussed, namely the E S -M 0 , logM 0 -M L, and M W -M L relations. Such analysis can be used to systemically examine the relationships among source parameters for earthquakes in Taiwan in the future.

Data
Velocity seismograms from the Broadband Array in Taiwan for Seismology (BATS; Institute of Earth Sciences (IES), Academia Sinica, Taiwan 1996) were adopted to analyze the radiated seismic energy (E S ) of the two moderate-sized earthquakes and their aftershocks occurring in the Shoufeng region, Hualien County, Taiwan. To reduce the effect of complicated structures on wave propagations, we only used data with epicentral distances of less than 50 km. Therefore, the seismic waves propagated approximately within a half-space material (Kanamori 1990). E S is mainly expressed in the S-wave, so E S is generally estimated through the S-wave (cf. Boatwright and Fletcher 1984;Kanamori et al. 1993;Venkataraman et al. 2006). However, it would be not easy to determine the S-wave train when using near-field seismograms. For this reason, we derived E S from only P-waves. Finally, the E S values from P-waves would be corrected to obtain the total energy by introducing energy partitioning between the P and S waves (Boatwright and Choy 1986;Newman and Okal 1998). Because of noise interference, for stably calculating E S , seismograms from various M L were filtered by different frequency bands. Here, a noise test was first made to determine the frequency band for filtering. Figure 2 provides an example of a noise test at station NACB for 2.5 < M L < 4.5 events. We summarized the noise test from each station. Finally, for events with M L ≥ 4.0, the frequency band for filtering was 0.1-10 Hz; for events with 3.5 ≤ M L < 4.0, we used 1-10 Hz to filter these seismograms; the events with 3.0 ≤ M L < 3.5 were filtered from 2 to 10 Hz; for events with 2.5 ≤ M L < 3.0, the band-pass filter with frequencies between 3 and 8 Hz was employed.

Radiated seismic energy (E S )
Radiated seismic energy from the P-wave, called E P S , is measured using the integration of the squared velocity records ( v(t) 2 ) from a given station (cf. Boatwright and Fletcher 1984;Kanamori et al. 1993), as follows: where α and ρ are the P-wave velocity and the density at the source area, respectively; r is the hypocentral distance, representing the geometrical spreading for seismic-wave propagation; t a and t b , determined manually, are the time band for integration; R is the radiation pattern depending on the focal mechanism; R is the average radiation pattern over a focal sphere and equates to 0.52 for the P-wave (Aki and Richards 2002). In Eq. (1), R 2 R 2 is assumed to be 1.0 by using the average radiation patterns for all stations, as proposed by Kanamori et al. (1993).
Before calculating E P S , several corrections to the P-waves needed to be made, including correcting the effects of the free surface and attenuation. The free surface effect was calculated using the incident angle of the P-wave to the receiver (cf. Okal 1992;Aki and Richards 2002), and t * = 0.029 was employed to correct the P-wave attenuation from source to receiver following the work of Chan et al. (2020) (also see 5.). Finally, the total radiated seismic energy (E S ), the sum of the P-wave energy ( E P S ) and the S-wave energy ( E S S ) was calculated as follows: where q is the ratio of the S-wave energy to the P-wave energy and is defined as 1.5 α β 5 and α and β are the P-wave and S-wave velocities at the source, respectively. For a Poisson material, q is equal to 23.4 due to α = √ 3β . Many observations also show q varying between 9 and 25 (cf. Shearer 2009). In this study, the travel times of the P-waves and S-waves were estimated from source to station for each earthquake to obtain the average α and β values, which were used to calculate q . Finally, for a given earthquake, we calculated the average of the available E S values from each station to be the E S of the event.
In theory, E S should be estimated over a frequency range of 0-∞ Hz; however, in practice, E S is generally calculated at a finite frequency bandwidth and is apt to be underestimated (Di Bona and Rovelli 1988;Ide and Beroza 2001;Wang 2004). Therefore, the finite frequency bandwidth limitation needs correcting, especially for high-frequency parts. Such correction is dependent on the corner frequency.

Corner frequency (f c ) and seismic moment (M 0 )
For correcting the finite frequency bandwidth limitation in the calculation of E S , we first derived the corner frequency (f c ) from the displacement and velocity P-waves. According to the ω −2 source model (Aki 1967;Brune 1970), Andrews (1986) derived f c as follows: (3) f c = 1 2π the observed velocity and displacement P-waves, respectively. t a and t b are the time band for integration. In the ω −2 source model, f c is defined as the crossover frequency between the low-and high-frequency spectra. For f < f c , the spectrum is approximately a constant; for f > f c , the spectrum decays with f −2 (also see 5.). Seismic moment (M 0 ) is expressed in the following form: is the low-frequency spectral level (Andrews 1986); �R� = 0.52 for the P-wave; α and ρ are the P-wave velocity and density at the source area, respectively. Before calculating 0 , the effects due to the free surface and attenuation on the P-waves were corrected. For a given earthquake, we averaged the available M 0 calculated from each station to be the M 0 of this event.

Finite frequency bandwidth correction
On the basis of Parseval's theorem, the integral of the square of a time function is equivalent to that of the square of its spectrum, that is, Because of the finite frequency bandwidth used for the observed P-waves, the real integrals are as follows: where f l and f u are the integral range with f l < f c < f u , and κ v (≤ 1.0) is the proportion of finite E S (observation) to full E S . Here, κ v is expressed in the following form (cf. Ide and Beroza 2001;Wang 2004): Ide and Beroza (2001) indicated that over 80% of the radiated seismic energy is created from high frequencies larger than the corner frequency. Because E S should be calculated over an infinite frequency, the observed E S from a finite frequency should be corrected by κ v , that is E SC = E S /κ v , where E SC is the corrected E S . In other words, through κ v , we can restructure E S for a given earthquake.
Page 5 of 10 Hwang et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:19 4 Results and discussion Figure 3 displays a plot of logE S versus M L for 2.5 < M L < 5.5 and the M L 5.8 and M L 6.2 events. The M L values were from the earthquake catalog of the CWB, which routinely reports earthquake parameters in Taiwan. From Fig. 3, E S increased with increasing M L , and a linear relationship between logE S and M L was present for 2.5 < M L < 5.5 events as logE S = 1.99 M L + 2.06 (E S in Nm). The empirical relation initially derived from 2.5 < M L < 5.5 events seemed to be able to predict E S for the M L 5.8 and M L 6.2 events. Therefore, by incorporating the two 2021 events into the regression, we obtained the logE S -M L relation as follows.

E S versus M L
Our results agreed well with logE S ∝ 2.0M L and are similar to the logE S -M L relation, that is, logE S = 1.96M L + 9.05 (Es in erg), proposed by Kanamori et al. (1993) for 1.5 < M L < 6.0 earthquakes in southern California. Likewise, Chan et al. (2020) analyzed E S for the 2019 Xiulin earthquake sequence to reveal logE S ∝ 2.0M L for 2.5 < M L < 4.3 events. Many previous studies have also indicated that logE S ∝ 2.0M L holds at different M L ranges, from 1 to 7 (Thatcher and Hanks 1973;Seidl and Berckhemer 1982;Kanamori et al. 1993;Dineva and Mereu 2009). M L is determined using seismograms recorded by Wood-Anderson seismographs (Richter 1935) and proportional to M 0 (Randall 1973;Deichmann 2006). Then, logE S will scale as 1.5M L following source self-similarity (Aki 1967; Kanamori and Anderson 1975). However, from numerical simulations for − 1 ≤ M L ≤ 6, Deichmann (2018a) further stated what conditions are valid (7) log E S = (1.98 ± 0.05) M L + (2.11 ± 0.17)(E S in Nm) for 2.5 < M L < 6.3.
for logE S ∝ 2.0M L . For small earthquakes (M L ≤ 3.0), logE S ∝ 2.0M L is inevitable because earthquake duration approaches constant; for larger earthquakes, it happens to be consistent with the logE S -2.0M L relation because of the effect of the Wood-Anderson response when calculating M L (Deichmann 2018a). Chan et al. (2020Chan et al. ( , 2021 investigated the logE S -M L relation for the 2019 Xiulin aftershocks (2.5 ≤ M L < 4.3) to obtain logE S ∝ 2.0M L as well as to find the constant-tending duration for those small events. Although Malagnini and Munafò (2018) noted a crossover magnitude at approximately M L = 4.3 for the M W -M L relation, where M W is the moment magnitude (Kanamori 1977;Hanks and Kanamori 1979), logE S was still proportional to 2.0M L for a wide M L range (1-7; Thatcher and Hanks 1973;Seidl and Berckhemer 1982;Kanamori et al. 1993;Dineva and Mereu 2009). Kanamori et al. (1993) concluded that M L < 6.5 events do not saturate with increasing E S because the derived relationship could not predict the E S of the 1992 M L 6.8 (Mw 7.3) Landers earthquake. Similarly, from the 2019 Xiulin aftershocks, the logE S -M L relationship, analyzed by Chan et al. (2020), also failed to predict the E S of the 2019 M L 6.3 Xiulin earthquake (mainshock). Figure 3 also illustrates the E S values for several M L > 6 earthquakes in Taiwan (Venkataraman and Kanamori 2004;Wang, 2004;Hwang 2012;Hwang et al. 2019Hwang et al. , 2022; Table 1). The derived E S -M L relation does not seem to satisfactorily predict the E S values for M L larger than 6.3 other than the E S of the 1999 Chi-Chi earthquake from Wang (2004). This also implies that M L saturates for events with M L > 6.3. Besides, whether a difference exists between the eastern and western fault systems for the E S -M L relation is worth investigating further.  (Kao et al. 1998;IES 1996). Such analysis has the advantage of calculating M 0 for small earthquakes. E S is a dynamic source parameter and can describe the stress state and friction drops during an earthquake (Houston 1990;Kanamori et al. 1993;Kanamori and Heaton 2000). A low E S /M 0 shows that the friction drops gradually during faulting; however, a high E S /M 0 denotes the fact that the friction drops rapidly. Then, it would result in various stress drops, including static and dynamic stress drops. Therefore, the E S /M 0 ratio provides a measurement of the average stress drop on the fault. Fig. 3 Relation between logE S and M L , where E S is the radiated seismic energy and M L is the local magnitude. The logE S -M L relation has a slope of 1.98; therefore, logE S is almost proportional to 2.0M L . For comparison, also included in this plot are the E S of several earthquakes in Taiwan (IRIS DMC 2013; Venkataraman and Kanamori 2004;Wang 2004;Hwang 2012;Hwang et al. 2019Hwang et al. , 2022; also see Table 1). The blue dashed lines denote the standard deviations Page 6 of 10 Hwang et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:19 Figure 5 illustrates the plot of E S vs. M 0 . As can be seen, E S increases with increasing M 0 . Additionally, E S /M 0 also varies with M L for the earthquake sequence, as displayed in the inset of Fig. 5. On average, for 2.5 < M L < 3.5, 3.5 < M L < 4.5, and M L > 4.5, E S /M 0 is 5.3 × 10 −6 , 6.3 × 10 −5 , and 2.7 × 10 −4 , respectively. Except that E S / M 0 = 6 × 10 -5 for 3.5 < M L < 4.5, the others differ from the global observations of 5 × 10 −5 (Kanamori 1977) and 3 × 10 −5 (Ide and Beroza 2001). Although the ratio E S /M 0 in this study reveals differences to the global values, variations in E S /M 0 with M L (or M W ) have also been observed in many other studies (e.g., Abercrombie 1995; Mayeda and Walte 1996; Kanamori and Heaton 2000). When ignoring the fracture energy, the corresponding Orowan stress drop, defined as �σ = 2µ E S M 0 , where µ is the rigidity (Orowan 1960;Houston 1990;Kanamori et al. 1993), varies from 3 (small events) to 40 and 160 (larger events) bar on average. In other words, the E S /M 0 ratio reveals the strength heterogeneities on the fault plane.

E
In theory, the difference in Es/M 0 between M L < 4.0 and M L > 4.0 events might be interpreted in a simple stressdrop model (cf. Kanamori and Heaton 2000). Following   Kanamori (1977) for large earthquakes, and E S /M 0 = 3 × 10 −5 is estimated by Ide and Beroza (2001) for a wide M 0 range. The insert illustrates the relation between E S /M 0 and M L . On average, E S /M 0 = 5.3 × 10 −6 , 6.3 × 10 −5 , and 2.7 × 10 −4 for 2.5 < M L < 3.5, 3.5 < M L < 4.5, and M L > 4.5, respectively. This reveals that the Orowan stress drops are approximately 3, 40, and 160 bar for different M L range, respectively Page 7 of 10 Hwang et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:19 Kanamori and Heaton (2000), the Es/M 0 ratio is proportional to (1 − Dc/D), where Dc is the critical distance and D is the total slip on the fault plane. For small earthquakes, Dc is comparable with D to have a large Dc/D, resulting in a low Es/M 0 ratio; by contrast, for large earthquakes, there is a small Dc/D due to the fact that frictional melting or fluid pressurization is likely to occur during an earthquake rupture. Then, it would have a high E S /M 0 . In other words, for a given M 0 , a high Dc/D has less E S released; reversely, a low Dc/D has more Es released. In fact, the value of Dc is not apt to be observed, especially for small earthquakes. The discrepancy in E S / M 0 shows different rupture dynamics for small and large earthquakes.

M 0 versus M L and M W versus M L
Both M 0 and M L represent the earthquake size. Therefore, logM 0 should scale as M L . That is, logM 0 increases with increasing M L . Figure 6A illustrates logM 0 varying with M L for the Shoufeng earthquake sequence. Another scale denoting earthquake size is the moment magnitude (M W ), defined as M W = (2/3)logM 0 − 6.07 (M 0 in Nm; Kanamori 1977;Hanks and Kanamori 1979). Figure 6B displays the relationship between M W and M L . M W is generally smaller than M L, as noted by Wu et al. (2005).
Because the M 0 -M L and M W -M L relations are two sides of a single entity, some interpretations mainly depend on the M W -M L relation. From Fig. 6, a magnitude turning point can be detected at M L = 4.0. Subsequently, the M W -M L relation was determined as follows: Equations (8) and (9)   Page 8 of 10 Hwang et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:19 of M L (Richter 1935), M L is proportional to logA; then, logM 0 ∝ (3/2)M L . Finally, M L ∝ M W . Subsequently, if T is assumed to be independent of M 0 , logM 0 will scale as logA, leading to logM 0 ∝ M L . Finally, M L ∝ (2/3)M W is obtained. Conversely, if M L ∝ M W is observed, it implies that M 0 ∝ T 3 ; if M L ∝ (2/3)M W is obtained, it implies that T is invariant. As a rule, the M 0 ∝ T 3 relation is based on the source self-similarity, in which the static stress drop (Δσ) is not related to M 0 (cf. Aki 1967; Kanamori and Rivera 2004). Kanamori and Rivera (2004) proposed that the M 0 ∝ T 3 relation may not necessarily obey the source self-similarity. Another possibility for M 0 ∝ T 3 is under the condition that ΔσVr 3 = constant, where Vr is the rupture velocity (Kanamori and Rivera 2004;Hwang et al. 2020). That is, both source self-similarity and ΔσVr 3 = constant can make M 0 ∝ T 3 hold. Of course, to verify this, it is necessary to probe into how T varies with M 0 from observations.

Conclusions
By systemically analyzing the source parameters for the 2021 Shoufeng earthquake sequence, we obtained several regression relations among these parameters, including the E S -M L , E S -M 0 , M 0 -M L , and M W -M L relations. A crossover magnitude, detected at M L = 4.0, divided the M 0 -M L and M W -M L relations into two parts, but logE S was still proportional to 2.0M L for 2.5 < M L < 6.3. Such results indicated variation between source duration and M 0 for the 2021 Shoufeng earthquake sequence. On average, a relatively large E S /M 0 ratio and high Orowan stress drop were present for M L > 4.5 events. That is, the stress state at the source area exhibited variation with earthquake size.