Plate‐Boundary Kinematics of the Afrera Linkage Zone (Afar) From InSAR and Seismicity

Studying the mechanisms of interaction between rift segments is key to understanding the kinematics of plate boundaries in continental rifts. However, the spatial and temporal evolution of deformation at rift linkage zones is rarely observed directly. Here, we combine InSAR data spanning 2005–2010 and 2014–2019 from ENVISAT and Sentinel‐1 satellites, respectively, with local seismicity from the Afar rift to investigate the plate‐boundary kinematics of the Afrera linkage zone, the junction between the Erta Ale and Tat Ali magmatic segments in Northern Afar (Ethiopia). We obtain time‐series of cumulative InSAR Line‐Of‐Sight (LOS) displacements that show deformation is accommodated by a series of active en‐echelon faults striking ∼NS and characterized by normal slip associated with a left‐lateral strike‐slip component. Additionally, we observe spatial variation in fault behavior with stick‐slip and creep. The faults in the center of the linkage zone behave primarily in a stick‐slip mode (with abrupt fault displacements up to ∼40 mm) and fault motions are associated with earthquakes of ML > 5. Conversely, faults at the edge of the linkage zone, near the magmatic segments, show creep and some stick‐slip behavior (with cumulative LOS displacement up to ∼30–40 mm over a ∼5‐year period) accompanied by low‐level seismicity. Some of the creeping faults are also spatially associated with hydrothermal springs. We interpret that the temporal behavior of the faults in the linkage zone is controlled by the interplay between tectonic extension, high heat flows, and fluid circulation near the magmatic segments where creeping of some faults is favored.

is distributed across a linkage zone and how much strain is accommodated by stick-slip faulting or by continuous creep is largely unknown. The complex fault patterns of linkage zones have prevented achieving a clear understanding of their kinematics. In recent years, deformation of different styles of oceanic transforms have also been measured in Iceland using geodesy, seismicity, and structural geology (e.g., Decriem & Árnadóttir, 2012;Einarsson, 2008;Geirsson et al., 2012;Metzger & Jónsson, 2014;Metzger et al., 2013;Opheim & Gudmundsson, 1989). However, little is known about active deformation in rift linkage zones of continental rifts with our knowledge is mainly based on theoretical models (e.g., Allken et al., 2012;Brune et al., 2017;Corti, 2012;Tapponier & Varet, 1974) and rare observations (e.g., Doubre, Déprez, et al., 2017;La Rosa et al. 2019;Pagli et al., 2019). Today, the growing body of SAR acquisitions on exposed continental rift systems offers the perfect opportunity to investigate, at high spatial and temporal resolution, the kinematics of rift linkage zones.
In this study, we combined InSAR time-series analysis and modeling with seismicity to probe the present-day kinematics of the Afrera Plain in Northern Afar (Ethiopia). Previous studies encompassing seismicity, numerical modeling, field observations, and geodetic measurements suggested that the Afrera Plain is an active linkage zone between the Erta Ale and Tat Ali rift segments, but whether this linkage zone may evolve into a transform fault is still debated (Bonatti et al., 2017;Illsley-Kemp et al., 2018;La Rosa et al., 2019). A recent study by La Rosa et al. (2019) used InSAR, seismicity, and structural geology to investigate a fault slip in Afrera associated with a M L 5.0 in 2007 and proposed a model of rift linkage consisting of a right-lateral, NW-SE-trending transfer zone where deformation is accommodated by dominant ∼NS-striking faults with oblique kinematics characterized by both normal and left-lateral components ( Figure 1). Here, we processed a large SAR data set from ENVISAT and Sentinel-1 (S1) catalogs spanning the time periods 2005-2010 and 2014-2019 to produce time-series of cumulative displacements and the related uncertainties. We identified not only deformation along different faults with both stick-slip faulting, in January 2018 but also fault creep throughout the observation period. Modeling of both types of faults showed motion along ∼NS-striking planes with oblique kinematics, similar to the 2007 event (La Rosa et al., 2019). We also investigated how seismicity accompany the observed deformation in space and LA ROSA ET AL.   Belachew et al. (2011), red dots are the seismicity during the 2007 seismic episode by La Rosa et al. (2019), and blue dots are the seismicity during 2011-2013 by Illsley-Kemp et al. (2018). The black lines are faults and fractures mapped in the area by La Rosa et al. (2019). The yellow diamonds are hydrothermal springs reported in the area by Brinckmann et al. (1970). (c) Simplified kinematic model for the Afrera Plain by La Rosa et al. (2019). time across the linkage zone. To this aim, we not only used the available catalogs 2005-2010 Illsley-Kemp et al., 2018;La Rosa et al., 2019) but also analyzed new data from two local seismic networks (Doubre, Leroy, et al., 2017;Keir et al., 2017;La Rosa et al., 2021) to study the seismic sequence in January 2018 and we identified two major earthquakes with M L > 5 (Figure 1a). This study presents one of the few direct observations of the fault kinematics in a rift linkage also detailing the temporal evolution of the strain. Our results show that Afrera is a linkage zone where the overall right-lateral motion between Erta Ale and Tat Ali is taken up by NS-striking faults with dominant left-lateral strike-slip and normal kinematics (Figure 1b). We document how creep and stick-slip are partitioned across the Afrera linkage zone and observe that creeping faults are located in proximity to the active magmatic segments and hydrothermal springs, suggesting that the creeping fault behavior could be controlled by high heat flows and hydrothermal fluids circulation.

Tectonic Setting
The Afar depression is the triple junction between the East African, Gulf of Aden, and Red Sea rifts accommodating the separation of the Nubian, Somalian, and Arabian Plates during the last 30 Ma (Barberi & Varet, 1970;Beyene & Abdelsalam, 2005;Manighetti et al., 1998). Starting from ∼11 Ma, the Red Sea rift stepped on land south of ∼16°N, creating the Danakil Depression in Northern Afar (e.g., Eagles et al., 2002). Extension in Afar was initially accommodated by large-scale border faults which today bounds the Afar rift floor from the Ethiopian Plateau to the west. Since the Quaternary, extension migrated to a series of en-echelon magmatic rift segments where the majority of magmatic activity and faulting now occur (e.g., Hayward & Ebinger, 1996;Keir et al., 2009).
The EA and TA segments strike ∼N330°E and are arranged en-echelon, with the EA segment left-stepping with respect to TA (La Rosa et al., 2019). The two segments partially overlap within a ∼20-km-wide area known as Afrera Plain (AP) (Figure 1). The AP is a depressed region reaching ∼100 m below sea level (b.s.l.) and hosting a salt water lake (Afrera Lake) fed by hydrothermal springs (Figure 1b) (Bonatti et al., 2017;Brinckmann et al., 1970). The AP is clearly bounded to the east by systems of west-dipping normal faults which are the northern-most tip of TA and also cut the eastern shore of the lake (Bonatti et al., 2017). Conversely, the western termination of the AP, as it links to the EA segment, is not as well defined. Evaporitic deposits and fresh basalts outcrop in the AP (Keir et al., 2013) with alignments of scoria cones and lava flows covering systems of tectonic fractures, suggesting that magmatic activity in the area is young, although these eruptive products have not been dated. High crustal V P /V s ratios (>2.0) were also measured in the area using receiver functions by Hammond et al. (2011) and were interpreted as due to the presence of magmatic fluids between EA and TA.
At the center, the AP is dissected by a system of ∼NS-striking faults (Figure 1b). Seismic catalogs from temporary networks in Afar Illsley-Kemp et al., 2018) showed seismicity along these fault systems indicating that the AP is tectonically active (Figure 1b). Seismicity at AP is characterized by low-magnitude earthquakes accompanying moderate earthquakes with M L > 5.0. The only previously analyzed earthquake is the M L 5.0 of October 2, 2007 (La Rosa et al. 2019). Fault modeling of InSAR data along with seismicity and structural data showed a combination of normal and left-lateral motions along ∼NS-striking faults, which were interpreted as evidence that the overall right-lateral shear between the EA and TA segments is accommodated by oblique slip on ∼NS-striking faults (Figure 1b

InSAR and Seismic Data
We analyzed a vast data set of SAR images by the European satellites ENVISAT and S1 covering the Afrera zone and spanning the time periods 2005-2010 and 2014-2019, respectively, to investigate how deformation is distributed across the linkage zone. We finally analyzed the local seismicity recorded by two temporary seismic networks active in Ethiopia during 2017-2018 ( Figure 1a) (Doubre, Leroy, et al., 2017;Keir et al., 2017;La Rosa et al., 2021).

InSAR Data Processing
We formed 92 interferograms from SAR images of the ENVISAT satellite in both ascending (track 028) and descending (track 049) geometries from 2005 to 2010. ENVISAT interferograms were generated using the ROI_PAC software developed by JPL/Caltech (Rosen et al., 2004). Topographic phase contributions were removed using an external 3-arc sec (∼90 m resolution) SRTM DEM (Farr et al., 2007). Interferograms were then filtered using a power spectrum filter (Goldstein & Werner, 1998) with strength of 0.6 and unwrapped using the ICU branch-cut algorithm (Goldstein et al., 1988). The unwrapped interferograms were then geocoded using the same SRTM DEM. We also produced 142 interferograms from S1 satellite images in both ascending and descending tracks (014 and 079) for the 2014-2019 period, using the JPL/Caltech/Stanford InSAR Scientific Computing Environment software package (Rosen et al., 2012). The SLCs coregistration and the topographic phase removal were performed using a 1-arc sec (∼30 m resolution) SRTM DEM (Farr et al., 2007). Residual noise and decorrelation were then filtered with an adaptive power-spectral filter with strength of 0.5 (Goldstein & Werner, 1998). Interferograms were then unwrapped with the ICU branch-cut method (Goldstein et al., 1988) and geocoded to the 1-arc sec SRTM DEM. We processed interferometric pairs with both small perpendicular and temporal baselines, which in the case of ENVISAT implies using all the available satellite acquisitions, with the exception of few interferograms with large baselines or that are noisy ( Figure S1). For S1, we processed almost all acquisitions between 2014 and 2016. Considering much more frequent acquisitions since 2016, we processed interferograms with small perpendicular baselines and spanning periods between 12 and 36 days, excluding some acquisitions with high level of noise ( Figure S2).

InSAR Time-Series Analysis
In order to minimize sources of errors in the interferograms, we estimated time-series of cumulative satellite Line-Of-Sight (LOS) displacements and their uncertainties for each InSAR track, using the Π-RATE software (Wang et al., 2012). We analyzed 51 interferograms for ENVISAT ascending track 028, 41 for EN-VISAT descending track 049, 55 for S1 ascending track 014, and 87 for S1 descending track 079 (Figures S1 and S2). Before the time-series analysis, we multilooked the S1 interferograms to a pixel size of 90 m × 90 m to reduce noise, this is also the same resolution as for the ENVISAT interferograms. Unwrapping errors were identified in ENVISAT data by adopting a phase closure method on Minimum Spanning Trees (Wang et al., 2012), while the excellent level of coherence in S1 data and lack of dense fringe patterns ensure that there are no unwrapping mistakes in this study. Using a network approach, we applied orbital filtering to the geocoded interferograms by fitting them with a linear-function (Biggs et al., 2007). We also removed topographically correlated atmospheric noise (Elliott et al., 2008) and applied an Atmospheric Phase Screen (APS) filter to minimize all other atmospheric disturbances (Ferretti et al., 2001). Abrupt displacements were also extracted from the time-series using the cross-correlation technique (Pagli et al., 2014) before applying the APS filter. To identify the time of abrupt displacements, we analyzed the series of interferograms and local seismicity catalogs, and we found that other than the previously reported M L 5.0 earthquake of Oc-  Figure S4 shows the results of cross correlation for the January 2018 episode. Another small sudden ground motion also occurred in February 2018 but only part of the signal is visible and some is covered by the Afrera Lake (Figures S3e and S3f). The signal is included in the time-series but cross correlation was not applied to it due to the lack multiple coseismic independent interferograms covering just this event. This is described with further detail in Section 3.  Xu et al., 2017), both of which occurred outside our study area. Dabbahu is ∼100 km away from the AP and our ENVISAT data set starts from October 2005, hence after the main intrusion occurred. The area affected by the Erta Ale volcano deformation is to the north of the AP but outside the selected study area. After abrupt displacements were identified and extracted from the interferograms ( Figure S4), we then applied APS using a combination of temporal high-pass and spatial low-pass filters with a temporal Gaussian filter with length of 0.5 year followed by a spatial Butterworth filter with low-pass cutoff estimated from the variance-covariance matrix of the spatially correlated noise. We then estimated time-series of incremental displacements and their uncertainties on a pixel-by-pixel basis using a weighted least-square approach and a Laplacian smoothing operator after APS removal (Wang et al., 2012). We selected a smoothing factor that minimizes the trade-off between the solution roughness and the Residual Sum of Squares of the displacements ( Figure S5). Our tests ( Figure S5) show the results are reliable because the inverted time-series are consistent with the raw observations with varying smoothing factors. Finally, the aforementioned estimated abrupt displacements of October 2007 and January 2018 were added back to give the resultant time-series.

Time-Series Results: 2005-2010
The cumulative time-series show three deforming areas of range increase (positive values) in both ascending and descending ENVISAT tracks (Figures 2 and 3) with deformation patterns corresponding to en-echelon ∼NS-striking faults mapped in the area. A series of profiles of cumulative InSAR displacement across different faults are displayed in Figures 4 and S6, showing that the InSAR deformation clearly match the faults. This demonstrates that a relationship between the InSAR deformation patterns and active faults exists. The range increase in both tracks is consistent with dominant normal faulting with downdip motion of the hanging-walls along 2-5-km-long faults.
The time-series of cumulative displacements also show that deformation is accommodated in different manners across the AP (Figures 2-4 and S6). We see sudden steps in ground motion indicating stick-slip faulting during the October 2007 earthquake (pixel 1 in Figure 3; Profiles A and B in Figures 4 and S6). However, we also see time-progressive continuous displacements indicating creep occurring along different subparallel faults during the observation period (pixels 2 and 3 in Figure 3; Profile C in Figures 4 and S6g-S6j). The total range increase along abruptly deforming faults is comparable to those deforming continuously at pixel 2 (∼30-40 mm), showing that stick-slip faulting accommodates as much strain as fault creep during [2005][2006][2007][2008][2009]. Conversely, less deformation is taken up by the fault at pixel 3 (∼10-30 mm) (pixel 3 in Figure 3; Profile D in Figure S6).
The InSAR observations show that faulting along ∼NS-striking faults is the dominant mode of deformation in the AP during 2005-2010. This is also consistent with InSAR and seismicity models of the M L 5.0 earthquake in 2007 (La Rosa et al., 2019). To further investigate the relationship between seismicity during 2005-2010 and the InSAR displacements, we analyzed the seismicity in the three deforming areas as identified by InSAR using local catalogs ( Figure 3) La Rosa et al., 2019). Figure 3 shows that the highest number of earthquakes occurred in the central fault system during the October 2007 seismic sequence when the M L 5.0 also happened. Few earthquakes are also observed in areas around pixels 2 and 3 indicating that fault creep is likely associated with minor seismicity.

Time-Series Results: 2014-2019
The S1 cumulative time-series show two main deforming areas of range increase in both the ascending and the descending tracks ( Figures 5 and 6). The deformation patterns are consistent with motions along short, ∼NS-striking fault planes. A deformation pattern observed during the 2005-2010 period also continued deforming in 2014-2019 (pixels 2 in Figures 2 and 3), and new abrupt displacements also occurred (pixel 4 in Figures 5 and 6).
The time-series at pixel 4 ( Figures 5 and 6) shows that the abrupt displacement occurred to the east of the October 2007 episode, confirming that extension in the center of Afrera is accommodated by stick-slip along a number of fault segments. The January 2018 displacement at pixel 4 also corresponds to a seismic sequence including two M L ∼ 5 earthquakes from January 10, 2018 to the January 31, 2018 (the seismicity is analyzed in Section 5). To the east of this event, the time-series at pixel 2 shows that fault creep is still ongoing similar to what was observed during 2005-2010 (Figures 5 and 6). The time-series in the descend-LA ROSA ET AL.
10.1029/2020JB021387 5 of 17 ing track show higher displacements in the eastern sector and near pixel 2 (∼40 mm) with respect to those measured in the same area in the ascending track (∼20 mm) (Figures 5 and 6). This is due to interaction with another abrupt displacement that occurred east of pixel 2 and north of the Afrera Lake in February 2018 ( Figures S3e and S3f), when a ∼19 mm of range increase along a short, ∼NNW-striking fault is observed. This sector of the study area hosts west-dipping faults bounding the eastern shore of the lake (Bonatti et al., 2017). The pattern is small and is observed just in the descending track likely because the normal and lateral component of the fault motion add constructively in the descending track but not in the ascending one. Yet the presence of the signal in two independent interferograms confirms that the signal is real. We tried to extract the February 2018 abrupt displacement before the APS filter using the cross correlation as done for the October 2007 and January 2018 episodes. However, this was not possible for a series of reasons, primarily due to the fact that we do not have two interferograms spanning the February 2018 event only (January 2018 event is included), the magnitude of the signal is small and its spatial pattern cannot be fully observed and it is partly covered by the lake (Figures S3e and S3f). We decided to include the February 2018 in the time-series analysis for completeness but without applying the cross correlation, hence the abrupt ground motion does not show as an offset in Figure 6d  Finally, the time-series from ascending track 014 also show some other patterns of range increase corresponding to mapped faults in the western sector of the AP (Figures 5a-5f), where time-progressive displacement signals were observed during 2005-2010 (pixel 3 in Figures 2-4 and S6). Such patterns are likely related to ongoing deformation along the same faults as observed in ENVISAT data. Indeed, similarly to ENVISAT, the displacements in S1 descending track are lower than those measured in the ascending track, if not completely zero. This could be due to the oblique fault kinematics and the fact that the strike-slip and normal-slip components add constructively in descending LOS but cancel out in ascending LOS.

Modeling of the January 2018 Earthquake
We modeled the displacement observed in January 2018 as fault slip by jointly inverting two S1 independent interferograms from both ascending (014) and descending (079) orbits assuming an Okada shear dislocation (Okada, 1985) with uniform slip across the fault plane in a conventional elastic half-space with a Poisson's ratio of 0.25 and a shear modulus (µ) of 3.2 × 10 10 Pa. For the modeling, we used a Monte Carlo simulated annealing algorithm, followed by a derivative based procedure, (Cervelli et al., 2001;Jonsson et al., 2002;Pagli et al., 2012). We inverted both the observed interferograms and those extracted with cross correlation. Initially, we inverted two 36-day long independent interferograms (this is referred as Model 1). The interferograms were subsampled using the quadtree partitioning algorithm from Jonsson et al. (2002). For Model 1, we set narrow bounds on the location and strike of the fault (between N330°E and N30°E) using the strike LA ROSA ET AL.  of mapped structures in the area as a priori knowledge, while we let the other parameters free to vary. The best fit Model 1 solution consists of a ∼NS-striking fault (N357°E) dipping to the east with an angle of 67° ( Figure S7). The fault is oblique with left-lateral (∼33 mm) and normal (∼45 mm) components ( Figure S7). Model 1 has relatively low Root-Mean-Square (RMS) residuals of ∼3 mm in both ascending and descending tracks (Table S1). However, the fault plane has an anomalous aspect ratio with a width of ∼10 km and a length of 3.5 km ( Figure S7 and Table S1). Furthermore, the best fit fault width corresponds to the upper bound of the search window. We then constrain the bounds on the fault width to be between 1 and 3 km (this is referred as Model 2) ( Figure S8). Model 2 shows that a solution with a 3-km fault width exists and the model is again an oblique left-lateral fault with geometry and kinematics comparable to Model 1. The fault strikes ∼N358°E and dips ∼76° to the east. The slip has 69 mm of left-lateral and ∼55 mm of normal components, resulting in a geodetic M w 5.0 ( Figure S8 and Table S2). Model 1 and Model 2 show similar RMS residuals of ∼3 mm in both ascending and descending tracks (Table S2), yet Model 2 has a geologically reasonable fault aspect ratio ( Figure S8).
Finally, we also inverted the January 2018 displacement signals as extracted from the time-series processing using the cross-correlation method (this is referred as Model 3). This allowed us to run the inversion minimizing the noise of the observed interferograms. The cross-correlated interferograms have been again quadtree partitioned (Jonsson et al., 2002). For the modeling, we again set relatively large bounds, except for the fault strike. Model 3 is again an oblique left-lateral fault striking N359°E and dipping to the east with an angle of ∼68°. The fault length is 3.6-km long and 4.8-km wide (Figure 7 and Table S3 (Table S3). The geodetic moment is 2.3 × 10 16 N m corresponding to a M w 4.9 earthquake (Table S3). We also calculated the uncertainties associated with each model parameter by adopting the Monte Carlo simulation of correlated noise Wright et al., 2003). We used the variancecovariance matrices of the input data to create 100 simulations of spatially correlated random noise. Such simulations have then been added to the input data and inverted ( Figure S9). The 90% confidence interval (CI) for each parameter has been calculated from the distribution of the 100 solutions (Table S4). The 90% CI for the fault length and width range between 3.5-3.7 km and 4.2-5.4 km, respectively. Normal slip is very well constrained with a 90% CI of 3.8-4.1 cm, as also the strike-slip component with a 90% CI of 1.1-1.5.
To summarize, all the models show similar fault kinematics characterized by oblique left-lateral slip along ∼NS-striking faults. The low variability in strike, kinematics, and magnitude suggests that these parameters are rather well constrained. Using the cross correlation (Model 3) has shown a great improvement in the quality of the input data, reducing the noise and the related local minima in the miss-fit space, and allowing the inversion for a better research of the global minimum. This is also shown by the very low variability in the fault parameters shown by the 100 Monte Carlo simulations ( Figure S9) and by the narrow 90% CI for the best fit solution (Table S4)

Modeling of Time-Progressive Displacement
We also modeled a creeping fault to understand the geometry and kinematics of these types of deformation. For the modeling, we selected the creeping fault just north of the Afrera Lake (pattern 2 in Figures 2 and 3) because it is far from other abrupt displacements. We jointly inverted the cumulative displacements between 2005 and 2010 of two InSAR ascending and descending tracks, using the same algorithms as for the stick-slip fault model. Our best fit model ( Figure S10 and Table S5) consists of an oblique left-lateral fault, striking ∼NS and dipping ∼64° to the east. The fault is 5.7-km long, 4-km wide and is located at a depth of 0.9 km, perfectly matching a mapped fault in the area (Figure 4). The fault accommodated ∼86 and ∼40 mm of left-lateral and dip-slip, respectively, during 2005-2010, which would correspond to a M w 5.2 earthquake. The creeping fault geometry, kinematics, and magnitude are similar to the stick-slip faults across the Afrera linkage zone.

Seismic Analysis
To better characterize the stick-slip events, we analyzed the seismicity accompanying the fault slip identified by InSAR in January 2018 and recorded by the seismic network deployed in Afar in 2017-2018(Doubre, Leroy, et al., 2017Keir et al., 2017;La Rosa et al., 2021). We inspected 1 month of continuous seismic recordings, from January 1 to January 31, 2018 to identify all the earthquakes during that period. We found a total of 499 earthquakes with the first P wave arrival at station N009 (located at Afrera) and manually picked both P and S waves for earthquakes recorded by four or more stations. We then located the earthquakes using the Oct-Tree search algorithm implemented into the NLLoc software (Lomax et al., 2000) and a 2.5D velocity model of Afar. We finally estimated local magnitudes (M L ) by measuring the zero-to-peak amplitude on simulated Wood Anderson seismometers and applying the distance correction for the Danakil region from Illsley  All the 499 located earthquakes cluster at the AP, along the main NNW-trending fault system, occur within the upper crust (1-10 km) (Figure 8 and Table S6) with average horizontal and vertical errors of ±2.8 and ±2.2 km, respectively (Table S6). The seismic sequence started the January 10, 2018 with a mainshocks of M L ∼ 5.2 ± 0.3 at depth of 9 ± 4 km ( Table S6). The mainshock was followed by several aftershocks with M L > 4.5 along with another shallower (∼2 km) M L 5.1 ± 0.4 on the January 11 ( Figure 8 and Table S6). The hypocentral locations are consistent with the InSAR time-series and modeling and indicate that two major earthquakes occurred. However, just a single fault was assumed in the InSAR modeling as the two events could not be separated temporally in both ascending and descending tracks. Since the depth of the M L 5.1 on January 11 was shallower (2.6 km) than the M L 5.2 mainshock on January 10 (∼9 km), it is likely that the displacement captured by InSAR is mainly due to the shallow earthquake, while the deeper M L 5.2 does not contribute much to the surface displacement. To test this hypothesis, we produced two simulated interferograms assuming an Okada shear dislocation source having the top edge located at a depth of 9 km. For the simulation, we used the same homogeneous elastic half-space, fault geometry, and kinematics of our preferred Model 3. We then set a fault slip of 134.5 mm, corresponding to a M w 5.2 and approximately comparable to the M L 5.2 calculated from seismic data. As can be seen in Figure S11, the simulations show negligible surface LOS displacements of 3.5 mm. This supports the hypothesis the only the M L 5.1 of January 2018 contributes to the observed InSAR displacement, as also indicated by the very good correspondence between LA ROSA ET AL.  the hypocentral depth estimated from both InSAR inversion (2.6 km) and seismic data (2.0 km) for the January 11 earthquake.
The cumulative seismic moment release (Figure 8b) shows earthquakes occurring in the 2 days (January 8 and 9, 2018) preceding the main earthquake. We attempted to process focal mechanisms for the earthquakes with M L > 4.5. However, due to the large azimuthal gaps of the seismic network in the study area, we have not been able to obtain unambiguous solutions. The majority of the aftershocks following the major earthquakes are shifted to the west (Figure 8a) likely suggesting that some small slip along nearby fault segments was triggered by the mainshock. We verified this by calculating the related Coulomb Stress Changes (CSCs) (see Section 6 and Figure S12).

Coulomb Stress Changes Calculation
Earthquakes modify the stress conditions on nearby faults by either promoting or inhibiting slip (e.g., Harris, 1998;King et al., 1994;Lin & Stein, 2004;Toda et al., 2005). These stress changes can be quantified with the Coulomb failure criterion, whereby positive CSCs caused by a fault slip (called the source fault) on a pre-stressed fault (called the receiver fault) can bring it to fail, while negative CSCs unload faults inhibiting their rupture. We explore whether some stress-triggering effect occurred between the three major earthquakes, the M L 5 in 2007 and the two M L ∼ 5 in 2018. For the CSCs calculation, we used the Coulomb 3 software (https://www.usgs.gov/software/coulomb-3), and faults are modeled as Okada shear dislocations within a uniform, elastic halfspace with the same elastic moduli as used in the InSAR modeling. The geometry and kinematics of both source and receiver faults were set based on the InSAR modeling results (La Rosa et al., 2019; this study).
We first analyzed the effect of the 2007 earthquake (source fault) on two receiver faults of January 2018. As the 2018 deeper event was not modeled by InSAR, we assumed that the fault has the same geometry and kinematics as for the shallow 2018 earthquake. The CSCs calculation (Figures S12a and S12b) shows that the 2007 earthquake caused negligible CSCs at depths greater than 6 km, therefore it is unlikely that it triggered the M L 5.2 seismic event of January 2018 as this occurred at ∼9 km depth. Furthermore, the 2007 earthquake caused a negative stress change (<−2 bar) at a depth of 2.6 km, where the M L 5.1 of January 2018 occurred (Figures S12a and S12b). The negative CSC means fault slip would have been inhibited. From this, in addition to the 2007 and 2018 seismic events being separated by a fair delay of 11 years, we conclude that the 2018 events were not triggered by the 2007 earthquake. Another factor such as tectonic strain accumulation is the likely cause of the 2018 seismic sequence.
We also calculated the CSCs imparted by the M L 5.2 earthquake on January 10, 2018 at ∼9 km depth (source fault) on the shortly following M L 5.1 earthquake on the January 11, 2018 at 2.6 km depth (receiver fault). Figures S12c and S12d shows that the M L 5.2 earthquake generated a positive stress change (>5 bar) at the base of the receiver fault with values of 0.5 bar at 2.6 km depth where the receiver fault hypocenter is positioned. Also, the area under positive CSCs at the surface fits well with the aftershock distribution. These observations indicate that the CSCs caused by the deeper earthquake on January 10, 2018 triggered the shallower earthquake on January 11, 2018.

Discussion
We used InSAR time-series and modeling, along with seismic data and CSCs calculations to document the present kinematics of the AP linkage zone in Afar. This data set, combined with the previous study on the AP ( linkage zone formed by the interaction of the EA and TA magmatic segments. We show that en-echelon, ∼NS-striking faults accommodate deformation at the AP. In the analyzed time period, different individual faults are observed to be active at different time periods suggesting that deformation is not limited to the central fault system (La Rosa et al., 2019) but instead it occurs on multiple faults with different behaviors across the linkage zone.
The major episodic displacements observed at the AP match the surface expression of mapped east-dipping faults. This is also supported by the east-dipping faults modeled from InSAR data for the major earthquakes that occurred in October 2007 (La Rosa et al., 2019) and January 2018. In addition, abrupt slip also occurs along west-dipping faults near the shore of the Afrera Lake in February 2018. This is also in agreement with field observations of dominant west-dipping faults along the eastern shore of the Afrera Lake (Bonatti et al., 2017). Overall, this indicates that extension at AP is mainly accommodated through dominant east-dipping faulting while active west-dipping faults characterize the eastern-most sector, where the AP connect to the TA segment.
The faults in Afrera also exhibit varying types of temporal behavior. Both the time-series and the individual interferograms show that abrupt fault slip events characterize the central fault system of the AP, suggesting that these fault segments have a dominant stick-slip behavior. Here, various fault segments accommodated deformation at different times with similar kinematics, fault size, and time scales of seismic sequences that included the three largest earthquakes of M L > 5.0. The most recent seismic sequence of January 2018 includes two M L > 5.0 along with several M L > 4.5 earthquakes. The best fit InSAR model shows that fault slip occurs along a shallow oblique fault striking in a ∼NS direction and dipping to the east (Figure 7). The two largest earthquakes perfectly correspond to the area with the highest displacement observed with In-SAR but just the second shallower (2.6 km) M L ∼ 5 of January 11 likely contributed to the observed surface displacement.
Stress triggering induced by the coseismic fault slips was explored using the CSCs modeling but the calculations show that either negative or negligible stress changes were imparted by the 2007 fault onto the 2018 faults (Figures S12a and S12b). These results, along with the long time-span between the seismic episodes, lead us to exclude any stress triggering from the 2007 earthquake. Conversely, the doublet ∼M w 5 earthquakes in January 2018 is a likely CSC triggering case (Figures S12c and S12d). Positive CSCs were caused by the earlier and deeper earthquake onto the shallower fault. Also, the area of positive CSC explains well the peculiar aftershocks distribution, all to the west of the faults.  Figures 1 and 3). It is thus possible that the fault creep is not aseismic, but instead associated with microseismicity. Conversely, a west-dipping fault at the eastern tip of the AP showed abrupt displacements in February 2018 (Figures S3e, S3f, and 5) indicating that some faults could creep and others stick-slip.
It is well known that several types of fault behavior can coexist along different portion of the same fault system or alternate in time along the same fault segment (Doubre & Peltzer, 2007;Harris, 2017). One of the best examples is the San Andreas fault where these phenomena have been well documented (e.g., de Michele et al., 2011;Harris, 2017;Rousset et al., 2019;Sammis et al., 2016). The fault behavior can be influenced by a wide range of factors encompassing temperature, presence of fluids, fault lithology, or a combination of these (e.g., Aochi et al., 2014;Byerlee, 1993;Byerlee & Brace, 1968Doubre & Peltzer, 2007;Harris, 2017;Vidale & Shearer, 2006). It has been shown that high temperatures in the deeper portions (>15 km) of a fault zone may change its rheological properties favoring fault creep (e.g., Brace & Byerlee, 1970;Harris, 2017). Similarly, shallower hydrothermal circulation resulting from the interplay between fluids and positive thermal anomalies alters the fault rocks and generates phyllosilicates that weaken the fault zone in the upper crust and reduce its shear strength, favoring fault creep (e.g., Moore & Rymer, 2007;Wintsch et al., 1995). Conversely, an increase in the pressure of fluids circulating within the fault zone may induce significant fault slip on pre-stressed faults and therefore cause seismicity (Aochi et al., 2014;Becken et al., 2011;Byerlee, 1993;Harris, 2017;Ross et al., 2020;Vidale & Shearer, 2006). Overall, the AP transfer zone links two en-echelon magmatic segments by oblique slip on multiple en-echelon faults. We primarily observed stick-slip fault behavior in the center of the linkage zone (Bonatti et al., 2017;Brinckmann et al., 1970). Conversely, we observe a combination of creep and stick-slip at the edges of the linkage zone, in the vicinity of the magmatic segments and where hot springs exist, close to the eastern faults. The observed variability of fault behavior, the results of the CSCs modeling, and the presence of hot springs suggest that the kinematics of the Afrera rift linkage zone is the result of the interaction between extensional tectonics, elevated heat flow, and presence of fluids near the magmatic segments which influences the fault behavior, either facilitating fault creep or promoting fault rupture by generating overpressures on the faults. This hypothesis is also supported by independent geophysical evidences which indicate the presence of magmatic fluids below Afrera .

Conclusion
In this study, we provided one of the few direct observations on how deformation is distributed across a linkage zone between two active magmatic segments. We combined a vast InSAR data set from ENVISAT and Sentinel-1 satellites, spanning 2005-2010 and 2014-2019, respectively, with local seismic recordings from the Afar rift. We show that deformation at the AP, between the EA and TA segments, is accommodated by several en-echelon, ∼NS-striking, oblique faults, which are arranged in two main structural architectures: dominant east-dipping faults characterize the center of the AP while dominant west-dipping faults are observed close to the eastern edge, where the AP meets the TA segment. Various fault segments are active at different time periods, showing great variability in their behavior. The fault segments at the center of the AP show a dominant stick-slip behavior characterized by episodic slip events associated with M L ≥ 5 earthquakes and related seismic sequences. Conversely, heterogeneous fault behavior encompassing creep, microseismicity, and minor episodic events characterizes the eastern tip. Here, a strong hydrothermal activity with hot springs and pools has been observed by Brinckmann et al. (1970) and Bonatti et al. (2017).
Our observations expand on our previous study showing the time-evolving behavior of fault segments at the AP (La Rosa et al., 2019). These new results support the kinematic model of rift linkage that we proposed for the AP where a right-lateral, NW-SE-trending transfer zone is accommodating deformation between the EA and TA segments. Within the middle of the transfer zone, ∼NS-trending oblique stick-slip faults accommodate most of deformation. In contrast, near the edges of the transfer zone, there is also deformation by fault creep, where we interpret that heat flow and hydrothermal activity influence the fault behavior and allow fault creep.

Data Availability Statement
The continuous seismic data used in this study are archived on the IRIS-DMC (http://www.fdsn.org/networks/detail/YQ_2017/; http://www.fdsn.org/networks/detail/YP_2017/). The data will become fully open access once the 3-year embargo period finishes at the end of 2021.

Acknowledgments
This study was developed in the framework of the PhD project of A.L.R. (XXXIII cycle of the Dottorato Regionale Pegaso in Earth Sciences) and is supported by the Ministero Università e Ricerca (MiUR) through PRIN grant 2017P9AT72. A.L.R. and C.P. acknowledge partial support by the University of Pisa grant PRA_2018_19. H.W. was supported by the NFSC project (41672225). ENVISAT ASAR and Sentinel-1 IW SLCs are provided by the ESA Online Dissemination (ENVI-SAT) (https://esar-ds.eo.esa.int/oads/ access) and by the Copernicus Open Access Hub (Sentinel-1) (https://scihub. copernicus.eu/). Seismic instruments were loaned by SEIS-UK. The facilities of SEIS-UK are supported by the Natural Environment Research Council (NERC) under agreement R8/H10/64. The temporary network deployment was part of an Actions-Marges and ISTeP project (#2016-82). We thank the regional authorities of the Afar, Tigray, and Amhara states for their administrative support. We also thank the wider staff body at the IGSSA of Addis Ababa University .