Subdaily Slow Fault Slip Dynamics Captured by Low‐Frequency Earthquakes

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.

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 ≤ 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 ≥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).

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 2005and 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. 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 , 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). 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.

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 = μAd, where μ (assumed μ = 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 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 LFE( ) into the low-frequency earthquake median moment rate: 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 3 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 * Δt, with Δt the average source duration of a low-frequency earthquake and N LFE the number of low-frequency earthquakes; we assume a Δ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 ̇0 ( ) over the slip duration as: 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 2 . 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) ⋅ 1 km 2 . 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: 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(Thomas et al., , 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 ̇ at both daily and subdaily sampling rates, represented in Figure 2c, to match the time evolution of the geodetic fault slip when building our fault slip model. and subdaily 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.

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 : where d is the modeled fault slip, t is time, and ζ 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: 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 ̇ 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 α to attempt to capture any potential non-linearity in the geodetic to seismic slip relationship: where C is simply a scalar that compensates for the physical units that change as a function of α so that d S (t) maintains units of slip (m) and ζ remains a unitless ratio of geodetic to seismic fault slip. Introducing this α ) to a powerlaw model. There are now two parameters of interest in our model: the powerlaw exponent of seismic fault slip α, 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 α (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 longterm convergence of >6.4 cm/yr to be non-physical. For every set of V L and α, we then estimate the aseismic to seismic slip ratio ζ (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 α 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  (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.

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. . 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 . 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 α 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 α = 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. 10.1029/2022AV000848 8 of 15 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 multiyear 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, 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 α = 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.
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 longterm continuous partial coupling are likely non-physical and are an artifact of estimating surface velocities over long time scales.

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 α. A power-law model with a α > 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 α 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 ζ that changes by a factor of 4 between the two phases (14 for "Inter-SSE" and 63 for "Co-SSE" period). To physically interpret ζ, we define the ratio of geodetic to seismic slip as follows: where 0 and A G are the geodetic (total) moment and slipping area associated with slow slip respectively, and 0 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 0 0 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 ζ 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 ζ 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 slipping patch of 300 × 150 km, while the M6.4 events rupture an 150 × 100 km patch (W. B. . 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 ζ 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;. 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 Only the 24-hr sampled fault slip time series are shown here for clarity; similar results were find for the 6-hr sampling rate (ζ 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 earthquakes 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 α ≈ 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. 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 α 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.

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).

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 ζ
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 ζ. We can define ζ as: which can also be rewritten and simplified as: where 0 and A G are the geodetic moment and slipping area respectively, and 0 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 2 assuming each low-frequency earthquake source has a rupture area A of 1 km 2 ( 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 2 for the large M7.5 slow slip event and 100 × 150 km 2 for the smaller M6.4 slow slip events (W. B. . We then constrain ζ 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 , 0 and ζ 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. . 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. , we conclude that our assumptions supporting our estimates of seismic and geodetic fault slip are physically reasonable.