Along‐Dip Segmentation of the Slip Behavior and Rheology of the Copiapó Ridge Subducted in North‐Central Chile

We studied the along‐dip influence of the Copiapó ridge subduction in the Atacama region, North‐Central Chile by building a new seismicity catalog, including similar events and non‐volcanic tremors (NVTs). We also obtained a 3‐D tomographic model for P and S‐waves velocity (and the implied Vp/Vs ratio). We identify along‐dip segmentation involving four distinct segments: a locked seismogenic zone hosting ordinary seismicity and clusters of similar events; a transition zone with NVTs and low seismicity; an aseismic zone with slow‐slip events; and a deep zone with abundant intraslab seismicity. The velocity models show differences among these zones, with low velocity anomalies of Vp and Vs coinciding with aseismic slip zones, indicating the possible presence of fluids. Due to the spatial distribution along‐strike and along‐dip of the aseismic zones, we propose that these differences in seismogenic behavior are generated by subduction of the heterogeneous seamounts associated with the Copiapó ridge.

remains challenging. In particular, we still need to understand why some regions are more prone to aseismic slip than others.
The Atacama North-Central region of the Chilean margin (26°-28.5°S), where the Copiapó Ridge (CoR) is subducting below the South American margin, is currently characterized by a seismic gap (Ruiz & Madariaga, 2018). The CoR is made of different seamounts of diverse geometries and considerable heights (∼2,000 m high) (Figure 1). The CoR subduction coincides with a zone of low interseismic coupling, separating two areas of high coupling that could potentially generate large subduction earthquakes (Mw > 8.5) (Klein, Metois, et al., 2018). In addition, seismic swarms have been reported in this area in 1973in , 1979in (Comte et al., 2002Holtkamp et al., 2011). The latter occurred between April and May 2006 with a total of 180 reported events ( Figure 1; Comte et al., 2006). A SSE in 2014 was detected deeper along the interface, at approximately 50 km depth (Klein, Duputel, et al., 2018) with an equivalent magnitude 6.9 (Figure 1). Continuous GPS data of the global network suggest that SSEs recurrently occur every 4-5 years in that area (Klein, Duputel, et al., 2018).
Large megathrust earthquakes also occurred in the past, the largest one in 1922 (Mw ∼ 8.5-8.8), and several moderate earthquakes have occurred more frequently for example, 1796, 1918, 1983(Beck et al., 1998Kanamori et al., 2019;Ruiz & Madariaga, 2018). The last event, with Mw 6.9, happened on 1 September 2020 and was part of a remarkable seismic and a-seismic sequence (Klein et al., 2021) (Figure 1). Most ruptures of moderate magnitude earthquakes (<7.5) that initiated north of the CoR stopped at the CoR, although the large megathrust earthquake of 1922 appears to ruptured through it.
We deployed a temporary seismological network to study the seismicity associated with the CoR and eventually to understand how it might control the along-dip seismic behavior of the plate interface. We conducted an effective search for similar events and NVTs, clues for potential aseismic slip. Using a non-linear tomography method, we built 3-D P-and S-wave velocity models and a Vp/Vs ratio model to investigate the rheology and the possible presence of fluids along the plate interface. Figure 1. Seismotectonic context. The pink lines and star indicate the historical and recent earthquakes. The green and lightblue triangles correspond to the two phases of the temporary broad-band network. The brown and yellow triangles correspond to semi-permanent broad-band stations of the Strasbourg university and the permanent stations of Centro Sismológico Nacional respectively. The white line indicates de CoR and its seamounts. The orange dots indicate the seismic swarm of 2006, and the red lines are slip contours of the 2014 deep SSE, each contour is 50 mm of slip (Klein, Duputel, et al., 2018). The segmented black lines correspond to Slab2.0 model (Hayes et al., 2018) and the white triangles indicate the trench.

Data, Methods and Results
Seismic monitoring is performed by five permanent stations of the Centro Sismológico Nacional (CSN) within 200 km around the study area. Three semi-permanent stations of the University of Strasbourg also monitor the region since 2019 (Zigone et al., 2019). In order to increase the detection capability for this study, we installed 10 temporary broad-band stations, in 2 different settings, during 4 months each (Figure 1). The first phase of acquisition ran from June 24th to October 8th of 2019 and the second phase from October 9th of 2019 to January 17th of 2020. In total, we acquired data for 8 months, with a total of 28 stations of which 8 were permanent during the entire experiment.

Seismic Catalog
During the 8 months of data acquisition, the network detected tens of thousands of local and regional events. In order to build a robust catalog, we implemented both manual and automatic phase-arrival picking. For the manual picking, earthquake detection was performed by a multistage approach. First, we roughly identified potential events by a conservative STA/LTA detection method on each station. Second, we filtered out the hundreds of very-low magnitude events and most of the outliers by considering only the potential wave-arrivals observed on for at least 5 stations. Finally, we gathered these potential-arrivals into potential-events and manually inspected them. We selected only events for which at least 1 station had an arrival time difference between P and S-wave smaller than 25 s, limiting the study area to approximately 200 km radius around the network. Picking was performed using Seisan (Havskov & Ottemöller, 2008). Due to the quality of the records obtained, we were able to manually picked a portion of the data between the months of August-September of the first phase and November-January of the second phase. We obtained 1,477 events made of 16,745 P and 14,653 S-wave arrival times.
To complete this catalog, we trained a neural network algorithm with the manual catalog to detect and pick arrival times. This algorithm was able to detect and pick accurately well identified events. To filter outliers, very small and quarry blasts events, we selected only events detected by at least 5 stations, with at least 7 arrival-times of which at least 1 corresponded to a P-wave and 1 to an S-wave. We located the resulting catalog to geographically select local events and uncover large arrival-time residuals potentially corresponding to misidentified data. We use this method for 35 days between September-October of the first phase and 45 days between November-December of the second phase, completing the rest of the registration that could not be done manually. This automatic approach identified 2,483 events with 18,914 P and 18,308 S-wave arrival-times. Therefore, our final catalog consists of 3,960 events between August 2019 and January 2020, with 35,659 P-and 32,961 S-wave arrival-times, with local magnitudes between 0.3 and 4.6, and the completeness magnitude of 1.45 (See catalogs in supplementary material).
In Figure 2b, we represent a projection of the seismicity on the profile A-A' identified on Figure 2a. The seismicity corresponds to all events between the two dashed black lines of Figure 2a. Three seismicity planes or double seismic zone can be recognized. This pattern has been observed in other regions along the Chilean margin (Bloch et al., 2014;Comte et al., 1999;Marot et al., 2013). The upper plane correspond to interplate seismicity and is bounded between 20 and 32 km depth. The intermediate plane, between 32 and 80 km depth, is apparently located in the oceanic crust (blue segmented line in Figure 2b). The lower plane (red segmented line in Figure 2b) is located in the upper oceanic mantle and presents most of the seismicity between 45 and 60 km depth, decaying toward 60 and 80 km depth, and increasing again between 80 and 110 km depth. The two deepest planes tend to merge between 100 and 120 km depth. The seismicity in the upper plate is scarce, with a marked eastward seismicity limit at 69.5°S. The clusters observed between 0 and 5 km depth correspond to mining activity ( Figure 2b).

Clusters of Similar Events
We investigated the possible presence of repeating earthquakes or similar events which could indicate if aseismic slip takes place, especially on the interface and near the region of the recurrent deep slow-slip observed by Klein, Duputel, et al. (2018). From the manually picked catalog we selected the events located at a depth of 60 km or above for a matched-filter search (Gibbons & Ringdal, 2006). This selection resulted in 908 templates with S-P times ranging from 6 to 25 s. By doing so the P-and S-waves won't be contained in the same window (e.g., Igarashi, 2020) however the moveouts, the different arrival times with respect to the first arrival across the network, will be constant in the matched-filter search. This allow us limit the decrease of similarity by using to large template waveforms while keeping the information of the S-P times intact.
To build the templates, we cut the waveform 1 s before the P-wave arrival on the vertical components and 4 s after the S-wave arrival on the two horizontal components to limit the overlap of the different seismic arrivals. Due to the magnitudes in our catalog, we conducted the scan using three components at each station, on data filtered between 2 and 8 Hz (Uchida & Matsuzawa, 2013) and 1-15 Hz (Uchida, 2019). P-and S-waves correlations were performed on the vertical and horizontal components respectively. We found repeating events with a frequency bands between 2 and 8 Hz, in more than three stations with a correlation coefficient greater than or equal to 0.95 (Figure 2c), however we did not find any in frequency bands between 1 and 15 Hz. Due to the difference in these results, we decided to define these events as clusters of similar events. We obtained a total of 27 similar events (Table S1) for the period 24 June 2019, to 17 January 2020. Twelve of these detections seem to be related to mining activity, with local magnitudes <1.0 and located outside our study area. The 15 remaining events represent 7 clusters of 2-3 events and are located around 27.5°S of latitude, at a depth of 18 and 29 km along the interface with one pair at a depth of 51 km, with a local magnitudes between 1.3 and 2.5 (Figure 2b). These similar events seem to be aligned with the CoR, potentially co-located with the deep SSE (Figures 1 and 2a) and in less coupled zones ( Figure S1 in Supporting Information S1).

Non-Volcanic Tremor Search
To detect non-volcanic tremors (NVTs) activity, we use the envelope correlation method (Ide, 2010(Ide, , 2012. Hypocenters of tremors are determined as follow: velocity data are band-pass filtered between 1 and 10 Hz. Filtered traces are squared, low-pass filtered at 0.1 Hz and resampled at 1 Hz. The distinct trace obtained is called the envelope data ( Figure S2 in Supporting Information S1). For all stations separated by less than 100 km, the envelopes of horizontal components are cross-correlated using a 5-min time window with a 150 s overlap. Following Sáez et al. (2019), a tremor is identified if the cross-correlation coefficient is greater than 0.6 for more than five pairs of stations.
Because of the elevated seismicity rate in the area, many local distant earthquakes were detected. We carried out a visual inspection to eliminate false detections. Figure 2a shows the spatial distribution of NVTs activity after visual inspection (Figures S3 to S8). This NVTs activity was identified uniquely on 20 September 2019 (Table  S2), with a clear lack of P-and S-wave arrivals (Figure 2d). The NVTs activity is located between 27°S and 27.5°S and 70.5°W-71°W, with dispersion in depth ( Figures S9 to S20), however, most of them are concentrated toward the downdip or transition zone (Figures 2a and 2b), which exhibits low coupling values ( Figure S1 in Supporting Information S1). The NVTs activity is updip of the SSE slip and, as the clusters of similar events, present a spatial correlation with the CoR subduction ( Figure 2a).

3-D Tomography Model
To build a 3-D tomographic model we use the INSIGHT code (Potin, 2016). The model consists of a set of Vp, Vs and Vp/Vs values at each node of a regularly spaced three-dimensional grid constituting the inversion grid. The inversion was carried out using a nonlinear minimization approach based on a stochastic description of the data and the model (see Text 1 in Supporting Information S1 for more details).
Using the arrival times from our catalog, we build local Vp and Vs velocity models and Vp/Vs ratio of the region from 26.5°S to 28°S and from 69°W to 72°W. These models were derived from the arrival times of 35,659 P-waves and 32,961 S-waves corresponding to 3,960 events on an inversion grid consisting of 1,175,878 cells with a longitudinal, latitudinal, and vertical resolution size of 3 km, 3 and 1.5 km, respectively (See Figure S21 and S22 for initial model used to make the Vp and Vs models and S23 to travel-time residues in the inversion in Supporting Information S1).
To visualize the velocity variations along-dip, we extracted a cross-section along profile A-A' (blue segmented lines in Figure 2a).   (Hayes et al., 2018), and the solid red line to SSE (Klein, Duputel, et al., 2018). The blue and red segmented lines indicate seismicity planes (black dots) in the oceanic crust and upper oceanic mantle. The less illuminated areas represent areas with lower resolution.

Discussion
Based on the seismicity distribution, clusters of similar events, NVTs activity and the velocity anomalies obtained from our Vp and Vs models, and Vp/Vs ratio (Figures 2b and 3) we inferred a clear and marked picture alongdip heterogeneity of the plate interface. From 20 to 35 km depths (zone A in Figure 4), interplate seismicity is concentrated in clusters up to a depth of approximately 32 km, coinciding precisely with high coupling values reported by geodetic studies (coupling >0.7, Figure S1 in Supporting Information S1). Intraslab seismicity presents some clusters potentially related to the re-activation of faults in the oceanic crust and also tends to form a seismicity plane toward the bottom of the oceanic crust. Clusters of similar events were identified in this zone at depths from 18 to 29 km (Figure 2b), close to the seismic swarms (Figures 1 and 2a) and in a zone with low coupling values ( Figure S1 in Supporting Information S1).
Similar and repeating events inside seismic swarms have been observed in the past (Pastén-Araya et al., 2018;Poli et al., 2017;Tréhu et al., 2015;Valenzuela-Malebrán et al., 2021). Several studies have proposed that these events might be caused by heterogeneities in the subduction zone due to the subduction of bathymetric features such as: (a) seamounts associated with oceanic ridges (Bilek & Lay, 2018;Valenzuela-Malebrán et al., 2021) and (b) more local structures like fractured zones, which are capable of transporting a greater quantity of fluids to the interface (Moreno et al., 2014;Poli et al., 2017). In this zone A, the Vp and Vs values are moderate and the Vp/Vs ratio is high, with values between 1.78 and 1.82 (Figure 3) suggesting the presence of fluids. According to Nishikawa and Ide (2018), fluid-rich zones are prone to exhibit aseismic slip and host seismic swarms with clusters of similar and repeating events. Therefore, our observations suggest that in zone A there are zones with different seismic behavior: zones with seismic slip characterized by elevated rates of ordinary seismicity with high coupling values (coupling >0.6; Figure S1 in Supporting Information S1) and zones with aseismic slip, characterized by low coupling (coupling <0.5; Figure S1 in Supporting Information S1), low seismicity rate, clusters of similar events, greater presence of fluids and regular seismic swarms (Figure 4). Subduction of seamounts or oceanic ridges can produce fracturing in the oceanic and upper crust and changes in frictional sliding or ductile deformation at the interface (Wang & Bilek, 2014). Between 20 and 35 km depths, clusters of similar events and swarms are spatially correlated with the subduction path of the CoR (Figures 1 and 2a). In addition, these similar events are located in less coupled zones ( Figure S1 in Supporting Information S1). Mochizuki et al. (2008) and recently Chesley et al. (2021) established the subduction of seamounts, with their different geometries, widths and altitudes can influence that the distribution of seismicity, produce low coupling above these structures and can trigger clustered repeating earthquakes and similar events. Plata-Martinez et al. (2021) related the distribution of slow earthquakes in the Mexico subduction zone to the difference in relief produced by the subduction of seamounts. Therefore, we propose that the geometric heterogeneities of the seamounts of CoR could influence the diverse seismic behaviors observed in this zone.
In zone B or transition zone (Figure 4), ordinary interplate seismicity is scarce, and NVTs activity is detected (Figures 2a and 2b). This boundary between interplate seismicity and the zone where the NVTs occurs is characterized as a transition between a zone of high coupling, unstable and rate-weakening friction and a zone of low coupling, stable and rate-strengthening friction ( Figure S1 in Supporting Information S1; Im et al., 2020). Under the assumption that zones hosting slow slip phenomena such as NVTs and SSEs may impede the propagation of seismic ruptures (Nishikawa et al., 2019;Rolandone et al., 2018), this boundary permits a rough estimation of the along-dip extent of potential future subduction earthquake ruptures that can be nucleated in the locked seismogenic zone. The low values of Vp and Vs and values of Vp/Vs ratio between 1.78 and 1.80 (Figure 3) indicate the presence of fluids in this zone and low normal stress, both factors favor the presence of NVTs (Shelly et al., 2006). Nishikawa et al. (2019) recognized NVTs activity in the southern segment of the Japan subduction zone associated to the subduction of seamounts. NVTs activity is mainly located updip of SSE and a little in zone A, similar to what has been observed recently in the Mexico subduction zone (Plata-Martinez et al., 2021). Moreover, NVTs activity is spatially correlated with CoR subduction (Figure 2a). Therefore we do not discard that the subduction of the CoR seamounts could affect the frictional properties in the zone hosting the NVTs activity.
Between 42 and 60 km depth (zone C in Figure 4), both Vp and Vs models potentially represent a large doming structure anomaly in the plate interface and in the oceanic crust ( Figure 3). This anomaly presents low Vp and Vs velocities and a Vp/Vs ratio between 1.76 and 1.78. Kato et al. (2010) recognized an anomaly with similar characteristics with low velocities in the Japan subduction zone, which they attributed to the subduction of an oceanic ridge. Therefore, we tentatively attribute this anomaly to the subduction of a large seamount associated with the CoR. This anomaly coincides precisely, along-strike and down-dip, with the slip of the SSE (Figures S33a-S33b in Supporting Information S1) and a notable lack of interplate seismicity (Figures 3 and 4). The low Vp and Vs velocities and Vp/Vs ratio suggest the presence of fluids, which could be brought through the ridge and released at these depths (Chesley et al., 2021). This idea is consistent with our observations, as fluid abundance is one of the possible causes of aseismic slip (Bürgmann, 2018;Nishikawa & Ide, 2018). These observations propose a change in the rheology of this zone, promoting an aseismic slip behavior at the interface without seismicity. Temperature variation could be another factor in the lack of interplate seismicity. Some studies suggest an increase in temperature due to crustal thickening caused by subduction of ocean ridges (Kato et al., 2010;Tsuru et al., 2002;Wang & Bilek, 2014). This increase in temperature would produce ductile deformation, which could also explain the lack of interplate seismicity.
Finally, the zone between 60 and 110 km depth (zone D in Figure 4) corresponds to the zone below the mantle wedge with stable sliding. There is an increase of the Vp and Vs velocities and the Vp/Vs ratio and a drastic increase of intraslab seismicity. We interpret this increase in seismicity as the onset of the eclogitization metamorphic reaction, due to the dehydration of hydrated minerals in the basaltic oceanic crust and serpentinized upper oceanic mantle (Behr & Bürgmann, 2021;Calvert et al., 2020;Ferrand et al., 2017). The fluids released by the eclogitization could migrate into the mantle wedge enhancing serpentinization in that zone (Hyndman & Peacock, 2003). Our Vp/Vs ratio values ∼1.80 ( Figure 3) and absolute velocity values of Vp ∼ 7.2 km/s and Vs ∼ 4.0 km/s ( Figure S34 in Supporting Information S1) support this interpretation.

Conclusion
Through seismicity analysis and a 3-D tomography model, we studied the along-dip influence of the CoR subduction in the Atacama region, North-Central Chile. The observed distribution of seismicity, the occurrence of clusters of similar events, NVTs and the Vp, Vs and Vp/Vs ratio anomalies allowed us to identify diverse behaviors with different zones along-dip hosting seismic and aseismic slip (Figure 4). These observations are spatially related with the subduction of the Copiapó ridge, a process that could promote along-strike and along-dip changes and favor different slip behavior. Our results present novel observations shedding light on how subduction of large bathymetric features can influence the distribution of seismicity, aseismic slip and the likely extent of along-dip megathrust seismic ruptures in the North-Central Chilean subduction zone.