Research ArticleGEOPHYSICS

Synchronized and asynchronous modulation of seismicity by hydrological loading: A case study in Taiwan

See allHide authors and affiliations

Science Advances  14 Apr 2021:
Vol. 7, no. 16, eabf7282
DOI: 10.1126/sciadv.abf7282


Delineation of physical factors that contribute to earthquake triggering is a challenging issue in seismology. We analyze hydrological modulation of seismicity in Taiwan using groundwater level data and GNSS time series. In western Taiwan, the seismicity rate reaches peak levels in February to April and drops to its lowest values in July to September, exhibiting a direct correlation with annual water unloading. The elastic hydrological load cycle may be the primary driving mechanism for the observed synchronized modulation of earthquakes, as also evidenced by deep earthquakes in eastern Taiwan. However, shallow earthquakes in eastern Taiwan (<18 km) are anticorrelated with water unloading, which is not well explained by either hydrological loading, fluid transport, or pore pressure changes and suggests other time-dependent processes. The moderate correlation between stacked monthly trends of large historic earthquakes and present-day seismicity implies a modestly higher seismic hazard during the time of low annual hydrological loading.


An unambiguous determination of the relationship between seismicity and hydrological loading cycles could provide valuable insights for improved regional hazard assessment. Previous studies have found the seasonal modulation of seismicity to be influenced by water forcing in a variety of tectonic settings, including volcanic areas, plate boundary zones, and plate interior regions that are subject to a wide range of interseismic tectonic strain rates (10−3 to 1 μstrain/year) and annual precipitation (0.3 to 3 m) (17). Often, seasonal modulation of seismicity is only resolvable when performing annual stacks of all events spanning a large number of years [e.g., (2)] or requires estimating Coulomb stress on representative fault planes and slip directions to document triggering [e.g., (4)]. However, statistical analyses are often hampered by the limited number of events associated with a low level of background seismicity, insufficient station coverage, or short catalog durations. To establish statistical significance, these analyses also require a careful consideration of catalog completeness and event clustering (810).

Taiwan is located on the convergent boundary where the Philippine Sea plate is colliding with the Eurasian plate at a rate of 85 to 90 mm/year (11, 12). This region has both frequent damaging earthquakes, with about one magnitude 6+ earthquake each year (13) and heavy seasonal precipitation of more than 2000 mm/year on average (14). The largest annual precipitation is close to 4 m, and the spatiotemporal distribution of rainfall is uneven (Fig. 1A). About 70% of the annual precipitation occurs from May to September and is primarily contributed from monsoons and typhoons. The seasonal water cycle causes annual fluctuations of 5 to 15 m in groundwater levels and 5 to 20 mm in Global Navigation Satellite System (GNSS) vertical displacement time series, reflecting the Earth’s elastic response to the seasonal water loading (14). Taking advantage of active seismicity (Fig. 1B) and high geodetic strain rates, as well as the densely distributed seismic, groundwater level, and geodetic networks in Taiwan, we quantify the spatiotemporal relationship between the hydrological cycle and earthquake seasonality by analyzing time series of seismicity, groundwater level, and the GNSS vertical time series.

Fig. 1 Rainfall, seismicity, depth distribution of earthquakes, and the magnitude of completeness (Mc) in 2002–2018.

(A) Difference of rainfall between the wet (month-month) and dry (month-month) seasons modified from figure 1 in (14). Red lines denote the boundaries of six physiographic regions in Taiwan: From west to east, these are the Coastal Plain (CP), the Western Foothills (WF), the Hsueshan Range (HR), the Central Range (CR), the Longitudinal Valley (LV), and the Coastal Range (CoR). (B) Seismicity with ML ≥ 2.5 and depth ≤ 30 km and large earthquakes with ML ≥ 6 (yellow stars). Black text shows the year, month, and magnitude of ML ≥ 6 events. Boxes show the selected regions for analyzing earthquake seasonality. (C to E) The depth distribution and Mc as a function of month for earthquakes located in WTW and ETW [locations shown in (B)], respectively. Black and colored lines in (C) and (D) show the depth distributions for earthquakes in the original CWB and various declustered catalogs (NNA, ETAS, and DL). Black arrows mark the focal depth above which 95% of all events in WTW and ETW occur.

We divide Taiwan into western and eastern subregions (Fig. 1B), in both of which GNSS-inferred dilatation rates reach up to 1 μstrain/year (15). Because of an insufficient number of events in the extensional stress regime (15) of the Central Range (CR in Fig. 1A), we exclude seismicity in this region. In western Taiwan (WTW), faults in a fold-and-thrust belt are mostly locked, and the deformation is taken up by crustal shortening above bedding-plane decollements (16). In eastern Taiwan (ETW), the Longitudinal Valley Fault is an east-dipping, plate boundary fault (17), and its southern portion is creeping at a rate of about 30 mm/year from the surface to a depth of about 3 km (18). In this study, we attempt to constrain the contribution of the hydrological cycle to earthquake modulation on the fault systems in WTW and ETW, respectively. The water forcing can involve an immediate effect of modulating seismic activity (synchronized mode) by static stress changes or exhibit a delay in earthquake triggering (asynchronous mode) through spatiotemporal variations of pore pressure, fluid transport, or other time-dependent processes. Time lags between the peak water loading and seismicity rate and the temporal characteristics of seismicity rates at different depth intervals provide critical constraints on triggering mechanisms and for evaluating the role of other processes in earthquake nucleation [e.g., (1921)].

We choose a constant magnitude threshold of ML ≧ 2.5 and a seismicity cutoff depth greater than the focal depths of 95% of all events in 2002–2018 (22). Detailed data processing procedures are described in Materials and Methods. The cutoff depths are about 28 km in WTW and 41 km in ETW, respectively (Fig. 1, C and D). We analyze seismicity using three different declustering methods: the epidemic type aftershock sequence (ETAS) model, the spatiotemporal double-link (DL) method (23), and the nearest neighbor algorithm (NNA) (24, 25) (fig. S1). Principal components analysis (PCA) is used to extract common hydrological signals in monthly time series of 40 groundwater level stations and 40 GNSS stations with continuous data in 2002–2018. The GNSS vertical position time series are converted to the effective water thickness (GNSS-EWT) using the method described by Hsu et al. (14). We next use empirical mode decomposition (EMD) to reconstruct groundwater level, GNSS vertical position, and seismicity rate time series into annual and interannual variations according to the frequency bandwidth of each intrinsic mode function (IMF). We specifically investigate earthquake occurrences with depth and time and explore possible physical mechanisms for synchronized and asynchronous modulation of earthquakes with seasonal and interannual hydrological load variations in Taiwan.


Seasonal variations in time series

The annual time series for the seismicity rate with ML ≧ 2.5 from the NNA declustered catalog, groundwater level, and GNSS annual vertical positions are shown in Fig. 2. Figure S2 shows the annual time series of seismicity rate from ETAS and DL declustering. The groundwater level and GNSS vertical motions show a clear annual periodicity, with low groundwater levels in winter during the peak GNSS annual uplift. The seismicity rate often reaches its peak in winter and spring around the time of low groundwater level and peak uplift. Peak earthquake occurrences occasionally appear in summer and autumn, for instance, in 2004 and 2009 (Fig. 2A). We find that the correlations between the seismicity rate and the other two time series are low (Fig. 2A) due to the inherent nonstationary and nonlinear nature of seismicity rate, thus the traditional cross-correlation analysis is not always the best tool to measure its similarity with the other datasets. Instead of showing a wiggle-to-wiggle match (i.e., a very high cross-correlation coefficient), the main purpose of Fig. 2A is to demonstrate the similar seasonal patterns (i.e., peak to peak and trough to trough) of seismicity rate, groundwater level, and GNSS vertical displacement. We plan to explore other methods to address the relations of nonstationary signals in the future.

Fig. 2 Annual time series of seismicity rate, groundwater level and GNSS vertical motion.

Time series of annual seismicity rate from NNA, groundwater level, and GNSS vertical motion (A) in WTW and (B) ETW. Red line indicates annual seismicity rate with ML ≧ 2.5 derived from the EMD. Blue and black lines represent the annual variations of groundwater level and GNSS vertical motions. Black vertical line indicates the time of ML ≧ 6 earthquakes labeled in Fig. 1B. Estimates of time lag and cross-correlation coefficient between data pairs are shown in the title of figure. The unit of time lag is in months. The time lag is the first data time series relative to the second data time series. Positive and negative lag values mean the first data time series is ahead of and behind the second data time series, respectively. The dashed lines in (B) show the predicted annual variations of groundwater level and GNSS vertical motions extrapolated using data from 2010–2018. Cross-correlation coefficients are estimated using data in 2010–2018.

In ETW, groundwater levels and GNSS vertical motions show annual periodicity with very similar temporal patterns in 2010–2018. Because of insufficient groundwater level and GNSS data before 2010, we use data in 2010–2018 to estimate the amplitude and phase of average annual motions and extrapolate the average seasonal pattern into 2002–2009. The temporal pattern of seismicity rate reveals a more complex variation in ETW compared to WTW. The peak seismicity rate (Fig. 2B and fig. S3) can occur both in winter and summer (e.g., 2002 and 2008–2012).

Annual stacked time series

To enhance the temporal characteristic of hydrological modulation, we produce annual stacked time series of seismicity rate, groundwater level, and GNSS-EWT in 2002–2018 (Fig. 3). Annual stacking is useful to resolve the average time of peak earthquake occurrence and estimate the associated time lags over the entire study period. We examine both magnitude and depth dependences of seasonality, as well as the periodicity of the subset of thrust-faulting events.

Fig. 3 Annual stacks of the EMD-extracted seismicity rate in 2002–2018 using different cutoff magnitudes.

(A) Results in WTW and (B) ETW. Numbers indicate the amplitude for normalization. Blue and black lines show the average of annual groundwater level and GNSS-EWT in the same time period. Seismic activity is high during lows of groundwater level in WTW and ETW.

We perform annual stacks for time series data extracted from the EMD (Fig 2 and figs. S2 and S3) and normalize stacking data to compare the different declustered seismicity catalogs with cutoff magnitudes ranging from 2.5 to 3.5 (Fig. 3). The annual stacks of raw monthly time series for seismicity rates with different cutoff magnitudes and depths (10 and 20 km) are shown in fig. S4. In WTW, the seismicity rate with ML ≧ 2.5 reaches peak levels in February to April, drops to a minimum in July to September, and gradually increases toward winter (Fig. 3A). The characteristics of earthquake periodicities with ML ≧ 2.5 and ML ≧ 3.0 are quite similar, with the lowest seismic activity in June to August (Fig. 3A and fig. S4). The seismicity rate is generally inversely correlated with the annual variations of groundwater level and GNSS-EWT. For ML ≧ 3.0 and ML ≧ 3.5 earthquakes, the seismic activity generally increases in the autumn, and the peak rate of ML ≧ 3.5 events occurs in November to December, about 3 months before the peak occurrence of ML ≧ 2.5 and ML ≧ 3.0 events. Note that the sample size for ML ≧ 3.5 is about five times and two times less than those for ML ≧ 2.5 and ML ≧ 3.0 events, respectively. We examine the annual stacks of raw monthly seismicity rates with different cutoff depths (10, 20, and 28 km in fig. S4) and find that the seismicity rates for ML ≧ 3.0 and ML ≧ 3.5 events are higher during the autumn compared to ML ≧ 2.5 earthquakes, in line with the EMD-extracted results (Fig. 3A). Moreover, the annual stacks of raw monthly seismicity rates in WTW show similar temporal features as those extracted by the EMD using different cutoff depths. We find peak GNSS-EWT occurs in August to September, about 2 months before the peak groundwater level in October (Fig. 3A and Table 1).

Table 1 Estimates of cross-correlation coefficient and time lag for annual stacked datasets.

We performed cross-correlation for each paired time series, including shallow/deep seismicity rate (EQ), groundwater level (GW), and GNSS-EWT. The depth thresholds for WTW and ETW are 11 and 18 km, respectively, which divide the respective catalogs into similarly sized numbers of events. The numbers on the left and right hand side of the slash are cross-correlation coefficient (cc) and time lag (t) with units of months, respectively. The time lag is the first data time series relative to the second data time series. Positive and negative lag values mean that the first data time series is ahead of and behind the second data time series, respectively.

View this table:

Figure 3B shows the normalized annual stacked data (2002–2018) in ETW with cutoff magnitudes ranging from ML 2.5 to 3.5. The annual stacks of raw monthly time series for seismicity rates with different cutoff depths (10 and 20 km) are shown in fig. S4. All declustered catalogs in ETW share some common features, with more frequent earthquake occurrences in April to June and a low seismicity rate in November to January. In general, the level of seismic activity is anticorrelated with the groundwater level. Both the peaks and lows in seismicity rate in ETW appear to have a delay of 2 to 3 months compared to those in WTW, which may reflect different triggering mechanisms. In ETW, the GNSS-EWT reaches its peak in July to August, 2 to 3 months before the peak groundwater level in October (Fig. 3B). The temporal characteristics of seismicity binned with different cutoff depths and magnitudes show more variability in ETW. The annual stack of ML ≧ 2.5 earthquakes is less dependent on the cutoff depths of seismicity (fig. S4), presumably due to a larger number of events compared to the ML ≧ 3.0 and ML ≧ 3.5 earthquakes.

To investigate the influence of hydrological modulation at different depths, we separate the declustered ML ≧ 2.5 seismicity catalogs into two groups (shallow and deep). The separation depths for WTW and ETW are about 11 and 18 km, respectively. We also examine the depth dependence using the annual stack of raw monthly seismicity rates and subsets of data, considering a fixed cutoff depth of 20 km (fig. S5). In WTW, the seismicity rate at shallow depths peaks in February to March, gradually decreases to its lowest level in August to September, and then rapidly increases in the winter and early spring (Fig. 4A). The annual stack of shallow seismicity rate is inversely correlated with the GNSS-EWT, which reaches its minimum in February to March and maximum in August to October. The cross-correlation coefficients (cc) between the shallow seismicity rate and GNSS-EWT range from −0.90 to −0.96, with a time shift of no more than 1 month (Table 1). The deep seismicity rate attains its highest value in March to April, drops to a minimum in August to September (Fig. 4A and fig. S5), and gradually increases toward the end of the year and into the next spring. This temporal pattern matches the variations of groundwater level and GNSS-EWT. The peak seismicity rate of deep earthquakes is delayed by about 1 month behind the annual peak of shallow earthquakes, with a cc of 0.78 to 0.93 (Table 1). Annual stacks of shallow and deep events in WTW extracted from the EMD are consistent with the annual stack of raw monthly seismicity rate, as well as with the temporal variability derived from a subset of seismicity when a cutoff depth of 20 km is used (fig. S5).

Fig. 4 Annual stacks of the EMD-extracted shallow and deep seismicity rate (M ≥ 2.5) and thrust events in 2002–2018.

(A) Data in WTW and (B) ETW. The depth thresholds for the separation of shallow and deep earthquakes are 11 km for WTW and 18 km for ETW, respectively. Numbers indicate the amplitude for normalization. Blue and black lines show the monthly average of annual groundwater level and GNSS-EWT over the same time period.

In ETW, the seismicity rates for shallow and deep earthquakes exhibit opposite temporal patterns, with more frequent shallow earthquakes in May to August and more frequent deep earthquakes in December to February (Fig. 4B). There, the temporal distribution of shallow earthquakes is generally positively correlated with the annual variation of groundwater level and GNSS-EWT (Table 1). In contrast, the temporal signature of deep earthquakes extracted from the different declustered catalogs shows a more coherent inverse correlation with the annual variation of GNSS-EWT. Additional annual stack results from raw monthly seismicity rate and the subset of seismicity with the fixed cutoff depth of 20 km are shown in fig. S5. We find an overall consistency between different approaches and cutoff depths (Fig. 4B and fig. S5). except for the seismicity rate derived from ETAS, which has a very limited number of events (~10) and shows a large variability (fig. S5).

To this point, we have included all seismic events, without considering different types of faulting. This is mainly because (i) it is difficult to resolve earthquake focal mechanisms of small-magnitude earthquakes and (ii) it is necessary to have a sufficient number of events to obtain statistical significance. If hydrological elastic loading plays a primary role in earthquake modulation, then we expect the seismicity rate of thrust-faulting events to be inversely correlated with the seasonal water cycle, as an increased overlying water load increases the normal stress and decreases the up-dip shear stress. Here, we choose only thrust-faulting events to perform the annual stacking. In WTW, the activity of thrust earthquakes is low in July to September and high in January to March (Fig. 4A), similar to the pattern of the overall seismicity and in line with elastic loading and unloading effects due to the hydrological cycle. In ETW, the annual stacks of ETAS and NNA declustered thrust events show a low seismicity rate in May to June and more frequent earthquakes in September to December, whereas the seismicity rate from DL exhibits an inverse temporal pattern (Fig. 4B). This discrepancy may suggest that the nucleation of shallow and deep earthquakes in ETW is driven by different physical mechanisms. We also examine the annual stack of thrust events using the subset of seismicity with a cutoff depth of 20 km (fig. S5), which results in a similar temporal distribution as shown in Fig. 4.


Physical mechanisms for earthquake modulation

In WTW, shallow earthquakes occur more frequently in late winter and early spring, coincident with peak crustal rebound in response to decreasing water loads indicated by a trough in the GNSS-EWT (Fig. 4A and fig. S5). The temporal distribution from the subset of thrust-faulting events also follows this pattern, suggesting hydrological modulation of seismicity as shown by the strong correlation among various datasets (Table 1). Deep earthquakes show peak occurrence in March to April, between the time periods with the lowest GNSS-EWT and groundwater level (Table 1). The GNSS-EWT may better reflect the temporal variation of total water storage, because peak groundwater storage likely lags behind the peak land water storage by about 2 months (14).

To better resolve earthquake-triggering mechanisms, we further examine the EMD-derived interannual time series of rainfall, groundwater level, GNSS-EWT, and declustered seismicity rate at all depths in WTW (Fig. 5A). For comparison, the interannual time series obtained using a zero-phase–shifted band-pass filter with periods of 2 to 7 years is also shown in Fig. 5B. If hydrological loading plays a dominant role, then we expect correlation between longer-term lows in water storage and peaks in seismicity rate. Notable time lags between minimum water storage and peak seismicity rate at multiannual time scales imply that hydrological loading may not be a resolvable driver of seismicity rate variations at this time scale (text S1 and table S1). An alternative possibility for a large time shift may result from mixing between annual and interannual signals due to imperfect separation. We find an increase of seismicity rate in 2010–2012, corresponding to a time period with relatively low water storage. A notable increase of shallow earthquakes (<15 km) for >2.5 years after the August 2009 typhoon Morakot was also reported by Steer et al. (26). They attribute the change of seismicity rate to crustal unloading due to rapid sediment export after typhoon Morakot (26), although the spatial extent of substantial seismicity rate changes extends well beyond the landslide zone. We speculate that, despite the typhoon-associated rainfall, the increase of shallow events lasting ~2.5 years after typhoon Morakot corresponds to relatively low water storage in Taiwan as indicated in Figure 5.

Fig. 5 Multiannual time series of rainfall, groundwater level, GNSS-EWT, and declustered seismicity rate.

(A) All parameters are extracted from the EMD and (B) with a zero-phase–shifted band-pass filter (2 to 7 years). The curves are offset for clarity. Blue shaded bands denote relatively wet periods. We use the average annual rainfall in the hydrological year, which is defined as the 12-month period from May 1st to April 30th of the following year. Wet periods are defined as the years with the average annual rainfall above one SD of the mean value in 2002–2018. Blue vertical lines indicate the time of typhoon Morakot.

In ETW, shallow and deep earthquakes show almost opposite temporal distributions (Fig. 4B). Frequent shallow events in May to August may suggest that fault zones are highly fractured, allowing fast fluid infiltration, as demonstrated by a high fault creep rate on the southern segment of the Longitudinal Valley Fault in the wet season from May to October and peak creep rate roughly between July and October (27, 28). Chang et al. (28) modeled the temporal variation of shallow fault creep (depth < 5 km) using a velocity-strengthening friction and one-dimensional (1D) groundwater diffusion model and reported a large hydraulic diffusivity of 0.05 to 0.8 m2/s in the fault zone (28). Highly permeable fracture zones can provide pathways to facilitate fluid transport and could promote triggering during the wet summer months (20, 29). However, most of the hypocentral depths of the shallow events range from 8 to 15 km in ETW (Fig. 1D), and it is unclear whether fluids can infiltrate to such great depths. A very high hydraulic diffusivity would be needed to produce a significant pore pressure change at depth within a short time (30). The insignificant time lag between the peak GNSS-EWT and the frequent occurrence of shallow earthquakes imply that pore pressure diffusion may not be responsible for triggering of the shallow earthquakes.

Stress changes due to poroelastic deformation (31) are generally small and rapidly decay with distance from the injection region or a fluid-pressurized fault zone (31, 32). Recent laboratory experiments and in situ observations have shown that aseismic slip in pressured regions of a fault can activate slip in nonpressurized fault zones (3335). Hydraulic fracturing-induced earthquakes distant from injection wells may be induced by pore pressure–driven aseismic slip propagating from reservoir depth to the seismogenic zone (36). We expect that the poroelastic stress variation due to pore pressure changes from precipitation and groundwater level changes will differ substantially from that caused by a localized fluid injection. A direct comparison of these two scenarios is not easy and is beyond the scope of this study. To provide a quantitative measure of the seasonal stress cycles, we use a 1D spherical layered Earth model and estimated Coulomb stress changes on a north-south trended, 30° east-dipping receiver fault geometry with a rake of 90° due to seasonal hydrological surface loading constrained by the GNSS-inferred effective water thickness change (14). The mean annual water thickness change across Taiwan is 0.53 ± 0.17 m, and the largest seasonal change of up to 0.91 m is found in southwest Taiwan (14). We estimate Coulomb stress changes ranging from 3 to 5 kPa (fig. S6), given a value of 0.4 for the apparent coefficient of friction, 40 GPa for the rigidity, and 0.25 for the Poisson’s ratio. Our estimates are comparable to seasonal stress change amplitudes of 2 and 10 kPa from hydrological loading found in California and southern Alaska, respectively [e.g., (4, 5)]. The Coulomb stress changes due to elastic loading are also similar to stress changes of 10 to 100 kPa at a few kilometer depths resulting from poroelastic effects of shallow hydraulic fracturing injections [e.g., (37)] and comparable to stress changes of a few kilopascal from pore pressure cycles associated with precipitation and groundwater variations, which feature phase lags that strongly depend on the hydraulic diffusivity of the upper crust [e.g., as illustrated for Parkfield, California, in figure S7 in (38)].

On the southern Longitudinal Valley Fault, the surface creep recorded by creepmeters shows clear seasonal fluctuations (27, 28), with 1-year periodicity inferred from slip rates of repeating earthquakes (39). However, time series of shallow and deep creep rates are not in phase nor do they have a constant time lag, suggesting that a single mechanism is insufficient to interpret these observations. We speculate that multiple mechanisms may act at the same time, leading to the complex behavior. For deep events in ETW, however, the temporal distributions of seismicity rate and GNSS-EWT show a strong inverse correlation, suggesting control of earthquake nucleation by elastic hydrological unloading as also observed in WTW.

Other periodic stress variations may modulate earthquakes as well. In Taiwan, effects from atmospheric, nontidal ocean, or global hydrological loads are likely to be modest, given the small amplitude (<3 mm) of seasonal vertical displacements from these sources, in contrast to 10 to 20 mm resulting from the water cycle (14). Thermoelastic strain due to seasonal surface temperature changes is another candidate mechanism (40, 41). However, Taiwan features a subtropical climate that does not show significant seasonal surface temperature variation. The amplitude of predicted annual vertical thermoelastic displacements is only ~0.1 to 0.2 mm, implying a very weak effect on seasonality (14).

Impact on seismic hazard

Our examination of earthquake modulation in Taiwan indicates a pronounced correlation between annual seismicity rate and the hydrological water cycle. The strong correlation implies that, at any given time, a large population of faults is critically stressed such that the small stress perturbations from the water cycle can lead to triggering of small events. In this way, the seasonal hydrological cycle may modulate the release of seismic energy for small-magnitude earthquakes. However, a higher background seismicity rate in winter may also increase the probability of higher-magnitude ruptures on large fault systems (3, 4).

We assess the seasonal preference of 63 Mw ≥ 6 events from 1604 to 2018 (42) and exclude large aftershocks of Mw ≥ 6 earthquakes by counting only the first event if multiple Mw ≥ 6 earthquakes occurred within a year. This method was also used by Heki (3) for earthquakes in Japan. In WTW, the occurrence of historic earthquakes peaks in March (Fig. 6A), in accordance with the high seismicity rate in spring shown in the modern records (Figs. 3 and 4 and fig. S4). Historic events in ETW show more frequent occurrences of earthquakes in January, April to May, and December (Fig. 6B), but there are only a total of 10 ETW events in our dataset. These results are similar to present-day background seismicity rate (Figs. 3 and 4 and fig. S4) despite the limited number of historic events. Since 1604, two Mw ≥ 7 earthquakes in WTW occurred in April and December, at the time of frequent earthquake occurrences in the modern-day record, but we note the exception of the 1999 Chi-Chi Mw 7.6 earthquake in September. The influence of seasonal loading on Mw ≥ 7 earthquakes remains obscure because there are too few large earthquakes to assess statistical significance. Nevertheless, by considering events with Mw ≥ 6, the overall temporal correlation between the peak occurrences of historic and present-day earthquakes imply an enhanced seismic hazard during times of low annual water storage.

Fig. 6 Annual stacks of Mw ≧ 6 to Mw ≧ 7 events in 1604–2018.

(A) Data from WTW and (B) ETW. The cutoff depths for seismicity are 28 km for WTW and 41 km for ETW. Connected dotted lines denote the number of events in each month. Solid and dashed horizontal lines represent the mean value of events per year and one SD above the mean value.

Seasonality of earthquake occurrence worldwide

Seasonal modulation of seismicity has been observed in various tectonic settings. Our observations for both shallow and deep events in WTW and deep earthquakes in ETW show an insignificant time lag (<1 month) between the lowest annual water storage in spring and annual peaks of ML ≧ 2.5 and ML ≧ 3.5 earthquakes. This synchronized temporal pattern is consistent with frequent occurrences of large intraplate earthquakes (M ≥ 7) across Japan (3) and shallow earthquakes (depth ≤ 20 km) in southwest Japan (21) during crustal unloading in the snow melting season. Moreover, studies of background seismicity rates in the Himalayas of Nepal (43), the New Madrid Seismic zone (2), and the western branch of the East African Rift System (6) also observed in-phase modulation of seismicity rates with the hydrological cycle. Ueda and Kato (21) found another peak of seismicity in autumn, which they attributed to an increase of pore pressure in fault zones due to heavy precipitation in summer. The two peaks of earthquake numbers in spring and autumn in southwest Japan (21) are similar to findings in the locked section of the San Andreas Fault (SAF) near Parkfield (19). However, the seismicity rate along the creeping section of the SAF shows a relatively broad peak extending from August to January, roughly following the rainy season (19).

In our study, the shallow events in ETW, about 56% of which are contributed by the southern creeping section of the Longitudinal Valley Fault, reach an annual peak during the summer rainy season. We find a 3- to 5-month time lag between the lowest annual water storage and peak occurrence of shallow earthquakes in ETW (Table 1). Therefore, hydrological loading is unlikely to be the primary driving mechanism for the observed asynchronous modulation of earthquakes. Pore pressure change is another candidate mechanism, proposed by previous studies to explain the 3-month time lag between the increase of annual seismicity rate and failure-encouraging annual stresses in southern Alaska (5). Strong modulation by pore pressure increases has also been invoked for shallow events (<4 km) in highly permeable host rock (20, 29). However, it is difficult to invoke pore pressure diffusion from near-surface groundwater to explain events with hypocentral depths well below 8 km in ETW (Fig. 1D). We speculate that an increase in small earthquake activity at shallower depth may be able to trigger subsequent ruptures along faults that are critically stressed and/or that aseismic slip may propagate from the pressured regions of fault zones to unstable stick-slip regions, as discussed previously. The bimodal distribution in seismicity rate observed on the SAF and in southwest Japan may be a superposition of two distinct peaks of shallow and deep earthquakes that correspond to different triggering mechanisms, similar to the modulation patterns in ETW.

Synchronized modulation of seismicity by water loading is observed in Taiwan. Frequent earthquake occurrences in late winter and early spring in WTW are coincident with the timing of the lowest groundwater and GNSS-EWT. In ETW, seasonal modulation of deep earthquakes is synchronous with annual water storage minima, but seasonal modulation of shallow earthquakes is not. Triggering of shallow earthquakes in ETW is difficult to explain by hydrological loading, fluid transport, and/or pore pressure changes. Alternative physical mechanisms may invoke subsequent triggering of secondary ruptures and aseismic slip propagating from the pressured region of fault zones to unstable stick-slip regions. Investigating depth dependence of earthquake modulation in various tectonic settings is important to evaluate the relative contributions of different underlying physical mechanisms. In WTW, the observed modulations of historic and present-day seismicity rates with annual hydrological cycles are similar, implying a modestly higher seismic hazard during the time of low annual hydrological loading. Our findings of both synchronized and asynchronous modulation between seismicity and hydrological loading cycles in Taiwan further contribute to the delineation of controlling factors of seismogenesis that, in turn, provides valuable insights into regional hazard assessment.


Seismicity catalog

We used seismicity between 2002 and 2018 from the Central Weather Bureau (CWB) in Taiwan (Fig. 1B) to minimize the influence of the 1999 Mw 7.6 Chi-Chi earthquake and to maximize the overall network stability, including GNSS and groundwater stations. Offshore earthquake epicenters were excluded from the analysis. We divided Taiwan into the western and eastern regions (WTW and ETW; Fig. 1B) according to the tectonic setting and fault slip characteristics. WTW covers the Coastal Plain, the Western Foothills, and the Hsueshan Range, whereas ETW consists of the Longitudinal Valley and the Coastal Range (Fig. 1A).

Estimates of the minimum magnitude of completeness (Mc) are crucial to quantifying the spatial and temporal evolutions of seismicity rate. Because of seasonal noise and the variability of the seismic network, the earthquake detection threshold varies with time and space. To ensure the statistical significance of the temporal distribution of earthquakes, we used the method proposed by Wiemer and Wyss (44) to estimate the Mc in WTW and ETW in each month of year (Fig. 1E). Mc varies from 1.8 to 2.0 in WTW and from 1.8 to 2.2 in ETW within a year. The multiyear variations of Mc in WTW and ETW were shown in fig. S7A. Mc is generally lower than 2 in WTW and lower than 2.2 in ETW (except for 2011). A significant drop of Mc after 2012 is associated with the upgrade of the CWB seismic network to 24-bit data recording and 100-Hz sampling rate (45). To examine the depth dependence of Mc, we also estimated the Mc using events in 2002–2018 in two depth bins: 0 to 10 km and 10 to 20 km in WTW and ETW, respectively. The values of Mc are 1.9 and 1.8 in WTW and 2.1 and 2 in ETW for events at shallow and deep depths, respectively. The Mc for shallow earthquakes (0 to 10 km) is slightly higher than that for deep earthquakes (10 to 20 km). This small difference does not warrant further interpretation but may be related to decreasing seismic attenuation at greater depth. Moreover, the temporal variations of Mc for shallow and deep earthquakes do not show a systematic difference (fig. S7B). We chose a constant minimum magnitude threshold of 2.5 for the entire study period to investigate earthquake seasonality in WTW and ETW. The seismicity cutoff depth was determined according to the focal depths of 95% of all events (22), and the threshold depths were about 28 km in WTW and 41 km in ETW, respectively (Fig. 1, C and D). The final catalog contains 5380 events in WTW and 3260 events in ETW.

Declustered seismicity

To remove events densely clustered in time and space, we analyzed seismicity using three different declustering methods, including the ETAS model, the spatiotemporal DL method (23), and the NNA (24, 25) (fig. S1). In the ETAS model, the earthquake sequence is represented by a point process, and probabilities are given for each earthquake being triggered by the preceding earthquake or by a background process (46, 47). We classified the background seismicity as the probability of being a background event greater than 50%. The time and spatial DL method aims to identify aftershocks using a space-time distance. Earthquakes closely related in space and time are considered aftershocks. Here, we removed aftershock sequences for the mainshocks with ML > 3.5 using 3 days and 5 km as the linking parameters. The NNA defines distances ηij between events in a space-time-magnitude domain asηij={tij(rij)d10bmi,tij>0,tij0(1)where tij = tjti and rij are the temporal and spatial distances between an event j and an earlier event i, mi is the magnitude of the event i, and b is the b value of the Gutenberg-Richter relation set to 1 in our case. d is the fractal dimension of the hypocenter distribution and is set to 2.5 (48). The calculated nearest neighbor distances ηij generally exhibit a bimodal distribution (24, 25), where the short-distance population defines the clustered events that are removed, and the long-distance population represents the independent background events that are used for the analysis.

After performing ETAS declustering, there were 3286 and 1994 events left in WTW and ETW, respectively, both account for about 61% of events in the nondeclustered catalog. The DL declustered catalog contains 4656 and 2443 events in WTW and ETW, respectively, about 87 and 75% of the original dataset. The NNA declustered catalog contains 4253 and 2049 events, respectively, in WTW and ETW, about 79 and 63% of the original CWB catalog (fig. S1).

Declustered focal mechanism catalog

Hydrological loading may favor triggering of events with different types of focal mechanisms in various seasons, as associated shear and normal stress changes depend strongly on the receiver fault geometry. We focus on investigating the seasonality of thrust-faulting events, because Taiwan is located at a convergent plate boundary with abundant thrust faulting earthquakes that can be used to test the influence of loading and unloading on earthquake modulation. The earthquake focal mechanism catalog based on first-motion polarities of P waves is available from previous studies (4951). Events with the desired focal mechanisms were selected from the declustered catalog for further analyses. For earthquakes with double-couple focal mechanisms in WTW (4951), about 45% are thrust, 3% normal, 43% strike-slip-faulting, and 9% oblique-slip mechanisms (fig. S7C), according to the classification method of Frohlich (52). In ETW, about 67% of events have thrust, 7% normal, 15% strike-slip-faulting, and 11% oblique-slip mechanisms (fig. S7C). The declustered focal mechanism catalogs following ETAS, DL, and NNA declustering contain 358 to 475 thrust events in WTW and 179 to 242 thrust-faulting earthquakes in ETW.

Groundwater level and GNSS data

The groundwater data consist of 1-hourly measurements from the shallowest wells in the groundwater level monitoring network operated by the Water Resources Agency, Ministry of Economic Affair (5355). In this study, we chose 40 stations with continuous data in 2002–2018 (Fig. 1A). Moreover, continuous GNSS data in the Taiwan GNSS array were used to study seasonal vertical deformation. We removed vertical motions associated with the atmospheric, nontidal ocean, and global hydrological loads from the GNSS time series to focus on studying the water storage variation in Taiwan. A least-squares regression was used to remove the long-term linear trends, earthquake-related signals, and instrumental offsets. We used the residual GNSS vertical time series to characterize the water cycle and estimate the temporal variation of EWT using Green’s functions (56) that relate the elastic vertical deformation to a point load based on the Preliminary Reference Earth Model (57). The detailed processing procedures for groundwater data, GNSS vertical time series, and GNSS-inferred EWT were described by Hsu et al. (14). We used 40 available GNSS stations in 2002–2018 (Fig. 1A) in this study.

Principal components analysis

The PCA is an effective tool to extract common signals among different stations characterized by a spatially correlated source. Applying PCA requires data to be continuous, and a linear interpolation was used to replace missing data values in each GNSS station’s time series. The vertical time series with 1-day sampling rate were assembled in matrix G with dimension of m × n where m is the number of stations and n is the number of data points. We decomposed G into three matrices using the singular value decomposition method asGm×n=Um×mSm×nVn×nT(2)

Here, U and V are matrices of eigenvectors for spatial and temporal variations, respectively. The diagonal elements of the diagonal matrix S represent the contribution of each spatial function in U and time function in V (58). The number of diagonal elements chosen in S, that is the number of selected principal components, determines how well the sum of the reconstructed time series can fit the original observations. The contribution of each element in S is arranged in descending order with the smaller contributions listed in the end. We determined the number of selected components for which the reconstructed time series can explain ≥95% of the variance of the data. We next used the empirical model decomposition to decompose the reconstructed time series into annual and interannual variations according to the frequency band width of each IMF.

Empirical mode decomposition

Most of the traditional data processing methods are developed under rigorous mathematical rules and may not be applicable to nonlinear and nonstationary processes in nature. Huang et al. (59) proposed the Hilbert-Huang transform consisting of EMD and Hilbert spectral analysis to better understand nonlinearity and nonstationary characteristics hidden in data. To separate signals with different time scales, we used the EMD to decompose time series of earthquakes into several IMFs through a shifting process. The goal of shifting is to eliminate riding waves on signals and to make signals more symmetric (59, 60). An IMF must satisfy two conditions: (i) in the whole dataset, the number of extrema and the number of zero-crossings must be equal or differ by one at most; (ii) at any point, the mean value of the upper and lower envelopes defined using local maxima and minima must be zero. The EMD process is terminated when there are no more local maxima and minima to generate envelops. The decomposition results in a number of IMFs with different frequencies from high to low. Figure S8 shows the decomposition results of declustered seismicity in WTW and ETW, respectively. The EMD-extracted annual groundwater level and GNSS vertical motions using PCA-derived time series are shown in fig. S9. We chose IMFs that contain annual and interannual signals to study seasonal modulation of earthquakes and long-term trends in seismicity rates, respectively. In addition, we compared annual seismicity rates extracted from the EMD and with a zero-phase–shifted band-pass filter (0.7 to 1.3 year) in WTW and ETW (fig. S10). The seasonal variations extracted by the EMD better represent annual signals compared to the filtered time series.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J.-P. Avouac and Z. Peng for helpful discussion of earthquake seasonality. We also thank the Central Weather Bureau (CWB) for provision of seismicity catalogs; CWB, the Central Geological Survey, and the Ministry of the Interior for continuous GNSS data; and the Water Resource agency, Taiwan for groundwater data. We are grateful to C.-H. Chen for the discussion of EMD and to S.-B. Yu, H.-H. Su, and Y.-C. Tsai for collecting and processing continuous GNSS data. We thank the associate editor, K. Furlong, and two reviewers for detailed and constructive comments. Funding: This research is supported by the Academia Sinica, AS-CDA-105-M05, and the Ministry of Science and Technology grant MOST 108-2116-M-001-021 –MY3. R.B. acknowledges support by the NASA Earth Surface and Interior program. Author contributions: Y.-J.H. designed the study, analyzed the data, and wrote the manuscript with contributions from all other coauthors. Y.-T.L., H.-H.H., Y.-F.H., and J.Z. performed earthquake declustering. All authors contributed to the discussion and interpretation of the results. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. This is IESAS contribution number 2389 and NRCan contribution number 20200734.

Stay Connected to Science Advances

Navigate This Article