GPS Time Series Analysis from Aboa the Finnish Antarctic Research Station

Abstract

Continuous Global Positioning System (GPS) observations have been logged at the Finnish Antarctic research station (Aboa) since February 2003. The station is located in Dronning Maud Land, East Antarctica. Almost 5000 daily observation files have been archived based on yearly scientific expeditions. These files have not been fully analysed until now. This study reports for the first time on the consistent and homogeneous data processing and analysis of the 15-year long time series. Daily coordinates are obtained using Precise Point Positioning (PPP) processing based on two approaches. The first approach is based on the Kalman filter and uses the RTKLIB open source library to produce daily solutions by unconventionally running the filter in the forward and backward direction. The second approach uses APPS web service and is based on GIPSY scientific processing engine. The two approaches show an excellent agreement with less than 3 mm rms error horizontally and 6 mm rms error vertically. The derived position time series is analysed in terms of trend, periodicity and noise characteristics. The noise of the time series was found to be power-law noise model with spectral index closer to flicker noise. In addition, several periodic signals were found at 5, 14, 183 and 362 days. Furthermore, most of the horizontal movement was found to be in the North direction at a rate of 11.23 ± 0.09 mm/y, whereas the rate in the East direction was estimated to be 1.46 ± 0.05 mm/y. Lastly, the 15-year long time series revealed a movement upwards at a rate of 0.79 ± 0.35 mm/y. Despite being an unattended station, Aboa provides one of the most continuous and longest GPS time series in Antarctica. Therefore, we believe that this research increases the awareness of local geophysical phenomena in a less reported area of the Antarctic continent.

1. Introduction

Antarctica is an important area due to its ongoing glaciation and contemporary climate-driven rapid changes. Geodetic and geophysical observations on Holocene changes in Antarctica will give us unique information when studying global climate-change related phenomena. Long time series of geodetic and geophysical stations provide data explaining parts of the interaction between solid earth, cryosphere, hydrosphere and atmosphere. They are also essential as ground truth when data from satellite missions are analysed.

Finland joined the Antarctic Treaty in 1984, became consultative party in 1989 and full member of the Scientific Committee on Antarctic Research (SCAR) in 1990 [1,2]. Five years later after signing the treaty, the Finnish research station, Aboa, was established in Dronning Maud Land (DML), East Antarctica. All Aboa activities were gathered under the general heading of the Finnish Antarctic Research Program (FINNARP). The first FINNARP expedition took place during the Antarctic summer of 89/90. Since then, a total number of 26 FINNARP expeditions have been arranged during the Antarctic summers (http://www.antarctica.fi/). There were no expeditions during the summers of 90/91, 94/95, 96/97 and 98/99. The research activities associated with the expeditions have been mainly financed by the Academy of Finland via regular calls for Antarctic research issued every fourth year. In addition, universities and research institutions invest financial resources of their own for Antarctic research [2].

In recent years, the Finnish Antarctic research has focused on geodesy [3,4,5], glaciology [6,7,8], geology and geophysics [9,10,11,12,13], air quality and climate change [14,15,16,17,18], meteorological and space physics [19,20,21,22], marine and structural technology [23,24], oceanography and marine biology [25,26]. Typically, the research work has been undertaken at Aboa and surrounding areas. These research activities have sometimes been complemented with trips and overnight stays to more distant locations, including the coast. In addition, foreign researchers have participated in the FINNARP expeditions and Finnish researchers have worked at foreign research stations and on board of foreign vessels. For example, Finnish researchers have performed ozone and ultraviolet (UV) research at the Argentine Marambio research station, cosmic rays research at the Italian-French Dome Concordia station, or marine and sea-ice research on South African R/V Agulhas II. Absolute and relative gravity measurements have been done at the Norwegian Troll station, Russian Novolazarevskaya station, Indian Maitri station, and South African SANAE IV station. SCAR GPS measuring campaigns have been conducted at several locations, including Swedish Wasa station neighbouring Aboa. Furthermore, logistic and research support has been provided to international collaborators, e.g., satellite controlled meteorological balloons on behalf of the Smith College, USA (http://www.antarctica.fi/).

Currently, the Aboa research station is equipped with several automated year-round measuring instruments, such as continuous GPS tracking equipment (2003), an automatic weather station (2004), and an automatic seismometre (2007). In addition, Finnish researchers have deployed other automated year-round measuring instruments outside Aboa, such as automated weather stations at the top of the Basen nunatak and near the coast, automated buoy at sea to measure the vertical temperature profile of sea ice, and cosmic ray detectors. All these measuring instruments ensure the collection of significant sample series and the continuation of long-term time series data in accordance with Finland’s Antarctic research strategy [2].

Numerous GPS and/or GNSS (Global Navigation Satellite System) continuous and episodic stations have been established in Antarctica in the last two decades (Figure 1). The development has lead to GPS becoming an essential component of geophysical and meteorological research for studying fundamental solid Earth and atmospheric processes that drive natural hazards, weather, and climate [27]. The first GPS measurements in Antarctica were conducted in the early 90s. The main objective was to gain detailed insights into the tectonic behaviour of the Antarctic plate [28,29]. The accumulation of the measurements allowed the definition of a precise reference frame [30]. Later, Thomas et al. [5] analysed GPS data from more than 50 stations across Antarctica. The authors estimated a more complete and precise vertical velocity field that allowed them to observe increase in uplift rates in the northern Antarctic peninsula after 2002–2003. Plate tectonics and vertical deformation studies also triggered the development of local or regional-scale projects [31,32]. Long time series of continuous GPS data have been used to constrain and test the reliability of glacial isostatic adjustment (GIA) models [33,34,35]. Furthermore, GPS measurements in Antarctica have allowed atmospheric studies, such as ionospheric scintillation [36,37] or ground-based water vapour retrieval [38,39].

Aboa is an important GPS station in Antarctica due to the sparsity of bedrock sites, especially in East Antarctica. As seen in Figure 1, most of the sites are located in West Antarctica, whereas the number of sites is much lower in East Antarctica. In addition, the distances between the sites are rather large and limited to coastal regions of the continent where bedrock can be found. Furthermore, only a limited number of sites provide continuous GPS observations. As of 26 October 2018 (http://www.igs.org/network), there are only 13 active permanent International GNSS Service (IGS) stations [40] on the Antarctic continent and only 10 of them have longer continuous time series than Aboa. This study covers a time span of 15 years. This is two times longer than in [5] who used data for 1959 days between 2003 and 2010, whereas [42,43] used reduced time periods. [42] covered a period of 1415 days between 2009 and 2013, while [43] analysed a time series derived for the 2010–2013 period. Furthermore, [33] computed a series of theoretical expected site velocity uncertainty for 5, 10, and 20-year long time series. Their computations show that the uncertainty decreases with the length of time series. Therefore, our results have the benefit of longer time series and thus higher reliability. Our study is the first to include an in-depth analysis of all available Aboa GPS data.

The rest of the paper is as follows. Section 2 describes the material and methods used to generate and analyse the coordinates time series. Section 3 reports the numerical results of our study. Section 4 discusses the obtained results and their relationship with other published data. Section 5 summarises the main achievements, the importance of our study and the way how our work will continue.

2. Materials and Methods

2.1. Site Description and Equipment

The Finnish Antarctic research station (Aboa) was established in 1989 [8]. The station is located in the western part of the Dronning Maud Land, approx. 130 km from the coast, on Basen nunatak in the northern part of the Kraul mountain range at 73.0437 S, 13.4071 W, and 485 m elevation above sea level. Aboa is occupied only during the Antarctic summers. The GPS permanent station was setup during the 2002/2003 Antarctic expedition. It became operational on 31 January 2003 at 10:10 UT [3].

The Aboa station tracks only GPS satellite signals. The station is equipped with an Ashtech choke ring antenna (ASH701945CM) covered with an Ashtech conical SNOW antenna dome. The antenna has remained unchanged since the beginning of the operation. It has been installed on bedrock on the top of a 1.5 m tall steel mast (Figure 2). Installation on the top of the mast is necessary since the snow environment changes constantly at the site. The height of snow has been close to the antenna height during several seasons, whereas during other seasons the bedrock below the mast appears visible. The GPS receiver is a Javad EURO80 GDA card (JPS HE GDA) custom-fitted into a metal enclosure. The choice of the receiver was made based on the low power usage (1.8–2.4 W) of the Javad receiver since the battery power is limited at Aboa during the winter months. For the same reason, the data is obtained only once per year during the FINNARP expeditions to Aboa. Actually, two identical receivers are used. The FINNARP expedition changes receivers at Aboa and returns the one from the previous year to FGI for data download and maintenance. This strategy allows for the FINNARP expedition to always have spare receiver in case of instrument failure during the previous winter season. Further technical information is given in [3].

2.2. Observation Data

Daily GPS observation files have been logged continuously at Aboa since 1 February 2003. The files have been retrieved on a yearly basis following the annual Antarctic expedition. The last observation day in the current archive is 8 December 2017. As a result, the covered period accounts for almost 15 years. The on-site observations were logged into the receiver manufacturer format and translated to RINEX (Receiver INdependent EXchange) format, version 2.11 using teqc utility [44] upon retrieval. The files include GPS-only observations at the epoch rate of 30 s. Epoch wise, the number of tracked satellites varied between 4 and 15, with a median number of 10 satellites per each epoch of observation.

Figure 3 gives the yearly distribution of the existing observation files. The overall daily availability is 91% (i.e., 4924 available days out of 5425 possible days). The largest gap of 381 days occurred after the 2006/2007 expedition. All files in 2006 and the first 17 days of 2007 were lost due to an unknown reason. Moreover, two other medium-length gaps were identified. One of the gaps occurred in 2008 and counted for 52 days between day-of-year 174 and 226. This was identified with the loss of power during the Antarctic winter. The other gap happened in 2010 and counted for 23 days between day-of-year 013 and 036. Furthermore, 13 days are missing in July 2004 as well as the last 8 days of December 2004. The other identified gaps count for less than 4 days. All together there are 23 gaps in the daily time series. Despite being an unattended GPS station, there is 100% daily availability for five years (2009, 2013, 2014, 2015 and 2016).

2.3. Data Processing

GPS observation data for high-precision applications are typically processed based on two basic approaches: relative positioning and point positioning.

Relative positioning allows estimating the position of a new point with respect to another known, reference point [45,46]. Baseline between a pair of points is observed and processed by forming double-differenced (DD) code and phase observations and solving for the integer phase ambiguities. DD approach (i.e., difference between-receivers and between-satellites) removes the observation errors common to both ends of the baseline, such as satellite orbital errors, satellite clock biases and receiver clock biases. In addition, errors due to troposphere and ionosphere reflection are eliminated or mitigated due to the correlation with the distance between the two points. The longer the distance, the higher the residual errors. Therefore, relative positioning may be extended to a local or regional network. However, the approach becomes computationally cumbersome if the number of stations within such a network increases significantly [27].

Precise Point Positioning (PPP) [47] uses un-differenced, dual-frequency, code and phase observations from a single user receiver along with precise satellite orbit and clock products. The dual-frequency observations are used to form ionosphere-free linear combinations in order to remove the first order ionospheric effects. Unlike in DD relative positioning, PPP has no baseline constraint and thus does not require simultaneous observations at two points. Although PPP itself uses data from a single station only, computation of the satellite orbits and clocks needed for its implementation require the use of a global tracking network [48]. Moreover, PPP requires more careful handling of all error sources affecting the satellite observations than DD relative positioning [49,50]. Furthermore, the physical models (e.g., Earth tides, antenna phase centre corrections) used in the PPP processing should be consistent with those used in the generation of the precise products. The parameters estimated with PPP processing are the 3D point position, the receiver clock, the zenith tropospheric delay and the non-integer ambiguity for each satellite. PPP positions are computed directly in the terrestrial frame in which the satellite orbits and clock parameters are given [27]. The current frame realisation is International Terrestrial Reference Frame (ITRF) 2014 [51].

A variety of tools for precise GPS/GNSS observation data processing have been developed at academic, research, and commercial institutions. They may be grouped into scientific, commercial and open-source software packages. In addition, numerous free online and commercial services exist.

Table 1. Precise GPS observation data processing software packages and services.

Tool	Type	Examples
Software	scientific	Bernese, GIPSY-OASIS II, GAMIT/GLOBK/TRACK, NAPEOS, PANDA
package	commercial	Trimble Business Centre, GrafNav, etc.
	open-source	RTKLIB , GPSTk, gLAB
Service	free on-line	AUSPOS, OPUS, SCOUT, APPS, CSRS-PPP, GAPS, magicPPP, CenterPoint RTX
	commercial	Atlas (Hemisphere), CORRECT-PPP (NOVATEL), OmniSTAR, RTX (Trimble), StarFire (Fugro), Ultra (Veripos), etc.

summarises the most common GPS data processing packages and services.

In this work, we chose the PPP approach instead of the DD to estimate station velocities independently of the other stations’ data. The DD solution would be strongly dependent of the reference velocities of the sparsely located ITRF stations and also affected by the data quality of the reference stations. Our PPP daily positions at Aboa were estimated using the following two tools: RTKLIB open-source package (http://www.rtklib.com/) and Automatic Precise Positioning Service (APPS) web-service (http://apps.gdgps.net), respectively. The former is a versatile tool that enables scientific processing of the data via numerous parameters, some with default values. Nevertheless, the tool requires a high level of understanding of the processing engine to exploit its potential and flexibility. On the other hand, APPS is an easy to use web service based on the GIPSY processing engine. The user uploads the file to the cloud server and receives back a summary processing report in text format. The service allows only a few options to be changed by the user. APPS service use was intended for validation purposes. We describe below how the observation data were processed using the above tools.

2.3.1. PPP Processing Settings (RTKLIB)

RTKLIB open-source package consists of a portable program library and several application libraries designed for GNSS standard and precise positioning. The package is distributed under the open source initiative. The current release is version 2.4.2, whereas the development version is 2.4.3b29. We used the latter one. In addition, RTKLIB supports various positioning modes: single, differential, baseline, static, kinematic, real-time, or post-processing. We exploited the PPP mode in this work. Furthermore, RTKLIB employs un-differenced measurement equations and Extended Kalman Filter (EKF) estimation in order to obtain the final solutions. The mathematical formulation is explained in detail in the Appendix E.8 of the package manual.

As mentioned, PPP requires precise products. In this work, PPP solutions were produced using the final satellite orbits, final clock product with a sampling of 5 s, and weekly final Earth rotation parameter products from Centre for Orbit Determination in Europe (CODE). The products were taken from the homogeneous CODE reprocessing, release 2017 [52]. They were available until the end of 2014. Afterwards, we used the operational products [53]. Further corrections were derived using several mathematical models and standard conventions. The elevation cutoff angle was set to 10 degrees. This is a best compromise between the quantity and quality of data for high-latitude regions [39]. An elevation-dependent weighting scheme was adopted with separate constant and elevation-dependent terms. The estimated parameters were station coordinates, receiver clock, troposphere parameters and float ambiguities. The estimated coordinates were obtained in ITRF2008 [54] before 29 January 2017 and ITRF2014 afterwards.

Table 2. Methods of correction for various observation errors in our PPP data processing strategy.

Source	Error	Method of Correction
satellite	orbit	CODE ${}^{\mathrm{a}}$ precise orbit products in standard product 3 (SP3) format with sampling rate of 15 min [52,53]
	clock offset	CODE precise clock product in extended RINEX clock format (CLK) with sampling rate of 5 s [52,53]
	phase center offset and variations	IGS ANTEX ${}^{\mathrm{b}}$ maps [55]
atmosphere	Ionophere	Ionosphere-free linear combination (1st order term only). Higher order corrections were not applied.
	Troposphere	Saastamoinen model [56] for a-priori dry tropospheric delay using standard atmosphere. The wet tropospheric delay together with horizontal gradients were estimated. Global Mapping Function (GMF) [57] to estimate the delays along the line of sight.
receiver	phase center offset and variations	IGS ANTEX maps [55]
	differential code biases	not applied as the receiver reported P1 observable
site	Solid Earth tide	IERS2010 conventions [58]
displacements	Pole tide	CODE weekly ERP files [52,53]; IERS2010 conventions [58]
	Ocean Loading tide	FES2004 model ${}^{\mathrm{c}}$ [59]
others	Phase windup	wavelength dependent [60]
	Phase ambiguity	estimated as float number

${}^{\mathrm{a}}$ftp://ftp.aiub.unibe.ch, ${}^{\mathrm{b}}$ftp://igs.org/pub/station/general, ${}^{\mathrm{c}}$http://holt.oso.chalmers.se/loading.

gives an overview of the various error components and corrections used in our PPP processing.

2.3.2. PPP Combined Daily Solutions (RTKLIB)

EKF is a recursive estimator that does not require all previous data to be saved and reprocessed when a new measurement is taken [61]. Thus, the estimated solution converges and its quality improves over time. PPP convergence time and positioning performance depend on several factors, such as the sampling interval of the precise products, the length of the observation session, the observation epoch rate, the observation quality, the receiver type, the multipath level, the geographical location, the number of tracked satellites and constellations, the positioning mode, and the ambiguity resolution [50,62,63]. For example, in static PPP, 20 min are needed for 95% of solutions to reach a horizontal accuracy of 20 cm [62]. Nevertheless, sub-centimetre accuracy levels are achievable for the 24-h data set [50,62]. Figure 4 depicts the convergence process for the forward and backward solutions at Aboa on 1 January 2017.

Typically, the filter is run in the forward direction and the solutions over the convergence period must be discarded. To cope with this, we used an unconventional way to run the filter. We run the filter twice, i.e. in both the forward and backward directions. This is possible only in post-processing mode. Next, the forward and backward solutions are joined into one combined solution using the inverse variance weighting approach. The following equation is used:

$$X_{cmb}=\frac{X_{fwd}·Q_{fwd}^{-1}+X_{bwd}·Q_{bwd}^{-1}}{Q_{fwd}^{-1}+Q_{bwd}^{-1}}$$

where the X and Q denote the position vector and its associated variance co-variance matrix from the forward ($fwd$), backward ($bwd$) and combined ($cmb$) solution, respectively.

As a result, the limitations coming from the convergence period are removed. The combined solution has a consistent standard deviation at each observation epoch. Figure 5 depicts the combined solution at Aboa on 1 January 2017.

Finally, we chose the median coordinates of the daily combined solution as the daily representative to form the coordinate time series.

2.3.3. Validation Solution Processing Settings (APPS, GIPSY-OASIS)

All available observation data at Aboa were also processed using the APPS web-service. The main intent was for validation purposes. APPS solution is based on the GIPSY-OASIS 6.4 software and uses Jet Propulsion Laboratory (JPL) final orbit and clock precise products. The web-interface allowed us to set two parameters: the elevation cut off angle was set to 10°, whereas the observation weighting scheme was chosen $srtq\left(sin\right)$. The station positions were estimated at 5 min intervals. All other processing parameters and correction models were defined automatically by the service and inaccessible to the user. According to the information on the service webpage, several signal and ground models were also applied by GIPSY, such as GPS yaw attitude, phase windup, GMT troposphere mapping function, a-priori hydrostatic delay as function of the station height above the ellipsoid, a-priori constant value for the wet delay of 10 cm, tropospheric gradients, relativity, pole and solid tides, ocean tides, and second order ionospheric delay. APPS produces several GIPSY and generic output files, including a summary of the PPP processing. APPS provides coordinates in ITRF2014 starting 5 June 2018.

2.4. Time-Series Analysis

PPP processing provides positions estimates (X, Y, Z) in an Earth-fixed, Earth-centred terrestrial reference frame in which the precise orbit are provided, i.e. ITRF2008 and/or ITRF2014. The Cartesian coordinates are then transformed into topocentric (or levelled) coordinates for a more intuitive and meaningful interpretation. Therefore, the coordinate time series are provided in terms of local displacements in the North, East and Up (N, E, U) directions with respect to the position at an initial epoch. We define the position of the topocentric origin in the geometrical centre of the time series.

Typically, time series analysis includes outlier detection, offsets and displacements identification, noise characterisation, trend and seasonal estimation, and residual analysis. The most popular analysis tools are GGMatlab (TSView) [64], FODITS [65], CATS [66], and Hector [67]. TSView is written in Matlab and complements the GAMIT/GLOBK package, FODITS is embedded into the Bernese GNSS software [68], whereas CATS and Hector are C/C++ written independent command line routines. Hector (http://segal.ubi.pt/hector/) is able to deal with missing data without the need for interpolation to fill in the gaps in the time series. Additionally, the parameter estimation is based on the Maximum Likelihood Estimation (MLE) method. Furthermore, Hector supports all the analyses mentioned earlier. Therefore, it was our choice for the coordinate time series analysis. Version 1.6 was used in this study.

Figure 6 depicts the entire data collection, processing and analysis work flow diagram.

2.5. Environmental Loading Factors

The loading of the crust due to environmental factors, such as varying atmospheric and hydrological conditions as well as non-tidal ocean variations are well known, but not implemented in the GNSS computation. The loading factors have generally amplitudes of a few millimetres up to a few centimetres in the Antarctic continent, but as the phenomena show a clear annual variability it may distort especially campaign measurements as well as long term trends.

We compared pre-computed loading time series from three service providers as shown in

Table 3. Environmental loading services.

Service	Reference	Link
International Mass Loading Service	[69]	http://massloading.net/
EOST loading service	[70]	http://loading.u-strasbg.fr/index.php
GFZ service	[71]	https://isdc.gfz-potsdam.de/esmdata/loading/

. These services provide pre-computed time series (EOST and IML) and gridded data (GFZ) of hydrological loading, non-tidal atmospheric loading and non-tidal ocean loading. All services provide data in both centre-of-the-mass (CM) and centre-of-the-solid-earth (CE/CF) frame. All three services have similarities and differences. They all provide data starting from approximately 1980 and update their products continuously. There are altogether 20 different products available, but only two non-tidal ocean and two atmospheric loading products cover the whole time period (2003–2017). We decided not to use the hydrological loading computed from hydrologic models as the models have severe problems on ice covered regions and cannot reliably present the annual variability.

3. Results

3.1. PPP-Derived Coordinate Time Series

This section reports on the quality of the PPP-derived coordinate time series at Aboa in terms of internal precision, accuracy and postfit observation residuals.

Figure 7 illustrates the topocentric coordinate time series at Aboa from the RTKLIB and GIPSY processing, respectively. The origin of the topocentric system is set to the mean coordinates over the entire period. Whereas the GIPSY solution shows no noticeable outlier, the RTKLIB solution exhibits a few obvious outliers. In addition, the RTKLIB solution appears somehow noisier. Nevertheless, both time series appear to follow a similar trend.

Figure 8 depicts the standard deviation of the estimated position as given by the processing engine. The median standard deviation along all the directions is below 5 mm. In the horizontal direction, the median precision is around 2 mm for RTKLIB and less than 1 mm for the GIPSY solution, respectively. In the vertical direction, both solutions exhibit about 4 mm median standard deviation. The GIPSY solution appears a little more robust for the entire interval period, whereas the RTKLIB solution exhibits larger standard deviations in the beginning of the time series. A significant improvement is visible in 2004. However, this may not be attributed to the processing engine but perhaps to the data quality or the difference in the precise products used by the two solutions. One may suspect a weaker performance of the products, perhaps in the Antarctic region, before April 2004. After that, the standard deviations appear rather consistent despite episodic mild jumps for a few days. Although these small differences, we can conclude that both time series are robust.

Figure 9 illustrates the accuracy of the RTKLIB solution with respect to the GIPSY solution as error distribution histograms of the local offsets. The accuracy validation requires that both series are expressed in the same frame. Thus, the RTKLIB coordinates were first transformed from ITRF2008 to ITRF2014. Then, the differences in Cartesian coordinates were calculated. These differences were converted to local NEU offsets for a more meaningful interpretation. Next, the outliers and trend were removed using a reiterative linear regression algorithm combined with the IQR (inter quartile range) robust statistic and a threshold of 2.2 for outlier removal. This threshold is roughly equivalent to the 3-sigma and it works well in our experiences. This process was carried out individually for each axis. The detrended offsets show peak-to-peak variations of less than ±8 mm for the North and East directions. The variations on the Up direction are twice as large and account for ±18 mm. Overall, 2.6 mm rms (root mean square) is obtained along the horizontal directions, whereas 6.0 mm rms is computed for the vertical direction.

Table 4. Conventional and robust statistics to assess the accuracy of RTKLIB solution with respect to GIPSY solution [unit: mm].

Conventional	North	East	Up	Robust Statistics	North	East	Up
Count	4886	4858	4886
Min	−7.9	−7.8	−17.9	Q1 (25%)	−1.8	−1.8	−4.1
Max	7.8	7.8	17.8	Q3 (75%)	1.8	1.8	4.1
Spread	15.8	15.7	35.6	IQR	3.6	3.6	8.2
Avg	0.0	0.0	0.0	MED	0.0	0.0	0.0
Std	2.6	2.6	6.0	MAD	1.8	1.8	4.1
Rms	2.6	2.6	6.0

summarises the conventional and robust statistics. These statistics confirm an excellent agreement between the two solutions.

Figure 10 shows the code and phase observation post-fit observation residuals as obtained from GIPSY processing. The code residuals indicate a dual behaviour over time. This comes from the usage of two receivers mentioned in Section 2.1. The receivers have been switched during the yearly FINNARP expedition, usually at year end/start.

3.2. Time Series Analysis

It is generally accepted that the GPS coordinate time series contain some form of temporally correlated noise along the expected white noise [33]. The effect of this noise has direct effect on the uncertainty of the estimated parameters. We used Hector [67] software that enables modelling the temporal correlation of the data and thus providing more realistic uncertainty estimates.

First, daily Cartesian coordinate time series were converted to local NEU coordinates with respect to the mean position, and the time series of each component was analysed separately. We found only one discontinuity in the time series analysis. It was linked to the switch from IGb08 to IGS14 in CODE’s operational solutions, which were used for the RTKLIB solution. Therefore, we solved an offset in the case of RTKLIB time series at 29 January 2017. The frame of the JPL products (used in the GIPSY solution) was the same for the entire period of the time series. Second, we removed the outliers using an IQR factor of 2.2. In addition, the offsets as well as annual and semiannual effects were taken into account in the outlier rejection. About 50 epochs (i.e., 1%) were rejected from the GIPSY solution and 100 epochs (i.e., 2%) from the RTKLIB solution, respectively. Third, we estimated the trends using three different noise models: Generalised Gauss-Markov (GGM), power-law (PL) and flicker noise implementation of the Hector (FlickerGGM). Typically, the white noise is added to the model, but in our time series the fraction for white noise was zero and, therefore, it was omitted from the model. The choice of the final model was based on the analysis of the power spectral analysis (PSD) and the Akaike and Bayesian information criteria (AIC/BIC). The model with lowest AIC or BIC values should fit best the time series.

Figure 11 depicts the PSD plots, whereas

Table 5. Estimated trend and seasonality parameters for different noise models for the GIPSY solution.

Component	Noise Model	Velocity (mm/y)	Annual (mm)	Semiann. (mm)	AIC	BIC
North	FlickerGGM	11.20 ± 0.03	0.23 ± 0.12	0.26 ± 0.12	17,712.2	17,757.7
	GGM	11.20 ± 0.03	0.19 ± 0.10	0.24 ± 0.11	17,694.5	17,752.9
	Powerlaw	11.23 ± 0.09	0.22 ± 0.12	0.24 ± 0.11	17,704.5	17,756.4
East	FlickerGGM	1.46 ± 0.05	0.38 ± 0.15	0.22 ± 0.10	15,048.7	15,094.2
	GGM	1.46 ± 0.03	0.36 ± 0.10	0.18 ± 0.08	14,970.1	15,028.6
	Powerlaw	1.46 ± 0.05	0.36 ± 0.10	0.18 ± 0.07	14972.6	15024.5
Up	FlickerGGM	0.71 ± 0.04	2.31 ± 0.23	2.07 ± 0.23	28,733.6	28,779.0
	GGM	0.71 ± 0.05	2.31 ± 0.28	2.07 ± 0.28	28,674.9	28,733.3
	Powerlaw	0.79 ± 0.35	2.43 ± 0.57	2.14 ± 0.42	28,736.1	28,788.0

and

Table 6. Estimated trend and seasonality parameters for different noise models for the RTKLIB solution.

Component	Noise Model	Velocity (mm/y)	Annual (mm)	Semiann. (mm)	AIC	BIC
North	FlickerGGM	11.10 ± 0.05	1.44 ± 0.23	0.37 ± 0.16	18,598.6	18650.4
	GGM	11.09 ± 0.04	1.43 ± 0.16	0.35 ± 0.12	18,553.9	18,618.7
	Powerlaw	11.14 ± 0.09	1.44 ± 0.17	0.33 ± 0.12	18,558.4	18,616.7
East	FlickerGGM	1.49 ± 0.06	0.85 ± 0.26	0.29 ± 0.15	19,952.4	20,004.2
	GGM	1.48 ± 0.04	0.83 ± 0.17	0.22 ± 0.11	19,850.1	19,914.9
	Powerlaw	1.50 ± 0.09	0.82 ± 0.17	0.22 ± 0.11	19,855.0	19,913.3
Up	FlickerGGM	0.38 ± 0.11	3.11 ± 0.49	1.12 ± 0.37	26,283.4	26,335.2
	GGM	0.39 ± 0.11	3.10 ± 0.48	1.11 ± 0.36	26,287.0	26,351.9
	Powerlaw	0.45 ± 0.34	3.11 ± 0.50	1.13 ± 0.35	26,297.3	26,355.6

summarise the AIC/BIC values. The results suggest that the GGM and PL models fit mostly equally best to the estimated spectrum. For the GIPSY Up component, the AIC/BIC analysis supports more GGM model, but the PSD plot shows that none of the model fits very well to the data. There are a few power peaks at roughly 27 cpy and 60–70 cpy corresponding roughly 14 and 5 day cycle, respectively. Those peaks probably cause the misfit of the models to the estimated spectrum.

Aboa coordinate time series contain temporally correlated coloured noise which spectral index is closer to flicker noise. The noise of the horizontal components is characterised by an average spectral index $\kappa =-0.8$ and average amplitude of 5 mm/y${}^{0.20}$, whereas the noise of the vertical component has similar nature $\kappa =-0.9$ but three times larger amplitude, i.e. 16 mm/y${}^{0.22}$.

Figure 12 depicts the coordinate residuals after removing the outliers, trend, annual and semiannual periodic signals. The rms error is reduced to less than 2 mm for the horizontal components and below 6 mm for the vertical component.

Figure 13 displays the horizontal movement. Since the beginning of operation on 1 February 2003, Aboa station has moved about 175.5 mm in the North direction at a rate of 11.23 ± 0.09 mm/y and about 29.5 mm in the East direction at a rate of 1.46 ± 0.05 mm/y (see Table 5 and Table 6).

Figure 14 illustrates the movement in the vertical direction and its signal decomposition. The vertical movement consists of a linear trend and seasonal signals. The trend was found to be upwards at a rate of 0.79 ± 0.35 mm/year. The seasonal signals include annual and semiannual components. The amplitude of the annual signal was estimated to be 2.43 ± 0.57 mm, whereas the amplitude of the semiannual signal was estimated to be 2.14 ± 0.42 mm. The maximum amplitude is reached on 23/24 August while the minimum in 13 December. However, these seasonal signals are modelled under the assumption of constant yearly amplitude. Whereas such assumption allows to determine a significant part of the variation, it does not capture the inter-annual variations.

3.3. Environmental Loading

The four time series that were used are depicted in Figure 15. The black curves are computed by the EOST. Their atmospheric loading is computed using the ECMWF operational model with 3 h resolution, averaged to daily values. The non-tidal loading is from the ECCO2 ocean bottom pressure model [72], which provides daily values. The red curves were interpolated from the GFZ loading grid using bilinear interpolation. The atmospheric loading is computed from the same 3 h ECMWF operational model, but for the non-tidal ocean loading the MPIOM model [73] is used with 3 h resolution. The atmospheric loading shows a variability of few millimetres in the horizontal components and a few centimetres in the vertical. The non-tidal ocean shows a smaller variation both in vertical and horizontal.

We studied the effect of the environmental loading on the time series rms and trends. We did the time series analysis before and after removing the atmospheric and non-tidal ocean loading. The effect was unremarkable, as removing the loading affected the trend below the 0.01 mm level and increased the residual rms in all three components. For example, using the EOST products, the trend of the up-component decreased from 0.79 to 0.74 mm/year and the uncertainty for the trend increased from 0.35 to 0.48 mm/year. The residual rms increased from 5.7 to 6.1 mm.

4. Discussion

4.1. What Do Aboa Results Tell to Us?

Two PPP-based solutions were produced, compared and validated. These solutions agree very well despite several differences in the processing engine and precise products in use. The achieved level of agreement (2.6 and 6.0 mm in rms error in horizontal and vertical, respectively) brings confidence in the quality of the PPP-derived coordinate time series and its analysis findings.

Although the receivers are identical, they have different firmware versions. For code residuals, median rms error accounts for 0.54 m in case of firmware version 2.3, whereas 0.48 m accounts for the firmware version 2.5. It is clear that the version 2.5 brings an improvement in the code tracking of a 10%. Moreover, the receiver with firmware version 2.3 exhibits a mild upwards trend in the code residual. This might be a sign of the equipment ageing. Furthermore, the Antarctic cold environment might also effect the tracking performance as the residuals over the Antarctic winter are slightly lower than the residuals over the Antarctic summer. For the phase residuals, there is little difference between the two firmware versions, with median rms around 7.6 mm. Despite several days exhibiting larger errors, more than 95% of the phase residuals are below 9 mm.

The annual and semiannual signals are well-known in GNSS time series and they are also visible in these time series. In addition, the vertical time series of the GIPSY solutions showed periodical signals at roughly 5, 14 and 27 day period at the Up component of the GIPSY solution. The 14 day signals have similar periods as the aliased periods resulted when combining unmodelled short-period tidal signals i.e., diurnal O1 tide and semidiurnal M2 tide [74,75]. The biweekly period might also be related to semi-monthly Mf tide. On the other hand, the time series reveals no clear draconitic signal at 351.2 days as reported in several studies [76,77]. According to the Rayleigh criterion, a time series longer than 25 years would be needed to be able to distinguish between the tropical year (365.25 y) and the GPS draconitic year (351.2 y) periodic signals.

Taking the environmental loading factors into account did not reduce the rms or the seasonality of the time series. Some components were maybe slightly affected, but no major improvements could be spotted. This is in accordance with our previous studies with Fennoscandian stations. We found that taking the environmental loading into account may or may not reduce the standard deviation in the time series [78]. Even when there is a clear visual correlation, the reduction of, for example non-tidal ocean loading, may not necessarily improve the results [79]. The effect on trend is not clear either, as the hydrologic or atmospheric models may not be that stable over decadal time scales. The usage of reprocessed atmospheric models is likely to reduce this uncertainty.

The choice of the noise model affects at a maximum 0.08 mm/y to the estimated trend (GIPSY Up component), but the estimated uncertainties are on average three times larger for the PL model. In addition, the RTKLIB solution includes larger annual variability compared to the GIPSY solution. The semiannual amplitudes are roughly on the same level. Furthermore, the GIPSY and RTKLIB velocity estimates differ 0.1 mm/y horizontally and 0.3–0.4 mm/y vertically. The difference in vertical velocity estimates also supports the choice of PL noise model for more realistic uncertainty estimation. The choice of PL noise is in agreement with other general studies [33,80,81,82].

The residuals may suggest a different trend for the 2003–2006 period. This can be supported by studies on the snow mass changes in Antarctica showing a significant increase in snow mass after the 2006 in the area of Aboa [83]. Increase of snow mass should be visible as decrease of the vertical rate. Basically, we could estimate a multitrend solution. However, in practice the uncertainty would be very large for the first solution (three first years) and non-usable. On the other hand, we validated our results by estimating trends using data only between 2007 and 2018 (not showed here). The estimated rates were very close to the estimate of the full time series. The uncertainties were slightly larger due to the shorter time span as expected.

4.2. How Do Our Results Relate to Other Published Data?

The horizontal motion describes mainly the dynamics of the tectonic plates. There are no published values on the horizontal rates for Aboa. On the other hand, Rülke et al. [84] reported on a campaign-based GPS horizontal rates at Wasa, a Swedish benchmark located 200 m from Aboa, used in several SCAR measuring campaigns. Using the DD approach and a 7.9 year time span, the authors estimated 10.4 ± 1 mm/y in the North direction and 0.5 ± 1 mm/y in the East direction that is within 1 mm/y from our results.

There are several global plate models including Antarctica based on either geological or geodetic data. For a comparison, we generated the plate motion velocities for Aboa as summarised in

Table 7. Tectonic plate motion at Aboa site using various plate motion models in the NNR* frame of the Antarctic plate. The velocity and speed are given in mm/y. Azimuth is measured clockwise from the North and expressed in degrees.

Model	Reference	N Vel. [mm/y]	E Vel. [mm/y]	Speed [mm/y]	Azimuth [deg]
GSRM v2.1	[86]	12.38	3.03	12.75	13.77
ITRF2014 [2016]	[51]	11.49	1.22	11.55	6.05
ITRF2008	[85]	10.86	0.66	10.88	3.46
NNR-MORVEL56	[87]	11.16	4.60	12.08	22.41
GEODVEL 2010	[88]	11.13	1.43	11.22	7.30
MORVEL 2010	[89]	11.16	4.61	12.08	22.42
This study		11.23	1.46	11.32	7.40

* NNR or no-net_rotation frame refers to the reference frame for a given model of plate motion that yields zero for the integral of the vector cross-product $vxr$ over the surface of the Earth.

.

Our GPS-derived horizontal rates at Aboa agree on 1 mm/y level for the north component, but varies up to 4 mm/y in the east component. The closest agreement is with the GEODVEL model. Also, the agreement with the ITRF2014 model is very good on the level of 0.3 mm/y. These models are based on geodetic data. The geological models (NNR-MORVEL56, MORVEL 2010) agree GEODVEL and our results very well in the North component but they differ 3 mm/year in the East component. This difference is partially explained in Altamimi et al. [85], who found a large global X-rotation rate between ITRF2008 and the NNR-MORVEL56 model equivalent to a systematic bias of 2.5 mm/y at the Earth surface.

The vertical rates from our two solutions are lower than those reported in several studies that have included Aboa in their results. Also, the uncertainty of our result is the lowest when comparing to other published values with reported uncertainties. Thomas et al. [5] reported 1.13 mm/y elastically adjusted GPS vertical rate for Aboa. Elastic response of the crust is caused by contemporaneous ice mass change and is interlaced in the vertical rate with a more long term trend. Their GPS-derived rate was 1.4 ± 0.84 mm/y, while the modelled elastic rate was 0.27 mm/y. In addition, they reported 0.09 mm/y ICESat elastic rate for Aboa. These findings are based on 1959 days of observations between 2003–031 and 2010–011. [42] assessed estimated GIA uplift rates with a set of elastic-corrected GPS vertical velocities at 67 Antarctic sites from 2009 to 2014 inclusive. However, Aboa data was limited to 2013-01-08. They reported an observation rate of 1.54 ± 0.84 mm/y. In addition, the authors reported predicted vertical rates at Aboa from several GIA models.

Table 8. Predicted vertical movement at Aboa from different GIA models [42].

Model	UP Vel. [mm/y]	Model	UP Vel. [mm/y]	Model	UP Vel. [mm/y]
W12	1.32	A13	0.65	G14	1.79
IJ05R2	0.24	AGE1b	0.19	RATES	1.89
ICE-6G	3.24	R09	1.93

summarises these values. Liu et al. [43] used similar GPS processing results as [42] but the start of the time series of beginning of 2010. After applying a spatio-temporal filtering, they found a vertical rate of 2.05 ± 0.50 for Aboa. Rülke et al. [84] did not include Aboa but reported 3.2 mm/y campaign-based GPS vertical rate at Wasa. Using ICESat observation, they report −1.2 mm/y elastic correction for Wasa. They also report the vertical rates predicted by two GIA models: W12a (1.3 mm/y) and IJ05R2 (0.7 mm/y). Their findings were based on 7.9 year time interval. An in-depth geophysical analysis of all the vertical rates (both from GPS analysis and GIA modelling) is beyond of the scope of this study. Such comparison needs extensive analysis and comprehension of the different assumptions behind each of the models above.

5. Conclusions

We estimated horizontal rates of 11.23 ± 0.09 mm/y for the East component and 1.46 ± 0.05 mm/y for the North component for Aboa in East Antarctica. These results reflect the dynamic properties of the Antarctic tectonic plate. They are in an excellent agreement with the motions derived from the GEODVEL and ITRF2014 velocity models. In vertical, the Aboa’s uplift rate is 0.79 ± 0.35 mm/y including all ongoing geophysical signals (e.g., GIA signals and current mass changes). The time series showed that the trends at Aboa have been constant since 2007 despite potential mass changes during the period of this study.

The location of the Aboa station is unique and the GPS results obtained provides valuable data for geophysical studies especially in East Antarctica. The unique location also brings challenges to the data processing as conventional models and methods have deficiencies for Antarctica. We produced rates with the highest accuracy and realistic uncertainties. Our vertical estimates are the most robust ones available for Aboa due to the thorough raw data analysis, the validation of processing and the amount of data. Our results will contributes to the ongoing discussion on better understanding the impact of Antarctica on our planet. For more in-depth geophysical interpretation, we will analyse the data history of absolute gravity measurements carried out at Aboa and include a larger scale GNSS analysis in the future.