The 2013 Slab‐Wide Kamchatka Earthquake Sequence

Studies of initiation of large earthquakes are usually focused on frictional instabilities occurring in the near vicinity of the future rupture. Possible contributions of long‐distance interactions with large‐scale tectonic instabilities remain unknown. Here we analyze seismic catalogs and geodetic time series during a few months preceding the 2013 M = 8.3 deep‐focus Okhotsk earthquake. This deep‐focus event is preceded by four intense seismic clusters in the seismogenic zone. GNSS time series in Kamchatka revealed a transient landward motion episode 1 month prior to the mainshock, consistent with an increase of seismogenic zone loading. This transient loading episode is accompanied by a doubling of the intermediate depth seismicity rate suggesting a transient slab pull as the origin. These observations question the constant subducting velocity hypotheses and may have implications in the understanding of the long‐distance along‐slab stress interactions and in their contribution to initiation of large deep‐focus earthquakes.

Supporting Information may be found in the online version of this article.
In this article, we report a detailed analysis of seismicity and ground displacement during the months preceding the M = 8.3 Okhotsk sea earthquake that occurred on 24 May 2013. The observations have been recorded by a regional seismic network (Chebrov et al., 2013;Chebrova et al., 2020) and continuous GNSS stations operated by the Kamchatka Branch of the Geophysical Survey of the Russian Academy of Science (KBGS RAS) (Shestakov et al., 2014). This deep-focus event originated in the Kamchatka subduction zone that accommodates 80 mm/yr of motion between the Pacific plate and the Okhotsk plate, a microplate part of the larger North America plate (Seno et al., 1996). The subduction terminates at 56°N, where the trench intersects the Aleutian arc (Levin et al., 2002. The slab interface has a dip angle of 55° south of 55°N down to depths of 600 km. North of 55°N, the broken slab stops at 300 km depth and the dip angle diminishes to 35° (Gorbatov et al., 1997;I. Y. Koulakov et al., 2011). The seismogenic section of the subduction has ruptured in five historical large earthquakes of M > 8.0 in the 20th century, the largest being a M = 9.0 earthquake in 1952 (Johnson & Satake, 1999;MacInnes et al., 2010). Geodetic measurements show a high interseismic locking of the seismogenic zone where the locked areas coincide spatially with historical earthquake asperities (Bürgmann et al., 2005).
The 2013 deep-focus M = 8.3 earthquake occurred at ∼620 km depth and produced horizontal and vertical surface displacements of up to 15 mm, recorded both by regional seismic and geodetic networks. The surface deformation can be modeled with a dislocation on a low eastward dipping plane compatible with the double couple focal mechanism, with slip amplitudes of a few meters (Shestakov et al., 2014;Steblov et al., 2014). The stress drop is 12-15 MPa, which is similar to the stress drops of shallow megathrust earthquakes, and an order of magnitude lower than the other historical deep-focus earthquake of M > 8 that happened in Bolivia in 2015 (Ye et al., 2013;Zhan et al., 2014). The mainshock has been followed by few aftershocks, which is common for large deep earthquakes. A M = 6.7 event is noticeable a few hours after and 200 km south of the mainshock, producing larger stress drops, possibly as a result of lateral stress heterogeneities within the slab.

Seismicity Observations
In this study, we use the earthquake catalog created with the routine analysis at the KBGS RAS. This analysis is based on records of the Kamchatka regional network (Chebrov et al., 2013). Both the network configuration and the earthquake determination approach (Droznin & Droznina, 2011) have not been significantly modified since 2011. KBGS RAS uses the modified regional magnitude scale (Gusev & Melnikova, 1990). The completeness magnitude over the Kamchatka -northern Kuril Islands regions is estimated to be M l = 3.5 (Chebrov et al., 2013;Levina et al., 2013) and even lower in some areas (Saltykov, 2019). Therefore, we only analyze earthquakes with M l ≥ 3.5 which is still much lower than available from global catalogs. Most of the stations of the network are installed in the Kamchatka peninsula and are located west of the Kuril-Kamchatka trench. With such one-sided station coverage, the hypocenter locations for earthquakes in the shallow part of the Benioff zone are subject to significant trade-offs between the lateral position and the depth with the errors in the letter reaching a few tens of kilometers (Droznin et al., 2019). KBGS RAS does not provide earthquake focal mechanisms in their catalog. Therefore, in our analysis we use Global Centroid Moment Tensor (GCMT) focal mechanisms (Ekström et al., 2012).
In this study, we focus on the months prior to the 2013 deep-focus M = 8.3 mainshock. During these months, the seismicity at large depths was very quiet. Two noticeable events of M = 7.7 at ∼630 km depth and M = 7.3 at ∼500 km depth happened in July and November of 2008, 5 years before. While seismicity was quiescent at large depths, an exceptionally intense seismic activity happened at shallow depths (Figures 1 and 2). Within the Kamchatka trench, the Pacific plate subducts below the Okhotsk plate at a rate of 8 cm/yr. In the months before the deep-focus M = 8.3 earthquake, the shallow seismicity is organized in spatial clusters, color-coded by their depth. The GCMT seismic moment tensors (Ekström et al., 2012) are represented for events of M l > 6. The red contours indicate the rupture areas of large historic earthquakes determined from aftershock distributions (Johnson & Satake, 1999) also represented in Bürgmann et al. (2005). The gray lines indicate the slab contours (Hayes et al., 2018 Figure 1). Half of the M l > 6 events from 2010 to 2018 (13 over 26) happened during these seismic clusters ( Figure 2). The four clusters are all located at the edges of the historical rupture area of the historical M = 9.0 earthquake of 1952 (Johnson & Satake, 1999;MacInnes et al., 2010) (Figure 1). The GCMT seismic moment tensors (Ekström et al., 2012) for M l > 6 events show reverse focal mechanisms for the clusters 1, 3, and 4 and a normal focal mechanism associated with the outer-rise events of the cluster 2. The last cluster with ∼6 times more events than the three other clusters, spanning depths from 0 to 100 km, occurred from 19 May to 22 May, two days before the deep-focus mainshock. In the GCMT catalog, 28 similar reverse focal mechanisms are computed for events with M l between 4.8 and 6.1 during the fourth cluster. Looking at the M l of individual events within the clusters ( Figure S2 in Supporting Information S1), the outer-rise cluster 2 seems dominated by a mainshock aftershock sequence. Clusters 1 and 3 still have events of large M l after the mainshocks, suggesting a mix of mainshock aftershock sequences together with swarms. The cluster 4 is however clearly a large swarm, without event of dominant M l . At intermediate depths, between 100 and 300 km, a steady seismicity rate is observed during these months with a transient doubling of the rate early April 2013, during ∼18 days (M 0 = 1.3 × 10 16 N⋅m). During this transient seismicity increase, events were not spatially clustered like in the shallow seismic clusters, but were spread over a large area ( Figure S1 in Supporting Information S1). This intense seismic activity at intermediate depths is unusual when looking at a 10-year long catalog of seismicity. Other than in April 2013, such seismicity rates only happened during aftershock sequences following large earthquakes in 2013 and 2016 ( Figure 2 and Figure  S3 in Supporting Information S1).

GNSS Observations
The GNSS data used in this study have been recorded by 10 continuous stations operated and maintained by the KBGS RAS. We have processed the time series for a 10-year period of time, from 2005 to 2015. Most of the sites have been active only a few years over the 10 years. Detailed GNSS analysis is described in Supporting Information S1. At first look, the 10-year long time series show transient oscillating signals with annual periods that are mainly seen by the vertical components but also on horizontal components.
When scrutinizing the months in early 2013, we could not detect any signal associated with the shallow seismic clusters. However, we identified a transient event early April 2013, with a duration estimated to be between 16 and 22 days depending on the GNSS sites, with average duration of 18 days, that coincides in time with the transient increase of intermediate depth seismicity rates ( Figure S4 in Supporting Information S1). This transient is mainly captured on the east position components. Transient signals in GNSS time series with duration of a few weeks can have various sources. They can either be common mode errors due to propagation of orbital or reference frame errors. Common mode errors have the same amplitude on all sites of a local network. The second possible source is related to surface water loading due to precipitations of rain or snow, with a dominant annual period and usually important local variations from site to site. The third possible source is tectonic and produces specific deformation patterns at the surface which depend on the physical processes at depth. We compared the transient event observed in 2013 with all other years available ( Figure 3). This event seems unusual compared to all other years, and the deviation in 2013 from other years is more significant at sites close to the coast (e.g., PETS, RADZ). Also, the spatial pattern of the transient signal ( Figure 4a) is specific, with larger amplitudes at the eastern coast (e.g., VODO, RADZ) then inland (e.g., ES1). This pattern is very similar to the decadal-estimated interseismic velocity field ( Figure S7 in Supporting Information S1). The transient is not recorded at sites far from the subduction like TILI located at 62°N (Figure 1), thus discarding the common mode hypothesis. We also looked at independent observations to evaluate the surface loading hypothesis. The transient snow melt episode recorded with satellite images show that the melting occurred from mid-May to mid-June ( Figure S5 in Supporting Information S1), more than a month after the observed GNSS transient. Given the larger amplitude at coastal sites, a possible ocean loading episode should be considered. We looked at the surface sea height anomalies along the Kamchatka coastline recorded with altimetry satellites ( Figure S6 in Supporting Information S1), and no unusual transient ocean loading event is recorded in April 2013.
The surface pattern of the transient event being similar to the interseismic loading pattern ( Figure S7 in Supporting Information S1) and the concomitance with the acceleration of intermediate depth seismicity rates (Figure S4 in Supporting Information S1) both suggest that a tectonic origin associated with the subduction zone is the most likely. Compared to the well-known transient slow slip events that are releasing part of the accumulated stress below the seismogenic zone at some subduction zones (Dragert et al., 2001), the displacements of the April 2013 transient deformation event are not trenchward, but landward. Figure S6 in Supporting Information S1 compares the interseismic velocities relative to fixed Eurasian plate and computed over a period of 8 years ( Figure S7 in Supporting Information S1) with the velocities during the transient event. The directions of the two velocity fields are very compatible, but the amplitudes of the velocities during the transient event are about 3.2 times larger than the ones during the interseismic period ( Figure S7b in Supporting Information S1), suggesting a transient increased loading event.

Kinematic Modeling
In order to obtain a tectonic explanation to the transient deformation observed with GNSS observations in April 2013, we consider dislocations compatible with the subduction interface model slab 2.0 (Hayes et al., 2018) at a range of depths from the trench to the bottom of the slab at 600 km depth as well as on the coseismic rupture plane published by Ye et al. (2013). For each model considered, we computed slip and slip deficit. The detailed modeling strategy is described in Supporting Information S1. The RMS error for both slip and slip deficit models are given in Table S1 in Supporting Information S1. To explain the landward motion of the transient event, all models located between the trench and 400 km depth have lower RMS errors with slip deficit, which suggest an increase of loading of the upper plate. Best models at deeper depths are compatible with actual slip, suggesting a release of accumulated stress. When comparing these different depth models ( Figures S8 and S9 in Supporting Information S1), in which the rake angle is fixed in the direction of the convergence, only the ones at seismogenic depths are able to fit the transient deformation polarities at all the stations of the GNSS network. Other models, and in particular the deepest ones, are not able to explain the northern sites polarities. One reason is the deep slab break-off for latitude higher than 55°N. The RMS error for the slip deficit model in the seismogenic zone is more than two times lower than the RMS errors of all the deeper models tested. We thus conclude that the surface transient deformation event is likely due to slip deficit on the seismogenic zone at depths lower than 60 km. Without the slab break-off north of 55°N, the signal could have been explained by a deep slip episode. Although the seismogenic zone slip deficit model is the most simple to explain the data, it is possible that deeper transient slip participates in the observed surface deformation. The best forward model with homogeneous slip deficit on the seismogenic zone has a root mean square (RMS) error of 0.93 mm for a slip deficit amplitude of 1.4 cm (Figures S8a and S8c in Supporting Information S1). The preferred inverse slip deficit model with a correlation length λ = 300 km has a RMS error of 0.78 mm. The fit to the east position time series with largest signal amplitudes is presented in Figure 4 and the fit to the north position time series is presented in Figure S10 in Supporting Information S1. The gain of RMS error between the inverse model with lateral slip deficit variations and the forward model with homogeneous slip deficit is 15%. This small gain shows that the main parameter explained by the GNSS transient offsets is the amplitude of the transient slip deficit, rather than the roughness of the lateral variations.
The interseismic locking C defined as the ratio between slip deficit rate δ and a reference velocity V r : = is a classical metric to define the degree of locking of the seismogenic zone (e.g., Avouac, 2015). The reference velocity V r corresponds to the interseismic subducting velocity of the subduction below the seismogenic zone. The long-term convergence velocity estimated at the trench is often considered as the reference velocity. In order to estimate the locking of the seismogenic zone during the 2013 transient event in Kamchatka, we consider an average slip deficit of 1.4 cm over a duration of 18 days. If we assume the long-term steady-state convergence velocity at the trench of 80 mm/yr, we obtain an average value of locking C = 3.6, which is physically impossible. Because the locking cannot be larger than 1, the reference velocity during this transient episode cannot be considered as the long-term convergence velocity at the trench. Assuming that the seismogenic zone locking keeps high values close to one as it was in average during the previous decades (Bürgmann et al., 2005), the reference velocity has to be ∼3.6 higher than the long-term convergence velocity to reach amplitudes of ∼28.8 cm/yr. Given that the slip deficit amplitude of 1.4 cm during the transient deformation episode corresponds to an upper bound since part of the signal could also be explained by deeper transient slip, the increase of the reference velocity by 3.6 also corresponds to an upper limit of acceleration. Since the subducting velocity at the seismogenic zone results from a force balance between the ridge push and the slab pull, one of these two forces has to transiently increase. No evidence tend to show an increase of the ridge push, however, the simultaneous increase of intermediate depth seismicity rate together with the subsequent deep-focus M = 8.3 event suggest an increased transient deformation at depth related to a transient slab pull acceleration.

Discussion
In the months before the M = 8.3 deep-focus Okhotsk earthquake, four intense seismic clusters happened in the seismogenic zone and outer-trench, 500-600 km above the large mainshock. The analysis of GNSS time series shows a transient signal which is best explained with a transient slip deficit on the seismogenic zone, compatible with a transient loading of the forearc. This transient event coincides in time with a doubling of intermediate depth earthquakes rate. Figure 5 shows an interpretation of the sequence of events that could explain the ensemble of observations. At steady-state, the subducting velocity is set by a balance between the ridge push and the slab pull forces that are both steady. The slab pull force results from the weight of the slab that pulls it inside the mantle and the viscous drag of the overriding mantle that resists the down-going motion. At shallow depths, the locked seismogenic zone compresses the forearc and produces landward motion of the GNSS often referred to as interseismic motion (Bürgmann et al., 2005). A disturbance from this steady-state system might be coming from a transient reduced viscous drag. Numerical simulations including a low viscosity zone on the subduction interface show that it could significantly increase the subducting velocity (Behr et al., 2022). The transient decrease of apparent viscosity of the slab interface could be due to increased fluid flow along the plate interface, knowing that slabs are thought to transport significant amounts of water to the deep mantle (Ohtani, 2005). Another possible explanation is partial melting at intermediate depths below the active volcanic arc, reducing the viscosity at the slab interface shear zone (Behr et al., 2021). The Kamchatka volcanic arc is one of the most active on earth with large melt volume production. Also, because of the proximity to the northern termination of the subduction, the toroidal mantle flow around the edge of the plate Yogodzinski et al., 2001;Peyton et al., 2001;I. Koulakov et al., 2020;Zhao et al., 2021) might modify melting conditions of the oceanic plate and have an influence on the subduction dynamics.
The transient reduced drag of the down-going plate induces a slab plunge acceleration that produces extension at intermediate and shallow depths resulting in increased intermediate depths seismicity rates and might also compress the bottom of the inner-slab. The increase of the slab pull force modifies the subducting velocity below the seismogenic zone. Given that the seismogenic zone is locked, the compression of the forearc is also increasing, as observed by the GNSS transient landward motion. We computed the coulomb stress change on the seismogenic zone associated with dip slip on the subduction interface below 60 km depth for an increased subducting velocity by 3.6 and it shows that the stress change on the seismogenic zone can get values up to 1 kPa. Such an amplitude is rather small, and likely cannot explain alone the triggering of the shallow swarms. This small amplitude, the time delay between the transient slab motion and the shallow seismic clusters and the swarm behavior of the seismic clusters indicate that the redistribution of stresses during this transient slab pull motion might have induced poroelastic fluid flow close to the seismogenic zone that triggered the intense swarm activity. Alternatively, the fact that shallow seismic swarms are happening before and after the transient deformation episode could also be due to a much longer duration of the transient slab pull, either continuous or as a cascade of several transient increased down-going motions, while only the fastest phase can be seen with GNSS observations.
On top of the intermediate depth extension, the transient slab pull could also have produced sub-vertical compression of the bottom inner-slab that might have played a role in the triggering of the deep-focus M = 8.3 earthquake that ruptured a near-horizontal eastward dipping plane (Ye et al., 2013) perpendicular to the compression a month later. Laboratory experiments show that acoustic emissions emerge for strain rates higher than 20% in the conditions of deep-focus earthquakes (Schubnel et al., 2013). The transient deep-slab compression might have participated to reach such a threshold.
Because GNSS and seismic signals are indirect consequences of the slab transient motion, we cannot quantify the exact volume displaced and the stresses induced at various depths of the subduction. Modeling of this large-scale subduction deformation event with realistic rheologies to reproduce an increase of subducting velocity at time scales of months could help quantifying these amplitudes. Such an event should have a viscoelastic response that might be measurable in GNSS time series in the following years. Modeling the viscoelastic response should help to constrain the associated transient rheologies. Finally, other observations, and in particular gravity measurements, might help obtaining constraints on the amplitude of the slab deformation (Bouih et al., 2022;Mikhailov et al., 2016;Panet et al., 2018).
This observation of a transient increase of the subducting velocity could have important implications for a better understanding of the deep-focus earthquake generation that might be the result of transient increased strain rate associated with slab plunge episodes. It also sheds light on slab-wide interactions between deep slab deformation and seismogenic zone activity. Finally, although this type of deviation from steady-state subducting motion might be rare, it potentially has implications for both large-scale geodynamic models and seismic cycle simulations in which the subducting velocities are often considered constant.

Data Availability Statement
The data used in the work were obtained with large-scale research facilities "Seismic infrasound array for monitoring Arctic cryolitozone and continuous seismic monitoring of the Russian Federation, neighboring territories and the world." The GNSS time series presented in this work can be found at: https://doi.org/10.5281/zenodo.6817009. The seismicity catalog used in our analysis is available from the website of the Kamchatkan Branch of Geophysical Survey of Russian Academy of Sciences: http://sdis.emsd.ru/info/earthquakes/catalogue.php.