Estimation of Lesser Antilles Vertical Velocity Fields Using a GNSS-PPP Software Comparison

Vertical land motion in insular areas is a crucial parameter to estimate the relative sea-level variations which impact coastal populations and activities. In subduction zones, it is also a relevant proxy to estimate the locking state of the plate interface. This motion can be measured using Global Navigation Satellite Systems (GNSS), such as the Global Positioning System (GPS). However, the inﬂuence of the processing software and the geodetic products (orbits and clock offsets) used for the solution remains barely considered for geophysics studies. In this study, we process GNSS observations of Guadeloupe and Martinique network (Lesser Antilles). It consists of 40 stations over a period of 18 years for the oldest site. We provide an updated vertical velocity ﬁeld determined with two different geodetic software, namely EPOS (Gendt et al, GFZ analysis center of IGS–Annual Report. IGS 1996 Annual Report, pp 169–181, 1998) and GINS (Marty et al, GINS: the CNES/GRGS GNSS scientiﬁc software. In: 3rd International colloquium scientiﬁc and fundamental aspects of the Galileo programme, ESA proceedings WPP326, vol 31, pp 8–10, 2011) using their Precise Point Positioning modes. We used the same input models and orbit and clock offset products to maintain a maximum of consistency, and then compared the obtained results to get an estimation of the time series accuracy and the software inﬂuence on the solutions. General consistency between the solutions is noted, but signiﬁcant velocity differences exist (at the mm/yr level) for some stations.


Introduction
can be done in other areas. Moreover, the islands of the volcanic arc are located too far from the trench, prohibiting most often the detection of significant velocity gradients with respect to the stable plate. Because of these reasons, the locking state deduced from GNSS observations, and thus megathrust risk is uncertain. Symithe et al. (2015) estimated with GNSS observations all along the volcanic arc a low coupling rate but did not exclude a megathrust possibility either. Since the estimation of a horizontal deformation rate is a difficult task for this area, vertical motion observations can then become a proxy to help the assessment of a potential strain accumulation (e.g. Mouslopoulou et al. 2016). Moreover, island areas are also threatened by the sea level rise, and extra subsidence can be an aggravating factor (Ballu et al. 2011). For these two reasons, measuring vertical motion in the Lesser Antilles is crucial. Some vertical movement assessments in this area were performed in the past. Paleo-geodesy based on coral reef growth tends to show a subsidence trend in Martinique and Les Saintes Islands (south of Guadeloupe archipelago) (Weil-Accardo et al. 2016;Leclerc and Feuillet 2019). This subsidence is corroborated by GNSS observations for a few stations within the vertical velocity ULR6 solution (Santamaría-Gómez et al. 2017). However, an uplift with decreasing rate for la Désirade Island (West of the Guadeloupe Archipelago) was measured (Léticée et al. 2019). Even though GNSS exploitation in this area is challenging, the islands are generally well instrumented, especially the two islands of Guadeloupe and Martinique. Therefore, this study has a double objective, technical and scientific: on one hand, we estimate a denser vertical velocity field using a maximum of GNSS data in the area. On the other hand, we compare and quantify the differences between the results (coordinates and velocities inferred from the time series) obtained with two GNSS processing software using homogeneous inputs.

Data
We considered the GNSS observations provided by three different station networks, namely the IPGP, the IGN/SONEL, and the ORPHEON networks. Maps of these networks on the two islands are represented in Fig. 1

Processing
The observation set described below was processed using two different GNSS processing software, namely EPOS (Gendt et al. 1998;Uhlemann et al. 2015) and GINS (Marty et al. 2011;Loyer et al. 2012), but using a similar strategy, equivalent models and identical product inputs. The underlying idea is to quantify the intrinsic differences due to different software. We used a Precise Point Positioning approach with float ambiguity resolution. A PPP processing is the most suitable one for this area because the reduced number of IGS reference stations prevents efficient differential processing. Moreover, those reference stations can be affected by the same tectonic processes as the geophysics stations. Thus, IGS stations in the area (namely ABMF and LMMF) are used here as "regular" stations. We considered only GPS observations since most of the stations recorded only this constellation during most of the period considered. The orbits and clock offset products have been generated beforehand by the GFZ Analysis Center in preparation of the IGS Repro3 campaign (Männel et al. 2020). These products are consistent with ITRF2014 (Altamimi et al. 2016).
Regarding the models used, we kept a consistency between the two software configurations. The same antenna eccentricities are used for both processings based on station site logs. Once the two daily coordinate sets are obtained, we select the intersection of both to get equivalent time series with the same daily coordinates. Indeed, some daily data were not properly computed by one or the other software. Stations STG0, SBL0 and PDB0 are completely excluded because of a lack of observations. For each time series, the corresponding velocities are determined using the trend estimation software HECTOR (Bos et al. 2013). We model the time series as combinations of a long term linear trend and an annual+semi-annual periodic signal, along with white and power-law noise. The term trend designates hereafter the linear component. The discontinuities introduced in the trend estimation are based on the material change site logs (an antenna change is systematically considered as a discontinuity) and on a supplementary visual detection . The same discontinuities are applied for both equivalent EPOS and GINS solutions.

Coordinate and Velocity Differences
To quantify the impact of the processings, we compute the differences between the two coordinate sets for the three topocentric components and the planimetric distance (Euclidean norm of East and North components). These differences are represented as a histogram in Fig. 3 and the statistical indicators are given in Table 1. We remark that the mean difference for the three components is not centered on zero but is shifted by some millimeters. We also remark that the Up difference doesn't respect a Normal distribution, which reflects the fact that the Up component remains the hardest component to estimate with GNSS technique, mainly because this geometrical parameter is highly correlated to the clock offsets and tropospheric delay parameters. For the planimetric distance, the mean difference is 12.33 mm  but respects a Gamma distribution. Regarding the standard deviations , North is three times smaller than the East , which can be explained by the better resolution of this component due to the general North-South trajectory of the satellites. We compute also the differences between the two sets of estimated velocities obtained at the end of the processing, represented in Fig. 4. For the absolute differences (Fig. 4a), we observe notable differences at the mm/yr level for the planimetric components, and for some stations a difference over˙1 mm/yr for the vertical component. Regarding the relative differences (Fig. 4b), we note differences on the order of˙10% for the planimetric component, but for the vertical component, these differences can vary by more than a factor of two (stations over 100%). Four stations have a negative relative variation, which means that they have an opposite velocity trend. Since the stations are sorted from the most complete set of data to the sparsest one, we observe no significative relation with the length of the time series.

Vertical Velocity Results
The vertical velocity values obtained for the two processings are given in Table 2 and represented in Fig. 5 for the Guadeloupe Archipelago, Fig. 6 for the specific area of the Soufrière Volcano, and Fig. 7 for the Martinique Island.
We note general subsidence for both islands. This tendency is consistent between the two solutions. For the Guadeloupe Archipelago, the subsidence is visible for most of the stations, nevertheless, a more complex behavior around the Soufrière area can be remarked, which might be related to local volcanic deformation but also the frequent hardware change due to the harsh conditions in the area (humidity, corrosion and frequent thunderstorms). On the Marie-Galante Island, South-West of the main island, the MGL0 station has a positive trend. Moreover, the four stations TDB0, ABD0 GOSI and DSD0, have opposite trend depending on the solution. For GOSI and DSD0, the velocities estimated in both cases are very close to zero, which make this opposite trend not significant. For GOSI, the different estimated velocities are still close to each other (with overlapping formal sigmas) but a clear velocity tendency for this station is also non-significant. The case of TDB0 is remarkable, since the time series is long and almost complete but the difference between the two solutions is important (4.6 mm/yr). A detailed view of the raw time series and the estimated tendencies are shown in Fig. 8. We observe that the different scatter for both time series lead to a completely different estimation of the trend (using the strategy we selected). Station DHS0 has a lot of corrupted raw data, which lead to a reduced amount of usable observation and an overestimated vertical velocity of almost 2 cm/yr for one solution. A similar statement can be made for FFE0 station on the main island, where several gaps in the time series along with several hardware changes might explain the positive trend estimated.
For the Martinique Island, the two solutions are also consistent and general subsidence is observed except for SAM0 station. The mean velocity rate measured is   The used days column refers to the common number of days correctly processed in both solutions, and used for the velocity estimation. The main purpose of each station in mentioned in the second column: we distinguish stations located in the vicinity of the volcano domes, and installed for volcanic deformation monitoring (V), stations for tectonic deformation monitoring (T), stations for reference frame definition and orbit determination (R), and stations for RTK surveying (S)

Comparison with Existing Solutions
To validate the consistency of our results, we compare the vertical velocities we determined with existing solutions for stations LMMF and ABMF (belonging to IGS network). We consider the velocities provided by SONEL for ULR6a (Santamaría-Gómez et al. 2017), NGL14 (Blewitt et al. 2018), JPL14 (Heflin et al. 2019) and ITRF14 (Altamimi et al. 2016) solutions. Values are given in Table 3. For LMMF, the estimated velocities are consistent with each other ( V D 0:44 mm), while for ABMF we remark more important differences between each solution ( V D 0:79 mm), just like we have also differences between the values of this work's solutions. Nevertheless, the negative trend remains significant.

Discussion
Using two different solutions but based on the same geodetic products and homogeneous models, we obtain significant disparities in terms of coordinates difference repeatability, especially on the East and Up components with a standard deviation at the centimeter level. Regarding the estimated vertical velocities using the same set of points and the same discontinuities, the differences are also significant. This result tends to motivate investigation on velocity combination strategies between different processing centers, as suggested and tested by Ballu et al. (2019) for instance, where a joint least square modeling is developed to combine equivalent time series from different Analysis Centers. A combination based on a maximum likelihood estimation would be also an relevant method. Nevertheless, for the studied area of the Guadeloupe and Martinique Islands, a negative velocity trend on the Up component is observed for most of the stations, which might suggest generalized subsidence of the area. This tendency is clear for the Martinique Island, but more complex trends for the Guadeloupe Archipelago can be observed, especially in the area around the Soufrière Volcano. This result can also be nuanced, since some stations have a positive trend, which might be due to local effects. A positive trend can also be due to an important number of discontinuities over the time series period, like the stations PAR1 (furthermore located inside the volcano area) and FFE0 (outside the volcano area). On another hand, a large number of discontinuities for the same station seem to lead also to an overestimated negative trend, like for instance the station LAM0, with a velocity estimated over 3 mm/yr for six discontinuities referenced. This statement reveals the necessity to maintain networks with a minimum of hardware discontinuities, i.e. by reducing the number of antenna changes. MGL0 station, located on the Marie-Galante Island, presents a singular behavior. It is the only station clearly uplifting, with a quasi-complete time-series of 3.5 years, without any visible discontinuity, while the other station on Marie-Galante (MAGA) presents a subsiding trend. Unfortunately, since this station belong to the commercial ORPHEON network, we have only a few metadata that prevent us to explain clearly this behavior.
We corroborate the paleo-geodesy studies carried out in the region. The coral reef records in Martinique (Weil-Accardo et al. 2016) and Les Saintes (Leclerc and Feuillet 2019) indicate also a subsidence but with a smaller order of magnitude of few tenths of a millimeter per year, which can be explained by the difference in the observation time spans (only a few years for GNSS, ca. one century for the coral records). According to those studies, long term subsidence can have multiple origins: volcanic activity, crustal faulting,  We used only one software for velocity estimation since we mainly focussed on the GNSS processing itself, but some other velocity estimator software are available (e.g. Blewitt et al. 2016;Santamaría-Gómez 2019). The impact of the velocity estimation software on solutions have been analysed for instance by Mazzotti et al. (2020).

Conclusion
This work brings a comparison of the coordinate time series obtained for the same dataset with two different software but using consistent parameters. New homogeneously calculated vertical velocity fields are made available for geophysical modeling, with unprecedented density for the two Martinique and Guadeloupe Islands. A general subsidence trend is observed for both islands.