Dominant periods and memory effect of the 2021 earthquake swarm in Hualien, Taiwan

An earthquake sequence that occurred in the Hualien area during April 7 to August 30, 2021 has been recognized as two swarms by the Central Weather Bureau. Its earthquakes with ML ≤ 6.2 (ML = local magnitude) and focal depths ≤ 25 km were located in an area from 23°46’ N to 24°04’ N and from 121°25’ E to 121°42’ E. The Morlet wavelet technique is applied to analyze the dominant periods of temporal variations in numbers of daily events for the earthquake sequence in two magnitude ranges, i.e., ML ≥ 3 and ML ≥ 4. Results show that the dominant periods are 30.8 and 38.0 days when ML ≥ 3; while the dominant period does not exist when ML ≥ 4. The fluctuation analysis technique in the natural time domain is used to study the memory effect in the swarm for two magnitude ranges, i.e., ML ≥ 3 and ML ≥ 4. Calculated results show that the memory effect is stronger for the time sequence of magnitudes than for that of inter-event times and higher for ML ≥ 3 earthquakes than for ML ≥ 4 events. Consequently, only the short-term corrected memory effect was operative in the earthquake sequence of the Hualien swarms. To explain the data of the 2021 earthquake swarm in Hualien To describe the Morlet wavelet and fluctuation methods Results show that the memory effect is stronger for the time sequence of magnitudes than for that of inter-event times and higher for ML ≥ 3 events than for ML ≥ 4 ones; and the short- term corrected memory effect was operative in the Hualien swarms To explain the data of the 2021 earthquake swarm in Hualien To describe the Morlet wavelet and fluctuation methods Results show that the memory effect is stronger for the time sequence of magnitudes than for that of inter-event times and higher for ML ≥ 3 events than for ML ≥ 4 ones; and the short- term corrected memory effect was operative in the Hualien swarms


Introduction
Taiwan is situated along the collision boundary between the Philippine Sea plate and the Eurasian plate (Tsai et al. 1977;Wu 1978;Tsai 1986). The former moves northwestward with a converging speed about 8 cm/year (Yu et al. 1997). The Philippine Sea plate has subducted underneath the Eurasian plate in northern Taiwan. This collision causes high seismicity in the region (Hsu 1961(Hsu , 1971Wang et al. 1983;Shin 1992;Wang 1988aWang , 1998bWang and Shin 1998). The Hualien area is located at the northern boundary of the two plates. Numerous large earthquakes have occurred in or near the area. Several main earthquake sequences are 1951 sequence (e.g., Chen et al. 2008;Lee et al. 2008), the 1986 sequence

Open Access
Terrestrial, Atmospheric and Oceanic Sciences *Correspondence: chenkc@earth.sinica.edu.tw (Liaw et al. 1986;Wang 1986, 1988;Wang 1988aWang , 1998bYeh et al. 1990), the 2002 sequence (e.g., Chen 2003;Chen et al. 2004), and the 2018 sequence (e.g., Hwang et al. 2018;Kou-Chen et al. 2018;Wen et al. 2018;Wu et al. 2019). The 2018 earthquake caused remarkable damage. Hence, it is important and significant to study the earthquakes occurring in the area for both scientific interest and social needs. The Central Weather Bureau (CWB) recognized two earthquake swarms that occurred in the Hualien area during April 7-August 30, 2021. Its earthquakes with M L ≤ 6.2 (M L = local magnitude) and focal depths ≤ 25 km occurred in an area from 23°46' N to 24°04' N and from 121°25' E to 121°42' E. The epicenters of the swarm seem to move from southwest to northeast and stopped almost on the south boundary of the 2018 M L 6.4 (M w 6.4) Hualien, Taiwan, earthquake sequence of February 6, 2018 (e.g., Hwang et al. 2018;Wen et al. 2018). Although the CWB recognized the whole earthquake sequence to include two swarms, we here only consider an earthquake sequence because the time difference between the so-called two swarms was very short. In order to explore the characteristics of the earthquake swarm, the evaluations of the dominant periods and the exploration of memory effect of the time sequence of the events will be made.
Although Fourier analysis is commonly applied to evaluate the dominant periods (or frequency) of a time series, it cannot provide temporal variation of the dominant periods. On the other hand, the wavelet transform can be used to analyze time series that contain non-stationary power at many different frequencies (Daubechies 1990). Hence, the wavelet analysis (e.g., Combes et al. 1989;Pyrak-Nolte and Nolte 1995) that is also known as the multi-resolution analysis is here taken into account. For this technique, a series of scaled and delayed oscillatory functions are used to decompose a time-varying signal into its non-stationary spectral components. Hence, the key advantage of wavelet analysis over traditional Fourier analysis is that the wavelet analysis provides information on how the spectral content varies with time delay. Wavelets also are advantageous over so-called windowed Fourier methods because with wavelets the relative accuracy of the delay and frequency remain constants over all of the delay-frequency parameter space. The application of wavelet analysis to geophysical problems can be seen in Torrence and Compo (1998). This method was applied to analyze the dominant periods of earthquake sequences in the Taipei Metropolitan Area (TMA) by Chen et al. (2015). Here, a non-orthonormal Morlet wavelet analysis (Morlet et al. 1982) is considered in this study.
In this study, we explore the correlation among events, in other words, we ask whether or not the earthquakes are correlated. Of course, it is necessary to further ask a question: Is the correlation between earthquakes long-term or short-term? Let n(t) be the number of earthquakes in an area at time t. When the changing rate of n(t) at time t, dn(t)/dt, is controlled only by n(t), the relationship between dn(t)/dt and n(t) can be represented by a linear 1-D difference equation: dn/dt = -λn(t). This equation gives a solution in the form of the exponential function, n(t) ~ exp(-t/λ), to show its temporal behavior. When dn(t)/dt is controlled not only by n(t) but also by the previous numbers, for example, n(t-δt), a memory effect exists in earthquakes. Hence, the relationship between dn(t)/dt and n(t) can be represented by a non-linear 1-D differential equation: dn/dt = -κn(t)n(t-δt). We have dn/dt = -κn 2 (t) as δt approaches zero. This gives a solution in the form of the power-law function, i.e., n(t) ~ κt −1 , to show its temporal behavior. Hence, power-law behavior of an earthquake sequence suggests the possible existence of a memory effect in earthquakes.
The radical problem now is in determining whether such a memory effect is long-term or short-term corrected in the earthquakes explored in this study. A longterm memory effect appears in several phenomena, for examples, in climate (e.g., Koscielny-Bunde et al. 1998), physiology (e.g., Peng et al. 1994), and the financial market (e.g., Liu et al. 1997). Lennartz et al. (2008Lennartz et al. ( , 2011 also assumed the existence of a long-term memory effect in earthquakes, especially for aftershocks. Wang et al. (2012a) applied three statistical models, the gamma, power-law, and exponential functions, to describe the single frequency distribution of inter-event times between two consecutive events for both shallow and deep earthquakes with M L ≥ 3 in the TMA from 1973 to 2010. Results show that the power-law function is the most appropriate one for describing the data points; while the exponential function is the least appropriate one. This suggests that the time series of earthquakes are not Poissonian, in other word, the correlation between two events exists. From the fluctuation analyses for the same data, Wang et al. (2012b) suggested that the M L ≥ 3 and M L ≥ 3.5 earthquakes are short-term corrected, while the M L ≥ 4 earthquakes are weakly corrected. Actually, it is difficult to exactly assign the time length for 'longterm' and that for 'short-term. ' The time lengths of the two terms can only be relatively defined. Considering the earthquake prediction as an example, the 'long-term' prediction means several tens or hundreds of years, while the 'short-term' means several months or few years. Since the total time length is about 146 days in this study, the 'short-term' may mean a time length of few days while the 'long-term' means the time length being longer than ten days.
The typical aftershock sequence has a strong memory effect, because it follows the Omori's power law. On the other hand, a swarm does not follow the Omori's power law. Meanwhile, unlike a typical earthquake sequence, a swarm is a set of many earthquakes with small differences in earthquake magnitudes. Hence, it is interesting to understand whether the memory effect exists in a swarm. Numerous authors (e.g., Kagan and Knopoff 1976;Wang 2013;Wang et al. 2017; and cited references therein) assume that the occurrences of earthquakes with M ≥ 7 are the pure Poisson process or the Poisson process with a weak memory effect. The degree of memory effect may increase from larger-sized earthquakes to smaller-sized events (e.g., Wang 2013; Wang et al. 2017). For the earthquake sequence of the Hualien swarms, the magnitude ranges from 3.0 to 6.2. Since the number of events with M L > 5 is small, we are interested in studying whether the memory effect exists in the earthquake sequences of two magnitude ranges, i.e., M L ≥ 3 and M L ≥ 4.
In order to study the existence of memory effect in the earthquake sequence, the fluctuation analysis (FA) technique (Koscielny- Bunde et al. 1998;Lennartz and Bunde 2009a) is the common choice. Essentially, those phenomena are assumed to be physically critical. Recently, the natural time is considered to be a good time system to represent critical phenomena (e.g., Varotsos et al. 2004Varotsos et al. , 2005Varotsos et al. , 2011Uyeda et al. 2009;Lennartz et al. 2008). Seismicity is also considered to be one of critical phenomena (e.g., Bak and Tang 1989;Main 1996;Turcotte 1997;Rundle et al. 2003). Hence, the temporal variation in earthquakes can be represented in natural time. Figure 1a shows the sequence of events (with magnitudes M i , i = 1, 2, 3 …, n + 1) in the conventional time domain. The interevent time (also denoted by inter-occurrence time in some articles) between events i and i + 1 is denoted by T i . In Fig. 1b, the earthquake sequence is represented in the natural time domain and denoted by the count, i, of an event. Hence, the inter-event time is just one unit, i.e., '1' , for all pairs of events in the natural time domain. Varotsos et al. (2011) introduced the natural time concept that is based on event counts as a measure of 'time' rather than the clock time. The conventional time is in the continuum of real numbers. On the other hand, natural time is not continuous and thus the values of natural time form a countable set of natural numbers. In the natural time domain each event is characterized by two terms, i.e., the natural time i that is the natural number of the i-th event and a physical quantity Q i . The natural time analysis has been applied to analyze complex time series and critical phenomena.
This study will focus on two parts: one for the evaluation of the dominant periods of time sequences in the numbers of daily events by using the Morlet wavelet analysis and the other for the exploration of memory effect of time sequences in the magnitudes and interevent times of the earthquake swarm by applying the Fluctuation analyses. As mentioned above, the Fluctuation analysis will be conducted in two magnitude ranges, i.e., M L ≥ 3 and M L ≥ 4. In order to study the possible difference in dominant periods for the earthquake sequences in the two magnitude ranges, the Morlet wavelet analysis will also be conducted in two magnitude ranges.

Data
Since 1991, the CWB has upgraded her old seismic network, by adding numerous new stations. This new network is named the CWB Seismic Network (CWBSN). In 1992 the TTSN was merged into the CWBSN. The earthquake magnitude of the earthquake catalogue has been unified to be the local magnitude. A detailed description about the CWBSN can see Shin (1992) and Shin and Chang (2005). At present, the CWBSN is composed of 72 stations, each equipped with three-component digital velocity seismometers. This network provides high-quality digital earthquake data to the seismological community. The local magnitude, M L , is used to quantify an earthquake by this network (Shin 1993). The earthquake data used in this study are directly retrieved from the CWB's data base (CWB 2021). Figure 2 displays the epicenters of 168 M L ≥ 3 earthquakes of the sequence that happened in the area 23°46' N to 24°04' N and from 121°25' E to 121°42' E. It seems that the events may be divided into two groups or two swarms. The first group occurred from April 7 to June 28 and the second one happened from July 1 to August  Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 30. For the first group, the earthquakes with M L ≤ 6.2 occurred in the area from 23°46' N to 23°57' N and from 121°25' E to 121°40' E; while for the second group, the events with M L ≤ 5.5 happened in the area from 23°50' N to 24°04' N and from 121°25' E to 121°42' E. The events of the second group were more or less located to the northeast of those of the first group. Of course, there was overlap of some events of the two groups. The migration of seismic activity from an area to the other demonstrates an interesting time history of readjustment of tectonic balance following the earthquakes. The phenomenon of earthquake migration can also be observed in the May 20, 1986 Hualien earthquake sequence Wang 1986, 1988). From the figure, we can see that most of the events are located around the Longitudinal Valley, some at the Costal Range, and some offshore. Figure 3 shows the depth distribution of number of events in a 5-km range along a selected longitude, because the error of focal depth is up to 5 km. A few events are located in the depth range 0-5 km. Most of the events are located in the depth range 5-10 km. Page 5 of 14 Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 This is similar to the conclusion reported by Wang et al. (1994) that inland earthquakes in Taiwan are located mainly in the depth range 0 − 12 km. The number of events remarkably decreases with increasing depth. The number of events with d > 25 km is only 3 that is very small. In the followings, only the events with d ≤ 25 km are taken into account. The crust-upper mantle boundary with v p = 7.5 km/s in the Taiwan region is mainly in the range 35 − 45 km as inferred from 3-D topographic inversed results by several authors (e.g., Rau and Wu 1995;Ma et al. 1996;Kim et al. 2005). Hence, an average depth of 40 km is taken as a boundary to classify the events: a crustal event with d ≤ 40 km and an uppermantle or subduction-zone event with d > 40 km. Hence, the events of the two earthquake swarms are mainly crustal earthquakes based on the average velocity structures of Taiwan. On the other hand, under the Costal Range in eastern Taiwan, Ma et al. (1996) suggested that the crust should be oceanic crustal due to the collision boundary between the Philippine Sea plate and the Eurasian plate (Tsai 1986;Tsai et al. 1977;Wu 1978) as mentioned above. Ma et al. (1996) also mentioned that the oceanic crustal thickness would be about 15 to 20 km inferred from the ascending of the 6.8 km/s contour around the longitudes of 121.25 − 121.50°E. The inversion results from the gravity survey by Yeh and Yen (1991) are based on the assumption that the oceanic crustal thickness near Taiwan is almost in the range of 10 − 20 km and that it increases from south to north. Hence, it is appropriate to take 20 km to be the lowerbound of such an oceanic crust. As mentioned above, the location error of focal depth can be up to 5 km for inland earthquakes and even higher for offshore events.

Depth distribution of earthquakes
Considering the location error, the selection of 25 km to be the lower-bound of focal depth seems acceptable. Figure 4a shows the time sequence of earthquake magnitudes for 165 M L ≥ 3 events with d ≤ 25 km. It is remarkable that the time difference between the two earthquake groups is shorter than 8 days. Although the largest event with M L = 6.2 occurred in the first group, there were five M L > 5 events in the second group. In addition, the average inter-event time interval is longer in the first group than in the second one. Figure 4b shows the time sequence of earthquake magnitudes for three M L ≥ 3 events with d > 25 km: two events with d = 30.2 km and d = 41.1 km, respectively, occurred during the time sequence of the first group and one event with 36.3 km occurred during the time sequence of the second group. The time intervals between any two events with d > 25 km are very long and these deeper events are not clearly associated with the shallow ones. Since the number of events with d > 25 km for each group is small (totally only 3 events), these deeper events are not further analyzed. Hence, totally 165 M L ≥ 3 events with d ≤ 25 km are taken into account in this study.

Methods
Since the methods applied in this study have been explained in details in several articles (e.g., Chen et al. 2015;Wang et al. 2017), only simple descriptions are given below. The detailed description about the technique can see Torrence and Compo (1998). In the followings, a time series is denoted by x i (i = 0, 1, 2, …, N-1).

Morlet wavelet analysis for dominant periods
The Morlet wavelet ψ(t) having a zero mean may be localized in both time and frequency space and it is composed  of a harmonic wave, with a constant k c subtracted from a plane wave, modulated by a Gaussian envelope (e.g., Farge 1992). The wavelet is: where The ω is the angular frequency (in Hz) and equals 2ℼ/T o where T o is the (characteristic) period of oscillations (Pyrak-Nolte and Nolte, 1995). In order to avoid some problems caused by small ω, ω is usually taken to be higher than 5 Hz, for instance, Farge (1992) took ω = 6 Hz that is related to T o = ~ 1 s. As ω > > 1 Hz, Eq. (1) becomes The continuous wavelet transform of x i is: where ψ*(t) is the complex conjugate of ψ(t), n is the localized time index, and s is the wavelet scale. A picture for different values of s and n can be plotted to show both the amplitude of any features versus scale (1) and how this amplitude varies with time. It is faster to calculate the wavelet transform by using Eq. (3) in the Fourier space. Since ψ(t) is usually complex, W n (s) is also complex and thus is the sum of its real part R[W n (s)] and its imaginary part. Hence, the amplitude, phase, and wavelet power spectrum are, respectively, |W n (s)|, θ = tan −1 {R[W n (s)]/I[W n (s)]}, and |W n (s)| 2 .
Let the discrete Fourier transform (FT) of x i be χ k : The FT of ψ(t/s) is F[ψ(sω)]. Based on the convolution theorem, W n (s) is: To ensure that the W n (s) at each s are directly comparable to each other and to the FTs of other time series, the FT of ψ(t/s) at each s is normalized to get unit energy, i.e., thus leading to ʃ|F[ψ o (ω)]| 2 dω = 1. Using Eq. (6) and referring to Eq. (5), the expectation value (EV) of |W n (s)| 2 (4)  Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 is equal to N times the EV of |χ k | 2 . For a white-noise, the EV is σ 2 /N, where σ 2 is the variance. For white-noises, the EV is |W n (s)| 2 = σ 2 for all n and s. Torrence and Compo (1998) defined the s as: s(j) = s o 2 jδj (j = 0, 1, 2, …, J) where s 0 is the smallest resolvable scale and J gives the largest one. In this study, s varies from s 0 to s 0 2 (J'δj) , where s 0 = 2δt with δt = 1 day and J' = 5/δj with δj = 0.1. This leads to 51 scales ranging from 2 to 64 days. The normalized Fourier power spectrum is N|χ k | 2 /2σ 2 . In order to explain the significance levels of the calculated values, the time series has a mean power spectrum. When a peak in the wavelet power spectrum is significantly above this background spectrum, it can be considered to be a true feature with a certain percentage of confidence. The 95% confidence level implies a test against a certain background level, and its interval refers to the range of confidence about a given value.
If x n is a normally distributed random variable, both the real and imaginary parts of χ k are normally distributed (Chatfield 1989). Since the square of a normally distributed variable is chi-square distributed with one degree of freedom (DOF), then |χ k | 2 is chi-square distributed with two DOFs, denoted by χ 2 2 (Jenkins and Watts 1968). To determine the 95% confidence level (significant at 5%), one multiplies the background spectrum by the 95th percentile value for χ 2 2 (Gilman et al. 1963). The 95% Fourier confidence spectrum will be displayed by a dashed curve in the following figures. Note that only at a few periods the power spectra will be above the 95% line. To meet the requirement of equal time interval, we consider the number of events occurring in a day and evaluate the dominant periods of time sequence by using the Morlet wavelet analysis.

Fluctuation analyses for memory effect
The auto-correlation of x i and x i+s of a time series x i that has been measured in an equal-time span will be made for different time lags or over different time scales s. In order to avoid a constant of set in the data, we take a new variable y i = x i -< x > in which < x > is the mean of x i . Hence, the auto-correlation function is C(s) = < y i y i+s > . where γ is a scaling exponent with 0 < γ< 1. Usually, it is not appropriate to directly calculate C(s) due to possible superposition of noises whose origins are unknown. Hence, we take an alternative way that is described below.
The power spectral density, S(f), of fractional noises is (Turcotte 1997): (7) C(s)µs −g The scaling exponent, β, is 0 for an uncorrelated white noise having a constant spectrum and 1 for a 1/f noise. For the signals with long-term correlations, we have 0 < β < 1.
When the long-term memory exists in an earthquake sequence, a larger-sized event is more likely to be followed by a larger-sized one and a smaller-sized event by a smaller-sized one. This clustering phenomenon leads to a 'mountain-valley' structure on every time scale (e.g., Bunde et al. 2005;Livina et al. 2005;Lennartz and Bunde 2009b). A 'mountain' that represents a group of a large number of events occurring in a short time interval is followed by a 'valley' that denotes a small number of events in the sequent time interval. In order to find the correct scaling law of fluctuations, non-stationarities that often exist in the data must be distinguished from the intrinsic fluctuations of the system. This task is not easy to perform due to numerous reasons. For example, the procedure to subtract some kind of moving average with a certain bin width would artificially induce the time scale into the data, thus destroying a possible scaling over a wider range of time scales. A convenient way to perform the task is the use of the fluctuation analyses (FA) technique. The fluctuation function F(s) is: Obviously, F(s) is associated with the auto-correlation function and power-spectral density. There are two steps to calculate the fluctuations. First, as mentioned above we calculate x i -< x > , where x i is the magnitude or inter-event time at natural i. Secondly, we sum up events within a window of length s, which divide the earthquake sequence into several segments. If the number of events is N, the number of segments is N s = N/s where N s is an integer. Since N is not necessary a multiple of s, the events in a small portion at the end of the sequence will be ignored. The squared fluctuation function is then the squared sum, averaged over all windows. Accordingly, F 2 (s) is, up to a factor s 2 , the variance of the mean values of s successive data points. For a long-term correlated time sequence, F(s) is with 1/2 ≤ α ≤ 1: α = 1/2 for white noise with β = 0 and α = 1 for 1/f noise with β = 1. As α < 1/2, the data are short-term corrected or uncorrelated; while as α > 1, the data are non-stationary, random-walk-like, and unbounded. The correlation, power-spectral density, and fluctuation function are related to each other by the following relationships (e.g., Kantelhardta et al. 2001;and (8) (10) F (s)µs α , Page 8 of 14 Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) Lennartz and Bunde 2009b): α = (2-γ)/2; α = (β + 1)/2; and γ = 1-β.

Morlet wavelet analysis for dominant periods
Results for the Morlet wavelet analysis are displayed in Fig. 5 for M L ≥ 3 events and Fig. 6 for M L ≥ 4 events. Each figure consists of three panels: (a) for the time sequence of number of monthly events; (b) for the wavelet power spectrum; and (c) for the global wavelet spectrum.
In panel (b), the logarithmic values of wavelet power spectrum for different periods (in month) at a certain time span are displayed by distinct kinds of color (from dark red to dark blue). The thick contour is the 95% confidence level, using a white-noise background spectrum. The black net region is the cone of influence, where zero padding has reduced the variance. The values of wavelet power spectrum inside the net have high uncertainties and thus cannot be taken into account. The period related to the local maximum is the local dominant period in a time span. In order to examine the dominant local maximum and related dominant period, it is The thick contour is the 95% confidence level, using a white-noise background spectrum. The black net is described in the text. The dashed line is the 95% confidence level for the global wavelet spectrum, using a white-noise background spectrum Page 9 of 14 Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 necessary to calculate the average of wavelet power spectra from panel (b) over time at a certain period. This average is named as the global wavelet spectrum. Results are demonstrated in panel (c), where the solid line represents the global wavelet spectrum and the dashed line denotes the 95% confidence level, using a white-noise background spectrum. The spectral intensity of white noise is uniform in the whole frequency range. In case the solid line does not intersect the latter, we cannot determine the period of an earthquake sequence. On the other hand, when the former intersects the latter at several periods, some segments of solid line are to the right of the dashed line, thus indicating that the global wavelet spectra in each segment between two intersection periods are significant. Hence, we may determine the periods of the earthquake sequence and the peak value of a period range for a segment is the dominant period of the period range.
The time sequences of number of daily events for M L ≥ 3 and M L ≥ 4 earthquakes are plotted in Figs. 5a and 6a, respectively. The two figures show that the number of events appeared in earlier April, then increased to the peak in middle July, and decreased from middle July The thick contour is the 95% confidence level, using a white-noise background spectrum. The black net is described in the text. The dashed line is the 95% confidence level for the global wavelet spectrum, using a white-noise background spectrum Page 10 of 14 Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 until the end of August. In addition, there were two small spikes in earlier June and earlier July. Figures 5b and 6b display the power spectra of the wavelet transform for the daily events for M L ≥ 3 and M L ≥ 4 earthquakes, respectively. The local maximums at several periods in different time spans can be seen. As shown in Figs. 5a and 6a, there was a large spike in middle July. The wavelet power spectra associated with this spike as shown in Figs. 5b and 6b are quite abnormal. In Fig. 5c, the solid line intersects the dashed line only during a few periods. Two relative maximum peaks above the 95% confidence level can be observed at 30.8 and 38.0 days in Fig. 5c. In between 30.8 and 38.0 days, the solid line is to the right of the dashed line, thus indicating that the global wavelet spectra in the two ranges are significant. There are several relative high peaks in Fig. 6c. However, the solid line to the left of the dashed line means that the global wavelet spectrum for each peak period is less significant. The dominant period is significant when its peak value is higher than the related 95% confidence level. The period associated with suck peak is taken to be the dominant period.

Fluctuation analyses for memory effect
In order to make fluctuation analyses, from Fig. 3a we plot the time sequences of earthquakes with d ≤ 25 km in the natural time domain for two cases, i.e., the earthquake magnitude and the inter-event time, in two magnitude ranges. Results are displayed in Fig. 7: (a) for magnitudes when M L ≥ 3; (b) for magnitudes when M L ≥ 4; (c) for inter-event times when M L ≥ 3; and (d) for inter-event times for M L ≥ 4 earthquakes. Clearly, for each case the plot for M L ≥ 3 earthquakes is different from that for M L ≥ 4 events. For the case of earthquake magnitude, the pattern is denser for M L ≥ 3 earthquakes than for M L ≥ 4 events. For the case of inter-event time, the vertical line segment is on the average higher for M L ≥ 4 earthquakes than for M L ≥ 3 events. The number of events and the maximum inter-event time for each magnitude range are given, respectively, in first parenthesis under the M-column and in that under the T-column of Table 1.
For the two magnitude ranges, the log-log plot of fluctuations F M (s) for magnitudes and F T (s) for interevent times versus time window, s, that varies from 0 to 146 days are, respectively, shown in Figs Chen et al. Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:24 (as listed in Table 1) Fig. 8, the plots for all cases slightly become flat or show roll-over at large log(s). (The "roll-over" means that the distribution of data points is changed into a shape of a part of a circle, in other word, the data points deflect from the regression line.) This might be attributed to the finite-size effect as suggested by Lennartz and Bunde (2009b). Figure 4a displays that for M L ≥ 3 earthquakes with d ≤ 25 km, in some time spans few line segments are close to one another. This suggests the possible existence of repeat of events with a short period, because the area where the Hualien swarm occurred is considered as a whole. Figure 4b displays that there are only three M L ≥ 3 earthquakes with d > 25 km, and the time interval between two events is longer than 50 days. The number of data of the time sequence is not enough to form a complete cycle to give a significant dominant period. Hence, it is difficult to evaluate the dominant period for deeper events due to a small number of data.

Morlet wavelet analysis for dominant periods
For M L ≥ 3 earthquakes, Fig. 5c shows that there are local peaks in the solid line at several periods. In between two dominant periods, i.e., 30.8 and 38.0 days, the solid line is to the right of the dashed line and there is a global peak of the solid line. The reasons to cause the dominant periods of 30 days are still unknown. It is interesting that such a value is similar to the period of lunar cycle around the Earth. Hence, one of the possible reasons is the effect caused by the lunar tide. Since we cannot find other swarms occurring the study area, it does not seem appropriate for us to confirm this reason at present.
For M L ≥ 4 earthquakes, Fig. 6c displays that the solid line is to the left of the dashed line and thus there is no peak above the associated 95% confidence level. This indicates that there is not any dominant period. Since the time period of the earthquake sequence is 146 days, the dominant period for M L ≥ 4 earthquake sequences might be longer than the total time period used in this study.   Table 1 Figure 7a, b show that the time variations of magnitudes in the natural time domain for two magnitude ranges are somewhat uniform. This indicates that abnormally large earthquake did not happen in the swarm during the time period. On the other hand, Fig. 7c, d show that the time variations of inter-event times in the natural time domain vary remarkably. Longer inter-event time appeared in the earlier and later time periods than in the middle time period. This indicates that the frequency of earthquake occurrences was lower in the former two time periods than in the latter one. This phenomenon can be seen in Fig. 4 as mentioned previously. In addition, the longest inter-event time is larger for M L ≥ 4 (~ 20. 375 days) than for M L ≥ 3 (~ 8.792 days). The two values are also listed Table 1.

Fluctuation analyses for memory effect
In Fig. 8, the data points are displayed by open circles for M L ≥ 3 and by crosses for M L ≥ 4. In Fig. 8a, the log-log plots of F M (s) versus s are well-distributed in two magnitude ranges and there is a significant reduction in log[F T (s)] at a particular value of log(s): log(s) = 1.4 (or s = 25) for M L ≥ 3 and log(s) = 1.14 (or s = 14) for M L ≥ 4. Such a particular value is smaller for M L ≥ 4 than for M L ≥ 3. This might be due to a fact that a smaller number of events for the former than for the latter. Figure 8b reveals that there is a significant reduction in log[F T (s)] at a particular value of log(s): log(s) = 1.1 (or s = 13) for M L ≥ 3 and log(s) = 0.7 (or s = 5.0) for M L ≥ 4. For the time sequences of both magnitudes and inter-event times, such a particular value is smaller for M L ≥ 4 than for M L ≥ 3. This might be due to a fact that segmentation of an earthquake sequence results in a decrease in the number of data points due to the exclusion of some data points in the later part of a sequence for practical calculations as mentioned above. Meanwhile, for those earthquake sequences as displayed in Figs. 7c, d longer inter-event times appear in the later part of a sequence. This can also result in a decrease in the calculated values of F T (s). In Figs. 8a and 8b, there is no data point at large log(s) > 1.5 when M L ≥ 4, because the number of events is small (see Table 1).
In Fig. 8a, a linear correlation exists between log[F M (s)] and log(s), i.e., log[F M (s)] = a M + α M log(s), for log(s) ≤ 1.40 or s ≤ 25 when M L ≥ 3 and for log(s) ≤ 1.14 or s ≤ 14 when M L ≥ 4. In Fig. 8b, a linear correlation exists between log[F T (s)] and log(s), i.e., log[F T (s)] = a T + α T log(s), for log(s) ≤ 1.1 or s ≤ 13 when M L ≥ 3 and for log(s) ≤ 0.7 or s ≤ 5 when M L ≥ 4. The ranges of log(s) for the existence of linear correlations are listed in Table 1. Hence, in the time sequences of magnitudes, an event can be influenced by several events that occurred before it. The number of such events is 25 when M L ≥ 3 and 14 when M L ≥ 4. While, in the time sequences of inter-event times, an event can be influenced by several events that occurred before it. The number of such events is 12 when M L ≥ 3 and 5 when M L ≥ 4. Clearly, the number is larger for the time sequence of magnitudes than for that of inter-event times. The linear regression equations inferred from the data points are displayed by thin solid lines in Fig. 8. The values of α M are 0.4784 for M L ≥ 3 and 0.4677 for M L ≥ 4 and those of α T are 0.4747 for M L ≥ 3 and 0.4849 for M L ≥ 4. These values of slope are given in Table 1. Clearly, the two values of α M and α T are smaller than but close to 0.5. This suggests that the memory effect in earthquake sequence of the swarm could be only short-term corrected or weakly corrected. For almost all cases, a linear correlation may be recognized for small s. Clearly, the memory effect is stronger for the time sequence of magnitudes than for that of inter-event times and higher for the M L ≥ 3 earthquakes than for the M L ≥ 4 ones because of the differences in the values of s.
The present observations are different to those obtained by Lennartz et al. (2008Lennartz et al. ( , 2011 for the earthquake sequences in northern and southern California, for which they assumed the existence of long-term memory effect. This might be due to a fact that many mainshock-aftershocks sequences with mainshocks of M > 6 were included in their study, while only the events of a swarm are included in the present data sets. The lasting time of aftershocks increases with the size of the mainshock. This is due to strong correlations between aftershocks and their mainshock and between an aftershock with others as described by Omori law (Omori 1896). Aftershocks are considered to be the result of stress alterations in the crust induced by mainshocks through timedependent processes, for example, the pore-fluid flow, viscous relaxation of the lower crust and upper mantle, and afterslip. Viscoelastic relaxation is a common mechanism for generating aftershocks (Scholz 1990;Chen et al. 2012). Therefore, the memory effect should be higher and longer for the mainshoch-aftershocks sequence with a larger mainshock than that with a smaller mainshock. This is the reason why Lennartz et al. (2011) studied the scaling law of F(s) versus s from the BASS model of aftershocks (Turcotte et al. 2007). Hence, it is not surprised that only the short-term memory effect was operative in the earthquake sequence of the Hualien swarm.

Conclusions
An earthquake sequence occurred in Hualien from April 7 to August 30, 2021. The Morlet wavelet technique is applied to analyze the dominant periods of temporal variations in numbers of daily events for the earthquake sequence in two magnitude ranges, i.e., M L ≥ 3 and M L ≥ 4. Results show that there are two dominant periods, i.e., 30.8 and 38.0 days, when M L ≥ 3. The dominant period cannot be observed for shallow M L ≥ 4 earthquakes. It might reflect the size effect, and thus cannot obtain the dominant period of the earthquake sequence. Hence, the dominant periods are 30.8 and 38.0 days only for the M L ≥ 3 earthquake sequence of the swarm.
The memory effect of earthquake sequences by applying the fluctuation analysis technique in the natural time domain. The earthquake sequences are represented in the temporal variations of magnitudes and inter-event times in the natural time domain. Two magnitude ranges, i.e., M L ≥ 3 and M L ≥ 4, are also taken into account. In the time sequences of magnitudes, an earthquake can be influenced by a few events occurred before it: 25 events when M L ≥ 3 and 14 events when M L ≥ 4. While, in the time sequences of inter-event times, an earthquake can be influenced by a few events occurred before it: 13 events when M L ≥ 3 and 5 events when M L ≥ 4. Clearly, the memory effect is stronger for the time sequence of magnitudes than for that of inter-event times and higher for M L ≥ 3 earthquakes than for M L ≥ 4 events. Consequently, the short-term corrected memory effect was operative in the earthquake sequence of the Hualien swarm. We assume that the results obtained in the present study will be the basic conditions for the construction of a physical model of a swarm.