Analysis of local surface displacement using repeated GPS measurements: a case study of the Guabirotuba area, Curitiba, Brazil

Abstract:

In recent years, the ability of GNSS systems to estimate high accuracy 3D positions has greatly contributed to the development of new monitoring techniques aimed at understanding the mechanisms by which the different Earth processes are generated. This paper presents the result of 9 GPS campaigns, between 2014 and 2019 for a group of 29 stations, held in the Guabirotuba urban area, Curitiba, Brazil, intending to estimate local crustal movements. The average magnitude and direction of horizontal velocities obtained for each site, allow demonstrating deformation of the southern zone of the area as a result of a local landslide process. The magnitude of the vectors varies between 1 and 16 mm/a horizontally. All sites had an absolute vertical movement, probably correlated with a local geological pattern or a physical site’s motion driven by environmental mass redistribution. The uncertainties in positions and velocities for a long-term survey (over 5 years), showed good consistency at most sites when compared with the coordinates and velocities precision of the latest SIRGAS solution. The typical precision for the station positions at the reference epoch was ±2 mm horizontally and ±5 mm vertically, and ±0.60 mm/a horizontally and ±1.37 mm/a vertically for the constant velocities.

Keywords:
GNSS; GPS; surface displacement; trajectory model; velocity field; GAMIT-GLOBK

1. Introduction

Recent advances in the development of the different satellite positioning systems, have led them to become a key element in the study of the different physical Earth processes. The ability of Global Navigation Satellite Systems (GNSS) to estimate high accuracy 3D positions has contributed greatly to the development of new monitoring techniques aimed at understanding how different tectonic and geological processes affect the Earth surface and fundamentally, to explain how these mechanisms impact the society and the environment (Bock and Melgar 2016). One of the main applications that currently make use of GNSS is the study and monitoring of geophysical phenomena (e.g., tectonic motion, fault zones, landslides, subduction processes, and volcanoes) (Luna et al. 2017, Duman and Sanli 2019, Sánchez and Drewes 2020). In this scenario, GNSS observations play an important role, since they provide estimates of crustal deformation with respect to a global terrestrial reference frame, while also provide constraints of the nature of the forces acting on the Earth surface (Lambeck 1988). Most of the innovative studies conducted in recent years, in which GNSS velocity fields are used to analyse tectonic and geodynamic processes, are based on the use of continuous GPS data (i.e., 24 h sessions) (Zhang et al. 1997, Mao, Harrison, and Dixon 1999, Dixon et al. 2000, Langbein 2004, Vernant et al. 2004, Serpelloni et al. 2013, Williams 2015, Choudhury et al. 2018, Duman and Sanli 2019). In all of them, their objective is the study of the stochastic properties of data in order to minimize the effect of long-range time-dependent correlation noise and establish a more realistic uncertainty of the velocities (Williams 2003, Pilapanta, Krueger, and Criollo 2018). Despite the good results achieved with this method, the need to study geophysical processes at local scales (i.e., landslide and subduction urban area processes) and the impossibility of carrying out campaigns with long duration (i.e., 12 to 24h session), has led to the use of so-called “repeated GPS measurements” (i.e., GPS campaign measurements with just a few hours of observation). An example of this is the study performed by Duman and Sanli (2019), who analysed the incidence of employed GPS campaign measurements with only 1-8 hours observation sessions for the calculation of velocity fields compared to the use of continuous data. The results revealed that through the definition of a robust local geodetic network formed by the combination of permanent (continuous measurements) and semi-permanent (repeated measurements) sites, 85% of the horizontal velocities were comparable to the accuracy derived from 24 h campaigns. Based on these studies and to analyse possible landslide processes, in this paper we examine the data obtained in nine survey-mode GPS campaigns carried out between 2014 and 2019 for a group of 29 GPS sites, divided into a survey (multi-period GPS sites) and a core (continuous GPS sites) local network and with sessions lasted from 1 to 24 hours.

2. Geographical and geological setting

The study area is located in the Guabirotuba urban area, in the central region of Curitiba city, Paraná, Brazil (Figure 1). Geologically, it is located in the Curitiba Sedimentary Basin, which is situated in the central-southern portion of the Curitiba plateau, Paraná (Giusti 1989). It comprises an area of approximately 3000 km² with a maximum thickness of 80 meters and altitudes ranging from 750 to 980 meters (Salamuni et al. 2003). The Basin is subdivided into four main morphological domains: 1) plains and alluvial terrains, 2) hilly zones with elongated and planar tops, 3) a mountain system and 4) mound systems and karst topography; and stratigraphically organized into three formations: 1) Guabirotuba Formation, 2) Tingui Formation and 3) Younger colluvial and alluvial sediments of the valley (Salamuni et al. 2003).

Figure 1:
Location map of the study area in the central region of Curitiba city, Paraná, Brazil.

As mentioned by Werlich (2017), several studies carried out in this area by the Interdisciplinary Centre for Environment and Development of the Federal University of Paraná (NIMAD / UFPR) show that due to its geomorphological composition and the presence of various geological processes in the vicinity of the Guabirotuba River, make this one of the zones with the highest incidence of landslides in the city.

Figure 2 shows three important points of the study area, where road and structure issues related to possible landslide processes can be observed. In Figures 2a and 2b, for instance, it is possible to observe cracks in the walls of two houses (GPS sites M014, M012 respectively) probably caused by earth and debris flows. Figure 2c, for its part, shows the collapse of the road that passes over the Guabirotuba River, as a result of a potential erosion process (Werlich 2017).

Figure 2:
Evidence of landslide process generated in the vicinity of Guabirotuba River.

3. GPS data analysis

3.1 Data sources

Nine survey-mode GPS campaigns were carried out between 2014 and 2019 for a group of 29 stations divided into two hierarchical networks: a survey local network (S-GPS) composed of a set of 26 multi-period GPS stations with sessions lasted from 1 up to 4 hours and a core local network (C-GPS) composed by a set of four semi-permanent stations, including an IGS Station (UFPR 41610M002), with sessions of 24 hours (Figure 1). The approximate distance between the survey local network (S-GPS) and the core local network (C-GPS) is 2 kilometres. In all cases, observations were conducted at a sampling rate of 30 seconds and an elevation mask of 10°. A more detailed explanation of different campaigns and sessions made are described in Table S1.

The main difference between S-GPS and C-GPS networks lies in the session characteristics (i.e., campaign and session duration), which have been defined based on the GPS site monumentation, type of receiver and antenna used, and the characteristics of the study area. Thus, for the C-GPS network, all stations have reinforced concrete monuments with forced-centring apparatus for antennas (Figure 3a and 3b). For the S-GPS network, being a non-continuous network, its stations are materialized through a tripod-mounted station placed directly on top of a survey benchmark (approximately 2 m above the Earth’s surface) (Figure 3c).

Figure 3:
Monumentation of sites: a) UFPR station; b) BCAL station and c) M018 site.

Regarding the type of GPS antennas used, in the case of the C-GPS network, a set of choke-ring ground plane antennas was implemented to mitigate multipath interference. For all other stations, a set of LEICA ATX900CS and ATX1230GG micro-strip Smart-Track antennas were used. A list of all receivers and antennas types used in the different campaigns is provided in Table S2. In the case of the LEICA ATX900CS model, no IGS calibration parameters exist (w.r.t. IGS14.atx), therefore, a type-specific calibration model has been calculated (Section S3). For all other antennas, an individual IGS calibration was used.

3.2 GPS data processing

All the data from the different GPS campaigns were processed using the GAMIT-GLOBK software, version 10.71 (Herring et al. 2018). First, daily loosely constrained solutions (D-LCS) were computed using a “baseline” strategy, in which, orbit parameters are considered fixed and only site coordinates, atmospheric zenith delays (ZTD), and Earth orientation parameters (EOP) are solved using IGS final orbits. The Global Pressure and Temperature 2 (GPT2) slant delay model (Boehm, Heinkelmann, and Schuh 2007, Lagler et al. 2013) and Global Mapping Function (Boehm et al. 2006) were used to mitigate the signal delay in the neutral atmosphere. Tropospheric parameters were estimated at each station every two hours. To ensure the accurate modelling of ionospheric effects during periods of high solar activity, the second and third-order effects were also modelled using the daily IGS TEC maps (Hernández-Pajares et al. 2009) and the approach described by Petrie et al. (2010). The ocean tide loading was estimated with the FES2004 model (Lyard et al. 2006) and the atmospheric tide loading caused by the semidiurnal constituents S1 and S2 were estimated following the model of Ponte and Ray (2002). Other effects, such as receiver/satellite antenna phase offsets (section 3.2) and non-loading effects, were modelled according to International Earth Rotation Systems Conventions (i.e., IERS Conventions 2010) (Petit and Luzum 2010). In a second stage, station positions and velocities for each site were estimated from the combination of the D-LCS solutions into a single one (C-LCS) (i.e., one per campaign) using a smoothing Kalman Filter (Herring et al. 2015).

To align the positions of the estimated sites to an ITRF realization, a set of 10 SIRGAS-CON stations were included (i.e. IMBT, NEIA, POLI, PRGU, PRMA, SCAQ, SCFL, SCLA, SPBO and SPPI stations). These are all part of the latest SIRGAS multi-year solution (i.e., SIR17P01), which is aligned to the IGS14 reference frame, epoch 2015.0 (Sánchez and Drewes 2020). The stabilization process of the individual solutions was performed by using generalized constraints. That is, a set of core sites with good a priori values (i.e., SIRGAS-CON stations) was selected and minimized their adjustments while estimating an overall translation and rotation (Helmert parameters) of the network in an iterative process (Herring et al. 2015). Finally, the internal consistency of solutions was evaluated by using the χ2 statistic values (i.e., post-fit observed minus computed observation) as established by Herring et al. (2018) and their reliability (external precision) with respect to the SIRGAS coordinates, through the generation of a minimum constraint solution for the network as established by Sánchez et al. (2015). The processing strategies described above were applied to each subset of sessions listed in Table S1.

4. TIME SERIES ANALYSIS

4.1 Trajectory model estimation

As established by Floyd and Herring (2020), rigorous estimation of the uncertainties of the different GNSS parameters (e.g., secular velocities, seasonal cycles, offsets, transient motions, and others) requires full knowledge of the error spectrum and both temporal and spatial correlations. Thus, for the description of site positions, the standard linear trajectory model (SLTM) defined by Bevis and Brown (2014) was adopted. For that, we assume that the positions of the sites taken at the time t _i , are represented by the sum of a deterministic model and stochastic noise, thus:

$x=\sum _{i=0}^{n_P}{P}_i{\left(t-t_R\right)}^i+\sum _{i=1}^{n_J}{b}_{j}H\left(t-t_j\right)+\sum _{i=1}^{n_F}{s}_i\text{sin}\left({\omega }_{i}t\right)+c_i\text{cos}\left({\omega }_{i}t\right)+{\epsilon }_x\left(t_i\right)$ (1)

where P _i are the coefficients of a n _P degree polynomial, H(t) is the Heaviside or unit step function, and vector b _j describes the direction and magnitude of the jump (e.g., the effect of the datum change) which occurs at the time tj (Bevis, Bedford, and Caccamise II 2020). To represent periodic signals, an annual and semi-annual term was also included s _i and c _i , defined as the Fourier coefficients for the harmonic with angular frequency ${\omega }_i=\frac{2\pi }{{\tau }_k}$. The term ( _x (t _i ) in the equation represents the error model (i.e., stochastic noise) assumed for data, which is represented as a linear combination of independent, identically distributed (IID) unit-variance variables, ((t _i ), and a sequence of temporally correlated random variables, ((t _i ):

${\epsilon }_x\left(t_i\right)=a\alpha \left(t_i\right)+b_{\kappa }\beta \left(t_i\right)$ (2)

The scale factors a and b _k are defined as the variance of white noise and coloured noise (power-law processes other than classical white noise) of spectral index k respectively (Pilapanta, Krueger, and Criollo 2018). The covariance matrix for the measurements x _i can be expressed as:

$C_x=a^{2}I+{b^2}_{\kappa }{J}_k$ (3)

where I is the identity matrix and J _k is the covariance matrix for the appropriate coloured noise. For the estimation of J _k , a maximum likelihood estimator (MLE) was used. In this case, it was necessary to maximize the probability function of the given observations x _i by adjusting the covariance matrix C _x (Bos et al. 2013, Pilapanta, Krueger, and Criollo 2018). Thus, the likelihood function for J _k , can be obtained from:

$\mathrm{l}\mathrm{n}\left(l\left(x,J_k\right)\right)=-\frac{1}{2}\left[N\mathrm{l}\mathrm{n}\left(2\pi \right)+\mathrm{l}\mathrm{n}\left(\mathrm{d}\mathrm{e}\mathrm{t}\left(J_k\right)\right)+{\widehat{v}}^T{J_k}^{-1}\widehat{v}\right]$ (4)

where N is the number of data, det is the determinant of a matrix and $\widehat{v}$ are the post-fit residuals. Finally, using a generalized least square estimator (LSE), the x ₀ and r values are obtained from:

$\widehat{y}={\left[A^T{C}_x^{-1}A\right]}^{-1}{A}^T{C}_x^{-1}x$ (5)

where x is the vector of observations and A is the matrix of regressors, also called the design matrix. The covariance matrix for the estimates is given by:

${\widehat{C}}_{\widehat{y}}={\left[A^T{C}_x^{-1}A\right]}^{-1}$ (6)

4.2 Treatment of discontinuities and outliers’ detection

To reduce the impact of unusually noisy measurements in the SLTM (i.e., atypical values), the 3σ-rule for outlier detection was used. Thus, a linear trend was fit to the data using an ordinary least-squares estimator and afterward, it was subtracted from the observations to create coordinate residuals. Then, the interquartile range (IQR) value was computed and any residual with an absolute value more than three times the IQR, below or above the median, was considered to be an outlier (Langbein and Bock 2004). If outliers appeared sporadically, the residual was removed from the series, but if outliers reflected a discontinuity, a new offset value was incorporated to term. An important aspect considered in this stage was the incidence in the series of the change of reference system in the IGS products (i.e., update from IGS08 to IGS14). Thus, according to Montecino et al. (2017), the effect of the datum change series can be considered to be a “jump” within the trajectory model. Finally, to account for time-correlated sources of error in the estimates of station position and to allow for more realistic velocity estimates, a random-walk noise at the level of 2 mm/a in the horizontal coordinates and a vertical down-weight by a variance factor of 10 in the determination of reference frame transformation parameters were applied.

5. RESULTS AND DISCUSSION

5.1 Analysis of solution reliability

As described in 3.2, to evaluate the reliability of final solutions, different comparisons were performed. Figure 4 summarize some statistics of this comparison in terms of the uncertainties obtained for each campaign.

Figure 4:
Repeatability and χ2 (chi-square) increments per session.

According to that, the repeatabilities obtained for each campaign are well below the minimum levels established for GPS survey sessions (i.e., ±2-4 mm for horizontal coordinates and ±10-15 mm for heights) and, except for the sixth campaign (Up component), in the same magnitude of the SIRGAS coordinate accuracy (i.e., ±2 mm for horizontal coordinates and ±5 mm for heights) (Sánchez et al., 2015). Therefore, the stabilization process can be considered robust, and consequently, all individual solutions properly aligned to the SIRGAS reference frame (i.e., SIR17P01 solution). Concerning χ² (chi-square) increments, typical values of adequately loosely-constrained solutions can be observed for all sessions, except for the sixth campaign again, where anomalous increment is presented (i.e., χ² > 1). As described by Herring et al. (2018), χ² values greater than one occurs when the sigma adopted for each solution are too small (i.e., sites are improperly weighted) or there are uncompensated outliers among the coordinate estimates. Since for the stabilization, no particular weighting factors have been used (i.e., all solutions have the same relative weight), the main reason for the high values, and consequently the degradation of its precision, should be related to the presence of a discontinuity (i.e., physical process) or inconsistency (i.e., systematic process) in the data, which can be verified in the time series.

5.2 Testing consistency of time series

Figure 5, shows the coordinate time series of 2 of the 8 sites belonging to the S-GPS network, that are located in the southern zone of the study area (Figure 1). The rest of the series can be found in the Section S4.

Figure 5:
Coordinate time series of M004 and M015 S-GPS sites.

In the series, a clear discontinuity corresponding to the sixth campaign (red arrow) is observed (mainly at the East component). By the way the discontinuity appears in the series, and mainly its magnitude, suggests that it may be correlated with a physical process as a result of an isolated local dynamics.

It is important to note that other stations located in the same area (e.g., M005, M011, M014, M18 sites) (Figure S3 and S4), also present a similar perturbation for the sixth campaign, but in many of them, the large uncertainties of data make that the discontinuity have a low weight in the inversion/Kalman-filter process (i.e., the scale of the outlier close to the error) and therefore it remains in the time series. This consistency of variations between nearby sites, as demonstrated by Herring et al. (2015), proves the effect of reference frame stability and suggests that apparent slight outlier may likely be an effect of poor reference frame stabilization or an irregular perturbation in the data. To evaluate this precept, Table 1 shows the statistics obtained for the final iteration performed for each campaign as part of the stabilization process.

Table 1:
Position system stabilization statistics per campaign.

Campaign	Sites used in stabilization	Pre RMS (m)	Post RMS (m)
1	9	0.02710	0.00166
2	9	0.04900	0.00130
3	9	0.03590	0.00169
4	9	0.02170	0.00124
5	9	0.03340	0.00159
6	8	0.05360	0.00165
7	8	0.01880	0.00176
8	7	0.04340	0.00193
9	7	0.03410	0.00191

As Table 1 showed, a minimum of seven stations was used in all iterations of the reference frame solution. It ensures the consistency of the stabilization process and allows the reduction of the “Post RMS” at the expected level for a local/regional adjustment (i.e., approximately 1 to 5 mm). Consequently, it proves the reference frame reliability of each solution and assumes, as the most possible explanation of the discontinuity in the time series, the presence of an irregular perturbation in the data as a result of an extremely physical event (e.g., non-linear motion, landslide process or seismological event) or a systematic error due to equipment or environment changes, the latter considered previously eliminated.

5.3 Velocities sites estimation and local movement analysis

First of all, and in order to evaluate the consistency of our velocities with reference to SIRGAS velocities, a direct comparison was performed between them and those obtained in the latest SIRGAS solution (i.e. SIR17P01) (Sánchez and Drewes 2020). Figure 6 and Tables 2, show the velocities and their absolute differences obtained for each of the 10 SIRGAS-CON stations used in the stabilization process and the UFPR station, which is also part of the SIR17P01 solution.

Figure 6:
Differences between velocities solutions for each of the 11 SIRGAS-CON stations.

According to that, mean absolute differences of 0.25 mm, 0.32 mm and 0.45 mm can be obtained for the East, North and Up component, respectively, which, in principle, can be considered statistically equivalent at the 0.05 significance level (Table S3).

Table 2:
SIRGAS-CON station velocities (Confidence interval: 95%).

Station	Estimated velocities						SIR17P01 solution
	East (mm/a)		North (mm/a)		Up (mm/a)		East (mm/a)		North (mm/a)		Up (mm/a)
	rate	sigma	rate	sigma	rate	sigma	rate	sigma	rate	sigma	rate	sigma
IMBT	-3.19	0.12	13.47	0.18	-1.59	0.43	-3.20	0.50	13.70	0.60	-2.00	0.70
NEIA	-3.45	0.82	12.48	0.74	-3.13	1.41	-3.30	0.70	12.20	0.70	-3.90	0.70
POLI	-3.70	1.10	12.70	0.22	1.72	0.46	-3.80	1.20	12.70	0.90	2.40	0.90
PRGU	-2.76	0.19	12.83	0.13	-0.67	0.45	-2.80	0.80	12.70	0.70	-0.20	0.80
PRMA	-2.79	0.82	13.75	0.24	-1.32	0.48	-3.00	1.30	13.80	0.80	-1.80	1.00
SCAQ	-3.04	0.90	12.32	0.80	0.87	0.35	-3.00	1.30	12.10	1.10	0.80	1.20
SCFL	-2.50	0.42	12.45	0.14	0.08	0.41	-3.10	1.30	12.20	0.90	-0.60	1.10
SCLA	-2.90	0.12	11.80	0.15	-1.27	0.39	-2.60	0.50	12.10	0.60	-0.80	0.70
SPBO	-2.43	0.16	11.97	0.27	-0.88	0.66	-2.10	1.00	11.30	0.80	-0.40	0.80
SPPI	-3.47	0.31	12.63	0.37	-1.32	1.06	-3.20	1.30	12.10	1.10	-1.20	1.10
UFPR	-2.52	0.25	13.92	0.24	-1.11	1.05	-3.20	0.50	13.10	0.60	-0.80	0.60

This can also be evidenced through the repeatability analysis. For that, Table 3 shows the repeatability of the final combined solution, which also are in accordance with minimum levels established for GPS survey sessions (i.e. ±2-4 mm for horizontal coordinates and ±10-15 mm for heights) and in the same magnitude of the SIRGAS coordinate and velocity accuracy (i.e. ±2 mm horizontally and ±5 mm vertically for the station positions at the reference epoch, and ±0.8 mm/a horizontally and ±1.2 mm/a vertically for the constant velocities), which ensures the reliability of the solutions.

Table 3:
Repeatability (coordinates + velocities) for the final combined solution.

Parameter	Core sites used in stabilization	Position			Velocities
		E	N	U	E	N	U
		(mm)	(mm)	(mm)	(mm/a)	(mm/a)	(mm/a)
WRMS	8	1.72	2.28	4.91	0.66	0.54	1.37
NRMS	8	3.82	4.85	3.49	3.23	3.06	3.11

Figure 7 shows the horizontal velocity vectors obtained for each of the S-GPS sites. The error ellipses of the absolute velocities (black arrows) are computed at 95 per cent confidence level taking into account the white and coloured noise components (Mao, Harrison, and Dixon 1999, Herring et al. 2015) as described in section 4.1. A table with all velocities (including vertical rates) obtained for each site are presented in Table S4.

Figure 7:
Velocity field map of S-GPS network.

The average magnitude and direction of horizontal velocities obtained for each site show a common pattern for most of them, except for some sites located in the southern part of the study area, where the magnitude of the deformation vectors varies gradually between 1 and 16 mm/a, as it moves to the NW away from the river zone, probably as a result of a local landslide process or the geomorphological composition of the area. It is important to point out that, due to the large uncertainties values of the vertical velocities of some of the sites (Table S4), which, in principle, would be correlated with the implementation of a down-weight variance factor in the determination of reference frame transformation parameters (section 4.2), these did not provide conclusive results in this analysis.

On the other hand, and in order to reveal the localized deformation of sites, the relative velocities (red arrows) and their uncertainties with respect to the UFPR station was calculated (i.e., UFPR velocity subtracted from each site). In principle, the uncertainties of relative velocities are smaller than the uncertainties of absolutes velocities, due to statistical correlation between deformation sites (M codes) and the reference site (UFPR), which allow us to establish a more realistic velocities value (i.e., with more reliable uncertainties). Table 4 shows the relative velocities of the sites with respect to the UFPR station. From this, a clear deformation of sites M12, M13, M14 and M16 can be observed, which allows confirming that stations near the river are significantly deforming but those far from the river are not. An important aspect to consider is that related to the magnitude and direction of the M18 site, which differ significantly from those values obtained for the other sites. Due to its high uncertainty and in order not to bias our analysis, this site will be considered a blunder error and therefore removed from our analysis.

Table 4:
Relative velocities with respect to UFPR station (Confidence interval: 95%).

Station	East		North		Up
	Rate	Sigma	rate	Sigma	rate	sigma
	(mm/a)	(mm/a)	(mm/a)	(mm/a)	(mm/a)	(mm/a)
A. C-GPS Network
BCAL	-0.04	0.27	0.02	0.28	1.90	1.16
HORT	2.84	2.23	17.45	2.10	-57.23	9.24
SANE	1.08	0.13	-2.08	0.15	1.06	0.60
B. S-GPS Network
M001	-2.53	0.94	2.46	0.91	-8.56	5.78
M002	0.49	1.12	-0.86	1.22	-1.71	6.60
M003	-0.34	1.14	-3.40	1.18	-1.27	6.12
M004	3.88	1.26	2.85	1.15	-32.25	10.11
M005	0.72	1.32	-0.94	1.31	-2.30	6.71
M006	1.86	2.46	-2.15	2.83	-11.70	16.42
M007	-2.40	0.96	2.64	1.01	-3.54	5.11
M008	-1.21	0.90	-1.60	0.89	-4.15	5.16
M009	-1.11	0.71	-1.82	0.80	-3.37	3.37
M010	-1.46	0.70	-4.51	0.70	-3.81	3.66
M011	-0.16	0.66	0.01	0.68	-2.63	3.39
M012	-11.06	3.98	4.21	3.76	-4.09	21.47
M013	-13.50	0.68	-5.23	0.76	-5.60	3.55
M014	-11.15	1.24	-2.63	1.34	-31.18	8.76
M015	-0.60	1.16	-5.41	1.17	-11.55	6.50
M016	-11.15	0.68	-7.04	0.68	-17.35	3.35
M017	0.58	0.73	-2.49	0.80	-0.67	3.68
M018	-50.63	3.30	-52.35	3.07	-8.93	18.47
M019	-1.75	0.78	-0.75	0.82	-1.83	3.99
M020	-7.25	18.09	2.85	17.22	10.16	97.41
M021	-1.08	1.38	3.28	1.53	4.99	7.39
M022	1.03	2.41	-4.37	2.48	-2.68	12.51
M023	0.62	20.36	11.85	21.48	-11.27	104.45
MAR1	5.87	3.47	10.52	3.42	-11.96	20.04
MAR3	1.68	2.60	-1.73	2.81	21.17	14.60

6. CONCLUSION

This work presented a simple, but a concise method for analysis of local surface displacement using GPS observations. In general, it was possible to show that local landslide monitoring can be carried out through the use of repeated GPS measurements with short sessions lasted from 1 up to 4 hours. In this way, when the final positional solutions were compared with the SIRGAS solutions, a high level of consistency between them could be verified, which in turn allowed to verify their reliability for the study of landslide processes.

The use of a type-specific calibration model for the antenna LEICA ATX900CS adhered to the requirements established by the IGS and the NGL for modelling antenna phase centre variations. This was reflected in the phase residuals obtained for each of 26 S-GPS sites that used the LEICA ATX900CS antenna model, which remained relatively constant over time (i.e., close to zero), with the higher oscillating phase residuals (up to 10 mm) at low elevations, probably, as a result of multipath effects, but generally without influence on station displacement estimates.

Regarding the velocity field obtained for the study area, a local displacement of the southern zone near the Guabirotuba River could be observed, probably as a result of a local landslide process or the geomorphological composition of the area. A clear discontinuity is also observed in the vertical component of the time series. This is probably correlated to the existence of a local physical process, which in turn produce relative velocities of several centimetres per year in the up component of some stations (that is, M004, M006, M015, M016) (Table 4).

Finally, it was shown that landslide monitoring can be carried out through repeated observations collected by the GNSS system, provided that the different sources of error present in positioning are minimized or correctly modelled (e.g., antenna phase centre model, atmospheric effects, and loading processes) and a rigorous analysis of the data (e.g., stochastic properties and outliers’ detection) is carried out to obtain reliable results.

ACKNOWLEDGMENT

The authors would like to thank the Coordination for the Improvement of Higher Education Personnel (CAPES) for their financial support. The authors also would like to thank Dr. Renata Werlich Cavalcanti for her help and contribution in the fieldwork as well as data collection and availability.

REFERENCES

• Bevis, M. Bedford, J. and Caccamise II, D.J. 2020. The Art and Science of Trajectory Modelling. pp.1-27. DOI:10.1007/978-3-030-21718-1_1.
  » https://doi.org/10.1007/978-3-030-21718-1_1
• Bevis, M . and Brown, A. 2014. Trajectory models and reference frames for crustal motion geodesy. Journal of Geodesy, 88 (3), pp.283-311. DOI:10.1007/s00190-013-0685-5.
  » https://doi.org/10.1007/s00190-013-0685-5
• Bock, Y. and Melgar, D. 2016. Physical applications of GPS geodesy: A review. Reports on Progress in Physics, 79 (10), pp.106801. DOI:10.1088/0034-4885/79/10/106801.
  » https://doi.org/10.1088/0034-4885/79/10/106801
• Boehm, J. Heinkelmann, R. and Schuh, H. 2007. Short note: A global model of pressure and temperature for geodetic applications. Journal of Geodesy , 81 (10), pp.679-683. DOI:10.1007/s00190-007-0135-3.
  » https://doi.org/10.1007/s00190-007-0135-3
• Boehm, J . Niell, A. Tregoning, P. andSchuh, H . 2006. Global Mapping Function (GMF): A new empirical mapping function based on numerical weather model data. Geophysical Research Letters, 33 (7), pp.3-6. DOI:10.1029/2005GL025546.
  » https://doi.org/10.1029/2005GL025546
• Bos, M.S. Fernandes, R.M.S. Williams, S.D.P. and Bastos, L. 2013. Fast error analysis of continuous GNSS observations with missing data. Journal of Geodesy , 87 (4), pp.351-360. DOI:10.1007/s00190-012-0605-0.
  » https://doi.org/10.1007/s00190-012-0605-0
• Choudhury, P. Gahalaut, K. Dumka, R. Gahalaut, V.K. Singh, A.K. and Kumar, S. 2018. GPS measurement of land subsidence in Gandhinagar, Gujarat (Western India), due to groundwater depletion. Environmental Earth Sciences, 77 (22), pp.0. DOI:10.1007/s12665-018-7966-5.
  » https://doi.org/10.1007/s12665-018-7966-5
• Dixon, T.H. Miller, M. Farina, F. Wang, H. and Johnson, D. 2000. Present-day motion of the Sierra Nevada block and some tectonic implications for the Basin and Range province, North American Cordillera. Tectonics, 19 (1), pp.1-24. DOI:10.1029/1998TC001088.
  » https://doi.org/10.1029/1998TC001088
• Duman, H. and Sanli, D.U. 2019. Assessment of geodetic velocities using GPS campaign measurements over long baseline lengths. Natural Hazards and Earth System Sciences, 19 (3), pp.571-582. DOI:10.5194/nhess-19-571-2019.
  » https://doi.org/10.5194/nhess-19-571-2019
• Floyd, M.A. and Herring, T.A. 2020. Fast Statistical Approaches to Geodetic Time Series Analysis. pp.157-183. DOI:10.1007/978-3-030-21718-1_5.
  » https://doi.org/10.1007/978-3-030-21718-1_5
• Giusti, D.A. 1989. Contribuição à geologia ambiental no município de Curitiba-PR. Universidade de São Paulo, São Paulo. DOI:10.11606/D.44.1989.tde-12082015-160804.
  » https://doi.org/10.11606/D.44.1989.tde-12082015-160804
• Hernández-Pajares, M. Juan, J.M. Sanz, J. Orus, R. Garcia-Rigo, A. Feltens, J. Komjathy, A. Schaer, S.C. and Krankowski, A. 2009. The IGS VTEC maps: A reliable source of ionospheric information since 1998. Journal of Geodesy , 83 (3-4), pp.263-275. DOI:10.1007/s00190-008-0266-1.
  » https://doi.org/10.1007/s00190-008-0266-1
• Herring, T.A. Floyd, M.A. King, R.W. and Mc Clusky, S.C. 2015. GLOBK Reference Manual
• Herring, T.A . King, R.W . Floyd, M.A . and Mcclusky, S.C. 2018. GAMIT. Reference Manual Massachusetts Institute of Technology.
• Lagler, K. Schindelegger, M. Böhm, J. Krásná, H. and Nilsson, T. 2013. GPT2: Empirical slant delay model for radio space geodetic techniques. Geophysical Research Letters , 40 (6), pp.1069-1073. DOI:10.1002/grl.50288.
  » https://doi.org/10.1002/grl.50288
• Lambeck, K. 1988. Geophysical Geodesy Oxford University Press.
• Langbein, J. 2004. Noise in two-color electronic distance meter measurements revisited. Journal of Geophysical Research: Solid Earth, 109 (4), pp.1-16. DOI:10.1029/2003JB002819.
  » https://doi.org/10.1029/2003JB002819
• Langbein, J. and Bock, Y. 2004. High-rate real-time GPS network at Parkfield: Utility for detecting fault slip and seismic displacements. Geophysical Research Letters , 31 (15), pp.2-5. DOI:10.1029/2003GL019408.
  » https://doi.org/10.1029/2003GL019408
• Luna, M.P. Staller, A. Toulkeridis, T. and Parra, H. 2017. Methodological approach for the estimation of a new velocity model for continental Ecuador. Open Geosciences, 9 (1), pp.719-734. DOI:10.1515/geo-2017-0054.
  » https://doi.org/10.1515/geo-2017-0054
• Lyard, F. Lefevre, F. Letellier, T. and Francis, O. 2006. Modelling the global ocean tides: Modern insights from FES2004. Ocean Dynamics, 56 (5-6), pp.394-415. DOI:10.1007/s10236-006-0086-x.
  » https://doi.org/10.1007/s10236-006-0086-x
• Mao, A. Harrison, C.G.A. andDixon, T.H . 1999. Noise in GPS coordinate time series. Journal of Geophysical Research: Solid Earth , 104 (B2), pp.2797-2816. DOI:10.1029/1998JB900033.
  » https://doi.org/10.1029/1998JB900033
• Montecino, H.D. de Freitas, S.R.C. Báez, J.C. and Ferreira, V.G. 2017. Effects on Chilean Vertical Reference Frame due to the Maule Earthquake co-seismic and post-seismic effects. Journal of Geodynamics, 112 (July), pp.22-30. DOI:10.1016/j.jog.2017.07.006.
  » https://doi.org/10.1016/j.jog.2017.07.006
• Petit, G. and Luzum, B. 2010. IERS Conventions (2010) Frankfurt am Main.
• Petrie, E.J. King, M.A. Moore, P. and Lavallée, D.A. 2010. Higher-order ionospheric effects on the GPS reference frame and velocities. Journal of Geophysical Research: Solid Earth , 115 (3), pp.1-8. DOI:10.1029/2009JB006677.
  » https://doi.org/10.1029/2009JB006677
• Pilapanta, C. Krueger, C.P. and Criollo, A.R.T. 2018. STOCHASTIC MODEL OF THE BRAZILIAN GPS NETWORK COORDINATES TIME SERIES. Boletim de Ciências Geodésicas, 24 (4), pp.545-563. DOI:10.1590/s1982-21702018000400033.
  » https://doi.org/10.1590/s1982-21702018000400033
• Ponte, R.M. and Ray, R.D. 2002. Atmospheric pressure corrections in geodesy and oceanography: A strategy for handling air tides. Geophysical Research Letters , 29 (24), pp.6-1-6-4. DOI:10.1029/2002gl016340.
  » https://doi.org/10.1029/2002gl016340
• Salamuni, E. Ebert, H.D. da Silva Borges, M. Hasui, Y. Costa, J.B.S. and Salamuni, R. 2003. Tectonics and sedimentation in the Curitiba Basin, south of Brazil. Journal of South American Earth Sciences, 15 (8), pp.901-910. DOI:10.1016/S0895-9811(03)00013-0.
  » https://doi.org/10.1016/S0895-9811(03)00013-0
• Sánchez, L. and Drewes, H. 2020. Geodetic Monitoring of the Variable Surface Deformation in Latin America. In: International Association of Geodesy Symposia DOI:10.1007/1345_2020_91.
  » https://doi.org/10.1007/1345_2020_91
• Sánchez, L. Drewes, H. Brunini, C. Mackern, M. V. and Martínez-Díaz, W. 2015. SIRGAS Core Network Stability. In: Chemistry of Heterocyclic Compounds pp.183-191. DOI:10.1007/1345_2015_143.
  » https://doi.org/10.1007/1345_2015_143
• Serpelloni, E. Faccenna, C. Spada, G. Dong, D. andWilliams, S.D.P . 2013. Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary. Journal of Geophysical Research: Solid Earth , 118 (11), pp.6003-6024. DOI:10.1002/2013JB010102.
  » https://doi.org/10.1002/2013JB010102
• Vernant, P. Nilforoushan, F. Hatzfeld, D. Abbassi, M.R. Vigny, C. Masson, F. Nankali, H. Martinod, J. Ashtiani, A. Bayer, R. Tavakoli, F. and Chéry, J. 2004. Present-day crustal deformation and plate kinematics in the Middle East constrained by GPS measurements in Iran and northern Oman. Geophysical Journal International, 157 (1), pp.381-398. DOI:10.1111/j.1365-246X.2004.02222.x.
  » https://doi.org/10.1111/j.1365-246X.2004.02222.x
• Werlich, R. 2017. Análise de alternativas de processamentos GPS para deteçao de deslocamentos verticais de Terra em áreas urbanizadas. Universidade Federal do Paraná.
• Williams, S.D.P . 2003. Offsets in Global Positioning System time series. Journal of Geophysical Research: Solid Earth , 108 (B6). DOI:10.1029/2002jb002156.
  » https://doi.org/10.1029/2002jb002156
• Williams, S.D.P . 2015. Description of GPS uncertainties within the Long Term Study on Anomalous Time - Dependent Subsidence. National Oceanographic Centre - Liverpool, pp.1-41.
• Zhang, J. Bock, Y. Johnson, H. Fang, P. Williams, S. Genrich, J. Wdowinski, S. and Behr, J. 1997. Southern California permanent GPS geodetic array: Error analysis of daily position estimates and site velocities. Journal of Geophysical Research: Solid Earth , 102 (B8), pp.18035-18055. DOI:10.1029/97JB01380.
  » https://doi.org/10.1029/97JB01380

SUPPLEMENTARY MATERIAL FOR ANALYSIS OF LOCAL SURFACE DISPLACEMENT USING REPEATED GPS MEASUREMENTS: A CASE STUDY OF THE GUABIROTUBA AREA, CURITIBA, BRAZIL

S1. Survey-mode GPS campaigns description

Table S1:
Observation time of 9 survey-mode GPS campaigns made in the study area.

Session number	Observation period	S-GPS network			C-GPS network
		# of stations	# of sessions	session duration	# of stations	# of sessions	session duration
1	Nov. 14 - Nov. 17, 2014	21	1	1 h	4	4	24 h
2	Jan. 28 - Feb. 02, 2015	23	1	1 h	4	6	24 h
3	Apr. 10 - Apr. 14, 2015	23	1	1h	3	5	24 h
4	Jun. 03 - Jun. 08, 2015	23	1	1h	3	6	24 h
5	Apr. 29 - May. 02, 2016	23	1	1h	3	4	24 h
6	Feb. 17 - Feb. 20, 2017	21	1	1h	3	4	24 h
7	Sep. 20 - Sep. 22, 2017	17	1	1h	1	3	24 h
8	Apr. 29 - May. 03, 2019	20	1	4h	3	5	24 h
9	Aug. 26 - Sep. 06, 2019	21	1	4h	3	9	24 h

S2. List of GNSS antenna models used

Table S2:
List of antenna combinations used in the campaigns.

Antenna Type	Calibration		Used at stations
	Type	IGS Model
TRM55971.00	ROBOT	igs14_2114.atx	UFPR
TRM115000.00	ROBOT	igs14_2114.atx	UFPR
LEIAR25	ROBOT	igs14_2114.atx	BCAL
LEIAX1202GG	ROBOT	igs14_2114.atx	HORT, SANE
LEIATX1230GG	ROBOT	igs14_2114.atx	MAR1, MAR2, M001, M002, M003, M004, M005, M006, M010, M016, M017, M018,
LEIATX900CS	ROBOT	Non-existing	MAR1, MAR2, MAR3, M001, M002, M003, M004, M005, M006, M007, M008, M009, M010, M011, M012, M013, M014, M015, M016, M017. M018. M019, M020, M021, M022, M023

S3. LEICA ATX900CS antenna calibration model

To calculate the type-specific calibration model of the LEICA ATX900CS antenna model, a relative field calibration approach was used as established by Krueger et al. (2015). The antenna was mounted on the Dresden Robotic Base (DRB) automatic calibration device (Figure S1.b), developed by the Technical University of Dresden - Germany (Huinca et al., 2015), which, in turn, is mounted on one of the three calibration pillars (Figure S1.c) that form the GNSS Antenna Calibration Basis of the Federal University of Paraná named BCAL/UFPR (Figure S1.a).

Figure S1:
a) BCAL/UFPR Calibration Basis; b) DRB calibration device; c) Calibration pillars.

While the DRB device rotated the GNSS antenna in a north-south and east-west direction, continuous measurements were made for 24 hours and at a sampling rate of 15s. Finally, the processing of data was made using the software Wasoft/Kalib, version 2.0 (Wanninger, 2021) and individual frequency-dependent GPS antenna parameters (i.e., Phase Centre Offsets, PCO and Phase Centre Variations, PCV values) were obtained. More details about the methodology used in this study to obtain the antenna calibration model can be obtained in (Huinca et al., 2015; Krueger et al., 2015).

As demonstrated by Herring et al. (2016), the elevation dependence of the phase residuals (averaged over all azimuths), provides a useful measure of the phase modelling (i.e., reliability of antenna calibration model) and multipath effects at each site (Herring et al., 2018). Thus, Figure S2. shows the mean ionosphere-free phase residuals (i.e., mean between 9 campaigns) obtained for each S-GPS site with the LEICA ATX900CS antenna model (Table S2) together with their total mean (i.e., mean of 26 S-GPS phase residuals).

Figure S2:
Mean LC phase residuals obtained for all C-GPS site with LEICA ATX900CS antenna.

For all sites, it is possible to see that the pattern of their phase residuals remains relatively constant over time, with the higher oscillating phase residuals (up to 10 mm) at low elevations, probably, as a result of multipath effects, but generally without influence on station displacement estimates.

As demonstrated by Herring et al. (2016), similar patterns to those obtained in the present study (i.e., individual residuals mean close to zero) are indicative of a good model for the antenna phase centre variations, which proves the consistency of the calibration model and validates its usability.

S4. Coordinate time series of GPS sites

Figure S3:
Coordinate time series of M004, M015, M016, and M017 S-GPS sites.

Figure S4:
Coordinate time series of M005, M011, M014, and M018 S-GPS sites.

S5. Two sample T-Test

Velocities analysis of SIRGAS-CON sites

Table S3:
Main statistics of Two Sample T-Test.

Parameter	East		North		Up
	Estimated Solution	SIRGAS Solution	Estimated Solution	SIRGAS Solution	Estimated Solution	SIRGAS Solution
Mean (mm)	-2.98	-3.03	12.76	12.55	-0.78	-0.77
Variance (mm)	0.19	0.19	0.48	0.56	1.68	2.56
Number of observations	11	11	11	11	11	11
t-test parameters
Degree of Freedom	20.00		20.00		19.00
T Stat	0.27		0.69		-0.02
P(T <= t) one-tail	0.39		0.25		0.49
t Critical one-tail	1.72		1.72		1.73
P(T <= t) two-tail	0.79		0.50		0.99
t Crítical two-tail	2.09		2.09		2.09

S6. Absolute velocities of sites

Table S4:
Summary velocity estimates.

Station	East		North		Up
	rate	Sigma	rate	Sigma	rate	sigma
	(mm/a)	(mm/a)	(mm/a)	(mm/a)	(mm/a)	(mm/a)
C-GPS Network
BCAL	-2.56	0.45	13.94	0.42	3.01	1.85
HORT	0.32	2.44	31.37	2.28	-56.12	10.08
SANE	-1.44	0.22	11.84	0.21	2.17	0.88
UFPR	-2.52	0.25	13.92	0.24	-1.11	1.05
S-GPS Network
M001	-5.05	1.14	16.38	1.08	-7.45	6.61
M002	-2.03	1.32	13.06	1.40	-0.60	7.43
M003	-2.86	1.34	10.52	1.36	-0.16	6.95
M004	1.36	1.47	16.77	1.33	-31.14	10.95
M005	-1.80	1.53	12.98	1.49	-1.19	7.54
M006	-0.66	2.67	11.77	3.01	-10.59	17.27
M007	-4.92	1.16	16.56	1.19	-2.43	5.94
M008	-3.73	1.10	12.32	1.06	-3.04	5.99
M009	-3.63	0.91	12.10	0.97	-2.26	4.18
M010	-3.98	0.90	9.41	0.87	-2.70	4.47
S-GPS Network
M011	-2.68	0.86	13.93	0.85	-1.52	4.20
M012	-13.58	4.19	18.13	3.94	-2.98	22.32
M013	-16.02	0.88	8.69	0.93	-4.49	4.36
M014	-13.67	1.45	11.29	1.52	-30.07	9.60
M015	-3.12	1.37	8.51	1.35	-10.44	7.33
M016	-13.67	0.88	6.88	0.85	-16.24	4.16
M017	-1.94	0.93	11.43	0.97	0.44	4.49
M018	-53.15	3.51	-38.43	3.25	-7.82	19.32
M019	-4.27	0.98	13.17	0.99	-0.72	4.81
M020	-9.77	18.30	16.77	17.41	11.27	98.27
M021	-3.60	1.59	17.20	1.71	6.10	8.23
M022	-1.49	2.62	9.55	2.66	-1.57	13.36
M023	-1.90	20.57	25.77	21.67	-10.16	105.31
MAR1	3.35	3.68	24.44	3.60	-10.85	20.89
MAR3	-0.84	2.81	12.19	2.99	22.28	15.45

Publication Dates

• Publication in this collection
  01 Apr 2022
• Date of issue
  2022

History

• Received
  31 Mar 2021
• Accepted
  16 Feb 2022