Geodetic Monitoring of the Variable Surface Deformation in Latin America

Based on 24 years of high-level GNSS data analysis, we present a sequence of crustal deformation models showing the varying surface kinematics in Latin America. The deformation models are inferred from GNSS station horizontal velocities using a least-squares collocation approach with empirically determined covariance functions. The main innovation of this study is the assumption of continuous surface deformation. We do not introduce rigid microplates, blocks or slivers which enforce constraints on the deformation model. Our results show that the only stable areas in Latin America are the Guiana, Brazilian and Atlantic shields; the other tectonic entities, like the Caribbean plate and the North Andes, Panama and Altiplano blocks are deforming. The present surface deformation is highly inﬂuenced by the effects of seven major earthquakes: Arequipa (Mw8.4, Jun 2001), Maule (Mw8.8, Feb 2010), Nicoya (Mw7.6, Sep 2012), Champerico (Mw7.4, Nov 2012), Pisagua (Mw8.2, Apr 2014), Illapel (Mw8.3, Sep 2015), and Pedernales (Mw7.8, Apr 2016). We see very signiﬁcant kinematic variations: while the earthquakes in Champerico and Nicoya have modiﬁed the aseismic deformation regime in Central America by up to 5 and 12 mm/a, respectively, the earthquakes in the Andes have resulted in changes of up to 35 mm/a. Before the earthquakes, the deformation vectors are roughly in the direction of plate subduction. After the earthquakes, the deformation vectors describe a rotation counter-clockwise south of the epicentres and clockwise north of the epicentres. The deformation model series reveals that this kinematic pattern slowly disappears with post-seismic relaxation. The numerical results of this study are available at https://doi.pangaea.de/10.1594/PANGAEA.912349 and https://doi.pangaea.de/10.1594/PANGAEA.912350.


Introduction
Geodetic reference frames comprise coordinates of station positions at a certain epoch and constant velocities describing a secular station motion. In active seismic regions, strong earthquakes cause large displacements of station positions and velocity changes disabling the use of such coordinates over any time periods. The continuous representation of station positions between different epochs requires the computation of reliable station velocity models. Whit these models, we can monitor the kinematics of reference frames, determine transformation parameters between pre-seismic and post-seismic (deformed) coordinates, and interpolate surface motions arising from plate tectonics or crustal deformations in areas where no geodetic stations are established. In the particular case of Latin America, the reference frame is called SIRGAS (Sistema de Referencia Geocéntrico para las Américas;cf. SIRGAS 1997) and it is a regional densification of the global International Terrestrial Reference Frame (ITRF; Petit and Luzum 2010). The first SIRGAS realisation was established by a GPS (Global Positioning System) observation campaign in May 1995 (SIRGAS 1997). It comprised 58 stations covering all South America. This network was measured again in May 2000 and it was extended to Central and North America including 184 stations . Since 2000, the Latin American geodetic reference network is materialised (and frequently extended) by continuously operating GNSS (GPSCGLONASS) stations (Brunini et al. 2012;Sánchez et al. 2013Sánchez et al. , 2015Cioce et al. 2018). As the western margin of Latin America is one of the seismically most active regions in the world, the maintenance of the SIRGAS Reference Frame implies the frequent computation of present-day (updated) surface deformation models. Such models were computed in 2003 (Drewes and Heidbach 2005), 2009 (Drewes and Heidbach 2012), 2015 (Sánchez and Drewes 2016), and 2017 (this paper). Here, we present the computation of the deformation model 2017 and its comparison with the previous models to show the very significant variations of the surface kinematics in Latin America during the past 15 years.

Surface-Kinematics Modelling Based on GNSS Multi-Year Solutions
Spatial continuous surface deformation may be inferred from pointwise velocities applying geophysical models or geodetic methods based on mathematical interpolation approaches. The approach used in the present study is the least-squares collocation (LSC, e.g., Moritz 1973;Drewes 1978). Previous studies applied also the finite element method used with geophysical models (e.g., Heidbach and Drewes 2003). It has been demonstrated that for the sole representation of the horizontal Earth surface kinematics, the results of both methods are very similar (e.g., Drewes and Heidbach 2005). The vertical deformation is not considered in this work because the station height variations are highly influenced by local effects, and the station distribution ( Fig. 1) is too sparse to apply correlations between neighbouring stations. In the modelling of the surface kinematics, we distinguish two components: the velocity field and the deformation model. In the latter one, a secular motion inferred from plate motion estimates ( Fig. 1) is removed from the station velocities, and the pointwise residual velocities are interpolated to a regular grid.
The least-squares collocation method is based on the analysis of the correlation of physical quantities between neighbouring points. The vector of the observations (in this case the station velocities) is divided into a systematic part (trend) and two independent random parts: the signal and the observation error (or noise). The parameters describing the systematic component and the stochastically correlated signals are estimated by minimising the noise. The spatial signal correlation is usually assumed as a function depending on the distance d and, presuming isotropy after trend removal, independent of the direction. The basic LSC formula is given by Heidbach 2005, 2012): v obs contains the station velocities obtained from the GNSS observations at the geodetic stations. v pred represents the velocities to be predicted at the grid points. C obs is the correlation matrix between the observed velocities. C new is the correlation matrix between predicted and observed velocities. C nn is the noise covariance matrix (it contains the uncertainty of the station velocities obtained within the multi-year solutions). The correlation between the observed velocities v i , v k at the (adjacent) geodetic stations i, k is determined under the stationarity condition over a defined domain by E is the statistical expectation and d ik is the distance between stations i and k. The C obs values are classified in d j class intervals and the respective cross-covariance C obs (d j ) and auto-covariance C obs (d D 0) D C 0 are determined using: n stands for the number of stations available at the defined domain, while n j represents the number of stations available at each class interval d j . After estimating the discrete empirical covariance values with Eq. (3), they are approximated by a continuous function C(d ik ), which here is the exponential function: The function parameters a and b are estimated by a leastsquares adjustment. C obs is symmetrical and its main diagonal (i D k) contains the values C 0 . Fulfilling the stationarity condition, the elements of C new are computed using the same Eq. (4) as a function of the distance between the grid node to be interpolated and the geodetic stations.
To satisfy the isotropy condition, we estimate a common rotation vector and remove this horizontal motion trend from all station velocities located in the same domain defined by d. Afterwards, we restore the removed trend to the velocities v pred predicted at the grid points (cf. Drewes 1982Drewes , 2009).

Existing Velocity Models for SIRGAS (VEMOS)
The first velocity model for SIRGAS (VEMOS) was released in 2003. It is based on the position differences between the two SIRGAS campaigns of 1995 and 2000, 48 velocities derived from the SIRGAS multi-year solution DGF01P01 (Seemüller et al. 2002), and 231 velocities from several geodynamic projects based on episodic GPS campaigns (cf. Drewes and Heidbach 2005). The different data sets were transformed to a common kinematic frame by deriving the rotation vector of the South American plate from the respective station motions located in the rigid part of South America, and reducing these plate motions from the particular data sets. The resulting residual motions were modelled to a continuous deformation field applying the finite element method and LSC approach as described in the previous section. The comparison of the results of both methods shows an agreement in the mm/a level. VEMOS2003 covers the South American area between the latitudes 45 ı S and 12 ı N. The second VEMOS model was released in 2009 (Drewes and Heidbach 2012). It considers 496 station velocities; 95 of them corresponding to the SIRGAS multi-year solution SIR09P01 (Seemüller et al. 2011) and the others derived from repetitive GPS campaigns. It covers the Latin American area between the latitudes 56 ı S and 20 ı N and the time-span from January 2, 2000 to June 30, 2009. The continuous surface velocity field was derived applying the same strategies as in VEMOS2003. The main advantages of VEMOS2009 with respect to VEMOS2003 are the increased number of input velocities, the better quality of measurements (due to an increase of continuously operating GNSS stations), and the extension of the velocity field to the Caribbean and the southernmost part of Chile and Argentina. The mean uncertainty of VEMOS2009 is abouṫ 1.5 mm/a. After the Maule earthquake in Feb 2010, the station velocities in the area between latitudes 30 ı S and 40 ı S changed dramatically (Sánchez and Drewes 2016). However, we could not compute a new VEMOS model immediately, because we required 5 years of observations after the earthquake in order to improve the modelling of the strong post-seismic decay signals detected at the affected SIRGAS stations. Consequently, a new VEMOS model was computed in 2015 using the LSC method with station velocities based on GNSS observations captured from March 2010 to March 2015 (VEMOS2015, Sánchez and Drewes 2016). VEMOS2015 is based on continuously operating GNSS stations only; it does not include episodic GPS campaigns. It covers the region from 110 ı W, 55 ı S to 35 ı W, 32 ı N with a spatial resolution of 1 ı 1 ı . The average prediction uncertainty is˙0.6 mm/a in the north-south direction and˙1.2 mm/a in the east-west direction. The maximum uncertainty (˙9 mm/a) occurs in the Maule deformation zone (Chile), while the minimum (˙0.1 mm/a) appears in the stable eastern part of the South American plate.

Present-Day Deformation Model and Velocity Field for Latin America (VEMOS2017)
The present study concentrates on the computation of a deformation model based on a set of 515 station velocities inferred from GNSS observations gained from January 2014 to January 2017 (Fig. 2). Station positions and velocities are defined at epoch 2015.0 and refer to the IGS14 Reference Frame (Rebischung 2016), which is based on the latest ITRF solution, the ITRF2014 (Altamimi et al. 2016). The estimated precision is˙1.2 mm (horizontal) and˙2.5 mm (vertical) for the station positions at the reference epoch, anḋ 0.7 mm/a (horizontal) and˙1.1 mm/a (vertical) for the velocities (Fig. 2). More details about the processing strategy for the determination of the station positions and velocities can be found in Sánchez and Drewes (2016) and Sánchez et al. (2015).
The complex on-going crustal deformation in the western margin of Latin America and the Caribbean has been studied intensively. Recent research concentrates on geophysical syntheses including geodetic constraints inferred from GNSS positioning to model tectonic evolution and associated geodynamic processes in this region. Most of these studies assume a segmentation of the Earth's crust and describe the surface kinematics by means of tectonic blocks or slivers rotating individually; see e.g., Brooks 2016), and references herein. This paper presents two main innovations with respect to the abovementioned publications: Firstly, we compute a deformation model for the entire Latin American and Caribbean region and not for isolated areas only. Secondly, we assume a continuous lithosphere deforming under certain kinematic boundary conditions (as suggested by Flesch et al. 2000;Vergnolle et al. 2007;or Copley 2008), without introducing small lithospheric blocks or slivers, which would enforce constraints on the kinematic model. For the collocation procedure, we consider the main tectonic plates South America (SA), Caribbean (CA), and North America (NA) (Fig. 1) according to the tectonic plate boundary model PB2002 (Bird 2003). Based on the velocities obtained in this study for the stations located on the stable part of the plates, we estimate plate rotation vectors following the strategy presented by Drewes (1982Drewes ( , 2009). These plate motions are removed from the pointwise velocities to get the residual velocities, which are interpolated to a continuous deformation model (1)-(4). The residual velocities with respect to the Caribbean plate are used for the LSC prediction in Mexico, Central America and the Caribbean (Fig. 3), while the residual velocities with respect to the South American plate are used in South America (Fig. 4). The collocation domain at every grid node is created by selecting the stations located up to a distance of 200 km. If no stations are available at this distance, the LSC is computed using the three nearest stations. In total 2,233 grid points are predicted. Once the LSC prediction is performed, the previously reduced trends (plate rotations) are restored to the interpolated residual velocities at the grid nodes to generate a continuous velocity field referring to the IGS14 (ITRF2014). The average prediction uncertainty is˙1.0 mm/a in the north-south direction and˙1.7 mm/a in the east-west direction. The maximum uncertainty values (up to˙15 mm/a) occur at the zones affected by recent strong earthquakes, not only in the Maule area but also in the northern part of Chile, Ecuador and Costa Rica. The best uncertainty values (about˙0.1 mm/a) are evident in the stable eastern part of the South American plate.

Discussion
The deformation model with respect to the Caribbean plate (Fig. 3a) shows an inhomogeneous surface kinematics. While the deformation vectors in Puerto Rico and the Lesser Antilles show small (less than 0.5 mm/a) relative motions, the direction of the deformation vectors in Hispaniola describes a southward rotation starting with an orientation of S70 ı W in the northern part and reaching a south orientation in the southernmost part of the island. The magnitude of the vectors also decreases with this rotation: the averaged deformation is about 12 mm/a in the North and less than 1 mm/a in the South. These deformation patterns are in agreement with the GPS results published in earlier studies, e.g.; Benford et al. (2012), Symithe et al. (2015), Calais et al. (2016). In the southern area of Central America (Panama block), we observe horizontal deformations in the range 5-15 mm/a relative to the Caribbean plate. These large magnitudes are dominated in the West by the north-eastward motion of the Cocos plate towards Central America (see Fig.  3a around longitude 84 ı W) and in the East by the eastern motion of the Nazca plate towards South America (see Fig. 3a around longitude 78 ı W). A progressive westward rotation of the deformation vectors toward the North American plate is detected over Nicaragua and Honduras (longitudes from 85 ı W to 90 ı W), where the very small magnitudes of the deformation vectors suggest that this region moves homogeneously with the Caribbean plate. Figure 3b presents the differences between this model (VEMOS2017) and the previous one (VEMOS2015). The largest differences in magnitude (about 12 mm/a) are a consequence of post-seismic displacement and station velocity changes caused by the strong earthquake of Nicoya (Mw 7.6, Sep 5, 2012), Costa Rica, (marked with B in Fig.  3b). This earthquake produced co-seismic displacements up to 30 cm at the GNSS stations located in the Peninsula Nicoya (see Fig. 9 in Sánchez and Drewes 2016). The postseismic relaxation process induces pre-and post-seismic station velocity differences up to 30 mm/a. Another relevant discrepancy between VEMOS2017 and VEMOS2015 is observed in Guatemala (marked with A in Fig. 4). In this case, the difference in the deformation magnitude (about 5 mm/a) is mainly caused by the Champerico earthquake (Mw 7.4, Nov 11, 2012).
The deformation model with respect to the South American plate (Fig. 4a) clearly defines the stable area belonging to the Guiana, Brazilian and Atlantic shields. Indeed, the present VEMOS2017 and the previous model VEMOS2015 are practically identical in this area (Fig. 4b). In contrast, the deformation vectors predicted in the Andean region are characterized by magnitudes up to 30 mm/a. These vectors are roughly parallel to the plate subduction direction and their magnitudes diminish with the distance from the subduction front as already stated by previous publications like Bevis et al. (2001), Brooks et al. (2011), Chlieh et al. (2011, Khazaradze and Klotz (2003) and references herein. However, we observe three zones with anomalous vector directions (oriented to the NW): the western part of Ecuador around latitude zero, the north of Chile around latitude 20 ı S, and the Maule region (around 38 ı S). As in the case of Central America, these abnormalities are also caused by recent strong earthquakes and post-seismic relaxations (Fig.  4b).
The surface deformation predicted for the North Andes (ND) block is characterized by two different kinematic patterns: a north-eastward motion with increasing magnitudes of about 9 mm/a in the southern part of Colombia (latitude 3 ı N) to 15 mm/a in the northern border area with Venezuela (72 ı W, 12 ı N); and opposite oriented deformation vectors in Ecuador (south of latitude 3 ı S). The latter is a consequence of the strong earthquake occurred in Pedernales (Mw 7.8) on Apr 16, 2016. This earthquake produced co-seismic station displacements up to 80 cm and station velocity changes of about 40 mm/a (see Fig. 4b, mark A). The differences between VEMOS2017 and VEMOS2015 in this area come up to 22 mm/a. South of this region, the poor station coverage in central Peru (latitudes 5 ı S to 12 ı S) prevents concluding statements about the deformation pattern in this area; however, our model agrees quite well with the findings published by Nocquet et al. (2014) and . Based on about 100 GNSS stations covering the area between latitudes 12 ı S and 4.6 ı N, they conclude that the southern Ecuadorian Andes and northern Peru (between latitudes 5 ı S and 10 ı S) move coherently 5-6 mm/a with an orientation of about S70 ı E. They also suggest that the internal deformation in this area is negligible (see Fig. 2a in Nocquet et al. 2014).
South of latitude 15 ı S the deformation model ( Fig. 4a) and its comparison with the previous one ( Fig. 4b) Fig. 2). This produces an apparent smaller deformation with respect to the South American plate and the differences between both VEMOS models reach magnitudes up to 20 mm/a (mark B in Fig. 4b). In the area Illapel (mark C in Fig. 4b), the postseismic effects of the 2015 earthquake superimpose the post-seismic effects of the 2010 Maule earthquake (mark D in Fig. 4b). Thus, it is not possible to distinguish their individual contributions to the deformation. As a matter of fact, the complex kinematic pattern south of latitude 25 ı S described by Sánchez and Drewes (2016, Fig. 18) persists. A large counter clockwise rotation around a point south of the 2010 epicentre (35.9 ı S, 72.7 ı W) and a clockwise rotation north of the epicentre are further observed (Fig.  4a). However, magnitude and direction of the deformation vectors considerably differ from those obtained in the previous model VEMOS2015. This is probably a consequence of the post-seismic relaxation process that is bringing the uppermost crust layer to the aseismic NE motion in this zone as suggested by e.g., Bedford et al. (2016), Klein et al. (2016) and Li et al. (2017). The surface kinematics shown in Fig. 4a again makes evident that the deformation regime imposed by the Maule earthquake reaches the Atlantic coast in Argentina. The comparison of the present deformation model with VEMOS2015 in the Maule surroundings presents discrepancies up to 25 mm/a (marks C and D in Fig. 4b).
To provide an integrated view of the changing surfacekinematics in the Andean Region, Fig. 5 presents an extract of the models VEMOS2003, VEMOS2009, VEMOS2015 and VEMOS2017.

Conclusions and Outlook
This paper presents the surface velocity and deformation models of the entire Latin American and Caribbean region over the time-span 2014-2017 and describes the evolution of the models from previous studies. The effects of the extreme changes in the surface kinematics complicate the long-term stability expected in any reference frame. Therefore, a major recommendation is to materialise the geodetic reference frames by means of a dense network of continuously operating stations and to repeat the velocity computations frequently. This ensures a permanent monitoring of possible reference frame deformations. Nevertheless, a reliable deformation modelling is not yet guaranteed. Some authors suggest the implementation of geodynamic models to predict the pointwise coordinate changes caused by coseismic and post-seismic effects (see e.g., Snay et al. 2013;Bevis and Brown 2014;Gómez et al. 2015). Since these models rely on hypotheses about the physical properties of the upper Earth crust, different hypotheses produce different results as demonstrated by e.g., Li et al. (2017). We based our analyses on the least-squares collocation as this approach respects the consistency of the geodetic observations and ensures a better agreement with the actual deformation. A problem in the geodetic use of pointwise velocities derived from multi-year solutions is their inconsistency after seismic events, i.e. their short-term validity. In the Andes region, like in any active seismic region of the Earth, there are large discontinuities in the station coordinate time series and considerable variations in the station velocities caused by strong earthquakes. The consequence is that the respective reference frames (e.g., ITRF) cannot be used or have to be frequently updated for geodetic purposes (like SIRGAS). An alternative of using multi-year solutions with station velocities is the release of frequent reference frames (e.g., every week or month). Our recommendation for the SIRGAS national reference frames in seismic active regions is to use the SIRGAS weekly coordinate solutions instead of velocities after seismic events. To consider discontinuities in the coordinates of non-permanently observed points, one has to interpolate them from the coordinate differences in reference stations. In any case, we shall continue the computation of short-period velocity and deformation models for the next future in order to enable the use of coordinates in close alignment with to the IGS reference frames.

Supplementary Data
In the preparation of the GNSS data solutions used in this study, we computed a new SIRGAS reference frame solution following the same procedure described in Sánchez and Drewes (2016