Co‐Occurrence of Low and Very Low Frequency Earthquakes Explained From Dynamic Modeling

Very low‐frequency earthquakes (VLFs) are characterized by longer source duration and smaller stress drop than regular earthquakes of similar magnitude. Recent studies have shown their frequent correlation with low‐frequency earthquakes (LFEs) on shared faults. The underlying source processes governing the occurrence of VLFs and their interaction with LFEs remain elusive. Here, we employ a slip‐weakening model for slow earthquakes. By comparing the source parameters of simulations and observations, it is suggested that VLFs are slow self‐arresting earthquakes that self‐terminate within the nucleation patch. Additionally, we adopt a composite model to reproduce the records of the simultaneous occurrences of a VLF and an LFE in the Nankai area. Our results present the possibility that VLFs, LFEs, and regular earthquakes can be distinguished using a unified dynamic framework.

While SSEs and nonvolcanic tremors/LFEs have been extensively studied (R. Ando et al., 2010;Ide, 2008;Kato, 2003;Liu & Rice, 2007;Rubin, 2008;Segall et al., 2010;Shibazaki & Iio, 2003), investigations into the source processes of VLFs remain comparatively scarce.Accumulating evidence shows that VLFs, like regular earthquakes, are caused by shear slip (Ghosh et al., 2015;Ide & Yabe, 2014).However, interpreting VLFs within the context of fast earthquakes poses a challenge, given that the characteristic duration of VLFs with magnitude 3.1-3.8are ∼20 s (Ito & Obara, 2006b;Matsuzawa et al., 2009), equivalenting to the source duration of a regular earthquake of magnitude 6-8.And the average stress drop and slip rate of VLFs are 2-3 orders of magnitude smaller than regular earthquakes of similar magnitude (Ito & Obara, 2006b).In addition to their inherent source characteristics, VLFs are frequently observed to spatiotemporal correlated with LFEs (Ghosh et al., 2015;Hutchison & Ghosh, 2019;Nakamura, 2017;Walter et al., 2013).Specifically, observations (Ito et al., 2007;Takemura, Matsuzawa, et al., 2019;Takeo et al., 2010;Tamaribuchi et al., 2019) show that the Nankai area has experienced synchronous VLF and LFE events up to 5 km apart, suggesting the potential for LFEs to occur within the VLF rupture ranges.These insights imply that both VLFs and LFEs may arise from shared processes along the plate interface during slow earthquakes.Nevertheless, the mechanism accountable for their concurrent manifestation remains elusive.In this scenario, developing a unified source model capable of generating VLFs, LFEs, and regular earthquakes may offer further insights into slow stress-releasing processes.
Here, we present a simple numerical dynamic source model to reproduce the source characteristics of VLFs and the co-occurrence of VLF and LFE.By simulating earthquakes under various frictional and stress conditions, it is indicated that VLFs might represent slow self-arresting ruptures (SSARs), stopping spontaneously within the nucleation zone.Comparing the source parameters of SSARs with those of observed VLFs reveals similarities in their source characteristics.Furthermore, we introduce a composite SSAR model that incorporates both a VLF and an LFE within a unified nucleation process.By comparing the theoretical records from this model with observations in the Nankai region, the results imply the capability to simulate two types of slow earthquakes with distinct corner frequencies.Our work suggests that SSARs at different scales could produce seismic events similar to LFEs and VLFs, with underlying source mechanisms linked to those of regular earthquakes.

Materials and Methods
Our starting hypothesis considers both VLFs, LFEs and regular earthquakes are individual events caused by shear sliding.Assuming that VLFs are generated from asperity patches on planar faults, we use a 3-D dynamic rupture model based on the slip-weakening friction law (Andrews, 1976a(Andrews, , 1976b;;Ida, 1972) (Figure 1a).We consider a homogeneous slip-weakening patch surrounded by a medium that cannot rupture (unbreakable medium) (Figure 1b).The slip-weakening friction law is as follows: where σ u is the peak shear strength, σ r is the residual stress, D is slip and D c is the critical slip distance.Once the rupture is triggered, the patch strength drops from peak strength σ u to a small residual value, denoted σ r .We denote the breakdown stress drop (σ u σ r ) as T u and the dynamic stress drop (σ 0 σ r ) as T e .We use the boundary integral equation method (Chen & Zhang, 2006;Zhang & Chen, 2006a, 2006b) to simulate the evolution of fault slip.Notably, the residual stress σ r , which appears in both the slip-weakening friction law and the boundary integral equation, will be eliminated in the final solution.Therefore, we assign σ r to 0 for simplicity without sacrificing generality (J.K. Xu et al., 2015).Initially, we consider that the nucleation patch is prestressed at a level of 90%T u (0.9T u ), while the background is prestressed with T e .By applying a stress perturbation (Text S1 in Supporting Information S1), we elevate the shear stress on the nucleation patch to 1.001T u to trigger the rupture (Andrews, 1976b;Day, 1982;J. K. Xu et al., 2015).The triggering starts from the center of the nucleation zone and is forced to propagate by a predefined η across the entire nucleation zone, ensuring asynchronous rupturing within the nucleation patch.After initiation, the subsequent process (including propagation and self-termination) is spontaneous and inherently governed by the slip-weakening law.Considering that simultaneous triggering of the nucleation zone over a large area is demanding in reality, we adopt this nucleation approach to better emulate real-world scenarios.
To provide a unified description of rupture under different slip-weakening parameters, based on previous studies (Wei et al., 2021;J. K. Xu et al., 2015), we adopt two dimensionless parameters T e , Dc ) : where R a is the effective radius of the nucleation asperity, and μ is the shear modulus of the medium.For a simulation, we fix the density, wave speeds within the nucleation patch unchanged so that the rupture style depends on the frictional parameters.

Phase Diagram for Simulating Dynamic Rupture
We set the radius of the slip-weakening friction patch to 10 km, the radius of the nucleation asperity to 3 km, T u to 1 MPa and the P-wave and S-wave velocities to 3 km/s and 1.25 km/s, respectively.By varying the normalized parameters T e , Dc ) from 0 to 1, we investigate the potential rupture styles and summarize them in the phase diagram.Each point in the phase diagram indicates a rupture simulation with the friction parameters at that point (Figure S1 in Supporting Information S1).Pioneering studies (Andrews, 1976b;Das & Aki, 1977;Madariaga & Olsen, 2000) have elucidated the influence of different combinations of T e , Dc ) on rupture conditions.Madariaga and Olsen (2000) delineated the bifurcation between supershear and sub-Rayleigh ruptures.Subsequent research (Wei et al., 2021; J. K. Xu et al., 2015) introduced the concepts of self-arresting rupture and SSAR.Currently, the phase diagram of strike-slip rupture in full space visualizes four different styles of rupture: supershear rupture, sub-Rayleigh rupture, self-arresting rupture and SSAR (Figure 1c).Supershear rupture only occurs in the in-plane slip direction, in which the rupture velocity is faster than the local S-wave velocity (Bouchon & Vallee, 2003), while sub-Rayleigh rupture denotes the rupture style of a regular earthquake in which the rupture velocity is slower than the local S-wave velocity (Figure 1d).Both types of rupture propagate across the whole slip-weakening patch until it encounters and is stopped by the unbreakable barrier.Self-arresting ruptures and SSARs (Figure 1e) refer to ruptures that self-terminate (i.e., without any external interference) (Wei et al., 2021;Wen et al., 2018;J. K. Xu et al., 2015;D. Y. Xu et al., 2023).The rupture front of the selfarresting earthquake continuously propagates for a short distance after exceeding the nucleation patch but spontaneously stops before encountering the unbreakable medium (Figure S2 in Supporting Information S1).More particularly, SSARs stop spontaneously within the nucleation asperity patch (Wei et al., 2021).SSARs emerge when Dc is relatively large, although T e changes across a rather wide range.Despite the rupture patch of an SSAR is only the nucleation zone, it has a longer source duration than a sub-Rayleigh rupture.Furthermore, the slip rate and released stress of an SSAR is 1-2 orders of magnitude smaller than that of the sub-Rayleigh and selfarresting ruptures (Figures S3 and S4 in Supporting Information S1).Building upon the slow characteristics of SSARs, which align with the source characteristics of VLFs, we proceed to investigate the potential of SSARs in elucidating the source processes of VLFs.

Comparison Between SSARs and VLFs
We quantitatively compare the source parameters between VLFs in the shallow offshore accretionary prism in the Nankai and Ryukyu area (Ide et al., 2007;Ito & Obara, 2006a, 2006b;Ito et al., 2009;Nakamura & Sunagawa, 2015;Passelègue et al., 2020;Sugioka et al., 2012;Takemura et al., 2022;Tamaribuchi et al., 2019) (Table S1 in Supporting Information S1) and SSARs here.We assume that sub-Rayleigh ruptures, self-arresting ruptures, and SSARs are potential models for VLFs.Using the phase diagram in the medium where VLFs occur (Figure 1c), we simulate these three types of ruptures with 149 possible frictional parameter combinations.To compare with the observed VLF, we set the rupture patch radius to 3 km and the dynamic stress drop T e to 1, 5, 10 MPa, respectively.As we employ a time-varying stress perturbation governed by η for rupture triggering (Text S1 in Supporting Information S1), the source duration in the simulation results comprises two components: the initiation time we have designated and the duration of propagation or self-arrest.We simulate sub-Rayleigh ruptures, self-arresting ruptures, and SSARs with η ranging from 200 to 1,200 (and ∞) (Text S2 in Supporting Information S1).The results indicate that decreasing η primarily extends the source duration of ruptures, while its impact on the simulated results of slip, stress drop, and seismic moment is limited (<1%) (Figures S5-S7 in Supporting Information S1).Based on observations of the rupture velocities and source durations of VLFs (Ito & Obara, 2006b), we set an average η of 1,000 (Table S2 in Supporting Information S1).
A comparison among all of the simulated rupture source parameters with those of VLFs suggests that SSARs have the potential to meet the source parameters of VLFs (Figure 2).For the slip, average slip rate, seismic moment magnitude and stress drop, we find that the simulated parameters of self-arresting and sub-Rayleigh earthquakes are generally 2-3 orders of magnitude larger than those of recorded events; in other words, in a rupture area with a radius of 3 km, these two kinds of ruptures cannot produce earthquakes as small as VLFs.An SSAR exhibits the ability to generate a relatively long source duration with limited slip and stress drop, owing to the decrease in shear strength along its rupture front that does not reach the minimum σ r .Thus, with the same rupture patch area, SSARs produce much smaller stress drops and slip rates.Notably, the VLFs in the Nankai area we compare here are shallow VLFs.We also simulate and compare the deep VLFs (Chaudhuri & Ghosh, 2022;Ghosh et al., 2015;Ide & Yabe, 2014;Ide et al., 2007;Ito et al., 2007Ito et al., , 2009;;Matsuzawa et al., 2009Matsuzawa et al., , 2015;;Maury et al., 2016;Passelègue et al., 2020) using a similar approach (Table S3 in Supporting Information S1).The source processes of SSARs in deep and shallow media are similar (Figures S8 and S9 in Supporting Information S1).The results suggest that SSARs also display characteristics akin to deep VLFs in terms of source parameters (Figure S10 in Supporting Information S1).Exhibiting source characteristics similar to VLFs based on the model parameter settings we have established, SSAR may serve as a potential source model for VLFs.

An SSAR Model for Synchronous Generation of a VLF and an LFE
LFEs and VLFs are both considered slow earthquakes due to their weak signals and low corner frequencies (Ito et al., 2007).Based on the comparison results in Section 3.2 and our previous work (Wei et al., 2021), the similarities among the source characteristics of SSARs, LFEs, and VLFs imply an interconnection in their rupture processes.
Here, we provide a potential composite SSAR model to interpret the relation between VLFs and LFEs (Figure 3a).(Ito & Obara, 2006a, 2006b;Ide et al., 2007Ide et al., , 2009;;Nakamura & Sunagawa, 2015;Passelègue et al., 2020;Sugioka et al., 2012;Takemura et al., 2022).Orange, blue, and green dots denote sub-Rayleigh earthquakes, self-arresting earthquakes, and slow self-arresting earthquakes, respectively.Different dot shapes represent ruptures with varying T e (circles for 1 MPa, triangles for 5 MPa, and squares for 10 MPa, respectively).The gray dashed rectangles show the ranges where the source parameters of the simulated earthquakes match those of VLFs.Initially, the simulation is configured with P-wave velocity at 4.9 km/s, S-wave velocity at 2.9 km/s, density at 2,500 kg/m 3 (Takemura, Kubo, et al., 2019), and T e at 1 MPa.We employ SSARs with radii of 5 km and 200 m for VLFs and LFEs, respectively (Ito & Obara, 2006b;Thomas et al., 2016).By adjusting the slip-weakening parameters and employing a trial-and-error method, we select the Dc and T e that best matches the observed source durations and recorded amplitude.Specifically, we use an SSAR with ( Dc , T e ) = (0.85, 0.9), representing (T e , T u , D c ) = (1.0MPa, 1.11 MPa, 0.20 m), to simulate the VLF.Similarly, we employ an SSAR with ( Dc , T e ) = (0.85, 0.1), denoting (T e , T u , D c ) = (1.0MPa, 10.0 MPa, 0.20 m) for the LFE (Figure 3b).The nucleation patch of the LFE is located within the nucleation patch of the VLF, with the model settings mainly refer to their relative positions.We still use the nucleation method described in Text S1 in Supporting Information S1.By setting the loading coefficient η to 1,000 m/s (Ito & Obara, 2006b), stress grows progressively from the center of the VLF, eventually enveloping the entire VLF and LFE nucleation patch.When the stress perturbation extends to the position of the LFE nucleation patch, stress that slightly exceeds its yielding threshold triggers its rupture.The variance in triggering times between LFEs and VLFs is also regulated by η (Figure S12 in Supporting Information S1).The influence of η is primarily manifested in the source duration and slip rate of both LFE and VLF.However, the impact of η on the simulated LFE and VLF in terms of stress drop, slip, and seismic moment is relatively minimal (Tables S4 and S5 in Supporting Information S1).With this composite model, we obtain rupture snapshots (Figure 3c) and the average slip rate history (Figure 3d) for this VLF and LFE.The rupture processes of the two SSARs exhibit mutual influence.Despite being influenced by stress disturbances from the larger SSAR, the smaller SSAR maintains its SSAR characteristics without transitioning into a sub-Rayleigh rupture.According to the simulation, the large SSAR (VLF) exhibits a source duration of 15.6 s with an M w of 3.94, while the small SSAR (LFE) that initiates within the large SSAR patch has a source duration of 0.64 s with an M w of 2.16, in agreement with the observations (Tables S4 and S5 in Supporting Information S1).The maximum slip rate of the LFE obtained from the simulation is approximately 10 times that of the VLF, which is consistent with the observed records.
To further confirm the similarity between the composite model and the VLF and LFE sources, we compare the waveforms recorded at each station with those obtained from the numerical simulations.For the observations, we remove the instrument response and filter the seismic records of F-net stations to the VLF frequency band, that is, 0.02-0.05Hz.For the synthetic seismogram, we use the slip rate of the fault plane obtained from our dynamic simulation as a point seismic source.Substituting it into the generalized reflectance transmission coefficient method (GRT) (Chen, 1993) simulation, we calculate the synthetic waveform at each F-net station, adopting a homogeneous layered model with wave velocities derived from Takemura, Kubo, et al. (2019).We set the fault strike of the VLF and LFE to 284.75°, dip to 14.84°, and rake to 146.81° (Takemura, Matsuzawa, et al., 2019;Tamaribuchi et al., 2019).To compare the VLF record, the theoretical waveforms are also filtered to 0.02-0.05Hz.The synthetic VLF waveforms are in agreement with the Nankai VLF records (Figure 4).The similarity of the two waveforms and the consistency of the source parameters suggest the potential capability of the SSAR to model two slow earthquakes of distinct scales and source durations.

Discussions and Conclusions
With the composite SSAR model, we generate LFEs and VLFs simultaneously.The D c for both the VLF and LFE are the same, indicating that the LFE and VLF share a nearly identical medium frictional condition.This property may help explain why LFEs and VLFs sometimes occur together in proximity.The T u of the LFE is greater than that of the VLF, implying that the LFE may occur in a medium with a higher yield strength and in the VLF patch with a more concentrated region of stress.Our research indicates that regular earthquakes and slow earthquakes represent distinct seismic source processes occurring under varying frictional conditions.Observations (Wang et al., 2023) further support this perspective by demonstrating that the low-frequency characteristics of LFEs primarily arise from an atypical rupture process rather than near-source attenuation.Intriguingly, in the Alaska region, seismic nucleation emits VLF-like signals before transitioning into/triggering regular earthquakes (Tape et al., 2018).Based on the phase diagram indicating that the critical conditions for SSARs and sub-Rayleigh ruptures are primarily determined by Dc , this suggests the potential for Dc to vary in nearby regions, hinting at the intricacies of sliding events in reality.
In the simulation, we introduce η to trigger points within the nucleation zone sequentially after a certain time delay.η primarily prolongs the duration of ruptures.Across ruptures with a 3 km effective radius, varying η from 200 to 1,200 (and ∞) potentially extends the source durations of SSARs by 31%, self-arresting ruptures by up to 26%, and sub-Rayleigh ruptures by 11% (Text S2 in Supporting Information S1).Consequently, when estimating the stress drop using the Brune model, the stress drop for SSARs could decrease by up to 55%, for self-arresting ruptures by 50%, and for sub-Rayleigh ruptures by 27% (Gallovic & Valentová, 2020), despite η having minimal impact on the actual stress drop, slip, and seismic moment of these ruptures (<1%).In the composite model, as η transitions from 1,000 to 2,000, the interval between the VLF and the LFE shortens from 1.26 to 0.63 s.This results in a 15% decrease in source duration for VLF and a 13% decrease for LFE (Tables S4 and S5 in Supporting Information S1), indicating the impact of the triggering in the nucleation zone on rupture behavior.
Additionally, various physical models have been developed to simulate slow earthquakes.For example, Kato (2003), Shibazaki and Iio (2003), and Shibazaki and Shimamoto (2007) employed friction laws with a cutoff velocity to model SSEs.Liu and Rice (2007) and Rubin (2008) investigated SSEs using the rate and state friction law.These pioneering studies portray slow earthquakes as essentially distinct-scale nucleation events without subsequent instability, a notion that resonates with SSAR.However, research on dynamic models of small-scale VLFs is relatively limited and entails additional complexity.Nakata et al. (2011) and Wu et al. (2019) examined fault zones with brittle-ductile rheology, generating VLF signals through a dynamic frictional-viscous model.These models prevalently adopt rate-dependent frictional laws to simulate slow earthquakes, while our study proposes the potential to interpret VLFs using a simple linear slip-weakening friction law.
Based on dynamic source simulations combining the slip-weakening friction law, we develop a potential source model for VLFs.Through a comparison between the source parameters derived from observations and those of numerically simulations, the results indicate that VLFs may be SSARs in nature, which self-terminate within the nucleation patch.We further simulate the VLFs and LFEs that occur in the Nankai Trough region using a composite SSARs model.The consistency between synthetic waveforms and observations suggests that the VLFs and LFEs are SSEs, just like regular earthquakes, are controlled by moment and corner frequency.Our results indicate that VLFs and regular earthquakes rupture under distinct frictional conditions, which can be comprehensively explained within the framework of the dynamic rupture phase diagram.The SSAR model may provide broad perspectives for the study of slow earthquakes and help us to infer the source processes of those seismic events.

Figure 1 .
Figure 1.Numerical rupture simulation model and dynamic rupture phase diagram.(a), Slip-weakening friction law.(b), Schematic diagram of the model geometry and frictional properties.The circular nucleation patch (yellow) is embedded in a slip-weakening region (gray) surrounded by an unbreakable fault domain (white).(c), Rupture phase diagram for strike-slip in full space.The open triangles (blue), open circles (red), asterisks (yellow), and stars (pink) represent supershear, sub-Rayleigh, self-arresting ruptures, and slow self-arresting ruptures (SSARs), respectively.(d), Sub-Rayleigh rupture.The blue background represents the simulated faults.White dashed circles outline the nucleation patch, while orange dashed circles enclose the slip-weakening region.(e), SSAR.

Figure 3 .
Figure 3. Source model of Nankai very low-frequency earthquakes (VLFs) and low-frequency earthquakes (LFEs).(a), Schematic of the seismic source distribution for Nankai VLFs and LFEs.VLFs (red circles) and LFEs (orange circles) occur in close proximity to each other.(b), Source models of the studied VLF and LFE in the Nankai region.The nucleation patch of the LFE (green dot) is located inside the VLF nucleation patch (red circle).The gray patch is the slip-weakening area.(c), Slip rate history for the VLF and LFE.Within the solid red line is the nucleation patch of the LFE.(d), Average slip rate history for the VLF and the LFE.

Figure 4 .
Figure 4. Comparison between the synthetic velocity seismogram of the composite model and the observed very low-frequency earthquake (VLF).(a), Comparison of the synthetic seismogram and VLF records.The red lines represent the synthetic velocity seismogram, and the blue lines represent the VLF recorded by F-net stations.(b), Map of the Nankai area, F-net stations, and associated topography.Low-frequency earthquakes (LFEs) are shown as orange dots, and the F-net stations are shown as purple triangles.The VLF, represented by the red star at coordinates (134.7°E,32.70°N), and the LFE, indicated by the blue dot at (134.64°E, 32.71°N), are used for comparison.