Seismic waves from large earthquakes induce ephemeral dynamic stress perturbations over large distances (Huang et al., 2004). Despite leaving no permanent stress changes once they pass, seismic waves are capable of triggering earthquakes at distances far beyondthe aftershock region (Hill, 2008; Brodsky and van der Elst, 2014; Miyazawa et al., 2021). Remote dynamic triggering is a rare ( < 2% of the time) (Pankow and Kilb, 2020) but ubiquitous phenomena (Velasco et al., 2008), and the number of evidence keeps growing since the first observations of entensive earthquake triggering in the United States following the 1992 moment magnitude (M _w ) 7.3 Landers earthquake (Hill et al., 1993) and the 2002 M _w 7.9 Denali Fault earthquake (Gomberg et al., 2004). In general, dynamically triggered earthquakes occur subsequent to the passage of the surface waves (Hill, 2008; Velasco et al., 2008) with periods of 15–20 s, since they dominate gound motion in the far field (waves with smaller periods tend to be scattered and attenuated) (van der Elst and Brodsky, 2010). However, the optimal conditions for triggering remain unknown and the latter represents likely a complex combination of wave characteristics (e.g., peak amplitude, duration) and local environment (fault types and geometry and wave incidence angle) (Parsons et al., 2014). If long-lasting (hundred of seconds) dynamic perturbations at moderate to high strain level ( > 100 nϵ or equivalently ∼ 3 kPa) may be an efficient combination for triggering (Pollitz et al., 2012; Johnson and Bürgmann, 2016), peak strain level of teleseismic waves represents a key contribution to long-range triggering (van der Elst and Brodsky, 2010) and critically stressed faults can also rupture under strain level that is near the background noise (∼ 2 nϵ) (van der Elst and Brodsky, 2010; Miyazawa et al., 2021) [Supplementary Figure S1 in the Electronic Supplement shows a case of remote triggering at very low regional strain in Central Greece (Bernard et al., 2006; Canitano et al., 2013)].

Establishing a temporal relationship between two events is simple because seismicity is produced continuously in active tectonic regions. More challenging is to establish temporal and spatial causation between two distant earthquakes. The longer the delay between the mainshock and the posited event, the more difficult to connect them. This delay, which represents a fundamental characteristics of remote triggering, may be difficult to explain by a simple Coulomb failure model only (e.g., Belardinelli et al., 2003). Conversely, a small, or the absence of delay between two events does not necessarily imply remote triggering, as still remains the possibility of coincidental occurrence of events. Indeed, Parsons et al. (2014) found that as many as five earthquakes with M _w ≥ 6 can occur purely by chance on any given day, for example. Since most of observed remote triggering cases concern microearthquakes and tectonic tremors (Peng et al., 2009; Peng et al., 2011; Gonzalez-Huizar et al., 2012), the magnitude range of the triggered earthquakes is still subject to debate. Many studies report the very rare occurrence of moderate to large earthquakes (M _w ≥ 5) remotely triggered hours to days following large teleseisms (M _w ≥ 7) (Parsons and Velasco, 2011; Johnson et al., 2015) while O’Malley et al. (2018) found that higher magnitude earthquakes occur more often than smaller events within 3 days following the passage of seismic waves. In Taiwan, detection of triggered microseismicity and tectonic tremors has also been reported (Chao et al., 2011; Sun et al., 2015), but whether larger events were triggered has not been demonstrated yet. Remote triggering is generally established using statistically significant variations in earthquake rate coincident with the passage of seismic waves (e.g., Johnson and Bürgmann, 2016; Yao et al., 2021). However, going beyond rate changes as the sole indicator for remote triggering (Pankow and Kilb, 2020) is essential for capturing and quantifying the timing of the triggered earthquakes and understanding the physical processes responsible for this delay whose remain major challenges in earthquake forecasting (Brodsky and van der Elst, 2014).

In this study, we propose a physical framework for the integration of borehole strainmeter time-series for the investigation of remote triggering of moderate to large earthquakes (M _w ≥ 4) in Taiwan. We first analyze the time-delay between any regional event with M _w ≥4 from 2000 to 2022 and the global earthquake with M _w ≥ 7 that precedes it (Supplementary Section S1; Supplementary Figure S2 in the Electronic Supplement) and then establish a selection of regional earthquakes based on a magnitude-dependent fluid-induced nucleation framework (Parsons et al., 2017). We then investigate pre-earthquake strain signals and propose an approach to compute far-field dilatation induced by an overpressure front diffusing in a fracture. Finally, we discuss the limitations faced for strain detection of the precursory phase related to hypothetical remote fluid-induced earthquakes.

Among the diverse mechanisms proposed to explain delayed remote triggering of earthquakes with all magnitudes (Parsons, 2005; Peng et al., 2011; Shelly et al., 2011), a mechanism in particular shows that delays observed in the case of triggered earthquakes are likely not randomly distributed but are rather proportional to the event magnitudes (Parsons et al., 2017). The mechanism relies on the fact that dynamic straining from seismic waves can break a fluid seal previously blocked into a fault zone (Brodsky et al., 2003) which then releases overpressurized fluid that progressively invades the fault region. Pore fluid infiltration reduces clamping normal stress (strength is reduced with increasing pore pressure) and thus helps promote failure according to the Coulomb failure criteria. Only critically stressed faults are hydraulically conductive (Barton et al., 1995) and thus prone to fluid-induced triggering. As such a process implies fluid flow through a fault zone, a delay is required before failure conditions are met. This fluid transit time increases with earthquake magnitude since the latter is directly proportional to the dimensions of the fault rupture (Kanamori, 1977). Besides, earthquake magnitude also scales with the critical nucleation dimension (2L _c , in meter) which represents the size of the smaller asperity where rupture may initiate (Ohnaka, 2000): 2 L c = 1 0 − 9 M 0 3 (1)

where M ₀ represents the scalar seismic moment (in N.m). The critical nucleation dimension represents the minimum transit distance for a fluid to induce an earthquake of a given magnitude, and measures 0.1 km to about 1.5 km for a magnitude range between 4 and 6.5, for example. Fluid transit from a highly pressurized source through a porous medium follows a diffusion behavior (Malagnini et al., 2012): r = 2.32 D t (2)

where r is the distance from the pressure source, t is time and D is the hydraulic diffusivity.

Since the time-delay between a transient stressing of the Earth’s crust and the eventual occurrence of remote earthquakes may be diagnostic of their nucleation process (Parsons et al., 2017), we establish a regional earthquake selection that may be compatible with a possible fluid-induced nucleation model. We retain the regional events for which estimated time-delays may be consistent with commonly observed fault zone diffusivity values (typically D ∼ 10^–2 to 10 m².s ⁻¹) (Parsons et al., 2017; Guo et al., 2021). We extend our selection to a delay of approximately 1 week to include some of the largest events (M _w ≥ 6.4) from the preliminary selection (Supplementary Section S1). We estimate maximum trough-to-peak dilatation of the Rayleigh wave using the waveform similarity between dilatation ϵ _v and the vertical component of the seismic acceleration a _z (Canitano, 2020): ϵ v = 0.12 T c π V r 1 − 2 γ γ a z (3)

where T _c is the Rayleigh wave main period (T _c = 15 s), γ denotes the Poisson’s ratio (γ = 0.25) and V _r the Rayleigh wave phase velocity (V _r = 4 km.s ⁻¹). A good waveform coherence is observed between signals (Supplementary Figure S3) while amplitude discrepancies (about 10%–15%) may be related to uncertainties in the borehole strainmeter calibration protocol (Canitano et al., 2018). We consider a Rayleigh wave trough-to-peak strain amplitude of 2 nϵ as limit and we adjust the delay by estimating the time-difference between the arrival time of the Rayleigh wave with maximum amplitude at the strong-motion station and the regional event onset (Figure 1). Maximum Rayleigh wave amplitude ϵ _m associated with the retained 22 regional events are ranging from 5 nϵ to about 900 nϵ with epicentral distances from about 1,400 km to 18,000 km (Supplementary Figure S3; Table 1). Regional events are showing delays ranging from about 20 min to 6 days and M _w from 4.12 to 6.80 (Table 2). Figure 2A presents the retained earthquakes integrated in a magnitude-dependent critical nucleation model.

We investigate pre-earthquake anomalies using Sacks-Evertson (Sacks et al., 1971) borehole strainmeter sensors. Given their very high sensitivity at short-to intermediate-period (minutes to weeks), strainmeters represent a powerful tool to search for crustal strain transients, including anomalies preceding earthquakes (Amoruso and Crescentini, 2010; Canitano et al., 2015). We process the 1-min sampling strain time-series [see Canitano et al. (2018); Canitano et al. (2021) for details] for stations located at a maximum distance of about 50 km from our target events (11 events). Other events are either too far from operating stations ( > 80 km) or occurred prior to the network deployment. Figure 3 and Supplementary Figure S4 present the temporal evolution of dilatation over a time-span of minutes to days preceding regional events.

In general, dilatation level after correction of external perturbations ranges from approximately 10^–1 to 5.10^–1 nϵ at the period of minutes to hours and then increases to about 1 nϵ at a daily period to a tens of nanostrain at the period of a few days. For events 1 and 2, which occurred during the passage of surface waves from global earthquakes, we observe no pre-rupture strain variations at remote strainmeters (40 km away). Event 3 (M _w 4.80) possibly represents the most favorable case for detecting subnanometric near-source strain variations (Figure 3C). It occurred at shallow depth (7 km), about 14 km away from SJNB station, and 2 h following the Rayleigh wave arrivals of a M _w 7.7 earthquake in Papua. Strain variations are strongly correlated with atmospheric pressure for about 45–50 min after the Rayleigh wave arrivals (15 min after the mainshock origin time), then suddenly the correlation breaks, suggesting the possible detection of an other source of deformation. We observe a gradual expansion starting about 1 h preceding the regional earthquake with a total deformation of about −7.10^–1 nϵ, well above the short-period strain noise of approximately 10^–1 nϵ. For event 11, large pore pressure variations ( > 25 kPa) induced by the intense Rayleigh wave straining in southern LV (about 1 μϵ) following the 2008 Eastern Sichuan earthquake may be expected. However, we find no evidence of pre-earthquake variations, which illustrates the complexity for detecting subnanometric deformation at periods larger than 12–24 h because signals are also impacted by environmental perturbations, especially rainfall (Hsu et al., 2015; Mouyen et al., 2017), and by low tidal strain noise remaining after correction (∼ 0.5-1 nϵ) (Supplementary Figure S4). Finally, strain contraction of −10 nϵ and −20 nϵ are observed hours to days preceding event 16 and event 21, respectively, and are analyzed in the next section.

To connect the strain observations to potential transient pressure anomalies, we follow a two-step approach. In a first step, we calculate the temporal variation of the theoretical pore pressure induced in the fracture by the passing seismic waves. In a second step, we compute the time-dependent evolution of dilatation generated by a remote diffusing pressure front at the sensor location.

The detection of crustal strain anomalies and their characterization as earthquake precursors is a complex problem (Bernard, 2001). Although the quality of the data and of the applied corrections play an important role in monitoring strain transients (e.g., the detection level), complexities are mainly related to the transient source characteristics and to the source-station distance (quasi-static deformation decreases with a factor of 1/R ²). In particular, the detection capability strongly depends of the strength of the source of deformation. Here, the latter is controlled by the size of the nucleation region (which depends of the earthquake magnitude) and by the pressure gradient, which is determined by the fracture hydraulic diffusivity and by the elastic wave pressure. We consider the volume change ΔV in the crack fracture resulting from a uniform pressure gradient ΔP as a proxy for the source strength. For an incompressible fluid, it is expressed as: Δ V = − 4 π 1 − ν Δ P G a 3 ∫ 0 1 t ϕ t d t (6)

where ϕ function is the solution of the Fredholm equation of the second kind (Fialko et al., 2001; Battaglia et al., 2013). Figure 6 presents the theoretical volume change in the fracture for earthquakes with magnitude ranging from 4 to 7. To detect the precursory phase of a remote fluid-induced earthquake at a subdaily period within a few kilometer radius from the source requires a source volume that is about 5–10 times larger than for the hypothetical cases considered in this study. Since strain measurement noise increases with period (Crescentini et al., 1997), this would require monitoring crustal deformation with a station located near a large regional earthquake that occurred shortly following a large teleseismic event; a configuration that is highly unlikely. Therefore, to improve our understanding of remote earthquake triggering in Taiwan, the protocol adopted here should be combined with a more conservative approach that relies on the search of earthquakes with a wide range of magnitudes following large global events.

Since it controls the dilatation level during computation, the incident-wave strain level may also represent a critical parameter for detection of fluid-induced seismicity. However, the effect of dynamic peak strain on triggered seismicity is not well understood (Brodsky and van der Elst, 2014), and a large strain perturbation doesn’t necessarily imply remote triggering, and vice versa (Parsons et al., 2014). The 2008 M _w 7.9 Eastern Sichuan earthquake, which generated strong regional ground motions due to pronounced rupture directivity effects (Kurahashi and Irikura, 2010), has induced dynamic stress over Taiwan about twice larger than any other great teleseismic earthquake (e.g., the 2011 M _w 9.0 Tohoku and the 2012 M _w 8.6 Sumatra events). Although, the event triggered seismicity in northern China (Peng et al., 2010), it had little to no impact on the seismicity in Taiwan (Chao et al., 2011), at least during the first days following the wave passing (Section 3). Recently, field (Bonini, 2020) and laboratory (Zheng, 2018; Jin et al., 2021) observations have shown that pore pressure induced in a fluid-filled fracture by a seismic wave can be amplified relative to the incident-wave pressure [up to three orders of magnitude (Zheng, 2018)]. This transient pore fluid amplification would suggest that even seismic waves with insignificant stress levels can induce substantial overpressure in fault zones. Besides, this also implies that fluid level can remain high after transiting in the fracture, therefore favoring nucleation in case of delayed triggering or inducing large normal stress unclamping in case of instantaneous (or nearly-instantaneous) triggering. Such a process may also favor detection at remote distances of transient events that generate little regional deformation, as analyzed here. For example, in the case of event 3 (Figure 5), we observe that the pore pressure rise-time and temporal evolution in the fracture is well predicted for station SJNB (nearby station TRKB was shutted off), but the strain level is underestimated by about three orders of magnitude. Although the analyze of a pressure amplification mechanism in the fault zone is beyond the scope of this study, we cannot rule out that such a process may be at play for some cases. On the other hand, moderate strain changes observed for events 16 and 21 are too large to reflect a transient deformation and likely represent the contraction response due to rainfall loading (Canitano et al., 2021), which illustrates the impracticability for detecting precursory strain during environmental disturbances.

Finally, we recognize that our analogue volcano deformation model is by no means complex enough to reflect the far-field deformation caused by fluid flow in active regions, but it provides a preliminary framework for the integration of borehole dilatation observations for the search of precursory signals associated with remote fluid-induced earthquakes (Parsons et al., 2017). A more elaborated model of fault valving (Zhu et al., 2020), integrating fluid-driven aseismic creeep (Shelly et al., 2011; Cebry et al., 2022) and seismicity (Shapiro et al., 2018), combined with a multidisciplinary geophysical monitoring is fundamental for enhancing detection and characterization of precursory phenomena in Taiwan (Fu et al., 2020) or in other active regions (Delorey et al., 2015).

We develop a protocol to integrate high-resolution borehole strain measurements for the investigation of remote triggering of moderate to large earthquakes (M _w ≥ 4) in Taiwan. The physical framework is compatible with a magnitude-dependent fluid-induced nucleation process where the nucleation region is simplified using the analogue volcano source represented by a horizontal circular crack. We observe that theoretical evolution of dilatation generated by a remote diffusing pressure front in a semi-infinite, elastic half-space is two to four orders of magnitude smaller than observations (∼ 10^–5 to 10^–3 nϵ). This suggests that detection of pore pressure-related deformation would have required change of volume in the nucleation region that is at least one order of magnitude larger than for the hypothetical cases considered here. The study represents the first attempt to analyze strain time-series for detecting pre-earthquake strain anomalies related to fluid-induced earthquakes and illustrates the challenge for detecting and characterizing intermediate-to far-field earthquake precursors caused by fluid flow in active regions.