InSAR-Based Mapping to Support Decision-Making after an Earthquake

It has long been recognized that earthquakes change the stress in the upper crust around the fault rupture and can influence the behaviour of neighbouring faults and volcanoes. Rapid estimates of these stress changes can provide the authorities managing the post-disaster situation with valuable data to identify and monitor potential threads and to update the estimates of seismic and volcanic hazard in a region. Here we propose a methodology to evaluate the potential influence of an earthquake on nearby faults and volcanoes and create easy-to-understand maps for decision-making support after large earthquakes. We apply this methodology to the Mw 7.8, 2016 Ecuador earthquake. Using Sentinel-1 Interferometric Synthetic Aperture Radar (InSAR) and continuous GPS data, we measure the coseismic ground deformation and estimate the distribution of slip over the fault rupture. We also build an alternative source model using the Global Centroid Moment Tensor (CMT) solution. Then we use these models to evaluate changes of static stress on the surrounding faults and volcanoes and produce maps of potentially activated faults and volcanoes. We found, in general, good agreement between our maps and the seismic and volcanic events that occurred after the Pedernales earthquake. We discuss the potential and limitations of the methodology.


Introduction
Earthquakes influence the crust around them, changing the stress state and influencing the faults and volcanoes in the vicinity [1,2].This can result in an increase in seismicity and volcanic activity in areas of stress build-up [3].Faults interactions by the static stress transfer have been demonstrated in numerous studies in the last decades, including the occurrence of aftershocks (e.g., [4][5][6][7]), the triggering of earthquakes in nearby faults or segments [8,9] or in the far field related to dynamic triggering [10,11], and the occurrence of volcanic eruptions or episodes of volcanic unrest after an earthquake [12,13].Similarly, volcanic processes can change the stress state in their vicinity, also interacting with active faults.As a consequence of these interactions and changes in the stress state of the crust, large earthquakes modify the seismic and volcanic hazard in their vicinity.
Space geodesy is now routinely used following an earthquake to image the displacement of the ground and estimate the rupture geometry and the distribution of slip over the fault rupture surface (e.g., [14][15][16][17]).Using the obtained source model, it is possible to infer the stress changes on nearby faults and volcanoes produced by the earthquake, which can be used to identify which faults and volcanoes are brought closer to failure or activation.
In recent years, several initiatives that use Earth Observation data to provide support to authorities managing the post-earthquake scenario have developed:

•
The International Charter Space and Major Disasters [18]: an international consortium including resources of fifteen space agencies that quickly provides satellite-derived imagery and supplemental information to any country that requires assistance in response to major disasters.

•
The UNAVCO consortium (http://www.unavco.org/) that supports geodetic networks and provides free and open geodetic data to help with preparedness and mitigation of hazards The ARIA project (https://aria.jpl.nasa.gov/)which is a collaboration between the Jet Propulsion Laboratory (JPL) and Caltech to exploit geodetic and seismic observations in support of hazard response.
The influence of an earthquake in nearby faults and volcanoes can be evaluated through different methods (see for example [19] and references therein).The Coulomb failure stress criterion is one of the most common analyses used to explain fault interactions and seismological observations after earthquakes, such as aftershocks distribution (e.g., [7,8]).
Although these procedures are commonly used today by the academic community, the transfer of these results to the authorities managing the post-earthquake situation is not straightforward and thus its usefulness is reduced in practice.Providing comprehensible maps to the final user is a pending task.
Here we propose a methodology to evaluate the potential influence of an earthquake on nearby faults and volcanoes and create easy-to-understand maps for decision-making support after an earthquake.We apply this methodology to the Mw 7.8, 2016 Ecuador earthquake that occurred on 16 April 2016 at 23:58:36.980UTC (https://earthquake.usgs.gov/earthquakes/eventpage/us20005j32#origin).Using Sentinel-1 SAR and continuous GPS data, we measure the coseismic ground deformation and estimate the distribution of slip.Then we use this model to calculate changes of stress on the surrounding faults and volcanoes.The results are compared with the seismic and volcanic events that have occurred after the earthquake.We discuss the potential and limitations of the methodology and the lessons learnt from discussion with local authorities (the Geological Survey of Ecuador, INIGEMM).

Study Area
The Pedernales earthquake occurred in the Ecuadorian subduction zone, where the Nazca and South American plates converge at 55 mm/yr [20,21].The earthquake of Mw 7.8 occurred at 20.6 ± 3.2 km depth (Figure 1) and the focal mechanism is consistent with the rupture of the subduction interface between the Nazca and South American plates.The studies published to date reveal that the Pedernales earthquake ruptured a 100-120 km long segment of the subduction interface and consisted of two main asperities [22,23].

Study Area
The Pedernales earthquake occurred in the Ecuadorian subduction zone, where the Nazca and South American plates converge at 55 mm/yr [20,21].The earthquake of Mw 7.8 occurred at 20.6 ± 3.2 km depth (Figure 1) and the focal mechanism is consistent with the rupture of the subduction interface between the Nazca and South American plates.The studies published to date reveal that the Pedernales earthquake ruptured a 100-120 km long segment of the subduction interface and consisted of two main asperities [22,23].Reference map of our study area in Ecuador (red box in the inset map).The relative Nazca-South American convergence rate and direction are shown with a black arrow [20] and the location of the trench is shown with a heavy black barbed line.The black dashed box shows the region in Figures 2 and 4. The location of the continuous GPS sites used in this study is shown as blue inverted triangles.The Global Centroid Moment Tensor (GCMT) focal mechanism of the 2016 Pedernales earthquake is shown.The topography and bathymetry are from the ETOPO1 1 arc-minute global relief model [24].The green segment indicates the approximate rupture area of the 1906 earthquake [25] and the red semi-transparent areas represent the rupture areas of the 1942, 1958, and 1979 earthquakes [26][27][28][29].The approximate rupture area of the Mw 7.8 2016 Pedernales earthquake from our variable-slip model is also shown (semi-transparent blue area).Reference map of our study area in Ecuador (red box in the inset map).The relative Nazca-South American convergence rate and direction are shown with a black arrow [20] and the location of the trench is shown with a heavy black barbed line.The black dashed box shows the region in Figures 2  and 4. The location of the continuous GPS sites used in this study is shown as blue inverted triangles.The Global Centroid Moment Tensor (GCMT) focal mechanism of the 2016 Pedernales earthquake is shown.The topography and bathymetry are from the ETOPO1 1 arc-minute global relief model [24].The green segment indicates the approximate rupture area of the 1906 earthquake [25] and the red semi-transparent areas represent the rupture areas of the 1942, 1958, and 1979 earthquakes [26][27][28][29].The approximate rupture area of the Mw 7.8 2016 Pedernales earthquake from our variable-slip model is also shown (semi-transparent blue area).
The relative movement of the plates produces M > 7.5 earthquakes (Figure 1) such as the Mw 8.8 1906 Ecuador-Colombia earthquake, which ruptured a large region of 500 km of the subduction interface [30], and the sequence of earthquakes that partially ruptured the same area between 1942 and 1979 (Figure 1).The 2016 Pedernales earthquake seems to have ruptured a similar area to the 1942 earthquake [29].Apart from megathrust earthquakes, the subduction process is responsible for the uplift of the Andean range [31], crustal faults on the upper plate, and Quaternary volcanism [32].
Continental Ecuador is divided into several morphostructural regions from west to east: the Coastal plain, the Cordillera Occidental, the Cordillera Oriental, and the upper Amazon basin.Active tectonics is dominated by large faults [33] that separate different tectonic regimes including the coastal region, where the subduction of the Carnegie Ridge and the oblique subduction control active structures; the North Andean Block, limited by a NNE dextral strike-slip fault system; and the Eastern Andean Frontal fault zone that limits the North Andean block along the Subandean region.These crustal faults are capable of producing earthquakes, such as the 1987 Ms 6.9 earthquake [34].

GPS Displacements
We processed daily RINEX (Receiver Independent Exchange) data from 28 cGNSS stations from the Continuous Monitoring Global Navigation Satellite System (GNSS) Network (REGME) in Ecuador.We estimated the site coordinates before and after the mainshock of the 16 April 2016, Mw7.8 Pedernales Earthquake by using Bernese GPS Software version 5.0 [35], and calculated coseismic displacements by taking the differences between the daily site coordinates.The data were processed on a daily basis using the Bernese Processing Engine (BPE), which applies a double-difference processing strategy.We applied International GNSS Service (IGS) precise orbits and earth orientation parameters, with an absolute antenna phase centre and a FES2004 ocean tide-loading model [36].Only observations made above an elevation cut-off angle of 5 • were used to estimate parameters, and an elevation-dependent weighting was applied.The tropospheric effect was modelled on a prior dry-Niell [37] model and completed by estimating zenith delay corrections for each site at 1-hour intervals using the wet-Niell mapping function.The ambiguity resolution is based on the Quasi-Ionosphere-Free (QIF) baseline-wise analysis.Baselines were constructed based on a strategy of maximum observations.
In order to establish a link with the global reference frame ITRF2008 (IGb08), we use 17 continuously recording GPS stations from the IGS network.Once we obtained the daily solutions and coordinates of REGME stations, we imposed "a posteriori" constraints on implementing the chosen reference frame.The coordinates of the IGS stations were constrained (NNT-No Net Translation option) to their ITRF2008 values.

InSAR Displacements
We processed Interferometric Synthetic Aperture Radar (InSAR) data from the European Space Agency (ESA) Sentinel-1 satellite to estimate the coseismic line-of-sight (LOS) ground motion.Using the Sentinel-1 processor developed at the Centre Tecnològic de les Telecomunicacions de Catalunya (http://geomatics.cttc.es/remote-sensing/),we combined 2 Sentinel-1 images from descending relative orbit 40, with dates 12 April 2016 and 24 April 2016.The 3-arcsec resolution NASA Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model [38] was used to remove the topographic phase contribution.The interferogram was filtered using a power spectrum filter [39] and unwrapped using the Minimum Cost Flow (MCF) method [40].The interferogram was geocoded and rectified with a 90 m resolution.
The interferogram was flattened using the GPS displacements as a reference (e.g., [14,41,42]).We first verified that GPS and InSAR data record essentially the same surface displacement, i.e., postseismic deformation in the interferogram is not significant.Although GPS time series covering the same postseismic period included in the interferogram (16-24 April 2016) were not available, we have GPS position time series during a period ranging between 2 and 6 days after the earthquake in most of the GPS network (Supplementary Table S1).In these stations, the difference between GPS coseismic displacements and GPS coseismic displacements plus postseismic displacements is less that 1 cm, which is even smaller when projected in the LOS direction.We can therefore assume that our InSAR and GPS data record essentially the same deformation field.We fitted a linear ramp to the difference between the GPS displacements recorded at the 7 GPS stations that lie inside the region covered by the interferogram (projected onto the local LOS vector) and the InSAR displacements at the same locations.The root mean square (RMS) difference between both datasets is 1.07 cm (Supplementary Table S2).

Modeling the Earthquake Source (Inversion Method)
We estimate the distribution of slip during the Pedernales earthquake from a linear inversion of the InSAR and GPS displacement field, modelling the earthquake as a dislocation in an elastic medium [43].
We first downsample the InSAR data using a quadtree approach (e.g., [44]) which samples the data according to their spatial variance, leaving more data in regions of high phase gradient near the source.After downsampling, the original ~235,000 data are reduced to ~1000 points.For the inversion, we use the 17 GPS stations closest to the epicentre (Figure 2 and Table 1).
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 24 most of the GPS network (Supplementary Table S1).In these stations, the difference between GPS coseismic displacements and GPS coseismic displacements plus postseismic displacements is less that 1 cm, which is even smaller when projected in the LOS direction.We can therefore assume that our InSAR and GPS data record essentially the same deformation field.We fitted a linear ramp to the difference between the GPS displacements recorded at the 7 GPS stations that lie inside the region covered by the interferogram (projected onto the local LOS vector) and the InSAR displacements at the same locations.The root mean square (RMS) difference between both datasets is 1.07 cm (Supplementary Table S2).

Modeling the Earthquake Source (Inversion Method)
We estimate the distribution of slip during the Pedernales earthquake from a linear inversion of the InSAR and GPS displacement field, modelling the earthquake as a dislocation in an elastic medium [43].
We first downsample the InSAR data using a quadtree approach (e.g., [44]) which samples the data according to their spatial variance, leaving more data in regions of high phase gradient near the source.After downsampling, the original ~235,000 data are reduced to ~1000 points.For the inversion, we use the 17 GPS stations closest to the epicentre (Figure 2 and Table 1).We build the geometry of the subduction interface based on the global subduction zone model of [45] at the latitude of the Pedernales earthquake.We divide it into an array of 450 elements, each measuring 10 km along strike and 10 to 14 km along dip, with a variable dip that increases with depth from 10 • to 30 • .We jointly invert the InSAR and GPS datasets to solve for the slip distribution along the 450 patches using a least-squares minimization with the non-negativity constraint on the slip.We impose the rake of 123 • (from the focal mechanism, http://earthquake.usgs.gov/earthquakes/eventpage/us20005j32#origin).To limit oscillations of the solution, we impose some smoothing on the solution, by minimizing the second-order derivative of the fault slip (e.g., [46]).The linear inversion includes also the weight matrix of each dataset: we tested different weighting to find a compromise between the RMS and the spatial density of the InSAR and GPS dataset.The optimal solution roughness for the final model is selected from the trade-off curve between misfit and solution roughness.

Estimation of Static Stress Changes
To estimate the static stress changes induced by the Pedernales earthquake we use Coulomb failure stress change (∆CFS) analysis (e.g., [8,47]).This takes into account the shear and normal stress changes (∆τ and ∆σ) induced by the earthquake on selected surfaces or receivers, such as faults or dikes, considering an apparent friction coefficient (µ ) to infer the total coulomb stress, using the equation: The Coulomb stress change was calculated using Coulomb 3.4 [48], assuming an elastic isotropic half-space characterized by a Shear modulus of 30 GPa and Poisson ratio of 0.25.The apparent (or effective) friction coefficient used for the receiver surfaces is 0.4.Two different source faults are used as input for our Coulomb stress calculations: on the one hand, the variable-slip source model obtained from the joint inversion of InSAR and GPS displacements; on the other hand, the source from the Global CMT focal mechanism with rupture dimensions obtained from [49] empirical relations.The interest of using these two sources is to compare the results using a faster but less well constrained solution (the Global CMT source) with the results obtained using a better-constrained source model from InSAR and GPS data (see Section 5.1).We use the coulomb criteria to estimate the variation of stress induced by the Pedernales earthquake on four different receivers:

•
Potential continental faults: we first select the faults from the catalogue of active faults in Ecuador of [50].We classify the main sets according to the type of fault: normal, reverse, and strike-slip faults.Additionally, we plot the frequency distribution of fault orientation in rose diagrams.In total, the database shows 278 faults, indicating the fault length, slip direction, plunge of the slip vector, and slip type, among others parameters.For each subset of faults we extract average values of strike, dip, and rake, which are used as receiver fault orientations.We compute stress changes on each fault subset at 10 km depth, which is the common depth for the maximum cortical elastic strength, and therefore the depth for the nucleation of the more relevant earthquakes.

•
Known fault traces: from the active faults database of [50] we extracted the fault traces and obtained points spaced 1 km apart along them.Then we estimate the Coulomb failure stress change induced by the Pedernales earthquake on those points with the correspondent potential rupture orientation characteristics (strike, dip, rake), at 5 km depth (the depth of the approximate centre of the modelled fault surfaces).

•
Subduction interface: Using the geometry of the subduction interface from [45] we extract the strike, dip, and rake of each node and estimate the Coulomb failure stress change induced by the Pedernales earthquake on each node.

•
Volcanoes: we used the Ecuadorian volcanoes with historic activity from The Global Volcanism Program database (http://volcano.si.edu/search_volcano.cfm).We select a set of five currently active volcanoes which have been extensively studied, and extract the orientation of their magma paths, using the criteria of [51] and references therein.We then estimate an average orientation that is used to calculate the induced normal stress change (∆σ) on the volcano magma pathway.
Positive values on a particular volcano indicate unclamping produced by the 2016 earthquake, which indicates dilatation in the crust at the depth of the magma chamber.In this case, the magma pathway is affected by a normal stress reduction or unclamping and this might promote new volcanic eruptions (e.g., [52]).We model the normal stress change at different depths, ranging from 0 to 10 km.

Final Maps
To produce the final maps we simplify the different maps obtained from the stress change analysis to produce more easy-to-understand maps.We use a traffic light colour coding (red, yellow, green) with three possible values to classify the volcanic and tectonic elements according to their level of potential activation by the Pedernales earthquake.We produce four different maps: • Map of potentially activated volcanoes: We classify the final ∆σ values obtained for the average orientation of the volcanoes magma paths in three levels of potential activation and colour the final map according to these levels: green (low level) indicates volcanoes with ∆σ < 0.01 bar, yellow (medium level) indicates volcanoes with 0.01 bar ≤ ∆σ ≤ 0.1 bar, and red (high level) indicates volcanoes with ∆σ > 0.1 bar.The difference in ∆σ values at depths of 0-10 km is almost identical due to the distance between the Pedernales earthquake source and the volcanic arc (~300 km), and therefore we only represent the value estimated at a depth of 10 km.

•
The values used to separate the three intervals of potential seismic or volcanic activation were chosen considering the values published in the literature: in the case of seismic triggering of faults, the lower value of 0.1 bar of ∆CFS is frequently used as threshold to show areas where aftershocks are triggered (e.g., [53][54][55]), while the higher value (1 bar) is typically related to clear earthquake triggering relations (e.g., [56,57]).In the case of volcanic triggering by unclamping of the magmatic chamber, variations of 0.1 bar have been demonstrated to be influential on the eruptive processes, while the lower value (0.01 bar) is a reasonable threshold to infer a moderate influence from a conservative hazard estimation point of view [58][59][60][61].Note that a global agreement about the exact values of stress separating triggering from not triggering does not exist and therefore the exact values for our intervals were determined by expert judgment based on our experience and bibliography.These values are subject to modification if the final users consider it appropriate from their experience using those maps.

GPS and InSAR Displacements
The coseismic interferogram (Figure 2) shows a lobe of ground deformation moving away from the satellite, which is located to the left of the image (white arrow labelled LOS in Figure 2).Maximum displacement is ~60 cm of movement away from the satellite.This interferogram covers a period of 12 days (from 12-24 April 2016.), and thus it can also include, apart from deformation associated with the Pedernales mainshock, deformation that occurred during the 8 days after the earthquake.
Position time series of the 17 cGPS sites used in this study are shown in Figure 3. Continuous GNSS data recorded displacements with almost 70 cm of horizontal trenchward motion and subsidence of more than 17 cm near the town of Pedernales (PEEC station, see Figure 2).The sites located near the coast also show subsidence, with areas of uplift found at inland stations (QVEC, ABEC, CXEC, QUEM).
The displacement field observed with both InSAR and GPS data is consistent with onshore ground deformation commonly observed in subduction earthquakes, characterized by subsidence and horizontal trenchward motion (ground experiencing uplift is closer to the trench and often under the ocean, see Supplementary Figure S1).

Earthquake Source from InSAR and GPS Data
According to our best model, the Pedernales earthquake ruptured a 100-km-long segment of the subduction interface and most of the slip (slip larger than 2 m) occurred between depths of 13 and 30 km (Figure 4).The equivalent geodetic moment of our model is 4.94 × 10 20 Nm, considering a shear modulus of 3.3 × 10 10 Pa, and 5.98 × 10 20 Nm, considering a shear modulus of 4 × 10 10 Pa, which is consistent with seismological estimates (Mw = 7.8).Slip maxima (asperities) locate in two main areas separated by a lower slip region (<2 m), in agreement with previously published models [28,29,62].Supplementary Figures S2 and S3 show GPS and InSAR residuals corresponding to this model.The optimal solution roughness for the final model is shown in Supplementary Figure S4.

Faults and Volcanoes Selected for This Study
The result of the analysis of active faults of Ecuador is shown in Figure 5.The 278 faults of the [50] catalogue can be divided into six families: two sets of normal faults, oriented N45°E and N320°E; two sets of reverse faults, oriented N350°E and N15°E; right-lateral strike slip faults, N60°E trending; and left-lateral strike slip faults that exhibit a conjugate set, N80°E and N330°E (Supplementary Table S3).Furthermore, the total length of each fault type is a relevant parameter to be considered in this analysis.In the case of reverse faulting, fault set with N15°E trending is the 85% of the reverse faulting for fault length larger than 10 km.The compilation of information on the selected volcanoes (Tungurahua, Guagua Pichincha, Reventardor, Sangay and Cotopaxi) revealed that most of them have vertical or sub-vertical magma paths, with NE-SW azimuths in most of them (see Supplementary Table S4).We use this average orientation to estimate the induced normal stress change (∆σ) on each volcano magma pathway.

Faults and Volcanoes Selected for This Study
The result of the analysis of active faults of Ecuador is shown in Figure 5.The 278 faults of the [50] catalogue can be divided into six families: two sets of normal faults, oriented N45 • E and N320 • E; two sets of reverse faults, oriented N350 • E and N15 • E; right-lateral strike slip faults, N60 • E trending; and left-lateral strike slip faults that exhibit a conjugate set, N80 • E and N330 • E (Supplementary Table S3).Furthermore, the total length of each fault type is a relevant parameter to be considered in this analysis.In the case of reverse faulting, fault set with N15 • E trending is the 85% of the reverse faulting for fault length larger than 10 km.

Faults and Volcanoes Selected for This Study
The result of the analysis of active faults of Ecuador is shown in Figure 5.The 278 faults of the [50] catalogue can be divided into six families: two sets of normal faults, oriented N45°E and N320°E; two sets of reverse faults, oriented N350°E and N15°E; right-lateral strike slip faults, N60°E trending; and left-lateral strike slip faults that exhibit a conjugate set, N80°E and N330°E (Supplementary Table S3).Furthermore, the total length of each fault type is a relevant parameter to be considered in this analysis.In the case of reverse faulting, fault set with N15°E trending is the 85% of the reverse faulting for fault length larger than 10 km.The compilation of information on the selected volcanoes (Tungurahua, Guagua Pichincha, Reventardor, Sangay and Cotopaxi) revealed that most of them have vertical or sub-vertical magma paths, with NE-SW azimuths in most of them (see Supplementary Table S4).We use this average orientation to estimate the induced normal stress change (∆σ) on each volcano magma pathway.The compilation of information on the selected volcanoes (Tungurahua, Guagua Pichincha, Reventardor, Sangay and Cotopaxi) revealed that most of them have vertical or sub-vertical magma paths, with NE-SW azimuths in most of them (see Supplementary Table S4).We use this average orientation to estimate the induced normal stress change (∆σ) on each volcano magma pathway.

Results of the Static Stress Analysis
We have produced 8 maps of static stress changes corresponding to the 6 crustal faults families, the subduction interface, the average stress increment taking into account all the fault families and the Ecuadorian volcanoes (Figure 6).We have produced the same set of maps using the Global CMT focal mechanism with rupture dimensions obtained from [49] empirical relations (Supplementary Figure S5).
The results of the stress change analysis on the active faults of Ecuador reveal different patterns of ∆CFS depending on the orientation of the receiver faults, the receiver fault orientation being one of the most influential parameters on the stress change results.The most extended lobes of ∆CFS correspond to faults with orientations N45-85, 180 (right hand rule orientation and slip rake), right lateral NE-SW strike slip faults (Figure 6b); and N230-75, -90, NE-SW normal faults (Figure 6d); while smaller lobes correspond to faults N15-35, 90 (Figure 6a) and its conjugate set N195-35, -90 (Figure 6c), NNE-SSW thrust faults.Towards the south of the rupture the families N110-75, -90 (E-W normal faults, Figure 6e) and N320-89, 0 (NW-SE left-lateral strike slip faults, Figure 6f) are stressed.Negative values of ∆CFS (i.e., inhibition of fault movement) occur mainly offshore and around Muisne.
The results of the stress change analysis on the subduction interface (Figure 6g) show positive values of ∆CFS downdip and updip of the rupture.Negative values of ∆CFS are located on areas that slipped during the earthquake.
The results of the stress analysis on the Ecuadorian volcanoes show that the 2016 earthquake increased the normal stress by ≥0.1 bar on 3 of the volcanoes (numbers1, 2, 3 in Figure 6h), with a maximum value of ∆σ of 0.162 bars at Tungurahua volcano (number 2 in Figure 6h).The volcanoes further north than Cayambe (number 9 in Figure 6h) show negative ∆σ, indicating that the 2016 earthquake induced a normal stress increase on the magma pathway of those volcanoes, which might hinder magmatic activity.

Results of the Static Stress Analysis
We have produced 8 maps of static stress changes corresponding to the 6 crustal faults families, the subduction interface, the average stress increment taking into account all the fault families and the Ecuadorian volcanoes (Figure 6).We have produced the same set of maps using the Global CMT focal mechanism with rupture dimensions obtained from [49] empirical relations (Supplementary Figure S5).
The results of the stress change analysis on the active faults of Ecuador reveal different patterns of ∆ depending on the orientation of the receiver faults, the receiver fault orientation being one of the most influential parameters on the stress change results.The most extended lobes of ∆ correspond to faults with orientations N45-85, 180 (right hand rule orientation and slip rake), right lateral NE-SW strike slip faults (Figure 6b); and N230-75, -90, NE-SW normal faults (Figure 6d); while smaller lobes correspond to faults N15-35, 90 (Figure 6a) and its conjugate set N195-35, -90 (Figure 6c), NNE-SSW thrust faults.Towards the south of the rupture the families N110-75, -90 (E-W normal faults, Figure 6e) and N320-89, 0 (NW-SE left-lateral strike slip faults, Figure 6f) are stressed.Negative values of ∆ (i.e., inhibition of fault movement) occur mainly offshore and around Muisne.
The results of the stress change analysis on the subduction interface (Figure 6g) show positive values of ∆ downdip and updip of the rupture.Negative values of ∆ are located on areas that slipped during the earthquake.
The results of the stress analysis on the Ecuadorian volcanoes show that the 2016 earthquake increased the normal stress by ≥0.1 bar on 3 of the volcanoes (numbers1, 2, 3 in Figure 6h), with a maximum value of ∆σ of 0.162 bars at Tungurahua volcano (number 2 in Figure 6h).The volcanoes further north than Cayambe (number 9 in Figure 6h) show negative ∆σ, indicating that the 2016 earthquake induced a normal stress increase on the magma pathway of those volcanoes, which might hinder magmatic activity.

Final Maps
Following the criteria defined in Section 3.5, we have produced four final maps: potentially activated mapped faults (Figure 7a), potentially activated continental faults (Figure 7b), potentially activated areas on the subduction interface (Figure 8), and potentially activated volcanoes (Figure 9).
The map of potentially activated mapped faults shows that faults with a high level of potential activation (red segments in Figure 7a) are located in the Coastal Plain (region of Manabi) between latitudes −1.5 and 0°N.These faults present values of ∆ > 1 bar induced by the Pedernales earthquake.
The map of potentially activated continental faults shows a main lobe of Coulomb Stress increment towards the southeast of the Pedernales rupture (Figure 7b), and a minor lobe towards the northeast, the later caused mainly by the stress increment on NNE-SSW thrust faulting.
The map of potentially activated areas on the subduction interface shows areas of the slab with a high level of potential activation (red areas in Figure 8) located around the patches of coseismic slip.The regions of the slab with a low level of potential activation (green areas in Figure 8) are found both on areas of coseismic slip and also away from the rupture.
The map of potentially activated volcanoes shows that volcanoes with a high level of potential activation (volcanoes inside the red area, i.e., numbers 2 and 3 in Figure 9) are located to the

Final Maps
Following the criteria defined in Section 3.5, we have produced four final maps: potentially activated mapped faults (Figure 7a), potentially activated continental faults (Figure 7b), potentially activated areas on the subduction interface (Figure 8), and potentially activated volcanoes (Figure 9).
The map of potentially activated mapped faults shows that faults with a high level of potential activation (red segments in Figure 7a) are located in the Coastal Plain (region of Manabi) between latitudes −1.5 and 0 • N.These faults present values of ∆CFS > 1 bar induced by the Pedernales earthquake.
The map of potentially activated continental faults shows a main lobe of Coulomb Stress increment towards the southeast of the Pedernales rupture (Figure 7b), and a minor lobe towards the northeast, the later caused mainly by the stress increment on NNE-SSW thrust faulting.
The map of potentially activated areas on the subduction interface shows areas of the slab with a high level of potential activation (red areas in Figure 8) located around the patches of coseismic slip.The regions of the slab with a low level of potential activation (green areas in Figure 8) are found both on areas of coseismic slip and also away from the rupture.
The map of potentially activated volcanoes shows that volcanoes with a high level of potential activation (volcanoes inside the red area, i.e., numbers 2 and 3 in Figure 9) are located to the southeast of the rupture while volcanoes with a low level of potential activation (volcanoes inside the green area, i.e., numbers 8-12 in Figure 9) are those located further north than latitude −0.17 • N.
These maps should be interpreted in terms of increased level of potential activation considering only the static stress changes produced by the Pedernales mainshock, but the occurrence of new seismic and volcanic events depends on the state of the system prior to the earthquake.Note that other phenomena not considered in this analysis may also influence the stress state of the crust and trigger or inhibit new seismic and volcanic events (see Section 5.2 and 5.3).

Models Using the Focal Mechanism
We performed our calculations of static Coulomb failure stress changes associated with the Pedernales earthquake slip using as an alternative source the Global CMT focal mechanism.The interest of using such source is the promptness with which focal mechanisms are available (even a few minutes after the earthquake) compared to the source derived from InSAR and GPS data (days to weeks).Figure 10a shows the map of potentially activated mapped faults of Ecuador calculated using the earthquake source from the Global CMT.Compared with the map created using our source (Figure 7), the level of potential activation on mapped faults is lower using the CMT source and fault segments with a high level of potential activation (red segments) are not the same in both maps.
We have computed the difference between both maps to identify their main similarities and differences (Figure 10b).The CMT model produces very similar results in faults far from the earthquake source (white segments in Figure 10b).For closer receivers the variation is, however, significant: for example, the fault labelled X 1 in Figures 7 and 10a shows a low level of potential activation in the map created using our source and a high level of potential activation in the map created using the CMT source.A similar case is the segment labelled X 2 in Figures 7 and 10a, with low values in the CMT map and high values in the map derived from our source.This effect is common on every modelled earthquake, in the near field the details and particularities of the slip distribution is one of the key parameters that define the stress change field [63]; while in the far field the slip details are less relevant and the approximate dimensions of the source and the earthquake magnitude is enough to define a valid stress change field.
Since the static stress change is highly dependent on the slip distribution, the more data we can add for the inversion of the source, the more reliable will be our final model.Using the focal mechanism to produce a preliminary version of the maps can allow fast results after the earthquake, but they should be updated using a better constrained source model when new observations, particularly from seismological and geodetic networks that accurately observe the earthquake process, are available.southeast of the rupture while volcanoes with a low level of potential activation (volcanoes inside the green area, i.e., numbers 8-12 in Figure 9) are those located further north than latitude −0.17°N.
These maps should be interpreted in terms of increased level of potential activation considering only the static stress changes produced by the Pedernales mainshock, but the occurrence of new seismic and volcanic events depends on the state of the system prior to the earthquake.Note that other phenomena not considered in this analysis may also influence the stress state of the crust and trigger or inhibit new seismic and volcanic events (see Sections 5.2 and 5.3).

Models Using the Focal Mechanism
We performed our calculations of static Coulomb failure stress changes associated with the Pedernales earthquake slip using as an alternative source the Global CMT focal mechanism.The interest of using such source is the promptness with which focal mechanisms are available (even a few minutes after the earthquake) compared to the source derived from InSAR and GPS data (days to weeks).Figure 10a shows the map of potentially activated mapped faults of Ecuador calculated using the earthquake source from the Global CMT.Compared with the map created using our source (Figure 7), the level of potential activation on mapped faults is lower using the CMT source and fault segments with a high level of potential activation (red segments) are not the same in both maps.
We have computed the difference between both maps to identify their main similarities and differences (Figure 10b).The CMT model produces very similar results in faults far from the earthquake source (white segments in Figure 10b).For closer receivers the variation is, however, significant: for example, the fault labelled X1 in Figures 7 and 10a shows a low level of potential activation in the map created using our source and a high level of potential activation in the map created using the CMT source.A similar case is the segment labelled X2 in Figures 7 and 10a, with low values in the CMT map and high values in the map derived from our source.This effect is common on every modelled earthquake, in the near field the details and particularities of the slip distribution is one of the key parameters that define the stress change field [63]; while in the far field the slip details are less relevant and the approximate dimensions of the source and the earthquake magnitude is enough to define a valid stress change field.
Since the static stress change is highly dependent on the slip distribution, the more data we can add for the inversion of the source, the more reliable will be our final model.Using the focal mechanism to produce a preliminary version of the maps can allow fast results after the earthquake, but they should be updated using a better constrained source model when new observations, particularly from seismological and geodetic networks that accurately observe the earthquake process, are available.S7 for the map calculated using the earthquake source from the Global CMT focal mechanism.

Comparison of Our Final Maps with Seismic and Volcanic Events
We have compared seismic events with magnitudes between 2 and 6.9 in the period 16 April 2016-18 February 2018 with our final maps.Using the seismic events from Instituto Geofísico Escuela Politécnica Nacional we have separated events that occurred in the subduction interface (selecting events at a distance ≤10 km from the slab) from events that occurred in continental faults (selecting events with distance to slab >10 km and depth <15 km). Figure 11a,b show the events that occurred in continental faults and in the subduction interface, respectively, coloured according to the traffic light colour coding in Figures 7b and 8. Red circles represent seismic events located in areas with a high level of potential activation (∆CFS > 1 bar), yellow circles represent seismic events in areas with a medium level of potential activation (0.1 bar ≤ ∆CFS ≤ 1), and green circles represent seismic events in areas with a low level of potential activation (∆CFS < 0.1 bar).We have quantified the number of events located on each region (histograms in Figure 11 c and d) for events that occurred on the subduction interface (Figure 11d), 50% of the events occurred in areas with a high or medium level of potential activation.For events that occurred at continental faults (Figure 11c), only 35% of events occurred in areas with a high or medium level of potential activation.This can be explained by different factors: on the one hand, the seismic events analysed are not relocated, and thus the final location of the events can vary from the one shown in Figure 11a,b.Moreover, our models only consider static stress changes induced by the Pedernales earthquake, but other significant seismic events have occurred in the period April 2016-February 2018 that have also changed the stress state of the crust.For instance, the three largest events that have occurred after the Pedernales earthquake (Supplementary Figure S8) with Mw 6.3, 6.7, and 6.9, located near Muisne, to  S7 for the map calculated using the earthquake source from the Global CMT focal mechanism.

Comparison of Our Final Maps with Seismic and Volcanic Events
We have compared seismic events with magnitudes between 2 and 6.9 in the period 16 April 2016-18 February 2018 with our final maps.Using the seismic events from Instituto Geofísico Escuela Politécnica Nacional we have separated events that occurred in the subduction interface (selecting events at a distance ≤10 km from the slab) from events that occurred in continental faults (selecting events with distance to slab >10 km and depth <15 km). Figure 11a,b show the events that occurred in continental faults and in the subduction interface, respectively, coloured according to the traffic light colour coding in Figures 7b and 8. Red circles represent seismic events located in areas with a high level of potential activation (∆CFS > 1 bar), yellow circles represent seismic events in areas with a medium level of potential activation (0.1 bar ≤ ∆CFS ≤ 1), and green circles represent seismic events in areas with a low level of potential activation (∆CFS < 0.1 bar).We have quantified the number of events located on each region (histograms in Figure 11 c,d) for events that occurred on the subduction interface (Figure 11d), 50% of the events occurred in areas with a high or medium level of potential activation.For events that occurred at continental faults (Figure 11c), only 35% of events occurred in areas with a high or medium level of potential activation.This can be explained by different factors: on the one hand, the seismic events analysed are not relocated, and thus the final location of the events can vary from the one shown in Figure 11a,b.Moreover, our models only consider static stress changes induced by the Pedernales earthquake, but other significant seismic events have occurred in the period April 2016-February 2018 that have also changed the stress state of the crust.For instance, the three largest events that have occurred after the Pedernales earthquake (Supplementary Figure S8) with Mw 6.3, 6.7, and 6.9, located near Muisne, to the north of the main shock rupture, may have triggered many of the seismic events located in that area in Figure 11a,b.Additionally, other processes apart from static stress changes, such as dynamic stresses or postseismic phenomena, may influence the stress of the crust.This will be discussed in Section 5.3.
the north of the main shock rupture, may have triggered many of the seismic events located in that area in Figure 11a,b.Additionally, other processes apart from static stress changes, such as dynamic stresses or postseismic phenomena, may influence the stress of the crust.This will be discussed in Section 5.3.If we only consider the largest events that occurred in the region after the Pedernales earthquake (Mw ≥ 6, Supplementary Figure S8) during the same period, the correlation with our maps is clearer.Since the location and focal mechanism of these events is consistent with the subduction interface between the Nazca and the South American plates, we have compared their location with the map of Coulomb failure stress change on the subduction interface (black circles in Figure 6g).Most of the aftershocks nucleated in areas of increased ∆CFS surrounding the mainshock slip patches.Overall, more than 87% of the Mw ≥ 6 aftershocks (7 of 8 events) occurred in areas of ∆CFS > 0.1 bar.We also compared the location of these events with the map of potentially activated areas on subduction interface (black circles in Figure 8).Fifty percent of the events occurred in areas with a high level of potential activation, 37.5% of the events occurred in areas with a medium level of potential activation, and 12.5% of the events occurred in areas with a low level of potential activation.Our map is therefore useful to constrain the areas where the largest aftershocks will nucleate after the mainshock.
Regarding volcanic activity in Ecuador, according to Instituto Geofisico (http://www.igepn.edu.ec/informes-volcanicos) at least two volcanoes in the Galapagos Islands have shown volcanic activity or periods of unrest after the Pedernales earthquake: the Fernandina volcano, in eruption since September 2017, and the Sierra Negra volcano, that shows seismic activity since July 2017 and an uplift rate of 70 cm/yr (http://www.igepn.edu.ec/servicios/noticias/1539-informe-especial-sierranegra-n-2-2017).However, Galápagos Islands are more than 1000 km away from the Pedernales earthquake source and thus the static stress changes induced by the Pedernales earthquake on these volcanoes are too small to trigger volcanic activity, although dynamic stresses associated to the elastic waves propagation could be responsible.In continental Ecuador, two volcanoes present a high level of potential activation (Cotopaxi and Tungurahua) and Sangay volcano is very close to the boundary of the high level too (Figure 9).Using data from https://volcano.si.edu/volcano.cfm?vn=352090 and https://www.volcanodiscovery.com,information about the volcanic activity of these volcanoes during the month after the earthquake (including ash emission, presence of thermal anomalies, and significant seismic activity) was collected.

•
Sangay volcano awoke one month before the earthquake (26 April 2016), according to reports about ash emission and thermal anomaly.Ash emission was ongoing during at least two months, and the thermal anomaly was recorded in the same period (12 May 2016).After that, both ash emissions and thermal anomalies were intermittent through July 2016.

•
In the vicinity of Cotopaxi volcano, two earthquakes of magnitude 3 and 3.1 occurred at a distance of 15 km from the volcano on 17 April 2016, the day after the Pedernales earthquake.

•
Near Tungurahua volcano, five earthquakes were recorded between 17 April 2016 and 19 May 2016, with magnitudes ranging between 2.8 and 3.2 and distances between 21 to 33 km away from the volcano.
Therefore, significant volcanic activity seems to have increased in these volcanoes around the date of the earthquake, suggesting a potential relationship between both phenomena.

Potential and Limits of the Methodology
Nowadays several SAR satellites such as Sentinel-1A/B, ALOS-2, CSK, TSK, allow the production of coseismic interferogramas relatively quickly (from hours to a few days) after an earthquake.The new European Space Agency (ESA) Sentinel-1 constellation is particularly well adapted to this task, since it acquires images globally every 6-12 days, that are distributed free of charge.In the coming years, several initiatives now in the testing phase to produce systematic interferograms should be fully operational (e.g., GEP, LicS, ARIA, EPOS).
Coseismic GPS displacements are not globally available and their accessibility depends on different factors, including the previous existence of a permanent and non-permanent GPS networks in the zone.Although more and more countries have a continuous GPS network, the data are rarely public.However, some international initiatives exist to give access to GPS data, such as UNAVCO consortium (http://www.unavco.org/) that provides Open Access GPS data with daily RINEX files suitable for processing, although they mostly cover the United States and Caribbean zones.Similarly, the European Plate Observing System will soon provide access to GNSS data in real time over Europe (https://www.epos-ip.org/tcs/gnss-data-and-products).The availability of real-time High-Rate Global Positioning System (HRGPS) data can be used to produce rapid magnitude estimation and produce rapid earthquake source models (e.g., [64]).

Potential and Limits of the Methodology
Nowadays several SAR satellites such as Sentinel-1A/B, ALOS-2, CSK, TSK, allow the production of coseismic interferogramas relatively quickly (from hours to a few days) after an earthquake.The new European Space Agency (ESA) Sentinel-1 constellation is particularly well adapted to this task, since it acquires images globally every 6-12 days, that are distributed free of charge.In the coming years, several initiatives now in the testing phase to produce systematic interferograms should be fully operational (e.g., GEP, LicS, ARIA, EPOS).
Coseismic GPS displacements are not globally available and their accessibility depends on different factors, including the previous existence of a permanent and non-permanent GPS networks in the zone.Although more and more countries have a continuous GPS network, the data are rarely public.However, some international initiatives exist to give access to GPS data, such as UNAVCO consortium (http://www.unavco.org/) that provides Open Access GPS data with daily RINEX files suitable for processing, although they mostly cover the United States and Caribbean zones.Similarly, the European Plate Observing System will soon provide access to GNSS data in real time over Europe (https://www.epos-ip.org/tcs/gnss-data-and-products).The availability of real-time High-Rate Global Positioning System (HRGPS) data can be used to produce rapid magnitude estimation and produce rapid earthquake source models (e.g., [64]).
One of the more time consuming steps when applying this methodology is the compilation of faults and volcanic elements to build the models.In the case of subduction earthquakes, the global model for subduction geometries of [45] can be downloaded from (https://earthquake.usgs.gov/data/slab/).Crustal faults are more difficult to obtain, and despite the important efforts towards the systematization and compilation of active sources (GEM, Global Active Faults: https://github.com/GEMScienceTools/gem-global-active-faults; SHARE, European Database of Seismogenic Faults: http://diss.rm.ingv.it/share-edsf/SHARE_WP3.2_Database.html)the completeness of the catalogue or the level of detail and the parameterization of the fault elements usually requires further analysis.The existence of national publicly available databases (e.g., the Quaternary Active Faults Database of Iberia, http://info.igme.es/qafi/, the New Zealand Active Faults Database http://data.gns.cri.nz/af/, the database of active faults in Greece [65]) is a great advantage in these cases.However, the main factor that determines the quality of the resulting maps is the quality, reliability, and completeness of the database (i.e., how reliable the faults and volcanoes parameters are).For instance, unmapped faults can result in unexpected triggered earthquakes, such as the case of the Christchurch earthquake [66].The reliability of the maps will increase as the knowledge and characterization of active faults and volcanic system improves.
The static Coulomb stress criterion, although a simplistic approach, is a powerful tool to explain many fault and volcanic interactions and is issued as a tool for earthquake forecasting (e.g., [7,19,67]).Different scientific software can be used to carry out these models (e.g., Coulomb 3.3 [48], RELAX [68]).
Apart from the fault-fault or fault-volcano interactions analysed in this study, other interactions can occur that could be taken into account, for example, aftershocks and postseismic deformation (e.g., [14]), volcanic interactions (e.g., [69]), or coseismic-induced landslides (e.g., [70]).Beyond the static stress changes estimated here, other mechanisms such dynamic stresses, pore-fluid diffusion, and viscoelastic rebound can also be important for different time frames after the earthquake.Additionally, the incorporation of relevant aftershocks into the models, as well as postseismic phenomena, can improve the hazard estimation and focus these tools towards an operational forecasting system.
Finally, a close collaboration with the final user is critical to create functional and useful products.The comprehension of our final maps was significantly improved thanks to the constructive criticisms and specific proposals of the Geological Survey of Ecuador staff (Instituto Nacional de Investigación Geológico Minero Metalúrgico, INIGEMM).

Conclusions
We present a methodology to evaluate earthquake-induced effects on neighbouring faults and volcanoes and create easily understandable maps addressed to decision-makers managing the post-earthquake situation.The first part of the method involves estimating the earthquake source model from space geodetic data, the second part consists in using this model as input to calculate static stress changes on tectonic and volcanic elements in the vicinity of the earthquake, and the third part is the creation of the final maps by synthetizing and simplifying the results for non-experts.
We apply the methodology to the Mw 7.8, 2016 Pedernales earthquake, Ecuador.Using Sentinel-1 InSAR and GPS data, we estimate the coseismic deformation field associated with the Pedernales earthquake and invert this field to model the source.Then, we calculate static stress changes on the subduction interface, on the active continental faults of Ecuador and on the volcanoes.The resulting maps show areas of high, medium, and low levels of potential activation of faults or volcanoes.When compared with seismic and volcanic event occurred after the Pedernales earthquake, we found in general a good agreement, with most of the largest aftershocks nucleated in areas of stress loading in our maps.Seismic and volcanic events located in areas of low level of potential activation in our maps are probably due to other phenomena not considered in our analysis, such as static stress changes from large aftershocks, dynamic stresses, or postseismic phenomena.
Our maps are useful to identify regions where new volcanic or seismic events are more likely to happen taking into account the static stress changes produced by the mainshock.But our models do not take into account the influence of other processes, such as dynamic stress changes produced by the earthquake, postseismic phenomena, or other tectonic and volcanic activity of the region.Hence, other seismic or volcanic events could occur out of the critical areas depicted in the proposed maps (red and yellow colours).
One of the main limitations of the methodology is the static stress approach, since other earthquake-related processes can also influence the behaviour of faults and volcanoes (dynamic stresses, pore-fluid diffusion, viscoelastic rebound).Also, an important limitation is the time necessary to have a new SAR acquisition over the earthquake area, which usually takes several days.This could be overcome using in the meantime the focal mechanism as a first-order earthquake source.
Considering the open data policy of Sentinel-1 SAR data, the availability of tools for the estimation of the earthquake source, and the evaluation of static stress changes, the methodology proposed here can be applied to earthquakes worldwide to quickly generate these maps after the earthquake and update them when additional data (e.g., GPS data, coseismic interferograms from other satellites, seismological and tsunami data) are available.

Supplementary Materials:
The following are available online at http://www.mdpi.com/s1, Figure S1: Ground deformation predicted by a uniform-slip model, Figure S2: Observed, modelled and residual GPS displacements, Figure S3: Observed, modelled and residual InSAR displacements, Figure S4: Smoothing coefficient for the optimal source model, Figure S5: Coulomb failure stress changes using the Global CMT focal mechanism, Figure S6: Map of potentially activated areas on subduction interface using as earthquake source the Global CMT focal mechanism, Figure S7: Map of potentially activated volcanoes on Ecuador using as earthquake source the Global CMT focal mechanism.Figure S8: Map with focal mechanism of seismic events with Mw ≥ 6 occurred in the period 16 April 2016-18 February 2018, Table S1: Differences between GPS coseismic displacements and GPS coseismic displacements plus + postseismic displacement, Table S2: Comparison between GPS and InSAR displacements, Table S3: strike/dip/rake values used for the Coulomb failure stress estimation on each family, Table S4: Dip and Azimuth of magma paths of the five selected Ecuadorian volcanoes used for the estimation of the stress changes.

Figure 1 .
Figure 1.Reference map of our study area in Ecuador (red box in the inset map).The relative Nazca-South American convergence rate and direction are shown with a black arrow[20] and the location of the trench is shown with a heavy black barbed line.The black dashed box shows the region in Figures2 and 4. The location of the continuous GPS sites used in this study is shown as blue inverted triangles.The Global Centroid Moment Tensor (GCMT) focal mechanism of the 2016 Pedernales earthquake is shown.The topography and bathymetry are from the ETOPO1 1 arc-minute global relief model[24].The green segment indicates the approximate rupture area of the 1906 earthquake[25] and the red semi-transparent areas represent the rupture areas of the 1942, 1958, and 1979 earthquakes[26][27][28][29].The approximate rupture area of the Mw 7.8 2016 Pedernales earthquake from our variable-slip model is also shown (semi-transparent blue area).

Figure 1 .
Figure 1.Reference map of our study area in Ecuador (red box in the inset map).The relative Nazca-South American convergence rate and direction are shown with a black arrow [20] and the location of the trench is shown with a heavy black barbed line.The black dashed box shows the region in Figures 2 and 4. The location of the continuous GPS sites used in this study is shown as blue inverted triangles.The Global Centroid Moment Tensor (GCMT) focal mechanism of the 2016 Pedernales earthquake is shown.The topography and bathymetry are from the ETOPO1 1 arc-minute global relief model[24].The green segment indicates the approximate rupture area of the 1906 earthquake[25] and the red semi-transparent areas represent the rupture areas of the 1942, 1958, and 1979 earthquakes[26][27][28][29].The approximate rupture area of the Mw 7.8 2016 Pedernales earthquake from our variable-slip model is also shown (semi-transparent blue area).

Figure 2 .
Figure 2. Coseismic displacement field from InSAR and GPS data.The coseismic interferogram (unwrapped) shows surface displacement in the line-of-sight direction.The white arrow indicates satellite Line-of-Sight direction (LOS) and satellite azimuth (Az).Black and blue arrows represent GNSS horizontal (black arrows) and vertical (blue arrows) displacements estimated from the Continuous Monitoring GNSS Network (REGME) in Ecuador, before and after the 16 April 2016, Mw7.8 Pedernales Earthquake.GNSS displacements and associated errors are shown in Table1.

Figure 2 .
Figure 2. Coseismic displacement field from InSAR and GPS data.The coseismic interferogram (unwrapped) shows surface displacement in the line-of-sight direction.The white arrow indicates satellite Line-of-Sight direction (LOS) and satellite azimuth (Az).Black and blue arrows represent GNSS horizontal (black arrows) and vertical (blue arrows) displacements estimated from the Continuous Monitoring GNSS Network (REGME) in Ecuador, before and after the 16 April 2016, Mw7.8 Pedernales Earthquake.GNSS displacements and associated errors are shown in Table1.

•
Map of potentially activated continental faults: We integrate the six different maps of Coulomb stress changes over the correspondent six fault subsets into one map of average positive stress change.We classify the final ∆CFS values obtained in three levels of potential activation and colour the final map according to those levels: green (low level) indicates areas with ∆CFS < 0.1 bar, yellow (medium level) indicates areas with 0.1 bar ≤ ∆CFS ≤ 1 bar, and red (high level) indicates areas with ∆CFS > 1 bar.• Map of potentially activated mapped faults: We classify the final ∆CFS values obtained over each fault trace in three levels of potential activation of the fault segment.Then fault traces are coloured according to those levels: green (low level) indicates faults with ∆CFS < 0.1 bar, yellow (medium level) indicates fault segments with 0.1 bar ≤ ∆CFS ≤ 1 bar, and red (high level) indicates faults with ∆CFS > 1 bar.• Map of potentially activated areas on the subduction interface: We classify the final ∆CFS values obtained over each node of the subduction interface in three levels of potential activation and colour the final map using those levels: green (low level) indicates areas of the slab with ∆CFS < 0.1 bar, yellow (medium level) indicates areas with 0.1 bar ≤ ∆CFS ≤ 1 bar, and red (high level) indicates areas with ∆CFS > 1 bar.

Figure 3 .
Figure 3. Position time series of the 17 cGPS sites used in this study.For each site, daily positions for components East, North and Up are shown in blue, red, and green, respectively.The vertical grey line represents the date of the Pedernales earthquake (16 April 2016).

Figure 3 .
Figure 3. Position time series of the 17 cGPS sites used in this study.For each site, daily positions for components East, North and Up are shown in blue, red, and green, respectively.The vertical grey line represents the date of the Pedernales earthquake (16 April 2016).

Figure 4 .
Figure 4. Coseismic slip distribution on the subduction interface from inversion of Sentinel-1 and GPS data.Colours show the magnitude of slip in metres.Depths on the fault plane are shown as green dashed lines.

Figure 5 .
Figure 5. Rose diagrams of the frequency of faults orientation according to the fault type, extracted from the Ecuador Active Faults database of [50].(a) Normal faults, two main sets with NE-and NW trending; (b) Reverse faults, showing a main orientation towards NNW; (c) Right Lateral strike-slip, with a main direction of NE-trending; (d) Left lateral strike slip faults with two conjugate sets, NW-SE and ENE-WSW trending faults.In (c) and (d), arrows indicate the sense of the movement for strike slips; (e) Fault segments (green lines) from the Ecuador Active Faults database of [50].Holocene volcanoes from http://volcano.si.edu/search_volcano.cfm are shown as red triangles.

Figure 4 .
Figure 4. Coseismic slip distribution on the subduction interface from inversion of Sentinel-1 and GPS data.Colours show the magnitude of slip in metres.Depths on the fault plane are shown as green dashed lines.

Figure 4 .
Figure 4. Coseismic slip distribution on the subduction interface from inversion of Sentinel-1 and GPS data.Colours show the magnitude of slip in metres.Depths on the fault plane are shown as green dashed lines.

Figure 5 .
Figure 5. Rose diagrams of the frequency of faults orientation according to the fault type, extracted from the Ecuador Active Faults database of [50].(a) Normal faults, two main sets with NE-and NW trending; (b) Reverse faults, showing a main orientation towards NNW; (c) Right Lateral strike-slip, with a main direction of NE-trending; (d) Left lateral strike slip faults with two conjugate sets, NW-SE and ENE-WSW trending faults.In (c) and (d), arrows indicate the sense of the movement for strike slips; (e) Fault segments (green lines) from the Ecuador Active Faults database of [50].Holocene volcanoes from http://volcano.si.edu/search_volcano.cfm are shown as red triangles.

Figure 5 .
Figure 5. Rose diagrams of the frequency of faults orientation according to the fault type, extracted from the Ecuador Active Faults database of [50].(a) Normal faults, two main sets with NE-and NW trending; (b) Reverse faults, showing a main orientation towards NNW; (c) Right Lateral strike-slip, with a main direction of NE-trending; (d) Left lateral strike slip faults with two conjugate sets, NW-SE and ENE-WSW trending faults.In (c) and (d), arrows indicate the sense of the movement for strike slips; (e) Fault segments (green lines) from the Ecuador Active Faults database of [50].Holocene volcanoes from http://volcano.si.edu/search_volcano.cfm are shown as red triangles.

Figure 7 .
Figure 7. (a) Map of potentially activated mapped faults of Ecuador.Coloured segments correspond to fault traces from [50].Colours green-yellow-red indicate the level of potential activation of the fault segment by the Pedernales earthquake.Faults traces labelled X1 and X2 are discussed in Section 5.1; (b) Map of potentially activated continental faults.

Figure 8 .
Figure 8. Map of potentially activated areas on subduction interface.The colour scale shows the level of potential activation of the subduction interface by the Pedernales earthquake.Green (low level) indicates areas with ∆ < 0.1 bar, yellow (medium level) indicates areas with 0.1 bar ≤ ∆ ≤ 1 bar and red (high level) indicates areas with ∆ > 1 bar.Contours of 25-km slab iso-depth from [45] are indicated by the black dashed lines.This map was calculated based on our variable-slip source model.See Supplementary Figure S6 for the map calculated using the earthquake source from the Global CMT focal mechanism.Black circles indicate the location of the Mw ≥ 6 events occurred after the Pedernales earthquake.

Figure 7 .
Figure 7. (a) Map of potentially activated mapped faults of Ecuador.Coloured segments correspond to fault traces from [50].Colours green-yellow-red indicate the level of potential activation of the fault segment by the Pedernales earthquake.Faults traces labelled X 1 and X 2 are discussed in Section 5.1; (b) Map of potentially activated continental faults.

Figure 7 .
Figure 7. (a) Map of potentially activated mapped faults of Ecuador.Coloured segments correspond to fault traces from [50].Colours green-yellow-red indicate the level of potential activation of the fault segment by the Pedernales earthquake.Faults traces labelled X1 and X2 are discussed in Section 5.1; (b) Map of potentially activated continental faults.

Figure 8 .
Figure 8. Map of potentially activated areas on subduction interface.The colour scale shows the level of potential activation of the subduction interface by the Pedernales earthquake.Green (low level) indicates areas with ∆ < 0.1 bar, yellow (medium level) indicates areas with 0.1 bar ≤ ∆ ≤ 1 bar and red (high level) indicates areas with ∆ > 1 bar.Contours of 25-km slab iso-depth from [45] are indicated by the black dashed lines.This map was calculated based on our variable-slip source model.See Supplementary Figure S6 for the map calculated using the earthquake source from the Global CMT focal mechanism.Black circles indicate the location of the Mw ≥ 6 events occurred after the Pedernales earthquake.

Figure 8 .
Figure 8. Map of potentially activated areas on subduction interface.The colour scale shows the level of potential activation of the subduction interface by the Pedernales earthquake.Green (low level) indicates areas with ∆CFS < 0.1 bar, yellow (medium level) indicates areas with 0.1 bar ≤ ∆CFS ≤ 1 bar and red (high level) indicates areas with ∆CFS > 1 bar.Contours of 25-km slab iso-depth from[45] are indicated by the black dashed lines.This map was calculated based on our variable-slip source model.See Supplementary FigureS6for the map calculated using the earthquake source from the Global CMT focal mechanism.Black circles indicate the location of the Mw ≥ 6 events occurred after the Pedernales earthquake.

Figure 9 .
Figure 9. Map of potentially activated Ecuadorian volcanoes.Circled numbers indicate the Holocene volcanoes from http://volcano.si.edu/search_volcano.cfm.Colours indicate the level of potential activation of the volcano by the Pedernales earthquake, which depends on the value of the induced normal stress change (∆ ) on the volcano magma pathway.Green (low level) indicates volcanoes with ∆ < 0.01 bar, yellow (medium level) indicates volcanoes with 0.01 bar ≤ ∆ ≤ 0.1 bar, and red (high level) indicates volcanoes with ∆ > 0.1 bar.This map was calculated based on our source model.See Supplementary Figure S7 for the map calculated using the earthquake source from the Global CMT focal mechanism.

Figure 9 .
Figure 9. Map of potentially activated Ecuadorian volcanoes.Circled numbers indicate the Holocene volcanoes from http://volcano.si.edu/search_volcano.cfm.Colours indicate the level of potential activation of the volcano by the Pedernales earthquake, which depends on the value of the induced normal stress change (∆σ) on the volcano magma pathway.Green (low level) indicates volcanoes with ∆σ < 0.01 bar, yellow (medium level) indicates volcanoes with 0.01 bar ≤ ∆σ ≤ 0.1 bar, and red (high level) indicates volcanoes with ∆σ > 0.1 bar.This map was calculated based on our source model.See Supplementary FigureS7for the map calculated using the earthquake source from the Global CMT focal mechanism.

Figure 10 .
Figure 10.(a) Map of potentially activated mapped faults calculated using the earthquake source from the Global CMT focal mechanism.The black rectangle delineates the source model.Faults labelled X1 and X2 are discussed in Section 5.1; (b) Difference of Coulomb stress change on mapped faults.This map results from the difference between map in Figure 7a and map in Figure 10a.

Figure 10 .
Figure 10.(a) Map of potentially activated mapped faults calculated using the earthquake source from the Global CMT focal mechanism.The black rectangle delineates the source model.Faults labelled X 1 and X 2 are discussed in Section 5.1; (b) Difference of Coulomb stress change on mapped faults.This map results from the difference between map in Figure 7a and map in Figure 10a.

Figure 11 .
Figure 11.(a) Coloured circles represent seismic events in the period 16/04/2016-18/02/2018 from Instituto Geofísico Escuela Politécnica Nacional filtered for continental faults (distance to slab >10 km and depth <15 km).Colours correspond to the traffic light colour code used in Figures 7 and 8: red circles represent seismic events located in areas with a high level of potential activation (∆CFS > 1 bar), yellow circles represent seismic events in medium level of potential activation (0.1 bar ≤ ∆CFS ≤ 1), and green circles represent seismic events in areas with a low level of potential activation (∆CFS < 0.1 bar; (b) Seismic events in the period 16/04/2016-18/02/2018 filtered for the subduction interface (distance to slab ≤10 km).Colours correspond to the traffic light colour code used in (a); (c,d) are the accumulated histograms of frequency (%) of seismic events after the mainshock on areas with the correspondent computed Coulomb Stress change for the continental faults (b) and the subduction interface (c).The colour of the bar indicates the level of potential activation on a traffic light colour code of Figure 8.

Figure 11 .
Figure 11.(a) Coloured circles represent seismic events in the period 16/04/2016-18/02/2018 from Instituto Geofísico Escuela Politécnica Nacional filtered for continental faults (distance to slab >10 km and depth <15 km).Colours correspond to the traffic light colour code used in Figures 7 and 8: red circles represent seismic events located in areas with a high level of potential activation (∆CFS > 1 bar), yellow circles represent seismic events in medium level of potential activation (0.1 bar ≤ ∆CFS ≤ 1), and green circles represent seismic events in areas with a low level of potential activation (∆CFS < 0.1 bar; (b) Seismic events in the period 16/04/2016-18/02/2018 filtered for the subduction interface (distance to slab ≤10 km).Colours correspond to the traffic light colour code used in (a); (c,d) are the accumulated histograms of frequency (%) of seismic events after the mainshock on areas with the correspondent computed Coulomb Stress change for the continental faults (b) and the subduction interface (c).The colour of the bar indicates the level of potential activation on a traffic light colour code of Figure 8.

Author
Contributions: M.B.-P.and J.A.A.G. designed the study.O.M. and M.B.-P.processed InSAR data.A.S. and M.P.L. processed and analysed GPS data.M.B.-P.modelled the earthquake source from InSAR and GPS data.J.A.A.G. performed the Coulomb failure calculation and created the final maps.R.P.-L., M.B.-P., and K.C. collected and analysed the catalogues of active faults and volcanoes.A.L., J.J.M.D., G.H., J.P.G., R.M.M. analysed the final maps and provided feedback to improve them.J.A.A.G. and M.B.-P.prepared the figures.M.B.-P.prepared the manuscript with advice and feedback from all team members.

Table 1 .
1. Coseismic displacements at the 17 cGPS sites used in this study.

Table 1 .
1. Coseismic displacements at the 17 cGPS sites used in this study.