# Subdaily Slow Fault Slip Dynamics Captured by Low-Frequency Earthquakes

Caroline Mouchon, William B Frank, Mathilde Radiguet, Piero Poli, Nathalie Cotte

► To cite this version:

Caroline Mouchon, William B Frank, Mathilde Radiguet, Piero Poli, Nathalie Cotte. Subdaily Slow Fault Slip Dynamics Captured by Low-Frequency Earthquakes. AGU Advances, 2023, 4 (4), 10.1029/2022AV000848. [hal-04201946]

HAL Id: hal-04201946

https://hal.science/hal-04201946v1

Submitted on 11 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
AGU Advances

RESEARCH ARTICLE

10.1029/2022AV000848

**Peer Review** The peer review history for this article is available as a PDF in the Supporting Information.

**Key Points:**

* Low-frequency earthquake activity can reproduce the geodetic record of the slow slip cycle

* Slow slip is the competition between subdaily dynamics and a fully coupled fault

* Low-frequency earthquakes are incidental to slow slip and do not play a role in mediating aseismic slip

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

C. Mouchon,
mouchonc@mit.edu

**Citation:**

Mouchon, C., Frank, W. B., Radiguet, M., Poli, P., & Cotte, N. (2023). Subdaily slow fault slip dynamics captured by low-frequency earthquakes. AGU Advances, 4, e2022AV000848. https://doi.org/10.1029/2022AV000848

Received 17 NOV 2022
Accepted 17 MAY 2023

**Author Contributions:**

**Conceptualization:** William B. Frank,
Mathilde Radiguet

**Data curation:** Caroline Mouchon

**Formal analysis:** Caroline Mouchon

**Funding acquisition:** William B. Frank,
Mathilde Radiguet, Piero Poli

**Project Administration:** William B.
Frank

**Resources:** Nathalie Cotte

**Supervision:** William B. Frank, Mathilde
Radiguet

**Validation:** Caroline Mouchon, William
B. Frank

**Visualization:** Caroline Mouchon

© 2023. The Authors.

This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

Subdaily Slow Fault Slip Dynamics Captured by
Low-Frequency Earthquakes

Caroline Mouchon<fnref n="1" /> ID, William B. Frank<fnref n="1" /> ID, Mathilde Radiguet<fnref n="2" /> ID, Piero Poli<fnref n="3" />, and Nathalie Cotte<fnref n="2" /> ID

**Abstract** Geodetic positioning is the geophysical record of reference for slow slip events, but typical daily solutions limit studies of the evolution of slow slip to its long-term dynamics. Accompanying seismic low-frequency earthquakes located precisely in time and space provide an opportunity to image slow slip dynamics at subdaily time scales. Here we show that a high-resolution time history of low-frequency earthquake fault slip alone can reproduce the geodetic record of slow slip that we observe to be dominated by subdaily fault slip dynamics. However, a simple linear model cannot accommodate the complex dynamics present throughout the slow slip cycle, and an analysis of different phases of the slow slip cycle shows that the ratio of geodetic to seismic fault slip varies as a function of time. This suggests that the low-frequency earthquake source region saturates as slow slip grows in moment and area. We propose that rheological heterogeneities at the plate boundary associated with low-frequency earthquakes do not play a significant role in the slow slip rupture process, thus implying that their activity is incidental to the driving aseismic slip.

**Plain Language Summary** Slow slip events can be observed in many subduction zones where they play an important role in the earthquake cycle. Decades after their discovery, slow slip events are now captured routinely in geodetic datasets with slip dynamics occurring over a broad range of time scales. Using high-time resolution seismological observations together with the geodetic record allows us to go beyond the coarse daily GNSS sampling rate to image slow slip dynamics at short time scales. Here we use the temporal evolution of seismic slip produced by low-frequency earthquakes to study the subdaily dynamics of a slow slip event cycle, reproducing the geodetic record of slow slip using only seismological observations. We develop a simple model where long-term loading is in competition with the intermittent release of stress tied to the seismic slip of low-frequency earthquakes. We show the full slow slip cycle is driven by bursts of slip at subdaily time scales that low-frequency earthquake events witness only in their immediate source region. This result implies that the low-frequency earthquake rupture process is incidental to slow fault slip and does not play a major role in the slow slip cycle.

1. Introduction

Slow slip is a major constituent of the seismic cycle (Avouac, 2015; Bürgmann, 2018), directly impacting the timing and extent of neighboring earthquakes (Obara & Kato, 2016). Predominantly present downdip of the seismogenic megathrust in subduction zones, slow slip can release as much built-up tectonic stress as major earthquakes (Beroza & Ide, 2011). While the fault motion associated with slow slip is too slow to radiate seismic waves itself, slow slip is often accompanied by spatiotemporally correlated seismic events (Rogers & Dragert, 2003), including tectonic tremor (Obara, 2002), low-frequency earthquakes (Shelly et al., 2006), and very low-frequency earthquakes (Ito et al., 2007). This seismicity is observed and detected over much shorter time scales than the geodetic record that directly captures slow slip, and thus tracks the evolution of aseismic fault slip. Recent work argues that these seismic signals and geodetically recorded slow slip events are distinct spectral signatures of the same broadband "slow earthquake" (Hawthorne & Bartlow, 2018; Kaneko et al., 2018; Masuda et al., 2020). Regardless of whether all of these geophysical signals represent the same phenomenon, the typical daily positioning derived from GNSS data captures the geodetic signature of slow slip at a low temporal resolution. To access slow slip dynamics at subdaily time scales, we use here seismic slow earthquakes as a guide to extract observational constraints from the geodetic record.

<footnote n="1">¹Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA,</footnote>

<footnote n="2">²Institut des Sciences de la Terre, Université Grenoble Alpes, Centre National de Recherche Scientifique, Grenoble, France,</footnote>

<footnote n="3">³Dipartimento di Geoscienze, Università di Padova, Padova, Italy</footnote>
**Writing – original draft:** Caroline Mouchon, William B. Frank

**Writing – review & editing:** Caroline Mouchon, William B. Frank, Mathilde Radiguet, Piero Poli, Nathalie Cotte

We focus here on low-frequency earthquakes, whose impulsive seismic arrivals and highly repetitive sources (W. B. Frank et al., 2014) act as a precise in situ monitor of slow slip activity (W. B. Frank & Brodsky, 2019). Tectonic tremors (noisy, long-duration seismic signals without any identifiable seismic phases) and very low-frequency earthquakes (observed beneath the microseismic frequency band between 0.05 and 0.2 Hz) are not as effective a proxy of slow slip at short time scales because of their emergent waveforms and long (>30 s) source durations. Low-frequency earthquakes are small ($M \le 2$) shear ruptures (Ide et al., 2007) thought to occur on small localized asperities within a predominantly aseismic fault rheology (Chestler & Creager, 2017; Thomas et al., 2018), with a characteristic deficit of high-frequency seismic radiation compared to ordinary earthquakes. Repeatedly driven to failure by the surrounding slowly slipping fault, low-frequency earthquakes only represent a fraction of the total geodetic moment release. As an indicator of underlying aseismic slip, the timing (W. B. Frank et al., 2018) and size (W. B. Frank & Brodsky, 2019) of low-frequency earthquakes can provide constraints on the dynamics of slow slip at higher resolution than geodesy alone can resolve.

The principal challenge in quantifying slow slip dynamics through the lens of low-frequency earthquake activity is bridging the nearly seven orders of time scales between seismic observables captured at $\ge 40$ Hz sampling rate and GNSS surface positioning typically estimated daily. One recent study bridged this gap to demonstrate that the seismic moment rate of the low-frequency earthquake source scales with the driving moment rate of slow slip (W. B. Frank & Brodsky, 2019). This empirical relationship implicitly redefines the rupture duration of slow slip to periods of low-frequency earthquake activity, ignoring the intermittent periods of loading hidden within the long-term geodetic signature of slow slip (Fujita et al., 2019; Rousset et al., 2019), to establish that slow slip moment scales with its duration cubed similar to ordinary earthquakes. This result demonstrates that observational constraints at the shortest time scales are needed to capture the physical mechanisms governing slow fault slip, the dominant faulting process just downdip of the seismogenic source region. We thus seek to push the limits of combining seismic and geodetic data to explore whether a catalog of low-frequency earthquakes can alone predict the geodetic record of a slow slip cycle.

Here, we develop a simple model to connect the history of low-frequency earthquake fault slip to the geodetic record of slow slip to access the dynamics of fault motion at daily and subdaily time scales throughout the slow slip cycle. We assess whether the cumulative slip of low-frequency earthquakes can reproduce slip observed geodetically during a slow slip cycle. The linear relationship (mediated by the Green's function coefficient) between surface displacements and fault slip at depth (Okada, 1985) presents the opportunity to directly project surface GNSS positioning at every sample onto the subduction interface with an assumed source model. We utilize different seismic observables that can be measured at the same sampling rate as the GNSS solutions to estimate low-frequency earthquake fault slip. Our model aims to capture the competition between long-term tectonic loading and transient aseismic fault slip. The hypothesis underlying our model is that the long term geodetic record of slow slip can be reproduced by the episodic low-frequency earthquake activity. We apply our model to the Guerrero segment of the Mexican subduction zone that hosts two cycles of slow slip: a M7.5 event every 4 years and M6.4 events every 75 days in average (Figure 1) (Husker et al., 2019; Kostoglodov et al., 2010; Radiguet et al., 2012; Rousset et al., 2017). We examine the seismic and geodetic records from January 2005 to April 2007, the temporal extent of a high-resolution catalog that captured low-frequency earthquake activity before, during, and after a large M7.5 slow slip event that dominated the GNSS position time series from April to November 2006 (Figure 1) (Radiguet et al., 2011).

## 2. Connecting Geodetic and Seismic Estimates of Fault Slip

### 2.1. Converting GNSS Surface Motion Into Fault Slip on the Subduction Interface

We characterize the Guerrero segment along a profile perpendicular to the trench (Figure 1a), where seismic and geodetic networks were deployed along a line perpendicular to the trench crossing Acapulco at the coast and Tempoal in northern Mexico between 2005 and 2008 (Perez-Campos, 2008). We compute daily and subdaily (24- and 6-hr) static GNSS solutions (Figure 1c) at three continuous GNSS stations, COYU, MEZC, and IGUA, positioned with increasing distance from the trench (Figure 1a). We choose these three stations: (a) for their relatively extended coverage of the area with at least one station (MEZC) located above the low-frequency earthquakes sources and close to the aseismic slip produced by the M7.5 slow slip and (b) to build a simple 2D model that runs through the middle of the slow slip source region (Figure 1a) to avoid the complexities of three dimensional effects in the surface displacement time series, such as along-strike migrations and edge effects (Radiguet et al., 2011). The static GNSS time series are processed with a Precise Point Positioning approach using the software GipsyX (Bertiger et al., 2020). Accurate satellite clock and orbit corrections from the Jet Propulsion Laboratory are used.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Geodetic record of slow slip on the Guerrero subduction interface. (a) Tectonic loading and slow slip is captured by three GNSS stations, COYU (blue), MEZC (yellow), and IGUA (orange). The M7.5 slow slip event that dominates the geodetic record generated more than 10 cm of slip (yellow patch) on the subduction interface. This aseismic faulting episode was accompanied by low-frequency earthquakes located at the downdip edge of the slow slipping area (red points; from W. B. Frank et al. (2014)). Depth contours of the subducting slab are shown as dashed lines (Hayes et al., 2018). The upper left inset shows three distant GNSS stations (green triangles) used to isolate the common mode noise of the position time series (see Materials and Methods for details). The black line indicates the location of the vertical profile in (b). (b) Dislocation model of the Guerrero subduction interface. We used a 2D elastic dislocation geometry (purple line) with the same spatial extent as the slow slip source region (yellow patch, same as in A) to map the surface displacements recorded at the three analyzed GNSS stations onto the subduction interface via elastic Green's functions ($GF_{coeff}$). (c) GNSS surface displacements at the three GNSS stations in (a). (d) Projected geodetic fault slip on the Guerrero subduction interface from the three analyzed GNSS stations in A using the elastic dislocation model in (b). The coherence of the projected fault slip across all three stations justifies our choice of dislocation geometry for both loading (backslip) and release (slip).</figcaption></figure>

Tropospheric delays are estimated using the Vienna Mapping Function every 30 s, and we use the International Reference Ionospheric model to correct for ionospheric effects. Ocean loading effects are corrected using the oceanic tide model FES2014b (Spiridonov & Vinogradova, 2020). The position time series are converted to displacement time series in a fixed North American plate reference frame to remove any displacement due to rigid plate rotation. We remove outliers and adapt our observations to the 2D geometry of our study to maximize the geodetic signature of tectonic motion by projecting the horizontal displacement time series along the direction of tectonic convergence (25° from North). To correct the position time series for non-tectonic signals such as seasonal deformation and other regional common modes, we subtract the stack of the detrended GNSS time series from three distant stations (INEG, TAMP, and UXAL, whose locations are shown in Figure S1a in Supporting Information S1) unaffected by tectonic deformation in the subduction zone (Márquez-Azúa & DeMets, 2003). The stacked common mode signal, as well as the corrected GNSS time series, are represented in Figure S1 in Supporting Information S1.

Subdaily GNSS positioning is often avoided due to a significant increase in the noise level, resulting from positions that are estimated from fewer samples and additional noise sources such as multipathing. This explains why the subdaily dynamics of slow slip, which intermittent low-frequency earthquake activity suggests (W. B. Frank et al., 2014), have yet to be observed even with subdaily GNSS solutions (Itoh
et al., 2022), which require significant spatiotemporal smoothing that likely hides the fine details of slow slip dynamics. The high temporal resolution of low-frequency earthquake activity is thus an opportunity to explore subdaily positioning, as the seismic catalog provides an independent verification of our 6-hr displacement time series.

To keep our model of fault slip simple, we estimate from both the low-frequency earthquake catalog and the surface geodetic records of slow slip the same physical quantity: fault slip at the subduction interface. For the geodetic estimate of fault slip, we project the surface displacements recorded by GNSS stations to slip on the fault interface, as shown in Figure 1, with a 2D elastic dislocation model (Okada, 1985). We choose our model geometry to match the extent of the M7.5 slow slip event as shown in Figure 1b (Radiguet et al., 2011); further details on our static dislocation model are discussed in Appendix A. This dislocation model collapses the geodetic records at the three analyzed surface stations into one fault slip time series (Figure 1d), justifying our choice of model geometry. We observe that 4 cm of surface displacement (Figure 1c) translates to about 10 cm of fault slip once projected onto the subduction interface. We note that the coherence of the fault slip time series before and after the 2006 slow slip event suggests this source model is also reasonable for tectonic loading. This simple geometry of the Guerrero slow slip source region allows us to map the surface signature of slow slip at the sampling rate of the GNSS positioning onto the plate boundary at depth.

## 2.2. Estimating the Seismic Fault Slip of Low-Frequency Earthquakes

With the geodetic record projected on to the fault interface, we now need the analogous seismic estimate of fault slip. Given their close relationship to slow slip (W. B. Frank, 2016), we define the seismic slip generated by low-frequency earthquakes as the seismic analog to the observed geodetic fault slip. We note that given significant difference in observed seismic and geodetic moment during slow slip (W. B. Frank & Brodsky, 2019), we expect our estimates of seismic fault slip to be orders of magnitude smaller than the geodetic estimates described above. Due to the low signal-to-noise nature of low-frequency earthquakes (Katsumata & Kamaya, 2003), it is difficult to robustly constrain the source parameters, including the amount of slip, of individual low-frequency earthquakes (Thomas et al., 2016). We therefore use a low-frequency catalog constructed from the continuous seismic record from January 2005 to April 2007 (W. B. Frank et al., 2014) to make time-averaged estimates of low-frequency earthquake source parameters that match the sampling rate of the GNSS positions (24 and 6 hr). This catalog classifies all low-frequency earthquakes into spatially distinct sources with each source emitting many repeating seismic events. Figure S2 in Supporting Information S1 shows that low-frequency earthquake sources are not present throughout the entire slow slip source region and the density of low-frequency earthquake sources increases with distance from the trench. This updip deficit of sources is not a detection artifact (W. B. Frank et al., 2014), and we suggest that the rheological conditions necessary to generate low-frequency earthquakes are not uniformly present on the subduction interface. We focus on the 58 sources most responsive to slow slip (W. B. Frank et al., 2018) that are within or close to the slow slip source region (red points in Figures 1a and 1b; sources contained within the yellow patch in Figure S2 in Supporting Information S1). Figure 2a shows the recurrence interval between consecutive events among these low-frequency earthquake sources, with intermittent bursts of activity indicating concomitant slow slip (W. B. Frank et al., 2016). Any accurate estimate of low-frequency earthquake fault slip must capture this episodic evolution that is characteristic of low-frequency seismicity (W. B. Frank et al., 2014; Shelly, 2017; Kato & Nakagawa, 2020; Poiata et al., 2021).

The definition of the seismic moment is $M_0 = \mu Ad$, where $\mu$ (assumed $\mu = 40$ GPa) is the shear modulus and $d$ is the seismic fault slip over the rupture area A. To estimate the seismic fault slip of low-frequency earthquakes $d_{LFE}(t)$ we must first estimate the seismic moment $M_0(t)$ and slipping area of low-frequency earthquakes $A^S(t)$. The seismic displacement amplitude is directly proportional to the moment rate (Aki & Richards, 2002). We thus use the low-frequency earthquake seismic amplitudes, estimated as the root-mean-square displacement amplitude $\bar{u}$ during every low-frequency earthquake S-wave (W. B. Frank & Brodsky, 2019); Figure S3a in Supporting Information S1 shows how the median daily and subdaily seismic amplitude evolves through time. We convert the time-averaged amplitudes $\bar{u}_{LFE}(t)$ into the low-frequency earthquake median moment rate:

$$ \dot{M}_0(t) = 4\pi r \rho \beta^3 \cdot \bar{u}_{LFE}(t), \quad (1) $$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Time history of low-frequency earthquake fault slip. (a) Recurrence time of consecutive low-frequency earthquakes. The black points represent the time difference between two consecutive seismic events from the low-frequency earthquake catalog of W. B. Frank et al. (2014). (b) Daily and subdailly low-frequency earthquake slip. We used three independent observables to estimate the low-frequency earthquake fault slip 24hr- (blue) and 6h- (orange) sampled (see Materials and Methods and Figure S3 in Supporting Information S1 for details): the seismic displacement amplitude, the number of low-frequency events, and the number of active sources. Bursts of slip associated with recurring M6.4 slow slip events are visible prior to the strong increase of activity during the large 2006 slow slip event. (c) Integrated increment of the low-frequency earthquake slip. We use the integrated (or cumulative) increment of seismic fault slip as the driving observable of our fault slip model. The 6-hr total fault slip produces about 97% of the 24-hr total fault slip over the study duration; in other words, the two cumulative seismic slip time series are similar with a difference of cumulative fault slip that is less than 3%. The similarity of the 24- and 6-hr cumulative slip time series confirms the consistency of our method to estimate the seismic slip.</figcaption></figure>

where *r* is the distance between the source and the seismic station, *ρ* is the density of the medium, and *β* is the shear wave velocity; aside from *r* = 40 km that is related to the source-receiver geometry (Figure 1; the seismic stations that recorded low-frequency earthquake activity were deployed along the same along-dip profile as the GNSS stations we use here to capture slow slip), we assume typical values of *ρ* = 2850 kg/m³ and *β* = 3750 m/s. We then need to estimate the total slip duration of the low-frequency earthquakes at every time *t* to integrate the seismic moment rate into seismic moment. The slip duration for a given time step is defined as $N_{LFE} * \Delta t$, with $\Delta t$ the average source duration of a low-frequency earthquake and $N_{LFE}$ the number of low-frequency earthquakes; we assume a $\Delta t$ of 0.5 s per low-frequency earthquake (Bostock et al., 2015; Thomas et al., 2016). We can then estimate the time-averaged low-frequency earthquakes moment $M_0(t)$ (Figure S3b in Supporting Information S1) by integrating the time-averaged seismic moment rate $\dot{M}_0(t)$ over the slip duration as:

$$M_0(t) = \dot{M}_0(t)N_{LFE}(t)\Delta t. \quad (2)$$

To estimate the seismic fault slip of low-frequency earthquakes, we must also evaluate the effective seismic rupture area $A^S(t)$ of the low-frequency earthquakes in a given time step *t*. Past estimates of an individual low-frequency earthquake source area are on the order of several hundreds of square meters to a kilometer squared (Bostock et al., 2015; Chestler & Creager, 2017; Thomas et al., 2016). We assume here a constant source area of all low-frequency earthquakes given their relatively moment-independent source duration (Bostock et al., 2017; Farge et al., 2020). For each time step, we consider a low-frequency earthquake source active if it generated any events during the time step. We can thus approximate the effective slipping area over one time step as the sum of all of the active low-frequency earthquake sources, considering that each low-frequency earthquake source has an area of 1 km². Practically, we compute this effective slipping area (Figure S3c in Supporting Information S1) as the sum of the total number of active low-frequency earthquake sources, $A^S(t) = N_{source}(t) \cdot 1$ km². At 6 hr, we observed that the source area $A^S(t)$ varied more significantly than the 24-hr area time series as expected, but occasionally became significantly smaller than the 24-hr source area time series (Figure S3c in Supporting Information S1), leading to excessive estimates of fault slip. We therefore preferred to use the 24-hr source area time series resampled every 6 hr to compute our 6-hr fault slip time series.

We then estimate the time-averaged low-frequency earthquake seismic slip at sampling rates of 24 and 6 hr, as illustrated in Figure 2b:

$$d_{LFE}(t) = \frac{M_0(t)}{\mu A^S(t)}. \quad (3)$$

This time history of low-frequency earthquake fault slip $d_{LFE}(t)$ represents the total slip across all active low-frequency earthquake sources in a given time step. The daily fault slip increments shown in Figure 2b are consistent with past (albeit indirect) measures of low-frequency earthquake slip on the order of microns (Thomas et al., 2016, 2018). The similarity of the 24-hr and 6-hr cumulative slip time series seen in Figure 2c confirms the consistency of our method to estimate the seismic slip. We now describe how we use the integrated increment of low-frequency earthquake fault slip $\dot{d}_{LFE}$ at both daily and subdailly sampling rates, represented in Figure 2c, to match the time evolution of the geodetic fault slip when building our fault slip model.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Linear fault slip models for two different loading velocities. The linear fault slip models at both 24-hr and 6-hr sampling rates are represented by the black and gray lines respectively, for two different loading velocities. The colored time series correspond to the geodetic fault slip from the three GNSS stations used in this study. (a) The two linear fault slip models (see Equation 4) using a loading velocity of 6.4 cm/yr (the convergence rate of the Cocos plate) don't match the geodetic fault slip at both 24-hr and 6-hr sampling rates. (b) To obtain well-fitted models, a physically non-reasonable loading velocity of 15 cm/yr (more than 200% of the long-term convergence rate) is necessary.</figcaption></figure>

### 3. Model of Fault Slip Driven by Low-Frequency Earthquakes

We use this time series of low-frequency earthquake slip to model the evolution of fault slip throughout the slow slip cycle. We start by defining our slip model $d$ as the simple competition between loading (backslip) $d_L$ and the seismic slip associated with low-frequency earthquakes $d_S$:

$$d(t) = d_L(t) - \zeta d_S(t) \quad (4)$$

where $d$ is the modeled fault slip, $t$ is time, and $\zeta$ is the unitless ratio between total (geodetic) and seismic slip. We define the backslip $d_L(t)$ and seismic slip $d_S(t)$ terms as:

$$d_L = V_L t + d_0 \qquad d_S = \int_0^t \dot{d}_{LFE}(\tau) d\tau \quad (5)$$

where the loading backslip $d_L$ is driven by a long-term loading velocity $V_L$ with some zero-time offset $d_0$, and the seismic slip $d_S$ is expressed as the integrated increment of low-frequency earthquake slip $\dot{d}_{LFE}$ over either a 24- or 6-hr time period (Figure 2c). Considering the time derivative of Equation 4, the difference between the first and the second term corresponds to the effective fault slip rate across our modeled fault dislocation (Figure 1b), where positive slip is tectonic loading or backslip, and negative slip is tectonic release or slip. $V_L$ functions as the maximum possible loading (backslip) rate when there are no low-frequency earthquakes and thus no modeled aseismic fault slip. Fitting our model to the geodetic observations (Figure 1d), we need a loading velocity $V_L = 15$ cm/yr for our linear model to reproduce the geodetic observations of fault slip over the full 2.5 years observational period as shown in Figure 3b. We consider such a linear model non-physical because this loading velocity is more than two times larger than the long-term convergence rate of the Cocos plate, of 6.4 cm/yr (DeMets et al., 2010). If we impose a long-term loading rate of 6.4 cm/yr, our model with a linear relationship between geodetic and seismic slip does not fit the observed fault slip record (Figure 3a). To achieve both a reasonable loading rate and a model that can fit the observations, we now complexify our model by making the relationship between geodetic and seismic fault slip non-linear.

Modifying the second term on the right-hand side of Equation 4, we raise the low-frequency earthquake slip $d_S(t)$ to the exponent $\alpha$ to attempt to capture any potential non-linearity in the geodetic to seismic slip relationship:

$$d_S(t) = \int_0^t C \dot{d}_{LFE}(\tau)^\alpha d\tau. \quad (6)$$

where $C$ is simply a scalar that compensates for the physical units that change as a function of $\alpha$ so that $d_S(t)$ maintains units of slip (m) and $\zeta$ remains a unitless ratio of geodetic to seismic fault slip. Introducing this $\alpha$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Distribution of best-fit non-linear models that can explain the fault slip time history throughout the slow slip cycle. Sets of loading velocities $V_L$ and power-law exponents $\alpha$ associated with the 5% best models for both the 24-hr and 6-hr fault slip models are shown by the color points. Only the overlap between the 24-hr and 6-hr models (black points) was considered when looking for the best model parameters (see text for details). The yellow star indicates the loading rate and exponent used for the two models presented in Figure 5a ($V_L$ = 6.4 cm/yr and $\alpha$ = 2). Gray shaded area illustrates the fault coupling as a percentage of the long-term tectonic convergence rate of 6.4 cm/yr (DeMets et al., 2010). The full misfit space is shown in Figure S4 in Supporting Information S1.</figcaption></figure>

parameter modifies our linear equation of fault slip (Equation 4) to a power-law model. There are now two parameters of interest in our model: the power-law exponent of seismic fault slip $\alpha$, and the loading rate $V_L$. We explore the model's parameters space with a grid search over various loading rates $V_L$ (0–13 cm/yr) and exponents $\alpha$ (0–5) for both 24- and 6-hr models; scanning a large range of loading velocities allows us to broadly image the shape of the parameter space, although we still consider loading rates larger than the long-term convergence of >6.4 cm/yr to be non-physical. For every set of $V_L$ and $\alpha$, we then estimate the aseismic to seismic slip ratio $\zeta$ (as well as the intercept $d_0$) via least squares. We finally calculate the L2-norm between each modeled time history of slip $d(t)$ and the projected GNSS time series (Figure 1d) as the model misfit to identify a set of best-fitting models shown in Figure 4; the full misfit space is shown in Figure S4 in Supporting Information S1.

To be consistent when comparing the 24- and 6-hr non-linear fault slip models, we focus on the overlapping area of the 5% best models shown in Figure 4. If we consider only models with realistic loading rates (<6.4 cm/yr) in this overlap (black points in Figure 4), the best-fit 24- and 6-hr fault slip models have a loading velocity $V_L$ of about 6.4 cm/yr and a power-law exponent $\alpha$ of 2 (yellow cross in Figure 4). We show in Figure 5 that these best-fitting fault slip models reproduce the GNSS position solutions at two different time scales throughout the slow slip cycle, leveraging only seismological observations of low-frequency earthquake activity. We observe modeled slip rates that vary over a broad range of time scales at both 24-hr and 6-hr sampling rates. The apparent smooth and monotonic slip rate of 10 cm/yr seen on the GNSS time series during the 2006 M7.5 slow slip event can be explained as the result of many episodic transient slip events separated by rapid periods of loading (W. B. Frank et al., 2018), as suggested by the

seismic fault slip time series in Figure 2b. A zoom on a M6.4 small slow slip event in November 2005 (Figure 5b) demonstrates that our non-linear models informed by low-frequency earthquake activity track small slow slip rates at both 24- and 6-hr sampling rates.

# 4. Discussion

## 4.1. Slow Slip Is the Competition Between Subdaily Dynamics and a Fully Coupled Fault

The competition between tectonic loading and release controls our model of fault slip (Equation 4) within the slow slip source region. We observe in Figure 4 that the best models consistently require a loading rate that is close to or above the tectonic convergence rate of 6.4 cm/yr (DeMets et al., 2010). A loading rate that is 100% of the convergence rate implies that the assumed source region of slow slip is fully locked (or coupled), with no seismic or aseismic faulting. If we normalize the time derivative of Equation 4 by $V_L$, we can compute the instantaneous fault coupling at every sample of our model, as shown in Figures 5c and 5d. In Figure 5c, positive fault coupling corresponds to loading, or backslip on the interface, while negative coupling (and negative slip rate values) represents release or slip on the fault; fault coupling <−1 means the fault slips faster than the convergence rate. We observe that tectonic release within the modeled slow slip source region is sporadic. Figure S5 in Supporting Information S1 shows that tectonic release occurs in rapid bursts that can reach slip rates of up to 1000 cm/yr, two orders of magnitude faster than the apparent long-term fault slip rate of the M7.5 slow slip event (Figure 5a).

Our models exhibit rapid periods of release before and after the 2006 M7.5 slow slip in both the 24- and 6-hr sampled data. This episodic release motion corresponds to the M6.4 slow slip events, which are geodetically observable only after stacking the position time series of multiple events (W. B. Frank et al., 2015). Zooming into one of these events in Figure 5b, we highlight that the 6-hr modeled fault slip is not uniformly distributed in a 24-hr time window. The consistency in cumulative slip between the 24-hr and 6-hr low-frequency earthquake fault slip (Figure 2c) suggests this subdaily variability is real. The instantaneous fault coupling during the same time period (Figure 5d) demonstrates that peak slip rates at 6 hr are nearly two times faster than at 24 hr; these
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Daily and subdaily geodetic fault slip models driven by low-frequency earthquake activity. (a) Modeled fault slip compared to geodetic fault slip (projected from GNSS data); we include an offset for the 24-hr sampled GNSS data and modeled fault slip for clarity. The modeled fault slip fits across a range of time scales the geodetic fault slip of the two plotted GNSS time series: IGUA (gray time series) and MEZC (blue for 24-hr-sampled data, orange for 6-hr-sampled data); the model was also simultaneously fitted to data from COYU (Figure 1) but its time series are not shown here for clarity. We show here the best-fit models with $V_L = 6.4$ cm/yr (regional convergence rate is 6.4 cm/yr) and $\alpha = 2$. The shaded gray box highlights the time period in (b). (b) Zoom on a M6.4 slow slip event that lasted several days in November 2005. Using only low-frequency earthquake activity, both our 24-hr- and 6-hr-sampled models (same colors as A) capture the slow slip event delineated by the black lines. For clarity, we show only the MEZC positions and an additional smoothed 6-hr-sampled GNSS time series (solid orange; dotted orange is the original 6-hr time series) prior to and after the slow slip event. (c) Instantaneous fault coupling across the slow slip source region at 24 and 6 hr (see text for details). Slip rates $<-6.4$ cm/yr are clipped for clarity; the full range of coupling and slip rate is shown in Figure S5 in Supporting Information S1, which highlights that the 6-hr slip rate is both more intermittent and faster than the 24-hr slip rate. (d) Instantaneous fault coupling (not clipped) during the same time period as (b). We observe intermittent and faster slip rates for the 6-hr model compared to the 24-hr model during the M6.4 slow slip event.</figcaption></figure>

rapid slip rates are likely still underestimated as they are spatially averaged over the entire slow slip source region. Although the long-period trend of the two models are very similar and reproduce well the geodetic records, small scale differences are observed between the 24-hr and 6-hr fault slip solutions.

This observation implies that the complex slow slip dynamics seen here at subdaily time scales, which are otherwise hidden within the noise of the geodetic record, dominate the tectonic release during slow slip. Our initial hypothesis that low-frequency earthquakes can capture the subdaily slow slip dynamics hidden in the geodetic record is validated by the fact that the rapidly evolving low-frequency earthquake activity reproduces the multi-year GNSS position time series. Past work has shown that subdaily variations in slow slip moment rate can occur (Hawthorne & Rubin, 2013). Our results generalize this observation to be fundamental to the slow slip rupture process and show how subdaily dynamics govern the overall moment release of slow slip. We speculate that future work on reducing the noise level of subdaily GNSS positioning will yield further details of slow slip dynamics at short time scales.

Between the observed rapid episodes of slow slip, we observe in Figure 5c that the plate interface is nearly fully coupled during the majority of the analyzed 2.5 years with an instantaneous fault coupling close to 1 (or a backslip rate of 6.4 cm/yr). Estimates of plate coupling are often focused on seismogenic faults, with the aim to quantify how much of the far field strain loading is being stored during the interseismic phase of the earthquake,
a rough estimate of the energy budget of a future earthquake (Avouac, 2015). These long-term measures of plate coupling are estimated by fitting linear trends over months or years of surface motion, glossing over any potential short-term dynamics present within the geodetic data. Our model allows us to go beyond this limitation to quantify subdaily fluctuations of fault coupling. If we estimate the linear trend of our modeled fault slip time series over a given time period (Figure 5a), we have the opportunity to estimate an analogous “long-term” estimate of plate coupling. If we do this over the time period prior to the 2006 M7.5 slow slip event, we obtain a plate coupling of 60% and 71% respectively for the 24- and 6-hr estimates. Both of these estimates are significantly lower than the near full coupling observed at daily and subdaily time scales (Figure 5c), but are compatible with previous estimates of plate coupling between major slow slip events (Radiguet et al., 2012). If we assume that plate coupling is purely a measure of temporal intermittence of slip (W. B. Frank, 2016), this implies that a partial coupling (≤1) directly tells just how often a given fault is slipping, being fully locked the rest of the time. While this is only one potential explanation of partial plate coupling, our results support the idea that estimates of plate coupling depend on the time scales over which they are calculated (W. B. Frank, 2016; Jolivet & Frank, 2020; Maubant et al., 2022). This bears out in our observations where we see that the instantaneous coupling at 6 hr is consistently higher than at 24 hr as shown in Figures 5c and 5d. This further implies that observations of long-term continuous partial coupling are likely non-physical and are an artifact of estimating surface velocities over long time scales.

## 4.2. Slow Slip Saturates the Low-Frequency Earthquake Source Region

We observe in Figure 3 that a linear model (Equation 4) cannot both have a physically reasonable loading rate and match the geodetic fault slip. In our non-linear version of this model, tectonic release is controlled by the observed low-frequency earthquake activity raised to some exponent $\alpha$. A power-law model with a $\alpha > 1$ means that during periods of intense seismic fault slip (or high low-frequency earthquake activity), the model must disproportionately produce more aseismic slip than during periods of low seismic fault slip, implying the ratio of seismic to aseismic slip is not constant. We observe that this power-law exponent $\alpha$ must be about 2 for our models to fit the GNSS time series for the full period of study (Figures 4 and 5a).

We explore this non-linearity between aseismic and seismic fault slip through the lens of our linear model (Equation 4) fitted to different phases of the slow slip cycle. We adjust our model to the GNSS data (a) during the Inter-SSE period, or the time period prior to the M7.5 slow slip event and (b) during the Co-SSE period, or the time period during the 2006 M7.5 slow slip event. Both geodetic time periods can separately be well-fitted using a linear model with a reasonable loading rate of 6.4 cm/yr (Figure 6); only the 24-hr linear models are shown for clarity, but the 6-hr model exhibits a similar behavior (see caption for details).

In this piece-wise model, the principal difference is the geodetic to seismic slip ratio $\zeta$ that changes by a factor of 4 between the two phases (14 for “Inter-SSE” and 63 for “Co-SSE” period). To physically interpret $\zeta$, we define the ratio of geodetic to seismic slip as follows:

$$ \zeta = \frac{d_G(t)}{d_S(t)} = \frac{M_0^G A^S}{M_0^S A^G}, \qquad (7) $$

where $M_0^G$ and $A^G$ are the geodetic (total) moment and slipping area associated with slow slip respectively, and $M_0^S$ and $A^S$ are the seismic moment and slipping area of low-frequency earthquakes. We note that the geodetically-estimated moment and rupture area overwhelmingly represent the aseismic rupture associated with slow slip, with seismic low-frequency earthquakes representing a negligible fraction of the total moment (W. B. Frank & Brodsky, 2019). We observe that during both the M6.4 slow slip events and the M7.5 slow slip event that respectively dominate the Inter- and Co-SSE periods: (a) the geodetic to seismic moment ratio $M_0^G/M_0^S$ is relatively constant over daily and subdaily timescales (Figure S3b in Supporting Information S1), and (b) the low-frequency earthquake rupture area is relatively constant (Figure S3c in Supporting Information S1), consistent with a source region that is structurally limited (Figure S2 in Supporting Information S1). Consequently, we consider that changes in $\zeta$ during the Inter- and Co-SSE periods are principally due to changes in the geodetic slipping area (see Appendix B for more details). The observed change in $\zeta$ thus implies that the slipping area of the M7.5 slow slip event in 2006 should be four times larger than the slipping area of the M6.4 slow slip events. This factor of four in slow slip rupture area is consistent with geodetic observations, which suggest the M7.5 event is hosted on a
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Linear models of different phases of the slow slip cycle and low-frequency earthquake source region saturation. (bottom) Best-fit linear models (Equation 4) when independently adjusted to two phases of the slow slip cycle. The red box ("Inter-SSE") corresponds to the best-fit model when adjusting to the period prior to the M7.5 slow slip event (red line, with a $\zeta$ of 14), from January 2005 to April 2007. The green box ("Co-SSE") represents the best-fit model when adjusting to the duration of the M7.5 2006 slow slip event (green line with a $\zeta$ of 63). Only the 24-hr sampled fault slip time series are shown here for clarity; similar results were found for the 6-hr sampling rate ($\zeta$ of 18 and 60 for "Inter-SSE" and "Co-SSE" respectively). (Top left) Schematic of seismic sources and aseismic rupture area for different slow slip magnitudes. The three cross-sections represent the subduction zone interface, and schematically illustrate the evolution of aseismic slip (yellow patches) and activated low-frequency earthquake sources (red points) for a daily transient event (1), a small M6.4 aseismic slip (2), and the large M7.5 slow slip event (3). (Top right) Schematic relationship between the aseismic to seismic slip ratio. The red and green tangent dashed lines represent the seismic to aseismic slip ratio of the "Inter-SSE" and "Co-SSE" slip models respectively, while a non-linear quadratic model with $\alpha \approx 2$ is required to fit the whole slow slip cycle (black line). The gray curve illustrates the saturation of low-frequency earthquake rupture area as aseismic moment increase toward M7.5 slow slip events.</figcaption></figure>

slipping patch of 300 × 150 km, while the M6.4 events rupture an 150 × 100 km patch (W. B. Frank et al., 2015). We suggest that only a model with a non-linear relationship between geodetic and seismic fault slip can reproduce the entire 2.5-year geodetic record because it is able to account for the observed transient changes in $\zeta$ to balance the rich dynamics of slow slip, from small, daily slip transients to major slow slip events.

Slow slip events are known to initiate at the deepest extent of their source region, where frequent low-frequency earthquake and tremor activity suggest near-continuous loading from further downdip (Obara et al., 2010; Wech & Creager, 2011; W. Frank et al., 2015). Once their rupture nucleates, slow slip can then only increase in moment by extending updip. The different sized events thus exhibit different updip extents to accommodate their different moments, schematically shown in Figure 6: (a) the frequent daily M5 transients are likely associated with slip downdip, along the flat part of the interface (W. B. Frank & Brodsky, 2019), driving a few low-frequency earthquake sources to fail (red circles); (b) the small M6.4 slow slip events rupture a larger area that extends further updip, activating more low-frequency earthquake families; (c) the slipping area of the large M7.5 slow slip event saturates the low-frequency earthquake source region, extending much further updip. Figure 6 illustrates that these slow slip events are systematically accompanied by low-frequency earthquakes in the same localized source region regardless of the spatial extent of the driving aseismic rupture.

Our results confirm that when present, low-frequency earthquakes act as an in situ gauge of neighboring aseismic slip, analogous to repeating earthquakes (Uchida & Bürgmann, 2019). As aseismic slip grows in moment, the low-frequency earthquake sources continue to respond to the increased aseismic slip rate, producing more
and higher amplitude low-frequency earthquakes (W. B. Frank & Brodsky, 2019), yet the spatial extent of the low-frequency earthquake source region cannot grow to match the aseismic rupture that extends updip (Figure S2 in Supporting Information S1). We observe such a saturation of low-frequency earthquake rupture area in Figure S3 in Supporting Information S1: while the low-frequency earthquake moment is larger during the M7.5 slow slip event than during the M6.4 slow slip events (Figure S3b in Supporting Information S1), the seismic rupture area of low-frequency earthquakes is the same during both sizes of slow slip events (Figure S3c in Supporting Information S1). Given the limited source area of low-frequency earthquakes that can be activated, a power-law exponent $\alpha$ is needed to accommodate wider range of slow slip rupture area, and thus moment (Figures 4 and 6). The schematic panel of Figure 6 (top right) illustrates how two linear models that individually fit well two phases of the slow slip cycle ("Inter-SSE" and "Co-SSE"), can be reproduced by a single quadratic non-linear model for the whole slow slip cycle. In other words, low-frequency earthquakes are only a powerful monitor of fault slip until the aseismic rupture grows beyond their immediate proximity.

### 4.3. Low-Frequency Earthquake Activity Is Incidental to Slow Slip

There are clear examples of the spatiotemporal correlation of slow slip and low-frequency earthquakes and tectonic tremor (Bartlow et al., 2011; Nishimura et al., 2013), yet this spatial overlap of aseismic and seismic source regions is often incomplete at different subduction zones. Notably in the Guerrero segment studied here, only a fraction of the slow slip source region hosts low-frequency earthquakes (Figure 1). Leveraging these natural laboratory controls, we observe that low-frequency earthquakes are only well correlated with neighboring slow slip within this spatial overlap, with decreasing sensitivity to aseismic slip moving away from the low-frequency earthquake source region. We thus suggest that low-frequency earthquake activity is incidental to slow slip. Because we observe slow slip's rupture area to first saturate and then extend beyond the low-frequency earthquake source region, it follows that the fault rheology responsible for slow slip must not necessarily account for low-frequency earthquakes. This implies that the physical mechanisms responsible for low-frequency earthquakes do not necessarily play a significant role in the aseismic rupture process. Recent work has reported a range of potential geological signatures of slow slip with mixed modes of shear failure (Barnes et al., 2020; Condit et al., 2022; Platt et al., 2018; Ujiie et al., 2018), with ductile mechanisms attributed to slow slip and brittle faulting linked to low-frequency earthquakes. If slow slip can occur without accompanying seismic faulting, we speculate that these rich images of mixed-mode deformation in subduction zones likely only represent a fraction of the geological possibilities that host slow slip. Slow slip that need not be accompanied by low-frequency earthquakes (or tectonic tremor) would also explain observations of slow slip that are associated with swarms of regular earthquakes, such as the slow slip events in the Boso channel in Japan (Hirose et al., 2014) or shallow slow slip in the Hikurangi subduction zone (Wallace et al., 2018). We speculate that the type of seismic event driven to failure by slow slip likely depends on the rheology that varies as a function of depth, as well as other varying fault conditions such as pore fluid pressure that are not present throughout the slow slip source region. We remark that this simplified conceptual model of slow slip, which does not depend on any interplay with seismic faulting, supports the observation that slow slip is a general feature of plate boundary faulting, occurring across a wide range of tectonic contexts (Jolivet & Frank, 2020).

## 5. Conclusion

Our results show that the long-term geodetic record of slow slip is driven by subdaily variations in aseismic slip rate that are highlighted by low-frequency earthquake activity. We observe that the ratio between the driving aseismic slip and the seismic slip of low-frequency earthquakes evolves throughout the entire slow slip cycle. This implies that low-frequency earthquakes are an incidental consequence of slow slip that dominates tectonic release downdip of the seismogenic plate boundary, and low-frequency earthquakes thus do not play a significant role in mediating aseismic slip. These observational constraints are only made possible by a seismic catalog with high temporal resolution, a necessary geophysical guide to extract the subdaily dynamics of slow slip from the geodetic record. We speculate that future constraints on slow slip dynamics at short time scales, especially those that combine the resolving power of seismic and geodetic observations, will shed further light on the physical mechanisms responsible for slow slip and the role that slow fault slip plays in the seismic cycle.
## Appendix A: Geodetic Dislocation Model of Fault Slip on the Plate Interface

Slow slip in Guerrero typically extends updip of the low-frequency earthquake source region, between 60 and 180 km away from the trench (Radiguet et al., 2016), with a majority of slip updip of the flat section of the subduction interface (Figure 1b). We thus model the slow slip source region as a sum of two fault dislocations in an elastic half-space that reproduces the geometry of the 2006 slow slip event (Radiguet et al., 2011): one along the dipping section and the other just downdip where the subduction interface goes flat (Figure 1b). The Green's function coefficient ($GF_{coeff}$) derived from unitary slip on a patch representing the spatial extent of the 2006 slow slip event relates the GNSS surface displacements onto the modeled dislocation. We assume that any slip allocated to this dislocation geometry is uniformly distributed in space, representing a spatial average of the slip on the subduction interface. We then simply divide the surface displacement time series by the $GF_{coeff}$ to project the GNSS position time series to fault slip on the interface. We can verify whether the assumed dislocation geometry is appropriate by examining the coherence of the resulting fault slip time series across the three stations (Figure 1d); because all three stations should produce the same fault slip time history after projection, the three time series should collapse together if the assumed geometry is reasonable. We note that the observed coherence of the fault slip time series before and after the 2006 slow slip event suggests this source model is also reasonable for tectonic loading or backslip. We compared whether one $GF_{coeff}$ per GNSS station worked better than using a single $GF_{coeff}$ associated with the reference station MEZC and fitting the two remaining stations to the MEZC fault slip time series. The latter option produced slightly more coherent fault slip time series than the former (Figure 1d), and the fitted coefficient for COYU and IGUA were within 15% of their original computed $GF_{coeff}$ using the dislocation model. We thus opted for the latter option using MEZC as a reference station. We choose MEZC as a reference to evaluate fault slip because this GNSS station is located just above the low-frequency earthquake sources and on the edge of the slipping area associated to the M7.5 2006 slow slip. We acknowledge that using a simplified 2D homogeneous model might induce bias in the fault slip estimations, but we focus here on the temporal evolution of slow slip rather than the absolute quantification of slip.

## Appendix B: Constraints on the Geodetic to Seismic Slip Ratio $\zeta$

The variation between the linear fault slip models fitted to the two different phases of the slow slip cycle is dominated by the change in the geodetic to seismic slip ratio $\zeta$. We can define $\zeta$ as:

$$ \zeta = \frac{d_G}{d_S}, \qquad (B1) $$

which can also be rewritten and simplified as:

$$ \zeta = \frac{M_0^G / A^G}{M_0^S / A^S} = \frac{M_0^G A^S}{M_0^S A^G}, \qquad (B2) $$

where $M_0^G$ and $A^G$ are the geodetic moment and slipping area respectively, and $M_0^S$ and $A^S$ are the seismic moment and slipping area of low-frequency earthquakes. We observe the time-averaged geodetic to seismic moment ratio is relatively constant (at daily and subdaily time scales) during different sized slow slip events (Figure S3b in Supporting Information S1). While there are 58 low-frequency earthquake sources, the maximum number of activated sources at any given time is 52, giving a maximum seismic source area $A^S$ of 52 km² assuming each low-frequency earthquake source has a rupture area $A$ of 1 km² (Figure S3c in Supporting Information S1). This maximum seismic slipping area is reached during both the large M7.5 slow slip event and the smaller M6.4 slow slip events. We thus consider the active seismic source region $A^S$ is the same during both small and large slow slip events. We estimate $A^G$ from geodetic slip inversions to be 150 × 300 km² for the large M7.5 slow slip event and 100 × 150 km² for the smaller M6.4 slow slip events (W. B. Frank et al., 2015). We then constrain $\zeta$ by fitting our models to the geodetic data and respectively estimate 14 and 63 for the Inter- and Co-SSE time periods for the 24-hr data set (we found 18 and 60 for the corresponding 6-hr models). Our static dislocation model assumes a homogeneous and constant geodetic slipping area, while our interpretation supposes a geodetic slipping area that evolves with increasing moment. We evaluated whether this poses an issue by changing the geometry of our dislocation model to match the extent of the low-frequency earthquake source region, and we find that the $GF_{coeff}$ of MEZC and IGUA do not change significantly. This implies that regardless whether we consider the full extent of the slow slip region or only limit our dislocation geometry to the low-frequency earthquake source region, the
fault slip models would be the same, with both geometries still requiring a power-law to accommodate evolving
differences in slipping area throughout the slow slip cycle.

As a consistency check of our approach and the assumptions underlying our model, we try to recover the moment magnitude of both the large M7.5 and smaller M6.4 slow slip events using the fitted parameters from our models and geodetic constraints. We use the estimates of $A^S$, $M_0^S$ and $\zeta$ from this study, as well as previous geodetic estimates of rupture area $A^G$ of the M7.5 and M6.4 slow slip events (W. B. Frank et al., 2015). Using these estimated values with Equation B2, we estimate moment magnitudes of M7.4 and M5.8 respectively for the large M7.5 slow slip event and the M6.4 slow slip events. Our estimated moment magnitudes are close to the geodetic estimates, with differences likely driven by rough estimates of the geodetic rupture area $A^G$, especially for the small M6.4 slow slip events. Given the consistency between our magnitude estimates and those independently constrained by geodesy (W. B. Frank et al., 2015), we conclude that our assumptions supporting our estimates of seismic and geodetic fault slip are physically reasonable.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

The Guerrero low-frequency earthquake catalog used to estimate the seismic fault slip in the study is available from the Slow Earthquake database (http://www-solid.eps.s.u-tokyo.ac.jp/sloweq/). The raw GNSS position time series used to estimate geodetic fault slip in the study are available here: https://doi.org/10.5281/zenodo.7921935. The raw GNSS time series are available for each three GNSS stations we focus on in this study (COYU, IGUA and MECZ), as well as the three distant stations use to correct from common mode (INEG, TAMP and UXAL), at both 24-hr and 6-hr sampling rates.

## References

Aki, K., & Richards, P. G. (2002). Quantitative seismology.

Avouac, J.-P. (2015). From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. *Annual Review of Earth and Planetary Sciences*, 43(1), 233–271. https://doi.org/10.1146/annurev-earth-060614-105302

Barnes, P. M., Wallace, L. M., Saffer, D. M., Bell, R. E., Underwood, M. B., Fagereng, A., et al. (2020). Slow slip source characterized by lithological and geometric heterogeneity. *Science Advances*, 6(13), eaay3314. https://doi.org/10.1126/sciadv.aay3314

Bartlow, N. M., Miyazaki, S., Bradley, A. M., & Segall, P. (2011). Space-time correlation of slip and tremor during the 2009 Cascadia slow slip event: Cascadia slow slip and tremor migration. *Geophysical Research Letters*, 38(18), L18309. https://doi.org/10.1029/2011GL048714

Beroza, G. C., & Ide, S. (2011). Slow earthquakes and nonvolcanic tremor. *Annual Review of Earth and Planetary Sciences*, 39(1), 271–296. https://doi.org/10.1146/annurev-earth-040809-152531

Bertiger, W., Bar-Sever, Y., Dorsey, A., Haines, B., Harvey, N., Hemberger, D., et al. (2020). GipsyX/RTGx, a new tool set for space geodetic operations and research. *Advances in Space Research*, 66(3), 469–489. https://doi.org/10.1016/j.asr.2020.04.015

Bostock, M. G., Thomas, A. M., Rubin, A. M., & Christensen, N. I. (2017). On corner frequencies, attenuation, and low-frequency earthquakes. *Journal of Geophysical Research: Solid Earth*, 122(1), 543–557. https://doi.org/10.1002/2016JB013405

Bostock, M. G., Thomas, A. M., Savard, G., Chuang, L., & Rubin, A. M. (2015). Magnitudes and moment-duration scaling of low-frequency earthquakes beneath southern Vancouver Island. *Journal of Geophysical Research: Solid Earth*, 120(9), 6329–6350. https://doi.org/10.1002/2015JB012195

Bürgmann, R. (2018). The geophysics, geology and mechanics of slow fault slip. *Earth and Planetary Science Letters*, 495, 112–134. https://doi.org/10.1016/j.epsl.2018.04.062

Chestler, S. R., & Creager, K. C. (2017). A model for low-frequency earthquake slip. *Geochemistry, Geophysics, Geosystems*, 18(12), 4690–4708. https://doi.org/10.1002/2017GC007253

Condit, C. B., French, M. E., Hayles, J. A., Yeung, L. Y., Chin, E. J., & Lee, C. A. (2022). Rheology of metasedimentary rocks at the base of the subduction seismogenic zone. *Geochemistry, Geophysics, Geosystems*, 23(2), e2021GC010194. https://doi.org/10.1029/2021GC010194

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. *Geophysical Journal International*, 181(1), 1–80. https://doi.org/10.1111/j.1365-246X.2009.04491.x

Farge, G., Shapiro, N. M., & Frank, W. B. (2020). Moment-duration scaling of low-frequency earthquakes in Guerrero, Mexico. *Journal of Geophysical Research: Solid Earth*, 125(8), e2019JB019099. https://doi.org/10.1029/2019JB019099

Frank, W., Shapiro, N., Husker, A., Kostoglakov, V., Bhat, H., & Campillo, M. (2015). Along-fault pore-pressure evolution during a slow-slip event in Guerrero, Mexico. *Earth and Planetary Science Letters*, 413, 135–143. https://doi.org/10.1016/j.epsl.2014.12.051

Frank, W. B. (2016). Slow slip hidden in the noise: The intermittence of tectonic release. *Geophysical Research Letters*, 43(19), 10125–10133. https://doi.org/10.1002/2016GL069537

Frank, W. B., & Brodsky, E. E. (2019). Daily measurement of slow slip from low-frequency earthquakes is consistent with ordinary earthquake scaling. *Science Advances*, 5(10). https://doi.org/10.1126/sciadv.aaw9386

Frank, W. B., Radiguet, M., Rousset, B., Shapiro, N. M., Husker, A. L., Kostoglakov, V., et al. (2015). Uncovering the geodetic signature of silent slip through repeating earthquakes. *Geophysical Research Letters*, 42(8), 2774–2779. https://doi.org/10.1002/2015GL063685

## Acknowledgments

We are grateful to everyone who partic-
ipates in the maintenance of the GNSS
network in Guerrero, including but not
limited to Vladimir Kostoglakov, Jorge
Real Perez, Andrea Walpersdorf, and
José Antonio Santiago. We thank Axel
Periollat, Michael Floyd, and Thomas
Herring for discussions related to GNSS
processing, and Louise Maubant and
Enrico Milanese for discussions that
improved the manuscript. We finally
thank Brent Minchew for keeping our
units honest. C.M. and W.B.F. acknowl-
edge support from the National Aeronaut-
ics and Space Administration under Grant
80NSSC22K0459 issued through the
Earth Sciences Division of the Science
Mission Directorate. C.M., W.B.F.,
and P.P. acknowledge support from the
Transatlantic Research Partnership Grant
93, a program of the FACE Foundation
and the French Embassy. C.M., M.R., and
P.P. acknowledge support from the Euro-
pean Research Council under the Horizon
2020 Grant 802777-MONIFAULTS.
C.M. acknowledges support from Math-
Works; the opinions and views expressed
in this publication are the authors' and do
not represent those of MathWorks. M.R.
acknowledges support from the Agence
Nationale de la Recherche under Grant
ANR-21-CE49-0023-SSDYN.
Frank, W. B., Rousset, B., Lasserre, C., & Campillo, M. (2018). Revealing the cluster of slow transients behind a large slow slip event. *Science Advances*, 4(5). https://doi.org/10.1126/sciadv.aat0661

Frank, W. B., Shapiro, N. M., Husker, A. L., Kostogllov, V., Gusev, A. A., & Campillo, M. (2016). The evolving interaction of low-frequency earthquakes during transient slip. *Science Advances*, 2(4), e1501616. https://doi.org/10.1126/sciadv.1501616

Frank, W. B., Shapiro, N. M., Husker, A. L., Kostogllov, V., Romanenko, A., & Campillo, M. (2014). Using systematically characterized low-frequency earthquakes as a fault probe in Guerrero, Mexico. *Journal of Geophysical Research: Solid Earth*, 119(10), 7686–7700. https://doi.org/10.1002/2014JB011457

Fujita, M., Nishimura, T., & Miyazaki, S. (2019). Detection of small crustal deformation caused by slow slip events in southwest Japan using GNSS and tremor data. *Earth Planets and Space*, 71(1), 96. https://doi.org/10.1186/s40623-019-1075-x

Hawthorne, J. C., & Bartlow, N. M. (2018). Observing and modeling the spectrum of a slow slip event. *Journal of Geophysical Research: Solid Earth*, 123(5), 4243–4265. https://doi.org/10.1029/2017JB015124

Hawthorne, J. C., & Rubin, A. M. (2013). Short-time scale correlation between slow slip and tremor in Cascadia. *Journal of Geophysical Research: Solid Earth*, 118(3), 1316–1329. https://doi.org/10.1002/jgrb.50103

Hayes, G. P., Moore, G. L., Portner, D. E., Hearne, M., Flamme, H., Furtney, M., & Smoczyk, G. M. (2018). Slab2, a comprehensive subduction zone geometry model. *Science*, 362(6410), 58–61. https://doi.org/10.1126/science.aat4723

Hirose, H., Matsuzawa, T., Kimura, T., & Kimura, H. (2014). The Boso slow slip events in 2007 and 2011 as a driving process for the accompanying earthquake swarm. *Geophysical Research Letters*, 41(8), 2778–2785. https://doi.org/10.1002/2014GL059791

Husker, A., Frank, W. B., Gonzalez, G., Avila, L., Kostogllov, V., & Kazachkina, E. (2019). Characteristic tectonic tremor activity observed over multiple slow slip cycles in the Mexican Subduction Zone. *Journal of Geophysical Research: Solid Earth*, 124(1), 599–608. https://doi.org/10.1029/2018JB016517

Ide, S., Shelly, D. R., & Beroza, G. C. (2007). Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface. *Geophysical Research Letters*, 34(3), L03308. https://doi.org/10.1029/2006GL028890

Ito, Y., Obara, K., Shiomi, K., Sekine, S., & Hirose, H. (2007). Slow earthquakes coincident with episodic tremors and slow slip events. *Science*, 315(5811), 503–506. https://doi.org/10.1126/science.1134454

Itoh, Y., Aoki, Y., & Fukuda, J. (2022). Imaging evolution of Cascadia slow-slip event using high-rate GPS. *Scientific Reports*, 12(1), 7179. https://doi.org/10.1038/s41598-022-10957-8

Jolivet, R., & Frank, W. B. (2020). The transient and intermittent nature of slow slip. *AGU Advances*, 1(1). https://doi.org/10.1029/2019AV000126

Kaneko, L., Ide, S., & Nakano, M. (2018). Slow earthquakes in the microseism frequency band (0.1-1.0 Hz) off Kii Peninsula, Japan. *Geophysical Research Letters*, 45(6), 2618–2624. https://doi.org/10.1002/2017GL076773

Kato, A., & Nakagawa, S. (2020). Detection of deep low-frequency earthquakes in the Nankai subduction zone over 11 years using a matched filter technique. *Earth Planets and Space*, 72(1), 128. https://doi.org/10.1186/s40623-020-01257-4

Katsumata, A., & Kamaya, N. (2003). Low-frequency continuous tremor around the Moho discontinuity away from volcanoes in the southwest Japan. *Geophysical Research Letters*, 30(1), 20–1–20–4. https://doi.org/10.1029/2002GL015981

Kostogllov, V., Husker, A., Shapiro, N. M., Payero, J. S., Campillo, M., Cotte, N., & Clayton, R. (2010). The 2006 slow slip event and nonvolcanic tremor in the Mexican subduction zone. *Geophysical Research Letters*, 37(24), L24301. https://doi.org/10.1029/2010GL045424

Márquez-Azúa, B., & DeMets, C. (2003). Crustal velocity field of Mexico from continuous GPS measurements, 1993 to June 2001: Implications for the neotectonics of Mexico. *Journal of Geophysical Research*, 108(B9), 2450. https://doi.org/10.1029/2002JB002241

Masuda, K., Ide, S., Ohta, K., & Matsuzawa, T. (2020). Bridging the gap between low-frequency and very-low-frequency earthquakes. *Earth Planets and Space*, 72(1), 47. https://doi.org/10.1186/s40623-020-01172-8

Maubant, L., Radiguet, M., Pathier, E., Doin, M.-P., Cotte, N., Kazachkina, E., & Kostogllov, V. (2022). Interseismic coupling along the Mexican subduction zone seen by InSAR and GNSS. *Earth and Planetary Science Letters*, 586, 117534. https://doi.org/10.1016/j.epsl.2022.117534

Nishimura, T., Matsuzawa, T., & Obara, K. (2013). Detection of short-term slow slip events along the Nankai Trough, southwest Japan, using GNSS data. *Journal of Geophysical Research: Solid Earth*, 118(6), 3112–3125. https://doi.org/10.1002/jgrb.50222

Obara, K. (2002). Nonvolcanic deep tremor associated with subduction in southwest Japan. *Science*, 296(5573), 1679–1681. https://doi.org/10.1126/science.1070378

Obara, K., & Kato, A. (2016). Connecting slow earthquakes to huge earthquakes. *Science*, 353(6296), 253–257. https://doi.org/10.1126/science.aaf1512

Obara, K., Tanaka, S., Maeda, T., & Matsuzawa, T. (2010). Depth-dependent activity of non-volcanic tremor in southwest Japan. *Geophysical Research Letters*, 37(13), L13306. https://doi.org/10.1029/2010GL043679

Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, 75(4), 128–1154. https://doi.org/10.1016/0148-9062(86)90674-1

Perez-Campos, X. (2008). MASE: Undergraduate research and outreach as part of a large project. *Seismological Research Letters*, 79(2), 232–236. https://doi.org/10.1785/gssrl.79.2.232

Platt, J. P., Xia, H., & Schmidt, W. L. (2018). Rheology and stress in subduction zones around the aseismic/seismic transition. *Progress in Earth and Planetary Science*, 5(1), 24. https://doi.org/10.1186/s40645-018-0183-8

Poiata, N., Vilotte, J., Shapiro, N. M., Supino, M., & Obara, K. (2021). Complexity of deep low-frequency earthquake activity in Shikoku (Japan) imaged from the analysis of continuous seismic data. *Journal of Geophysical Research: Solid Earth*, 126(11), e2021JB022138. https://doi.org/10.1029/2021JB022138

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostogllov, V., & Cotte, N. (2011). Spatial and temporal evolution of a long term slow slip event: The 2006 Guerrero Slow Slip Event. *Geophysical Journal International*, 184(2), 816–828. https://doi.org/10.1111/j.1365-246X.2010.04866.x

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Walpersdorf, A., Cotte, N., & Kostogllov, V. (2012). Slow slip events and strain accumulation in the Guerrero gap, Mexico. *Journal of Geophysical Research*, 117(B4), B04305. https://doi.org/10.1029/2011JB008801

Radiguet, M., Perfettini, H., Cotte, N., Gualandi, A., Valette, B., Kostogllov, V., et al. (2016). Triggering of the 2014 M w 7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico. *Nature Geoscience*, 9, 829–833. https://doi.org/10.1038/ngeo2817

Rogers, G., & Dragert, H. (2003). Episodic tremor and slip on the Cascadia subduction zone: The chatter of silent slip. *Science*, 300(5627), 1942–1943. https://doi.org/10.1126/science.1084783

Rousset, B., Campillo, M., Lasserre, C., Frank, W. B., Cotte, N., Walpersdorf, A., et al. (2017). A geodetic matched filter search for slow slip with application to the Mexico subduction zone. *Journal of Geophysical Research: Solid Earth*, 122(12), 10498–10514. https://doi.org/10.1002/2017JB014448

Rousset, B., Fu, Y., Bartlow, N., & Bürgmann, R. (2019). Weeks-long and years-long slow slip and tectonic tremor episodes on the south central Alaska megathrust. *Journal of Geophysical Research: Solid Earth*, 124(12), 13392–13403. https://doi.org/10.1029/2019JB018724
Shelly, D. R. (2017). A 15 year catalog of more than 1 million low-frequency earthquakes: Tracking tremor and slip along the deep San Andreas Fault. *Journal of Geophysical Research: Solid Earth*, 122(5), 3739–3753. https://doi.org/10.1002/2017JB014047

Shelly, D. R., Beroza, G. C., Ide, S., & Nakamula, S. (2006). Low-frequency earthquakes in Shikoku, Japan, and their relationship to episodic tremor and slip. *Nature*, 442, 188–191. https://doi.org/10.1038/nature04931

Spiridonov, E. A., & Vinogradova, O. Y. (2020). Oceanic tide model FES2014b: Comparison with gravity measurements. *Izvestiya - Atmospheric and Oceanic Physics*, 56(11), 1432–1446. https://doi.org/10.1134/S0001433820110092

Thomas, A. M., Beeler, N. M., Bletery, Q., Burgmann, R., & Shelly, D. R. (2018). Using low-frequency earthquake families on the San Andreas Fault as deep creepmeters. *Journal of Geophysical Research: Solid Earth*, 123(1), 457–475. https://doi.org/10.1002/2017JB014404

Thomas, A. M., Beroza, G. C., & Shelly, D. R. (2016). Constraints on the source parameters of low-frequency earthquakes on the San Andreas Fault. *Geophysical Research Letters*, 43(4), 1464–1471. https://doi.org/10.1002/2015GL067173

Uchida, N., & Bürgmann, R. (2019). Repeating earthquakes. *Annual Review of Earth and Planetary Sciences*, 47(1), 305–332. https://doi.org/10.1146/annurev-earth-053018-060119

Ujiie, K., Saishu, H., Fagereng, A., Nishiyama, N., Otsubo, M., Masuyama, H., & Kagi, H. (2018). An explanation of episodic tremor and slow slip constrained by crack-seal veins and viscous shear in subduction Mélange. *Geophysical Research Letters*, 45(11), 5371–5379. https://doi.org/10.1029/2018GL078374

Wallace, L. M., Hreinsdottir, S., Ellis, S., Hamling, I., D'Anastasio, E., & Denys, P. (2018). Triggered slow slip and afterslip on the south-ern Hikurangi subduction zone following the Kaikōura earthquake. *Geophysical Research Letters*, 45(10), 4710–4718. https://doi.org/10.1002/2018GL077385

Wech, A. G., & Creager, K. C. (2011). A continuum of stress, strength and slip in the Cascadia subduction zone. *Nature Geoscience*, 4, 624–628. https://doi.org/10.1038/ngeo1215