Strain Pattern and Kinematics of the Canary Islands from GNSS Time Series Analysis

: Following the 2004 seismic unrest at Tenerife and the 2011–2012 submarine eruption at El Hierro, the number of Global Navigation Satellite System (GNSS) observation sites in the Canary Islands (Spain) has increased, o ﬀ ering scientists a useful tool with which to infer the kinematics and present-day surface deformation of the Canary sector of the Atlantic Ocean. We take advantage of the common-mode component ﬁltering technique to improve the signal-to-noise ratio of the velocities retrieved from the daily solutions of 18 permanent GNSS stations distributed in the Canaries. The analysis of GNSS time series spanning the period 2011–2017 enabled us to characterize major regions of deformation along the archipelago through the mapping of the 2D inﬁnitesimal strain ﬁeld. By applying the triangular segmentation approach to GNSS velocities, we unveil a variable kinematic behaviour within the islands. The retrieved extension pattern shows areas of maximum deformation west of Tenerife, Gran Canaria and Fuerteventura. For the submarine main seismogenic fault between Tenerife and Gran Canaria, we simulated the horizontal deformation and strain due to one of the strongest (m bLg 5.2) earthquakes of the region. The seismic areas between islands, mainly o ﬀ shore Tenerife and Gran Canaria, seem mainly inﬂuenced by the regional tectonic stress, not the local volcanic activity. In addition, the analysis of the maximum shear strain conﬁrms that the regional stress ﬁeld inﬂuences the E–W and NE–SW tectonic lineaments, which, in accordance with the extensional and compressional tectonic regimes identiﬁed, might favour episodes of volcanism in the Canary Islands. in a framework of both extensional and compressional tectonic regimes, and our present results seem to conﬁrm this hypothesis.


Introduction
The use of permanent GNSS (Global Navigation Satellite System) stations in active geodynamic areas is increasingly gaining a pivotal role in investigating the kinematics related to tectonic plate movements, as well as crustal strain associated with tectonics, seismic cycles and/or volcanic processes (e.g., [1][2][3][4]). Presently, GNSS databases are routinely analysed by research groups and analysis centres of scientific organizations, such as the EPN (EUREF Permanent Network) from EUREF (Regional Reference Frame for Europe) or the Nevada Geodetic Laboratory [5][6][7]. In this study, we select an optimal period of approximately 7 years of continuous daily GNSS solutions from the MAGNET database [6] to identify tectonic and/or geodynamic features along the Canary Islands archipelago (Spain). To improve the signal-to-noise ratio, the common-mode component [8] has been subtracted from the time series of each observation site to retrieve high-precision geodetic velocities. Then, the strain has been computed using triplets of GNSS stations in triangular elements across the entire archipelago. The study led us to obtain a reasonable description of the displacement and velocity distributions to characterize the mechanism of the current deformation in that region.
The Canary Islands is a volcanic archipelago, where space geodetic techniques, such as GNSS, might shed light on the present deformation state, which consequently is useful to manage seismic and volcanic hazards and even to improve hazard warning systems in that region. There has been limited work, by [9,10], based on GNSS data applied to study tectonics and/or crustal deformation across the Canary Islands. Ref. [9] established a GPS network that was used by different surveys and concluded that any significant deformation was not detectable throughout the archipelago due to the poor sensitivity of the geodetic network at that time. Ref. [10] used seismic data and four permanent GPS stations located on El Hierro, Tenerife, La Palma and Gran Canaria islands to depict the geodynamics and tectonic activity prior to the 2011-2012 submarine eruption at El Hierro. They detected significant surface deformation in those islands and modelled a tectonic extension between the sites in the E-W orientation from 2003 onwards. Additionally, they strongly recommended long-term monitoring of tectonic activity (seismicity and deformation) as a useful tool for the volcanic hazard assessment of future eruptive events in the Canaries. A number of papers have focused on geodynamic, tectonic and/or volcanic studies on one island of the archipelago using GNSS (e.g., [11][12][13][14]). The presence of the Teide-Pico Viejo active volcanic complex in Tenerife, for instance, makes this type of study necessary due to the enormous volcanic hazard underlying this densely populated island. The interesting papers [15,16] highlight once more the need for this kind of study. Ref. [16] reported on the deformation state of Tenerife and claimed the need for a better definition and development of CGPS (Continuous GPS) networks in that area to improve the understanding of the Tenerife kinematics and the models on rest/unrest dynamics of the Teide-Pico Viejo volcanic complex. More recently, [15] analysed GNSS data from 2005-2015, highlighting a movement of Tenerife towards Gran Canaria that is worth studying in the context of the entire Canary Islands. This would make it possible to establish whether the seismic zones between the islands affect this movement or whether it is due to the continuous slow deformation of Earth's crust.
We present a new study of the deformation pattern along the Canary archipelago through the analysis of continuous GNSS data. We use data for the 2011-2017 period to explore this active area and provide precise horizontal displacement measurements and a broad calculation of the strain, which have not been done in previous studies. With our results, we analyse the deformation obtained using the space-geodetic GNSS technique in connection with the regional tectonics, which is proven to play a key role in the volcanic eruptions that take place in the archipelago.

Geodynamic and Tectonic Context
The Canary Islands are an intraplate volcanic archipelago lying on oceanic lithosphere close to the Moroccan passive continental margin ( Figure 1). The archipelago is located on the Nubian (African) plate, which is a very slow-moving tectonic plate. Each island forms an independent edifice, separated by narrow channels, except for Lanzarote and Fuerteventura, which share a NNE-trending submarine ridge parallel to the African coast (e.g., [17][18][19]). Each island in the Canary archipelago thus has an independent history and starting activity at different times and is evolving differently, making it difficult to establish a general evolutionary model [20][21][22]. Although the eastern islands (Lanzarote, Fuerteventura and Gran Canaria) are in a post-erosional stage of growth, the western islands (Tenerife, La Gomera, La Palma and El Hierro) are in a shield stage of growth, although La Gomera, which is a volcanic edifice with a complex and extensive evolution and a mid-level erosional stage, represents a post-shield gap [23].  [24]. (c) Some geological features of Tenerife: rifts (solid yellow lines); massifs (yellow areas); the Teide-Pico Viejo volcanic complex; the Las Cañadas caldera wall (dashed yellow line); the volcano-tectonic direction (dashed white line), which is an inferred submarine alignment based on the volcanic seamount distribution (brown triangles) and a slope break in the SW flank of Gran Canaria proposed by [25]; and the sinistral strike-slip fault (solid white line) located between Tenerife and Gran Canaria.
The origin of the Canary Islands is still controversial: a mantle plume, regional tectonics or a combination of both are still under debate, although the unifying model of [26] is the most accepted model. Several authors (e.g., [26][27][28]) have proposed a link between the Atlas Mountains in Africa and the genesis of the Canaries based on magmatism and tectonic analogies. Ref. [26] hypothesized the existence of a residual plume under North Africa, the Canary Islands and the Western Mediterranean and a fracture system that allows the magma to ascend. Recent geodynamical models propose intraplate tectonics driven by asthenospheric upwelling (mantle plumes), while others are based on plate boundary processes or lithospheric thinning that favour asthenospheric upwelling [29,30]. Ref. [10] showed that anomalous tectonic activity at different spatiotemporal scales in the Canary Islands region since 2003 preceded the El Hierro magmatic activity that gave rise to the recent submarine eruption in 2011. In addition, [31] emphasized that there is no conflict between the mantle flow model and the involvement of lithospheric processes in magma generation and ascent in the Canary Islands. Additionally, they aimed to determine the distribution of volcanism through the interaction among the delamination of subcontinental lithospheric mantle material, inflow of mantle plume material and reactivation of geological structures driven by the African-Europen convergence.
The Eurasia-Africa kinematics, the Mid-Atlantic Ridge and other regional structures control the plate-scale stress field in the Canary archipelago [32]. Likewise, the regional fault systems influence the tectonics in the region. In particular, [24] suggested the connection between the fault systems related to the Mid-Atlantic Ridge (NE-SW orientation) and the Trans-Alboran Fault System (NE-SW orientation), which runs along the Atlas Mountain range (northwestern Africa), to explain the kinematics of the Africa-Eurasia plate. These authors proposed that the Moroccan microplate passed through the Canary archipelago at its southern limit, limiting the extension of the Trans-Alboran Fault System to the Atlantis Fracture Zone (see Figure 1). The Moroccan microplate is decoupled from the Eurasian plate by tectonic activity in the Gibraltar-Azores belt, the Gorringe thrust, the Gloria fault and the Terceira ridge. The strike-slip motion at the Canary-Transalboran fault system results in decoupling from the Nubian plate, and the boundary between Moroccan microplate and North American plate undergoes an E-W stretching [24]. A submarine sinistral strike-slip fault, with reverse components, is located between Tenerife and Gran Canaria, and its motion results in the underthrusting of the western block (Tenerife) and a southwest relative displacement of Tenerife with respect to Gran Canaria and Africa ( Figure 1). Refs. [28,33] proposed the existence of this crustal fracture based on seismic and gravity data analysis. In addition, [34] studied the focal mechanism of the Mw 4.0 earthquake preceding the 2011 eruption at El Hierro. They found NNW-SSE horizontal compression and E-W quasi-horizontal extension, suggesting that the rupture process was related to tectonic stress. In addition, [35], after analysing the seismicity that occurred at El Hierro during 2011-2014, proposed that the movement of a regional fault triggered the 2011 eruption and speculated that local tectonics, not magma overpressure, facilitated magma ascent. The NNW-SSE orientation matches the orientation of stresses obtained by [33] using data from the moderate Mw 5.0 earthquake that occurred between Tenerife and Gran Canaria in 1989 and agrees with the orientation of the regional compression model described by [36,37] for the archipelago. Recent aeromagnetic studies by [38][39][40] in the central islands of the archipelago revealed elongated anomalies with E-W (Tenerife and La Gomera) and ENE-WSW (Gran Canaria) orientations. Furthermore, [39] accounted for the intrusive body found in the submarine edifice of Gran Canaria, which has an ENE-WSW orientation that coincides with the aforementioned fracture suggested by [28]. All these facts shed light on the role of regional tectonics in the origin of the Canary Islands and the magma upwelling from the mantle.

Analysis of GNSS Data
GNSS-derived strain rates have countless applications in geosciences, including, for instance, in seismotectonic and geodynamic studies aimed at comparing stress fields and rates to explain earthquake occurrence or interpret peculiar plate motion features in diffuse plate boundary areas [41]. Currently, there are two possible approaches to estimate strains from GNSS observations. A gridded approach, called least-squares collocation (LSC), consists of inverting a uniform velocity field to the strain rate field [42,43] or calculating strain rates in triangular or more complex segments [44]. A very detailed description of the pros and cons of the different methods of strain retrieval can be found in [45]. Briefly, in the LSC method, the velocities of n stations obtained in t 1 and t 2 GNSS survey campaigns or in a time interval T = t 2 − t 1 from permanent stations are used to reconstruct a continuous velocity field. However, to move from a discrete collection of measurements in time/space to a continuous field, a generalized interpolation problem must be solved. Ref. [45] found that the strain rate field obtained in such a way leads to clearer deformation patterns than that retrieved from the discrete triangle or "segmentation" method. Although many recent works (e.g., [46][47][48]) are still based on the segment (triangular) approach, the main motivation for this choice is that it allows avoidance of any distortion of the GNSS network caused by the datum itself. Moreover, this method has proven to be more suitable for not only highlighting discontinuities of the strain field, such as faults, which make the strain field highly nonuniform, but also strain analysis over small areas with a limited number of measurement points. This is why geodesists still use similar methods for local strain analysis related to, for instance, mine exploitation or volcano monitoring [47,49]. For these reasons, we decide to use the same approach. With this method, the area is segmented into several geometrical units, the local strain and rotation rates are computed for each unit, and finally, all the individual results are integrated. Although the segmentation method, in principle, is not exclusively limited to triangulation (it can have any geometrical shape), Delaunay triangulation of the measurement points remains the most used for this approach.
The analysed GNSS data set comprises a time series of 3D (North, East, Up) site positions, which are available as daily solutions produced by the Nevada Geodetic Laboratory [6]. The solid earth and pole tide contributions are reduced according to the "IERS 2010 Conventions" [50], while ocean tide loading (diurnal, semidiurnal, Mf and Mm waves) accounts for FES2004 [51]. The data set encompasses 6.7 years between 2011 and 2017 and refers to GNSS sites spanning the Canary archipelago (Table 1, Figure 2). The daily GNSS solutions refer to the IGS14 realization of the International Terrestrial Reference Frame 2014 (ITRF2014) [52].  It is common knowledge that antenna motion, not ground displacement, is directly observed; thus, to observe crustal movements, all site self-motions should be eventually excluded so that the observed displacement of antennas could be interpreted as a real tectonic motion. This is mainly true in regions where crustal deformations are small; consequently, noise in calculated geodetic strain rates, which in general are smaller and less reliable, could be considerable and even larger than the tectonic signal.
Therefore, in our analysis, after removing outliers due to nearby earthquakes and known equipment changes from the GNSS time series, the signal-to-noise ratio in the daily solutions is improved by using the common-mode error according to the formulation proposed by [8]. This method is a stacking technique resulting in spatial filtering aimed at removing random noise and extracting spatially correlated transients, the so-called common-mode component (CMC), that deviates from a zero mean over the span of the position time series of a continuous GNSS network. A well-known kind of CMC in continuous GNSS position time series, especially in the vertical direction, is the seasonal motion, with mostly annual and semi-annual periodicities. Such an approach [8] consists of a three-step procedure: (1) Residuals-Stacking-Filtering: for each daily solution (d) and GNSS site (s), the residual (r) [observed (O)-predicted value (C)] is computed: (2) Computation of the common mode by means of stacking, by averaging the residuals from all the sites (S) in the network: Equation (2) represents the common-mode bias, which, in principle, should improve with the number of sites. (3) Filtering: for each site, the common-mode bias is subtracted from the observed position, resulting in a filtered position:Ô Thanks to this filtering, we are hopefully able to remove any signal common to a subset of GNSS stations from position time series, regardless of the origin.
As a result, we obtain a lower standard error on the retrieved linear trend (ITRF absolute velocity). For such improved velocities (Table 1), we compute, for each geometrical unit (triangle), the bi-dimensional horizontal infinitesimal strain from 23 sets of three non-co-linear GNSS stations each, forming triangles, as formulated by [53] for a flat earth.
Under the assumption that the partial derivatives of the displacement, denoted in vector notation by u of components u x , u y are small, several approximations can be made in the calculations, which lead to the following results about the normal strains corresponding to the x and y directions at point (x,y): Customarily, negative values of ε indicate contraction, while positive values represent stretching. In addition to stretching and contraction, the change in the angle between the positions before and after deformation is known as the shear strain. In particular, the shear strain ε xy is defined as half the increase in the angle formed by two infinitesimal line segments that initially are parallel to the x and y axes.
Following this definition and considering symmetry, ε yx = ε xy . These components are the four elements of the strain matrix.
Therefore, as discussed above, the strain is a function of the normal strains, shear strains and rotation (Θ).
= f ε xx , ε yy , ε xy , Θ By setting i = 1, 2, ε 1 and ε 2 become the two principal horizontal strains, and the maximum shear strain is In our analysis, the following first-order assumptions are implicitly made: (1) As the study area is small, the sphericity of Earth is neglected, and the Earth's surface is approximated by a flat surface model. Hence, plane-strain deformation remains in the horizontal plane. (2) Considering the principles of mechanics, we assume that the deformation of networks reflects the real movements of Earth's crust and that the area between geodetic points could be considered as a continuous medium in which the strain is homogeneously distributed across the triangular area between the three GNSS sites. (3) Vertical velocities do not significantly affect the physical interpretation of the strain. Even if the latter assumption above is naive in areas displaying significant uplift or subsidence, in our case, stations with negligible vertical displacement (<3 mm/year) have been selected; hence, we are confident that 2D strain can provide useful insights into ground deformation patterns. The parameters of the strain were calculated from displacements (u x , u y ) estimated by a constant velocity for the whole study period (6.7 years between 2011 and 2017). The velocities measured at the three GNSS sites are the result of the three components of crustal velocity: translation, rotation and distortion. For each of the GNSS sites, we know the horizontal coordinates of the initial site location (x o , y o ) as well as the east-west and north-south instantaneous velocities (V x , V y ). We solve for the "3-site triangle strain" problem, which is a well-known perfectly constrained problem [43] (p. 157). The principal strain axes are found by computing the eigenvectors of the 2-D strain tensor . The eigenvectors are unit vectors in the directions of the principal strain axes. The eigenvalues of are the principal extensions in the principal directions. The larger of the two eigenvalues are the greater principal extension, ε 1 , and the length of the semi-major axis of the horizontal strain ellipse is equal to the greater principal strain, S 1 = ε 1 + 1. The semi-minor axis is S 2 = ε 2 + 1, where ε 2 is the lesser principal extension. For each triangle, we compute the strain in terms of elongations (strain unit) and resulting dilatation (mean strain), the latter of which is referred to as the triangle centroid. The main strain parameters for each triangle are listed in Table 2. Included in the Supplementary Materials is Figure S1, which illustrates the geometric configuration of each triangle.  Figure 2 displays the velocities (horizontal component) and the velocity space calculated after the CMC analysis with respect to the ITRF2014 frame for all GNSS stations used in this study. The retrieved velocity field is relatively dense compared to the individual size of each island. The velocities are mostly grouped, except for three sites around the central islands (IZAN, TERR, ALAJ in Tenerife, Gran Canaria and La Gomera, respectively). For all the sites, the mean value of the velocity module is 23.8 ± 0.38 mm/year with an azimuth of 46.7 ± 0.95 • E, which agrees with the Nubia plate velocity from the NUVEL-1A model [54]. This is in accordance with the mean value of approximately 23.2 ± 1.03 mm/year in the NNE direction calculated by [9] for three permanent stations located on Lanzarote, Gran Canaria and La Palma during 2002-2009.

Results
The horizontal displacements were estimated from the improved velocities after CMC filtering (see Section 3), which were considered constant for the time interval of 6.7 years used in this study. Then, to interpret the crustal motions along the Canary archipelago, we computed the strain parameters for each triplet of GNSS stations following the segmentation (triangle) method. Table 2 lists all the numerical results in terms of baseline length and triangle area change. To facilitate the reading of the results obtained on the triangular mesh network covering the entire archipelago, the largest elongation and shortening values are collected in Table 3. These should reportedly reflect the main tectonic changes unveiled by our strain analysis. The pattern that emerges from this analysis clearly shows that the baselines encompassing the known fault line intersecting Tenerife and Gran Canaria undergo the maximum elongations (Figure 3).   It is worth noting the particular case of Tenerife. The mean velocity computed from five GNSS stations (Table 1) is 23.6 ± 0.3 mm/year, which is in agreement with the value reported by [16], 23.3 ± 0.2 mm/year, determined by analysing daily GNSS solutions of four sites distributed across the island for the period 2008-2015. This value is also close to the velocity estimated by [15], 23.1 ± 0.5 mm/year, using a set of nine GNSS stations between 2005 and 2015.
Regarding Lanzarote, we found that the baseline YAIZ-HRIA has a remarkable shortening of approximately 10 −8 . This is consistent with a recent detailed geodetic study by [14] of the Timanfaya volcanic area. The authors, relying on three years (2014-2017) of continuous GNSS data and approximately twenty years (1997 to 2017) of continuous tilt and strain measurements, detected a subsidence and a peculiar behaviour in the southwest sector of Lanzarote. This result is also in agreement with that achieved by [55] using InSAR data pointing to 6-7 mm/year subsidence in the Timanfaya volcanic area. Other baselines along the archipelago show negligible elongations or shortenings on the order of 10 −9 ( Table 2). Tables 2 and 3, the most relevant changes on the order of 10 −7 are encountered on Tenerife. The largest elongation occurs along the baseline IZAN-GRAF (1.530 × 10 −7 ), and the largest shortenings occur along the baselines STEI-IZAN (−2.529 × 10 −7 ) and TN03-IZAN (−1.009 × 10 −7 ). This indicates a differential kinematic behaviour of the central and western areas with regard to the northeast Tenerife. Indeed, the northeast area comprises the stable sector of the Anaga massif, which is a portion of the old basaltic shield structure of the island together with the Teno and Roque del Conde massifs [56], and part of the Dorsal Rift fissure system (see Figure 1). Ref. [16] retrieved different deformation patterns for Tenerife, such as compression in the central-western part and extension to the northeast part (parallel to the Dorsal rift). Likewise, [15] pointed out the subdivision of Tenerife into two areas with different kinematic behaviours: extension in the northeast and central areas of the island and compression in the west. According to these authors, the different behaviours observed in the two deformation patterns along the island could be produced by the effect of the submarine fault ( Figure 1). Ref. [57] described the large-scale deformation pattern found along Tenerife as due to predominant vertical kinematics and proposed gravitational sinking of the island's edifice as a mechanism. Later, [16] even reported a decrease in elevation, although they explained the observed subsidence by the combined effect of the lateral spreading of the elastic lithosphere and the descent of the water table level in the island. To understand these findings within the general framework of tectonics of the entire archipelago, our result suggests that the behaviour of the northeast sector of Tenerife follows the predominant extensional regime in the archipelago (see Figures 4 and 5) and could be a consequence of the effects of the submarine fault in the tectonic regime in this area.    Figure 4a shows the horizontal strain centroids determined from the principal axes of the 2D strain tensor calculated for the 23 triangles (see Table 2) defined by the GNSS stations. On the one hand, the magnitude of the tensile strains ranges from −4 × 10 −9 to 203 × 10 −9 , and that of the compressive strains from −254 × 10 −9 to 4 × 10 −9 (respectively, blue and red colour in Figure 4a). The general pattern of the principal strain axes reveals that the largest magnitudes are localized around the central islands of the archipelago. There are significant changes in the strain orientation pattern as well. The extension (tensile strain) of the eastern islands is generally trending NE, whereas that of the central and western islands is trending NNW. There is evidence of a wide area of crustal compression around the central and western sectors of Tenerife, from a dominant WNW direction. This finding is in accordance with the stress field modelling of [32] for the Canary Islands, which indicates that the stress distribution and orientation vary considerably after the Atlantis Fracture Zone is included (Figure 1), enabling the detection of areas of extension and strike-slip and not just a compressive regime through the whole archipelago. Table 4 lists the largest elongation and contraction values found in this triangle strain analysis. On the other hand, the strain tensor calculation allows us to draw a dilatation map of the archipelago, which gives us an idea of the areas of horizontal strain accumulation that might be related to surface deformation produced by tectonic stress or fault activity. Thus, the negative values in Figure 4b indicate compression, and the positive values indicate extension. The figure also shows the seismic activity recorded by the National Geographic Institute (IGN) of Spain localized in the surroundings of the central islands of the archipelago for the period of study (2011)(2012)(2013)(2014)(2015)(2016)(2017). A clear extension pattern along the archipelago might be identified, and the areas of maximum deformation are located to the west of Tenerife, Gran Canaria and Fuerteventura. Additionally, the distribution of areal strain results in the predominant tensile behaviour in terms of the principal strain above the major fracture between Tenerife and Gran Canaria. This is in agreement with the extension in the E-W direction since 2003, as suggested by [10] after analysing GPS data from two stations situated on La Palma and Gran Canaria. Table 4. Main tectonic change, as derived from the strain analysis, expressed as dilatation or contraction of the triangles (triangle numbers as in Table 2). Units are 10 −9 .

Triangle
Largest Dilatation Triangle Largest Contraction Generally, areas undergoing compression in the archipelago overlap the location of most of the epicentres during this period of study. Most of the seismicity is concentrated in the area between Tenerife and Gran Canaria (Figure 4b). Moderate-to-large (Mw > 6.0) tectonic earthquakes are likely to occur due to the presence of the northeast-southwest trending submarine fault, which is the main tectonic feature of the archipelago [58]. Compression occurs in the central area of Tenerife and to the south of La Gomera and Tenerife. In fact, there is a considerable change in the strain pattern along Tenerife, where maximum compression appears in central Tenerife and moderate extension occurs to the northeast, which is a clear consequence of areas with differential kinematics.
The aforementioned local shortening should cause compression in central Tenerife and may correspond to the ground subsidence in the Teide volcanic edifice. Indeed, the 2004-2005 seismic crisis of Teide (e.g., [59]) was a consequence of magma injection in the northwestern rift of Tenerife that later produced a fluid pressure increase during magma migration and was responsible for the gravity increase detected by [60]. As [15] pointed out, a decrease or stabilization of pressure after such a fluid migration might trigger compressional processes in the central area of Tenerife. However, our analysis spans the years 2011-2017, years after that seismic crisis, and the seismicity (see Figure 1) and monitoring techniques during that period did not show other symptoms of Teide volcano reawakening. In addition, previous studies based on GNSS and InSAR that cover different periods after the 2004 seismic crisis in Tenerife (e.g., [15,16,57]) all refer to subsidence in the central area of the island. The change from the dilation to compression pattern, from north to south along La Gomera Island, exhibits a well-defined E-W orientation. This is consistent with the regional fractures in the strike-slip tectonic regime suggested by [28] and later supported by other authors [10,32,33,61]. More recently, refs. [39,40] suggested that such an orientation could be responsible for the ascent and emplacement of magmas in the early formation phases of the central islands of the archipelago. Moreover, such a deformation pattern with an east-west orientation in the central islands is consistent with the existence of the major tectonic feature that connects the Trans-Alboran Fault System with the Mid-Atlantic Ridge, as suggested by [24] (see Figure 1).
To deepen into the spatial distribution of the seismicity, we compute and map the magnitude of the 2D maximum shear strain for most of the archipelago ( Figure 5). An abrupt change in shear strain (of approximately 100 × 10 −9 ) clearly appears between Tenerife and Gran Canaria along the NE-SW trending fault. This is an important seismogenic area where a submarine reverse fault separates two insular blocks [33]. This area, during the period of observation, was affected by a cluster of earthquake epicentres, ranging from low to moderate magnitude (from m bLG 1.2 to m bLG 3.4).
We used Coulomb software (release 3.4) [62,63] to simulate the horizontal deformation and strain caused by the 1989 m bLg 5.2 earthquake that occurred offshore between Tenerife and Gran Canaria at a focal depth of approximately 15 km. The study of the focal mechanism of the main shock by [33] reported a sinistral strike-slip fault oriented NE-SW, matching the strike of the aforementioned submarine fault [28]. According to Okada's formulae [64], Coulomb 3.4 calculates the 3D component of the displacement field due to a finite rectangular source dislocation buried in an elastic half-space. This formulation approximates the surface movements produced by earthquake faulting, making it appropriate for modelling crustal deformation measured by geodetic techniques. Dilatation is computed from components ε 1 + ε 2 + ε 3 ; that is, by summing the normal strains in the three directions (X,Y,Z), supposing a 3D representation. In our computation, we used the 2D plane strain approximation. Figure 6 maps the dilatation and maximum shear strain (Equation (10)) produced by a rectangular source whose top and bottom range from the surface down to 15 km. The calculation depth is fixed at 4 km b.s.l., but in our simulation, we tested different depths (up to 0 km, seafloor) and observed negligible differences, so Z is constant in the resulting map. Okada's source is embedded in an elastic half-space homogeneously parametrized as having a given Young's modulus (E = 80 GPa), typical of oceanic crust, and Poisson's ratio (σ = 0.25). Table 5 lists the values of the main parameters used in our computation. The resulting fields of both dilatation and maximum shear strain are in good agreement with the triangle strain fields (compare Figure 6a and Figure 4b, Figure 6b and Figure 5). The negative values of maximum shear strain (Figure 6b), which accordingly indicate a sinistral strike-slip, define a patch of the field oriented towards Gran Canaria, trending in the NW-SE direction, similar to the result of the GNSS triangle strain ( Figure 5) that shows two maximum shear strain highs located to the North of Tenerife and to the West of Gran Canaria. This demonstrates that earthquakes cause static stress changes within the lithosphere able to trigger measurable transient deformation of the Earth's surface. GNSS can sense both the spatial and temporal characteristics of pre-and post-seismic deformation. In the far field, geodetic data have been used to characterize both coseismic deformation and viscoelastic relaxation following large strike-slip earthquakes [65], i.e., the long-term seismic cycle. This is often called the stick-slip behaviour of a fault. During the earthquake, the elastic strain built-up in the surrounding crust during the interseismic interval is transformed into permanent slip. Deformation following a large earthquake can continue aseismically for months or even years in the form of viscoelastic relaxation detectable with GNSS networks, even at a large distance from the fault.  The map of maximum shear strain ( Figure 5) displays large shear onshore and offshore Gran Canaria in NE-SW orientation, accounting for a strain of 130 × 10 −9 , with increasing strain generally from North to South. Such a pattern seems to be influenced by the regional stress field (e.g., [32]). The seismicity in that strip is concentrated offshore, on the one hand to the west, very near the seismic swarm related to the submarine fault. On the other hand, seismicity near the east coast of Gran Canaria appears even more scattered, preferentially distributed along the volcano-tectonic line described by [28]. Other areas within the squared zone of Figure 5 exhibit a considerable reduction in the maximum shear strain.
The findings for Tenerife are particularly interesting because the areas of maximum shear strain correlate well with the areas of localized strain (see Figure 4a). The central area of Tenerife mainly includes the volcanic edifice called Las Cañadas, whose upper part is a caldera structure. Las Cañadas was formed as a result of intense explosive volcanic activity, with the Teide-Pico Viejo (TPV) twin stratovolcanoes in the northern sector of the caldera [66][67][68][69]. TPV is one of the highest volcanic complexes on Earth (7200 m from the seafloor) and worthy of study due to its proximity to densely populated areas and proneness to future eruptions. TPV suffered an episode of unrest during 2004-2005, which was interpreted as deep magma injection beneath the volcano in the northwest area that later migrated to the southeast flank. Although no significant ground deformation has been reported, the observed gravity changes correlate well with seismicity [59,60,70,71]. Ref. [72], after reviewing the seismic database from the Spanish National Geographic Institute (http://www.ign.es) from 2003 to 2016, reported that except for the reawakening of 2004-2005 and the isolated seismic swarm of October 2016, the seismicity in Tenerife remained within the background levels. According to Figure 5, a relative maximum of shear strain appears in the central-western area of Tenerife, and another high of greater magnitude (up to 160 × 10 −9 ) is located to the northeast of the island. The relatively high values of shear strain situated in the central-western zone form an elongated area with an E-W orientation, with a cluster of epicentres concentrated along its orthogonal direction. This higher strain area encompasses a high-density crustal structure detected by [73,74] in the southwestern part of the Las Cañadas caldera, which produces a high residual gravity anomaly linked with the Boca Tauce central volcanic edifice. The area between these two highs of shear strain along Tenerife defines a fringe with a slight reduction in shear strain trending NW-SE (approximately 100 × 10 −9 ). This fringe overlaps an important number of epicentres situated in the Dorsal (NE) rift slope and along the Güimar valley. This area is affected by important tectonic activity, and the SE part of this area could be linked to the movements of the active faults situated between the islands of Tenerife and Gran Canaria. Indeed, the authors of [12] have already detected the presence of horizontal ground displacements along this NW-SE axis. In the context of the geodynamics of Tenerife, the SE sector is acknowledged to be an area of intense tectonic activity and anomalous kinematics. The kinematic behaviour of the NW sector of this NW-SE axis was already associated with volcanic activity during the 2004-2005 crisis. In light of the strain regime retrieved from GPS data, the location of the epicentres, the low-density materials and the low P-wave velocities found in that zone [73][74][75]. Ref. [16] proposed a geodynamic scheme that confirms NW Tenerife as more susceptible to volcanic activity than surrounding areas.
The higher values of shear strain in central west Tenerife continue offshore, connecting with another high situated to the north of La Gomera Island, also with an E-W orientation ( Figure 5).
The epicentres follow the volcano-tectonic line reported by [28] that crosses Tenerife and La Gomera at the areas of maximum shear strain, mostly to the north of Tenerife, although seismicity rarely occurs to the northeast of La Gomera.
The high shear strain area with an elongated form is remarkable, and the E-W orientation is observed northeast of Tenerife, extending inland and offshore on both sides of the coast. There are also epicentres falling within this high shear strain area, mostly offshore. Many of these epicentres follow the NW-SE line crossing the island that is coincident with volcanic seamounts and the volcano-tectonic lineation reported by [25]. This line continues west of Gran Canaria, where other high shear strain appears. Indeed, this last high occupies most of the island territory and extends offshore (to the East) with an E-W orientation. The lower limits of the contour map we produced prevent a clearer definition of this high shear strain area that could continue to the south of Gran Canaria. However, the seismicity reported during the period of analysis is located in the area of high shear strain situated to the East, offshore of Gran Canaria, and follows another volcano-tectonic line oriented in the NE-SW orientation proposed by [28].
The analysis and modelling of aeromagnetic anomalies in Tenerife [38] and La Gomera [40] revealed that E-W tectonic features controlled the emplacement of magmas during the early formation of these islands since the most important intrusive complexes formed at the initial stages of their volcanic history are elongated in that direction. This implies that the E-W tectonic lineaments inferred from the present GNSS study represent long-lived tectonic features that have been active during the entire history of the Canary archipelago, from the Miocene to the present. This statement can also be extended to tectonic features with orientations from NE-SW to ENE-WSW between Tenerife and Gran Canaria, first suggested by [28] and later confirmed by other geophysical studies (e.g., [33,39]) and the present work.

Conclusions
From approximately 7 years of GNSS daily solutions (between 2011 and 2017), we reconstructed the present-day pattern of the active deformation in the Canary Archipelago through strain mapping. The 2D infinitesimal strain field was retrieved by computing high-precision geodetic horizontal velocities for 18 GNSS permanent stations, spanning the Canary Islands, and applying the triangular segmentation approach. This has enabled us to characterize major regions of deformation along the archipelago. The extension pattern along the archipelago is clear, with the sectors of maximum deformation found west of Tenerife, Gran Canaria and Fuerteventura. The analysis of the maximum shear strain mostly suggests E-W and NE-SW orientation patterns, suggesting that the regional stress field likely influences such orientations. The sharp change in shear strain between Tenerife and Gran Canaria identifies a sector of intense seismicity, which is mostly associated with the major submarine fault that separates the two insular edifices. On this submarine main tectonic structure, we have produced a tentative simulation of the horizontal deformation and strain caused by one of the strongest (m bLg 5.2) earthquakes of the region. The retrieved dilatation, maximum shear strain and displacement agree with the strain field determined from the GNSS data. In addition, the inferred submarine alignment that crosses Tenerife in the NW-SE direction also concentrates offshore seismicity, which is distributed along areas of higher shear strain elongated in the E-W direction.
We observe a varying kinematic behaviour within the islands, with the Tenerife area being the most active with regard to the regional geodynamics of the entire archipelago. The absence of eruptive activity during the period analysed in this work, except for the El Hierro submarine eruption (October 2011-March 2012) (e.g., [70]), which falls outside the calculated strain maps drawn in this study, the seismic areas between islands, mainly offshore Tenerife and Gran Canaria, seem predominantly influenced by regional tectonic stress rather than by volcanic activity. As suggested by [39], the different episodes of volcanism in the Canary Islands could have occurred in a framework of both extensional and compressional tectonic regimes, and our present results seem to confirm this hypothesis.
Finally, the present work might be useful to improve hazard warnings in active volcanic areas, such as the Canary Islands region; certainly, GNSS networks must be an inherent part of earthquake and/or volcanic eruption warning systems. However, more geophysical and geodetic data of the Canary Archipelago are necessary for producing a complete map of its kinematics and tectonics.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/20/3297/s1, The material shows the configuration of the 23 triangles formed for each triplet of GNSS stations. Figure S1: Triangles formed from 3 non-co-linear stations each, for the 23 GNSS stations used in this work (see Table 1). The horizontal strain was subsequently calculated for each triangle's centroid.