Over the last 3 decades, GPS measurements have been instrumental in quantifying tectonic plates current motion and deformation. Complex patterns of deformation along the plate boundaries revealed heterogeneous coupling on the plates interface and imaged seismic segments at different stages of their seismic cycle. Along the South-American trench in Chile, where large earthquakes occur frequently, continuous GPS observations (cGPS) captured both the long-term plate motion and the transient deformations associated with the seismic cycle. Over the years, a network of hundreds of cGPS stations has been deployed all across the South-American continent by many different institutions for many different sorts of purposes ranging from geographic reference to tsunami early warning. We report here on the processing of 20 years (2000–2020) worth of data over a selection of cGPS stations, devoted to the quantification and analysis of the deformation along the Chilean subduction zone between 18°S and 40°S. We use all available data near the trench in Chile and a less dense network inside the continent where the gradient of deformation is lesser. Our database, named SOAM_GNSS_solENS, provides time series of precise daily station position, obtained from double difference (DD) processing and expressed in the International Terrestrial Reference Frame 2014 (ITRF14). These time series allow to quantify, with sub-millimetric precision, any kind of ongoing deformation process, either from tectonic origin such as interseismic deformation, co- and post-seismic displacements associated with earthquakes, transient deformation associated to seismic swarms and/or a-seismic slow-slip events, or of other origin such as hydrological loading (for ex., the Amazonian basin load) or any other type of loading affecting the surface of the Earth (tides, atmosphere, etc.). We also provide a database of coseismic displacements associated with close to 60 earthquakes of Mw larger than 6.5 that occurred over the last 20 years within the observation area. All time series are directly accessible through a deposit and we plan to make them available through a web interface that will allow any user to perform elementary operations like estimating offsets, detecting outliers, detrending, filtering and stacking. That database will evolve with time, aggregating more data. In the future, we also plan to complement that database with a rapid solution in quasi real time processed in Precise Point Positioning (PPP), and with hourly atmospheric delays associated with water vapor content of the lower layer of the atmosphere.

Au cours des trois dernières décennies, les mesures GPS ont permis de quantifier le mouvement et la déformation actuelle des plaques tectoniques. Des modèles complexes de déformation le long des frontières des plaques ont révélé un couplage hétérogène à l’interface des plaques, et ont permis d’imager des segments sismiques à différents stades de leur cycle sismique. Le long de la fosse sud-américaine au Chili, où de grands tremblements de terre se produisent fréquemment, des observations GPS continues (cGPS) ont capturé à la fois le mouvement des plaques à long terme et les déformations transitoires associées au cycle sismique. Au fil des ans, un réseau de plusieurs centaines de stations cGPS a été déployé sur tout le continent sud-américain, par de nombreuses institutions différentes, pour toutes sortes d’applications allant de la référence géographique à l’alerte précoce aux tsunamis. Nous présentons ici le traitement de 20 années (2000–2020) de données sur une sélection de stations cGPS, consacrées à la quantification et à l’analyse de la déformation le long de la zone de subduction chilienne entre 18°S et 40°S. Nous utilisons toutes les données disponibles au Chili (proche de la fosse de subduction) et un réseau moins dense à l’intérieur du continent où le gradient de déformation est moindre. Notre base de données, nommée SOAM_GNSS_solENS, fournit des séries temporelles de position quotidienne précise des stations, obtenues par traitement en double différence (DD) et exprimées dans le système international de référence terrestre 2014 (ITRF14). Ces séries temporelles permettent de quantifier, avec une précision sub-millimétrique, tout type de processus de déformation en cours, qu’il soit d’origine tectonique comme la déformation intersismique, les déplacements co- et post-sismiques associés aux tremblements de terre, la déformation transitoire associée aux essaims sismiques et/ou aux événements de glissement lent a-sismiques, ou d’origine différentes comme les charges hydrologiques (par exemple, la charge du bassin amazonien) ou tout autre type de charge affectant la surface de la Terre (marées, atmosphère, etc.). Nous fournissons également une base de données des déplacements cosismiques associés à près de 60 séismes de Mw supérieure à 6.5, survenus au cours des 20 dernières années dans la zone d’observation et détectés par GPS. Toutes les séries temporelles sont directement accessibles via un dépôt et une interface web actuellement en développement permettra à tout utilisateur d’effectuer des opérations élémentaires telles que l’estimation des sauts, la détection des valeurs aberrantes, la correction de tendance linéaire, le filtrage et la superposition. Cette base de données évoluera avec le temps, en agrégeant davantage de données. Dans le futur, nous prévoyons également de la compléter avec une solution rapide en temps quasi réel traitée en Precise Point Positioning (PPP), et avec des séries temporelles de retards atmosphériques horaires associés au contenu en vapeur d’eau de la couche inférieure de l’atmosphère.

The Chilean subduction zone is extremely seismically active with on average one Mw 8 earthquake every ten years and at least one mega earthquake of magnitude 9 and over per century (Beck et al., 1998, Lomnitz, 2004, Ruiz and Madariaga, 2018). The largest earthquake ever recorded by seismometers, the 22nd May 1960 Valdivia earthquake of magnitude 9.5, occurred along the southern part of the Chilean subduction zone and generated a giant trans-Pacific tsunami that caused catastrophic damage along the coasts of Hawaii and Japan (Cifuentes, 1989).

Since the early 90's, GPS was instrumental in quantifying details of the tectonic plates motion and deformation in South-America (e.g., Larson et al., 1997, Angermann et al., 1999, Norabuena et al., 1999, Klotz et al., 2001, Kendrick et al., 2001, Ruegg et al., 2002, Kendrick et al., 2003, Brooks et al., 2003, Altamimi et al., 2007). As in many subduction zones worldwide, GPS measurements in Chile showed that the velocity field derived from several years of observation are spatially heterogeneous, drawing large patches where velocity gradient or deformation is intense separated by narrow zones where it is low. Those variations of deformation are evidence of variable coupling along the plate interface, with patches where deformation is high corresponding to long segments of high coupling (e.g., Moreno et al., 2011, Métois et al., 2012, 2016). These segment accumulate deficit of slip and store elastic energy in the surrounding medium, making them prone to a seismic rupture if enough time has elapsed, leading to the concept of seismic gap (e.g., McCann et al., 1979).

In Chile, several of these large seismic gaps ruptured during the last ten years. The Maule earthquake of 27th February 2010 (Mw 8.8), the Iquique earthquake of 1st April 2014 (Mw 8.2) and the Illapel earthquake of 16th September 2015 (Mw 8.4), all occurred in one of these previously identified seismic gaps, where GPS had revealed that deformation was accumulating prior to the earthquake (e.g., Moreno et al., 2010, Vigny et al., 2011, Ruiz et al., 2014, Schurr et al., 2014, Tilmann et al., 2016, Klein et al., 2017). Other identified seismic gaps pause an acute seismic hazard today. The “Great North gap”, locus of the megathrust earthquake of 1877 (Comte and Pardo, 1991), the “Atacama gap”, locus of the megathrust earthquakes of 1819 and 1922 (Willis, 1929), and the “Valparaiso gap” locus of the 1985 earthquake (Comte et al., 1986), are regions most likely to be approaching the end of their seismic cycle.

GPS, sometimes associated with InSAR, also revealed transient deformation at different time and space scales, associated with seismic swarms and/or a-seismic deformation (Pritchard and Simons, 2006, Klein et al., 2018, Klein et al., 2021) and the post-seismic deformation that immediately follows an earthquake which can reach several millimeters per year, as far as thousands of kilometers away from the epicenter and last for decades or even longer (Khazaradze et al., 2002, Suito and Freymueller, 2009, Trubienko et al., 2013, Klein et al., 2016, Hu et al., 2016). The latter results from the relaxation of stress in the underlying viscous mantle and takes place at the continent scale. Detecting and quantifying such small deformation is a challenge, but sheds light on fundamental properties of the Earth mantle like the viscosity of its different layers (Trubienko et al., 2014). In turn, a more precise quantification of these numbers allows to improve the prediction of long term deformation induced by megathrust earthquakes, glacial-melt isostatic rebound, and relative ground sea-level variations at coastal areas.

We report here on the processing of a large cGPS network deployed by many different institutions across the South American continent over the years. This network has the primary goal of measuring crustal deformation associated with the physical process of plate tectonics and earthquake cycle. Combining data in Chile with publicly available data from other national and international cGPS networks in South America, we provide a precise and consistent solution spanning 20 years at the continental scale. This solution is provided in the form of a database of time series at ~ 220 selected GPS stations, starting as early as 2000 when possible and currently until the end of 2020. These time series depict precise daily coordinates in the global reference frame ITRF2014 (Altamimi et al., 2017). Many different signals are present in those time series, such as instantaneous jumps attributed to either instrumental changes or earthquakes, and transient deformation attributed to different processes of either tectonic origin (e.g., slow-slip earthquakes) or non-tectonic origin (e.g., seasonal variations of hydrological loads). The precision of the processing is such that it allows the detection and quantification of very small co-seismic offsets (as small as 1 mm) at large distances from moderate size earthquakes (~ Mw 6.5) that frequently occur in Chile. Thus, we provide tables of co-seismic displacements for 55 identified earthquakes. That database will evolve with time, aggregating more data. In the future, we also plan to complement that database with a rapid solution in quasi real time processed in PPP, and with hourly atmospheric delays associated to water vapor content of the lower layer of the atmosphere.

Data description

As of today, we only process observations from the GPS constellation, but futures updates of this solution will include observations from other GNSS constellations, and in particular from the European Galileo. We process most of the GNSS data from the national Chilean network (Báez et al., 2018), between 18°S and 40°S, i.e., in the area of Chile where the subduction is rather linear and “simple”. North of 18°S, in the Arica bend, the trench curvature induces a rapid change of obliquity and increasing distance of the coastline from the trench. South of 40°S, plate tectonics becomes more complex with partitioning of the oblique convergence along the Liquiñe-Ofqui strike-slip fault (Cembrano et al., 1996). Additionally, towards the South, the continent becomes narrower, making more difficult the definition of a stable plate interior reference frame to study the tectonic deformation. Also, logistics become increasingly difficult at higher latitudes, which renders the cGPS network less dense and more difficult to maintain over a long period of time. Throughout Chile, the national data set is also complemented by local dense networks aiming at studying specific questions in targeted areas (e.g.Klein et al., 2021). In addition and in order to realise a continental scale processing, we include data of a selection of stations from the Argentine RAMSAC network (Piñón et al., 2018) and of the Brazil RBMC network. This selection is made in order to reasonably spatially sample the whole continent. We also include all the IGS stations available on the South American continent. This whole data-set is composed of up to 220 stations providing data between 01/01/2000 and currently 31/12/2020.

Processing strategy

We processed this data-set using the GAMIT/GLOBK software following the classical MIT methodology (Herring et al., 2010a, Herring et al., 2010b). GAMIT uses a so-called double-difference of GPS satellite phase measurements among all pair of sites. Because of the tremendous networks expansion through time, the double-difference approach cannot be applied to the whole data-set in a single network over the whole period. Before 2008, all stations are processed simultaneously in a single network (net CHIL). From 2008 onwards, we divided the processing into 3 sub-networks, based on geographical criteria (Fig. 1). First, we design a continental-scale sub-network, with a total of 75 stations spanning the whole South-America (subnet SAM1). The second group, includes about 70 stations spanning Central Chile and Central Argentina (subnet CHL1). Finally, the third sub-network includes slightly more than 60 stations spanning the northern regions of Chile, Argentina and Brazil (subnet CHL2). Up to 30 stations are common to the 3 sub-networks to ensure a reliable and stable combination of the sub-networks over the years. The numbers of stations in the sub-networks and in total varies quite significantly over time, due to stations showing up late or being decommissioned at some stage. The actual number of stations per day in each subnetworks is represented in Figure 2.

In a first step, we reduce the 30 s sampling data of 24-hour sessions to daily estimates of station positions using the GAMIT software (Herring et al., 2010a), using the ionosphere-free combination, and fixing the phase ambiguities to integer values. Overall, we follow the conventions recommended by IGS in the frame of the third reanalysis of the full history of GPS data collected by the IGS global network since 1994. We use precise orbits from the International GNSS Service for Geodynamics (IGS) (Dow et al., 2009) and IGS tables to describe the phase centers variations of the antennas. We use the Vienna Mapping function (VMF1) derived from numerical weather models for the mapping of the hydrostatic and wet tropospheric delays (Boehm et al., 2006). We estimate one tropospheric vertical delay parameter per station every 2h and 2 spatial gradient parameters on E-W and N-S directions per day. The complete parameters file used in GAMIT is available on the database repository. We also provide our station.info file listing all the known equipment changes used in our processing. We update this file over the years, checking and cross-checking metadata made available to us from different local sources (station descriptions, logsheets, rinex headers, field reports, stations visits, …).

In a second step, daily h-files that include adjusted parameters with their associated variance/covariance matrix for each sub-network are combined into one single continental h-file per day, using GLOBK. Such a solution is called “loosely” constrained, because only relative coordinates are precisely determined. We then use the PYACS toolbox to map the loosely constrained solutions onto the ITRF2014 (Altamimi et al., 2017). This is done by estimating the 7 parameters (3 translations, 3 rotation and a scale factor) of the Helmert transformation between our solution and an up-to-date cumulative IGS reference solution (here IGS21P01.ssc). PYACS uses the IGS discontinuity file regularly updated by the IGS analysis combination center (Rebischung, 2021), together with the ITRF/IGS parametric model for IGS sites experiencing non-linear post-seismic deformation. This procedure can be implemented either with PYACS or GLOBK. In theory, both should allow to reach similar results and GAMIT/GLOBK is undoubtedly the most advanced package. However, we choose to use PYACS, rather than GLOBK for several practical reasons. First, Helmert parameter estimation is performed by minimizing a L1 norm on coordinates in the topocentric local frame. This approach allows to separate the horizontal from the vertical component and is particularly efficient in identifying the frequent outliers in the vertical component. PYACS also provides a fast mode by neglecting the full variance/covariance matrices for rapid solutions. This mode eventually provides equivalent coordinates precision. This approach, in addition to the choice of using continental stations rather than the full global combination, results in considerable time saving, considering our current computing facilities.

We provide, in the Supplementary Material of this article and on the database repository, a complete list of processed stations, along with their geographic and Cartesian coordinates and their period of record.

Quality of the solution

Starting with the stabilization part performed using the PYACS software, we note that because of the network density evolution with time, the amount of reference stations also evolves with time, increasing from 10 to 25 over the first decade and eventually stabilizing to 30 during the second decade (Fig. 3). Figure 4 shows the stabilisation network. It is to note that, using such a network aperture, part of the real signal associated to hydrological continental loading (the largest being the Amazonian Basin), can be absorbed in the scale parameter of the Helmert transformation. But only part of it (and hopefully small), since this large scale deformation is affected by gradients that allow to partly separate its contribution from reference frame rotation, translation and scale factor (Fu et al., 2013, Chanard et al., 2018). In any case, such a stabilisation is optimal for the analysis of tectonic signal of the plate boundary, which is our primary concern, as well as the analysis of local sources of any kind (tectonic, hydrological, etc. See examples of deformation processes in section 4) but is probably less appropriate for the analysis of continental scale loading. We show in Figure 5 residuals time series after stabilisation at 3 selected IGS sites, highlighting the epoch used (in red) and the rejected one (in gray) during the Helmert estimation. Stations that show a transient deformation due to the post-seismic processes can still be used as stabilisation stations, because these processes are accounted for in the prediction of the station position in the ITRF2014. For example, at Antuco, at the foothill of the Andes (station ANTC), clearly affected by an intense post-seismic transient deformation, the residuals between the predicted station position and the measured station position are of the order of several millimeters only (Fig. 5A). On the other hand, the case of station URUS in Bolivia highlights the handling of discrepancies between a priori and actual velocities by Pyacs (Fig. 5C). That station was only integrated in the ITRF at the end of 2013 and anterior a priori coordinates are extrapolated. The definition of the ITRF is therefore not to blame. The important aspect that we want to underline here, is that PYACS automatically and efficiently rejects the station from the stabilisation process before 2014, because there is some inconsistency between ITRF2014 a priori values and actual observations, and re-integrates it after 2014 when ITRF2014 a priori values and actual observations are in agreement. The tipping date of 2014 corresponds to the occurrence of the 2014 Iquique earthquake, for which new coordinates (position and velocity) were derived after the earthquake. Indeed, although the station has a relatively small co-seismic offset, it exhibits a large change of velocity between before and after the earthquake, especially on the North component (cf.section 4.1.1).

We show in Figure 3 the evolution of the number of sites used to compute the Helmert parameters together with the time series of median and the weighted-root-mean-square (WRMS) of residuals (final coordinates minus IGS). Over the 20 years processed, we obtain a quite constant L2-norm postfit WRMS, averaging to less than 1.5 mm in both planimetric components, although variations appears quite larger over the first 3–5 years (Figs. 3B and 3C). Unsurprisingly, the average of L2-norm postfit WRMS on the vertical component is higher, reaching 5 mm, there again with significantly larger variations over the first 3–5 years. Note that we find in the postfit WRMS the seasonal variations that are not taken into account in the current parameterization of the ITRF2014. They become even more visible as the precision gets better with time. We expect to obtain even much flatter postfit WRMS curves when these seasonal oscillations will be accounted for in the ITRF, expected with the next release ITRF2020. However, the mean WRMS are already very small, less than 1.5 mm on both horizontal components for both the L1-norm and the L2-norm.

The daily scatter, defined here as the median of the differentiated time series, is a good indicator of the high frequency noise level of a time series. It is not affected by transient signals such as postseismic deformation. We can therefore compute it on all time series and over the complete period of record. Daily scatter as function of latitude for the 3 components, highlight a geographical pattern (Fig. 6). It appears significantly smaller for stations located between 18°S and 50°S, with sub-millimetric values for the horizontal components and less than 3 mm on the vertical. It is larger for stations located north of 18°S plus south of 40°S (with median slightly larger than 1 mm on the horizontal components and almost 5 mm on the vertical). This seems directly correlated with the density of stations. A lower density implies longer and scarcer baselines, inherently less precise because phase-ambiguity resolution becomes more difficult as baseline length increases. Recent study also showed that non-tidal atmospheric and oceanic loading deformations affect the stochastic properties of vertical motions in a latitude-depend way, which could explain part of the larger daily scatter observed at the extremes latitude of our network (Gobron et al., 2021). In any case, this analysis of daily scatter suggests that we can actually extract sub-millimetric short-term signals from horizontal time-series, like for instance tiny co-seismic offsets.

Comparison with existing solutions

We compare our solution with two existing ones: i) the third IGS reprocessing campaign forming the IGS contribution to forthcoming ITRF2020, Repro3 (Rebischung et al., 2019), ii) the Nevada Geodetic Laboratory solution (NGL, Blewitt et al., 2018). Both Repro3 and NGL solutions are global ones. But these solutions differ on several aspects. Repro3 is a combination of solutions from 10 different Analysis Center (ACs, which are COD, ESA, GFZ, GRG, JPL, MIT, NGS, TUG, ULR, WHU). This combination therefore includes solutions that differ regarding the included constellations (some ACs process only GPS, other include GLONASS and/or Galileo). The NGL solution uses a PPP approach including more than 15 000 stations globally available.

Here, we compare time series of specific interest from a selection of stations across the whole continent: SANT (near Santiago de Chile) is an IGS station located near the central part of the Chilean subduction zone (Fig. 7-I.A), CHPI (between Sao Paolo and Rio de Janeiro) is an RBMC-IGS station located on the so-called Brazilian craton (Fig. 7-II.A), and RIO2 (in Patagonia) is a RAMSAC-IGS station located at the extreme south of the continent (Fig. 7-III.A). All 3 have long lasting and complete times series and contribute to our stabilization network. The 3 time series are simply detrended but otherwise unfiltered. SANT is showed only between 2000 and 2010, period over which its motion is linear with time and unaffected by any large earthquake. CHPI and RIO2 are showed over the time period over which the 3 solutions provide a coordinate. Our solution is overall smoother, with a lower level of high frequency noise. It also depicts smaller WRMS on the 3 components and in particular on the vertical component. This WRMS reduction comes from low frequency noise reduction, rather than from the reduction of the high frequency noise, since the latter is only a small contribution to the overall series variance (see Fig. 7B). It has been showed that regional GPS solutions show lower noise than global solutions, which is the case of the NGL and IGS solutions (Williams et al., 2004). Therefore, part of the noise reduction described by our solution can result from continental stabilisation, which aperture remains significantly larger than expected from a regional solution. At SANT, low frequency “oscillations” are observed (Fig. 7-I.A). However, their amplitude are significantly lesser in our solution. Therefore, they are most likely to be coming from the stabilization process. Most of the reference frame oscillation, which appears as continental scale common mode noise, is absorbed by the continental stabilisation made with PYACS. RIO2 is an exception on which we observe strong oscillations on the north component of our solution (Fig. 7-III.A). We also observe, on all 3 solutions, offsets in 2012 and a following change of trend until 2020 when equipment was changed again. The 3 solutions disagree on the amplitude of the documented instrumental changes, which comes from how it is handled in the processing (for ex. the weight applied on observations or the quality of the wet troposphere corrections).

For these 3 stations, we compute the Lomb-Scargle Power Spectral Density (Fig. 7B) of the 3 detrended components. All spectrum show similar pattern with white noise dominating at high frequency corresponding to periods of a few days to weeks, and power increasing with lower frequency attesting for colored noise (Williams et al., 2004).

In the details, the north component of SANT time series shows an interesting behaviour (Fig. 7-IA): Recurrent small incremental southward motion (i.e., in 2003.5, 2004.5, 2005.5). These signals are visible on the 3 solutions but with different characteristics. On the NGL solution, which has a significantly larger daily scatter, they show as ramps, which can be approximated by a sinusoidal function and may be partly absorbed by a seasonal filter. In our own solution, with lesser daily scatter, they show as offsets. On the Repro3 solution, with intermediate daily scatter, they show as a mixture of ramps and offsets. Finally, on our solution with the smaller daily scatter, some structures start to emerge when they are drowned in the higher noise of the other solutions. This is discussed further in section 4.3.

One of the advantage of a post-processed solution is to ensure, through frequent databases synchronisation and reprocessing, that time series are complete over time. In some regions, this is a critical aspect. For example, most stations located in the region of the Illapel earthquake of 2015 in Chile, show a 2–3 years data gap in the NGL process between ~ 2010 and 2013 (cf. example of station CNBA, Fig. S1, in Supplementary Material). This gap has been pointed out in recent articles (e.g., Yuzariyadi and Heki, 2021), but is only a processing gap, not a data gap since these stations were operational and their data collected, but possibly not available to NGL at the time of their processing. Such large gaps can challenge the interpretation of observed signals. For example, a vertical offset shows up during this time window at all these stations in the NGL solution. These offsets can be interpreted as coseismic offsets, or anything else happening to the stations during this time frame, but are actually not real since they do not appear on any of our time series.

In the following, we qualitatively describe several examples of the signals that are present in our time series and that shed light on the different types of studies that can be performed using them. In addition, we also point out a few signals that still need further investigation to be fully understood.

Seismo-tectonic processes

Surface deformation associated with the seismic cycle along the Chilean subduction zone

Because of the occurrence of 3 major earthquakes along the subduction zone (Mw 8.8 in 2010, Mw 8.1 in 2014 and Mw 8.3 in 2015), several years after the beginning of the deployment of cGPS networks in Chile (end of the 90's), the 3 phases of the seismic cycle (pre-, co- and post-) could be observed in great details. To illustrate what is actually observed, we show stations affected by the 2010 Maule earthquake (ANTC & MAUL, Fig. 8A), the 2014 Iquique earthquake (URUS, Fig. 8B) and the 2015 Illapel earthquake (PFRJ, Fig. 8C).

Convergence of the Nazca plate towards the South American plate results in deformation along the Chilean margin induced by locking at portions of the subduction interface. This slow and continuous deformation can last centuries. This is the interseismic phase, of which we observed only the last decade or so, and which appears as a linear trend before 2010 at MAUL and ANTC, and before 2014 at URUS. When enough deformation has accumulated, an earthquake occurs, which generates a deformation of the surface in only a few seconds to minutes. In daily time series, these coseismic offsets show as instantaneous jumps between the positions of the days before and after the earthquake. Depending on whether the earthquake occur early or late in the day may shift the jump from the difference between the day before and the day of the earthquake to the day of the earthquake and the day after. Following the largest earthquakes (Mw larger than 7), a significant transient deformation, called post-seismic deformation, occurs and can be observed. Often approximated by a logarithmic or exponential functions, this deformation can last for decades but decreases with time (see for ex. the 1960 Mw 9.5 Valdivia earthquake, Barrientos et al., [1992]). Because of the geometry of the Chilean subduction striking North-South and the almost East-West motion of the Nazca plate relative to South America, most of the surface deformation is usually observed on the East component of GNSS time series. This is what we observe at ANTC and MAUL in February 2010, where the westward coseismic offset reached close to 1 m against only about 20 cm on the North-South component. The post-seismic deformation observed at these stations also follow that rule (see a detailled analysis of the postseismic deformation following the Maule earthquake in Klein et al. [2016]).

Station URUS, located at about 200 km in the North-East of the epicenter of the 2014 Iquique earthquake is particularly interesting because it does not follow that rule (Fig. 8B). The coseismic offset on the East component is about twice as large as the one on the North component (1.6 cm westward vs. a 0.8 cm southward, see table 2014-04-01.dat from the coseismic data base hereafter described). However, four years after the earthquake, the trend on the East-West component is almost back to the pre-earthquake trend. This is not the case on the North-South component, where the postseismic trend still remains very different and has not returned yet to the pre-seismic trend. This is due to the different geometry: the location of the station with respect to the earthquake, as far north as east from the epicenter.

Station PFRJ, located very close to the 2015 Illapel earthquake epicenter shows more complexity (Fig. 8C). As expected, the co-seismic displacement is dominated by the East-West component. However, the lesser displacement of the North-South component reveals that if the co-seismic jump is southward, the post-seismic deformation that follows is northward, and affected by short term transients. A similar observation can be made on the vertical component: the co-seismic displacement is downward, but the displacement that follows is upward and rapidly changing with time. These post-seismic displacements reveal a combination of jumps due to aftershocks and rapid sliding due to after-slip on the subduction interface. Both phenomenons occur on the subduction interface, which is close to the coastline where the station is located.

Catalog of coseismic deformation due to the major Chilean earthquakes

Having a coherent reanalysis of 20 years time series, we built a data base of coseismic offsets estimated for all events of Mw larger than 6.5 that occurred over the period along the Chilean subduction zone. We found 55 earthquakes on the USGS catalog, that are represented in Figure 9 with their focal mechanisms extracted from Global CMT (Dziewonski et al., 1981, Ekstrom et al., 2012). About a third of the events belongs to aftershock sequence of one of the 3 megathrust earthquakes of magnitude above 8: 2010 (Maule), 2014 (Iquique) and 2015 (Illapel). For events that are occurring on the same day (for ex. a Mw 6.8 and a Mw 6.5 occurred on 1 September 2020 in the Atacama region) or on two consecutive days (for ex. a Mw 7.7 occurred on 14 November 2007, the Tocopilla earthquake, its largest aftershock of Mw 6.8 occurred on the following day, Fig. 9-inset), we estimated cumulative coseismic offsets induced by both events. Eventually, this leaves a list of only 44 events.

In order to provide an homogeneous data base, we applied an automatic approach identical for all events. We use the raw time series without correcting for outliers or trends nor applying any filtering. A detailed description of our algorithm can be found in the Supplementary Material. Eventually, 36 events trigger a non-Null offset at at least one station (cf.Table in the Supplementary Material). For each of these 36 events, we provide an ASCII table (and a regional map of the coseismic offsets (for ex. Fig. 9-inset). All tables include a comprehensive header detailing its content (cf.Supplementary Material). All figures and tables are available on the ftp deposit.

Considering the methodology applied, care should be taken when considering one particular event. Our goal here is to provide a database that is homogeneous and includes small earthquakes at the limit of detectability of the network. The fact that coseismic displacements corresponding to these small earthquakes are detectable in the time series of several stations is often ignored. However, the quantification of a very small offset is often difficult and some tuning of the algorithm would be required depending on the specific conditions of the affected stations (data gaps, noise level, presence of simultaneous transients, etc.). Therefore, when interested in a particular event, it is best to re-evaluate the offsets, using a more thorough and adaptive methodology on filtered time series. For that same reason, the coseismic tables extracted here for the 3 megathrust earthquakes can hardly be compared with previously published tables (such as, Vigny et al., 2011, Ruiz et al., 2014, Klein et al., 2017) for which additional filtering and corrections, including removal of after-slip, were applied.

Seasonal deformation

GPS time series are notoriously affected by seasonal oscillations. Typically, these seasonal oscillations are often described by annual and semi-annual terms and thus removed by fitting a pair of sinusoidal terms to each time series (Nikolaidis, 2002). Although powerful, this method has some drawbacks. First, seasonal signals do not necessarily have constant amplitude nor phase cycles, they can be shifted in time (either because the meteorology itself is delayed, or because of delays in the Earth response to a changing seasonal load). Thus, the fit by a pair of sinusoidal terms of constant phase and amplitude is not perfect and can leave a significant misfit. Worse, it can introduce a bias (through aliasing) on the long term trend determination.

Despite the fact that most of Northern-Chile is a desert, there is a seasonal cycle visible in GPS data in this area. This cycle certainly originates from the hydrological load of the Amazonian basin and possibly from local atmospheric loading. The Amazonian basin hydrological load is huge: a vertical cycle of up to 4 cm and an horizontal cycle of up to 1 cm peak-to-peak is clearly observed at PortoVelho (POVE, Fig. 10a) which is located in a meander of river Rio Madeira, the biggest tributary of the Amazon, no more than 1500 km North-East of Iquique in North Chile. It is a very long wavelength signal, so stations in a given area of Chile will have almost the same cycle and most of this signal will cancel out when differenciating. However, small but significant differences of several millimetres (both in horizontal and vertical components) can be observed at stations several hundreds of km apart (e.g., Fu et al., 2013, Chanard et al., 2018). Finally, the periodogram of POVE (Fig. 10b) reveals a richer frequency content than the simple annual/semi-annual standard estimation. We do not detect any other clear period, such as the draconitic period (351.6 days), which cannot be separated from the annual solar period with less than 26 years of data. However, it reveals that the 6-month peak amplitude is in reality much smaller (two order of magnitude, on the 3 components) than the 12-month peak. This is very different from the classical simple estimation of only an annual and semi-annual terms, which usually finds the two terms to be of comparable amplitude. This finding suggests that the amplitude of this semi-annual term can be significantly but artificially inflated through aliasing with nearby frequencies when estimated alone. Some periodograms indeed show several pics around but not exactly at 6 months (Fig. 7B).

Transient deformation

GPS time series also reveal a lot of transient signals, lasting between several days to several years, potentially repeating with more or less regular periods. Transient deformation can be generated by a whole range of different sources such as volcanic deformation, oceanic or atmospheric loading effects, and again hydrological loads, depending on the ground composition and/or seasonal and multiyear groundwater variations (e.g., Silverii et al., 2016).

However, many transients observed in subduction zones are associated with faulting processes, such as postseismic deformation (following major earthquakes, explained in section 4.1.1) or Slow Slip Event. These can have variable duration and potentially be recurrent, as in the region of Copiapo (Klein et al., 2018). In some cases, they can precede earthquakes and contribute to the initiation of the main shock rupture, as for the 2014 Iquique earthquake (Ruiz et al., 2014, Socquet et al., 2017, Boudin et al., 2021) or the 2017 Valparaiso earthquakes (Ruiz et al., 2017).

Some observed transient signal remain hard to relate to any known sources. As an example, we show the time series of SANT and DGF1, both in the Santiago de Chile metropolitan area, separated by about ~ 35 km (Fig. 11). SANT time series exhibits regular events on the North component (previously mentioned section 3.2), characterized by a transient northward sliding, interrupted by a sudden southward jump, followed by the resumption of the northward sliding (for ex. mid 2003, mid 2004 and mid 2005). The station DGF1 was only installed in 2004, but rapidly exhibits a northward transient simultaneously with SANT. Santiago is located in a sedimentary basin, which is plausibly affected by hydrological processes which could have variable impact depending on the location and the monumentation of the stations (outcrop or building rooftop, depth of anchoring, etc.).

Such transient signals of currently unknown origin are difficult to deal with, because they are small and show differently on different solutions: jumps and offsets in our solution, continuous ramps of changing slope in NGL solution (Fig. 7-I). Because it has a low daily scatter, our solution appears appropriate to extract the details of these kind of events, but extra-care should be taken to avoid over-interpret signal that could be related to the monumentation of the stations or, in some extent, to the data processing. Typically, a careful consideration of offsets and different noise types is necessary for a correct modeling of time series, which is most of the time very challenging.

All products (time series, coseismic tables) can be downloaded on a ftp deposit (ftp://ftp.geologie.ens.fr/incoming/tempo/SOAM_GNSS_solENS/). They are all named SOAM_GNSS_solENS, followed by the type of product (currently _Pos or _Coseismic) and the release version. The present solution is the release v1.0 and will evolve with time: a reprocessing due to parameter changes will be reported by an increase of the first index (for ex. v2.0), the addition of new data will be reported by an increase of the second index (v1.1). For coseismic offset tables (and all future products), the release version refers to the GNSS solution used. Station’s time series come as pos file, a standard format, explained in detail in the file header. We knowingly choose to disseminate raw time series, ie unfiltered and unmodeled, in order to allow users of any community (seismotectonics, atmosphere, hydrology, …), to use them without any misconception about the data content or filtering.

Our data base is meant to evolve with time. First, and as previously explained, the current solution will evolve, including more recent data (a solution up to the end of 2021 will be made available early 2022), or integrating new processing parameters. The coseismic offset data base shall be updated simultaneously. Second, we also plan to propose more products, first with a solution processed in quasi real time following a PPP strategy (Zumberge et al., 1997). This solution will be significantly less complete since some data are arriving with up to several months of delay. But it would allow the release of a very quick partial solution, for example after future major earthquakes. Similar quick solutions could be provided for high rate GNSS as well, at least at all stations where 1 Hz data are acquired.

We also plan to provide hourly atmospheric zenithal delays at every station. Tropospheric zenithal delay time series can be used as a proxy of the water vapor contained in the troposphere above the stations and be used in meteorology.

Finally, the diffusion channel of our products will evolve. We are currently developing a web application, composed of i) a web mapping showing the processed network and allowing to access to stations metadata, ii) a web application that will allow to visualize, manipulate and analyse position time series, based on the Pyacs suite. ASCII files of time series will therefore be downloadable through both pages, raw or analysed.

For more than 20 years, our research group has been using GNSS data to work on the Chilean subduction. To do so, a long and thorough work of processing of the data produced by more than 200 stations located across the whole South-American continent has been led over time. Today, and after a complete and thorough reprocessing of data acquired between 2000 and 2020, we aim at sharing that consequent database with all research communities interested by GNSS data in South-America. Even though we are not an observatory and neither mandated nor funded to do so, we will do our best effort to maintain, upgrade and share this solution.

We would like to warmly thank the numerous people who have, over the years, installed and maintained the stations, and also made possible the free access to data and metadata. This work, and any of the great scientific contributions based on these GNSS data, could not have been done without their collective work. This includes (but not only), all the CSN staff, in particular J.C. Baez, all the RAMSAC staff, in particular Demian Gomez and Diego Piñón and all the RBMC staff. These acknowledgements are also directed to all the French and Chilean participants who contributed to the beginning of permanent GNSS development in Chile among whom (but not limited to, in alphabetical order) Carlos Aranda, Carolina Bermejo, Sylvain Bonvalot, Jaime Campos, Daniel Carrizo, Jean-Bernard de Chabalier, Germinal Gabalda, Raúl Madariaga, Sylvain Morvan, Manuel Olcay, Ismaël Ortega, Alain Rudloff, Jean-Claude Ruegg, Anne Socquet, as well the Caltech and OSU north-American teams, and the German GFZ teams. The authors also want to sincerely thank Mike Floyd for his very precious and constant help regarding the GAMIT/GLOBK process, as well as the whole MIT team developing the software. We also want to thank Paul Rebischung for facilitating access to the Repro3 solution and J.B. Besnier, T. Jezequel, T. Mesure and C. Peutin for their significant contribution in the development of the web applications. Finally, we sincerely thank our two reviewers Nicola D’Agostino and Alvaro Santamaria-Gomez for their very constructive reviews.

Maps are made with Generic Mapping Tools GMT (Wessel et al., 2013). Python tool-boxes used: PYACS+ PYEQ (https://github.com/JMNocquet/pyacs36), ASTROPY.

Cite this article as: Emilie Klein, Christophe Vigny, Jean-Mathieu Nocquet, Hugo Boulze. 2022. A 20 year-long GNSS solution across South-America with focus in Chile, BSGF - Earth Sciences Bulletin 193: 5.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.