Geopositioning time series from offshore platforms in the Adriatic Sea

We provide a dataset of 3D coordinate time series of 37 continuous GNSS stations installed for stability monitoring purposes on onshore and offshore industrial settlements along a NW-SE-oriented and ~100-km-wide belt encompassing the eastern Italian coast and the Adriatic Sea. The dataset results from the analysis performed by using different geodetic software (Bernese, GAMIT/GLOBK and GIPSY) and consists of six raw position time series solutions, referred to IGb08 and IGS14 reference frames. Time series analyses and comparisons evidence that the different solutions are consistent between them, despite the use of different software, models, strategy processing and frame realizations. We observe that the offshore stations are subject to significant seasonal oscillations probably due to seasonal environmental loads, seasonal temperature-induced platform deformation and hydrostatic pressure variations. Many stations are characterized by non-linear time series, suggesting a complex interplay between regional (long-term tectonic stress) and local sources of deformation (e.g. reservoirs depletion, sediment compaction). Computed raw time series, logs files, phasor diagrams and time series comparison plots are distributed via PANGAEA (https://www.pangaea.de).


Background & Summary
We present the results from the analysis of a new dataset of GNSS stations installed on onshore and offshore industrial settlements provided by Eni S.p.A. (https://www.eni.com). This is one of the research activities started by INGV in the framework of the program "CLYPEA: Innovation network for future energy", as part of the "subsoil deformations" project 1 . Such a dataset has been collected at continuous GNSS stations installed, respectively closely to 13 onshore infrastructures (e.g. storage centers) located along the Italian Adriatic coast and 24 offshore hydrocarbon production platforms anchored to the seabed in central and northern Adriatic Sea (Fig. 1a).
This dataset improves the current density of GNSS stations along the Adriatic coastal belt, allowing also to capture local and/or regional crustal deformation on a large sector of the northern Adriatic Sea, therefore representing a potentially invaluable dataset, although it involves GNSS stations not realized for geophysical purposes. Moreover, due to its offshore extension, the dataset must be considered as unique since until now, with the exclusion of very few GNSS stations (such as the HARV station installed on the Harvest platform, approximately located 10 km off the coast of central California 2 ), not so extensive offshore datasets have been publicly released.
The raw GNSS dataset covers a period ranging from 1998 to the end of 2017, with time series covering time spans from 0.71 to 19.08 years (Table 1 and Fig. 1b). Raw data were provided by Eni S.p.A. in the RINEX format

Methods
This section illustrates the analyses carried out by different data processing centers, which analyzed the same dataset by means of different geodetic software and, within the same software adopted different procedures for the reference frames realization. With the aim of producing the most accurate final computation and, at the same time, obtaining a comparable analysis strategy, each data analysis center processed the raw GNSS data by using (i) a double difference approach, applied using the Bernese software, (ii) a double-difference distributed approach, applied by means of the GAMIT/GLOBK software package and (iii) a Precise Point Positioning (PPP) approach, applied by means of GIPSY.
The GNSS raw dataset has been processed at four different analysis centers: • Istituto Nazionale di Geofisica e Vulcanologia -Sezione di Bologna (hereinafter INGV-BO); this analysis center provided 1 solution processed with GAMIT/GLOBK and QOCA (named IBO_GAMIT).  Table 1 for additional details).
In the following, a brief description of each adopted procedure for obtaining daily positions is provided. Table 2 resumes some relevant characteristics of each processing scheme, including specific models applied in the analysis. raw data processing by using the Bernese software. The Bernese 8 GNSS software is a scientific, high-precision, multi-GNSS data processing software developed at the Astronomical Institute of the University of Bern. Data processing with version 5.0 of this software package was performed at INGV-ONT. Daily solutions of station positions were estimated in loosely constraint reference systems. The raw observations were processed forming Ionosphere Free linear combinations and solving for the troposphere biases and phase ambiguities using the Quasi Ionosphere-free approach 9 . The troposphere modeling consisted in an a priori dry-Niell model fulfilled by the estimation of zenith delay corrections at 1-hour intervals at each site using the wet-Niell mapping function (see Table 2). In addition, one horizontal gradient parameter per day at each site was estimated. Ocean loading was computed using the FES2004 tidal model coefficients as provided by the Ocean Tide Loading provider (http:// holt.oso.chalmers.se/loading). The GPS orbits and the Earth's orientation parameters were fixed to the final IGS www.nature.com/scientificdata www.nature.com/scientificdata/ products and the site coordinates were constrained to an apriori sigma of 10 m, thus the daily coordinates were estimated in a loosely constrained, unknown reference frame. In order to express the ONT_BERNESE solution in a unique reference frame, the daily covariance was first projected imposing tight internal constraints (at millimeter level), and then the coordinates were transformed into the ITRF2008 reference frame by a 4-parameter Helmert transformation (translations and scale factor). The regional reference frame transformation used 45 anchor sites among the IGb08 stations, located on the Eurasian Plate (see Table 2). raw data processing by using the GAMIT/GLOBK software. GAMIT/GLOBK 10 is a GNSS analysis package, designed to run under any UNIX operating system and developed at the Massachusetts Institute of Technology (MIT), the Harvard-Smithsonian Center for Astrophysics, Scripps Institution of Oceanography and Australian National University. This software was adopted at DICAM-UniBO (version 10.61) and at INGV-BO and INGV-CT (version 10.7). All the three analysis centers adopted the IGS "Repro2 campaign" 11 standards during the raw dataset processing. In this step, INGV-CT and DICAM-UniBO included into the cluster processing some high-quality IGS 12 stations (Table 2) in order to improve the overall configuration of the network and to tie the regional stations to an external global reference frame. The INGV-BO solution is part of a continental-scale geodetic analysis, including >3000 continuous GNSS stations 13,14 . Because of this large number of sites, the GAMIT analysis was performed independently for several sub-networks, each made by <50 stations, with each sub-network sharing a set of high-quality IGS stations, which are used as tie-stations in the combination step.
During the processing step, the GAMIT software uses an ionosphere-free linear combination of GNSS phase observables by applying a double differencing technique to eliminate phase biases related to drifts in the satellite and receiver clock oscillators. GPS phase data were weighted according to an elevation-angle-dependent error model 10 using an iterative analysis procedure whereby the elevation dependence was determined by the observed scatter of phase residuals. In this analysis the parameters of satellites orbit were fixed to the IGS final products. IGS absolute antenna phase center models (igs08.atx and igs14_wwww.atx available at ftp://ftp.igs.org/pub/station/general/; see also Table 2) for both satellite and ground-based antennas were adopted in order to improve the accuracy of vertical site position component estimations 15,16 . The first-order ionospheric delay was eliminated by using the ionosphere-free linear combination, while a second-order ionospheric corrections 17 were applied using the IONEX files from the Center for Orbit Determination in Europe (CODE). The tropospheric delay was modeled as a piecewise linear model and estimated using the Vienna Mapping Function 1 (VMF1 18 ) with a 10° cutoff. The Earth Orientation Parameters (EOP) were tightly constrained to priori values obtained from IERS Bulletin B. The ocean tidal loading was corrected using the FES2004 19 model. The International Earth Rotation Service (IERS) 2003 model for diurnal and semidiurnal solid Earth tides was also adopted.
In a successive step, the GAMIT solutions, in the form of loosely-constrained H-files, were aligned to the IGb08 reference frame. Such an alignment was performed by INGV-CT and DICAM-UniBO through the GLOBK/GLORG 10 software package by minimizing the deviations between horizontal positions and velocities of achieved solutions and those available in the IGb08.snx 20 reference solutions (Table 2). INGV-BO, instead, performed the alignment to the IGb08 reference frame through the ST_FILTER program of the QOCA 21 www.nature.com/scientificdata www.nature.com/scientificdata/ by applying generalized constraints 21 . Specifically, the reference frame was defined by minimizing the velocities of the IGS core 23 stations, while estimating a seven-parameter transformation with respect to the GNSS realization of the ITRF2008 6 frame, i.e., the IGb08 3 reference frame. raw data processing by using the GIPSY software. GIPSY is a GNSS-Inferred Positioning System and Orbit Analysis Simulation software package (https://gipsy-oasis.jpl.nasa.gov/) developed by the Jet Propulsion Laboratory. Data processing with this software package was performed at DICAM-UniBO (which used GIPSYX, version rc0.4) and INGV-ONT (which used GIPSY-OASIS II, version 6.3).
Raw observations were processed by both analysis centers in a precise point positioning mode applied to ionospheric-free carrier phase and pseudorange data 24 and using JPL's final fiducial-free GNSS orbit products. Ocean loading tidal loading and companion tides, computed using the FES2004 19 tidal model coefficients (http:// holt.oso.chalmers.se/loading) and were applied as a station motion model. The wet zenith troposphere and two gradient parameters were estimated every 5 minutes as a random walk process 25 by using the VMF1 18 with a 10° cutoff. A second-order ionospheric correction 26 was applied during the data processing. Station clock errors was treated as a white-noise process. The ambiguity resolution was performed by using the wide lane and phase bias (WLPB) method 24 , which phase-connects individual stations to IGS stations in common view. Resolving ambiguities reduced significantly the scatter mostly in the east component of time series. Satellite orbits and clock parameters were provided by JPL who determine them in a global fiducial free analysis using a subset of the available IGS core stations as tracking sites.
In this step, the ONT_GIPSY solution was aligned to IGS14 4 by applying a daily seven-parameter Helmert transformation (three rotations, three translations, and a scale component obtained from JPL) specifically calculated using the IGS14 Cartesian coordinates and velocities of 132 stations selected by specific quality criteria (e.g. long observation time, continuity through time, position and velocity constrained with sub-millimeter level accuracy, etc.) and located in and around the Eurasian plate and in the Mediterranean area 27 . This reference frame realization method applies a continental-scale spatial filter to the station coordinates, leading to a reduction of the common-mode errors 28 and to an increase of the signal-to-noise ratio.
The UBO_GIPSY solution was aligned to the IGb08 3 reference frame by applying a daily seven-parameter Helmert transformation which was ad hoc calculated using as reference a regional subset of IGS tracking network (see Table 2). The reference coordinates of this regional subset can be found in the IGb08.snx 20 file provided by IGS.
Each analysis center, in the end, provided position time-series in a common format (e.g. the UNAVCO PBO "pos" ASCII file format 7 ); the ONT_GIPSY solution was aligned to IGS14 while the other solutions were aligned to IGb08. It must be noted that some few stations show time series with different temporal length because each analysis center requested the RINEX dataset to Eni S.p.A. at different times. Moreover, the ONT_GIPSY solution spans the 2000.0-2017.99 time interval since the limited number of continuous stations doesn't allow to adequately constraint the adopted continental-scale spatial filter prior to 2000.

Data records
All the computed raw time series are distributed via the PANGEA repository 29 so that the interested researchers can use them for future studies. The raw time series (in ASCII pos format) were stored in different folders and accordingly renamed, based on the solutions described in paragraph "Raw data processing".
We included in the same repository: 29 • A file, named "ENI_Offsets.json", containing all manually picked offsets; as mentioned above it can be used as input for later use in the TSAnalyzer software 24 . • All logs files as resulting from the analysis performed with the TSAnalyzer software; the logs files were renamed accordingly to the associated solution. These logs contain all the parameters discussed in section 3. • The phasor diagrams in png format.
• The time series plots (observed, modelled and residual components) in png format; the plots were renamed accordingly to the associated solution. • The time series comparison plots in png format.
The original raw RINEX data are property of Eni S.p.A. company. Interested users can obtain free access to data under a specific confidential agreement by contacting Eni S.p.A. company (see the Supplementary Information for details).

Technical Validation
As mentioned above, the position time series computed by each analysis center as well as the estimated velocity fields were compared in order to highlight possible differences arising from the use of the different software and from the specific models that were selected and applied in the analysis.
Position time series comparison. To perform a simple comparison between the position time series computed by each analysis center we used the TSAnalyzer 30,31 open-source software which reads GNSS position time series with different formats and allows to remove outliers as well as to simultaneously estimate linear and seasonal signals.
Roughly speaking, position time series can be modelled as the superposition of three main types of signals: (i) linear long-term deformation, (ii) seasonal signals, (iii) instrumental (e.g., GNSS antenna changes) and tectonic (co-seismic displacement) offsets. Therefore, site motion y(t) for each component can be modelled as: www.nature.com/scientificdata www.nature.com/scientificdata/ ∑ π π π π = + + + + where t i for i = 1…N are the daily solution epochs (in units of years). The terms a and b are the site position and the linear rate; coefficients c and d describe the annual periodic motion while e and f stand for the semi-annual motion. The summation term corrects for any number (n g ) of offsets, with magnitudes g and epoch T g using the Heaviside 32 step function H and the last term, ν i , is the residual.
A first visual inspection of all raw coordinates time series reveals that they contain both offsets and outliers. Outliers showing large uncertainties (e.g. blunders) were removed by setting in TSAnalyzer an arbitrary sigma criterion (we set 15, 15 and 20 mm, for North, East and Up components, respectively). Offsets were manually picked and recorded in the "ENI_Offsets.json" file (available on the repository) for later use (see also Table 3). Offsets classified as "tectonic" correspond to co-seismic deformation related to the M >6 earthquakes striking central Italy on 24 August 2016 and on 26 October 2016 33 . Offsets classified as "equipment" correspond to change in antennas while offsets classified as "unknown" are characterized by suspicious motion whose nature should be related to service operations on platform.
Finally, estimated parameters for each solution are reported as ASCII log files in the online repository. The following analysis is based on the results achieved with this time-series analysis.
The Weighted Root Mean Squares (WRMS) values of all time-series were also estimated by using the TSAnalyzer software (see log files on the repository), after filtering outliers and estimating parameters according to Eq. (1). In such a step, each single daily solution was weighted with respect to its formal sigma without taking into account the correlation between the North, East and Up components. Frequency histograms of WRMS values for horizontal and vertical components of each different solution analyzed in this study are reported in Fig. 2, while range, mean and median values are reported in Table 4. All solutions are characterized by very similar values, highlighting similar patterns between the solutions as for instance, (i) the mean and median WRMS values for vertical component are about three times larger than the horizontal ones and (ii) the higher WRMS values are observed always at the same stations (CERA and AMEB; see the related logs file in the repository for additional details). Moreover, computed WRMS values are in the same range of the ones estimated for on-shore stations with comparable observing time windows 13,34 .
All the examined time series contain significant seasonal signals in both horizontal and vertical directions. Values estimated for each solution show a large agreement between them; annual amplitudes range between 0.5-5.3 mm horizontally and 3-8 mm vertically with higher values observed at sites installed on the offshore  Table 3. Offsets manually picked and recorded in the "ENI_Offsets.json" file which is available on the repository.
www.nature.com/scientificdata www.nature.com/scientificdata/ platforms (see plots reported in the online repository). Semi-annual amplitudes are usually <0.5 mm horizontally and ~1 mm vertically. Additional features on seasonal signals are provided by considering phasor diagrams of amplitudes and phase signals (see the example reported in Fig. 3). Regarding the annual signals, most of the sites installed on the offshore platforms show maximum amplitudes on January, September and August for North, East and Up components, respectively. Moreover, most of these sites (especially those closely located to the coastal area) are subjected to horizontal oscillations with a prevailing NNW-SSE attitude and to vertical oscillation associated with uplift during summer and subsidence during winter. Sites installed along the coastal www.nature.com/scientificdata www.nature.com/scientificdata/ area are characterized by maximum amplitudes largely scattered and no clear oscillation patterns can be recognized. Concerning the semiannual signals, most of the sites show maximum amplitudes concentrated during July-September for North component and during April-June for both East and Up components.
Some position time series (mainly recorded at the stations located on the offshore platforms) exhibit a curved shape, while others show weak post-seismic decay (recorded at few stations located along the coastal area, eastward of the 2016 seismic sequence epicentral area). Both features result in deviations from the linear trend; to quantify such a derivation for each time series, we defined a simple non-linearity index (I NL ) as the ratio between the standard deviation of the smoothed time series and the standard deviation of its intrinsic noise:  Table 5). Estimations of I NL index on a set of high-quality IGS stations ( Table 2) highlight that values <0.5 reflect time series with pure linear behavior, therefore in the following we arbitrarily set as "gentle" and "moderately" non-linear all the time series with 0.5 ≤ I NL ≤ 1 and I NL > 1, respectively (Table 5). Based on this simple definition and excluding FRAT station because of its short time series, 11   Table 1). Amplitudes of the estimated sine and cosine parameters (see Eq. (1)) are plotted for the North, East and Up components. The upper-right plot represents the key to correlate the maxima phase direction with the day year period. Phases are referred to January 1 and time increases clockwise.
www.nature.com/scientificdata www.nature.com/scientificdata/ stations are generally characterized by linear time series (for all the three components) while the other are characterized by gentle (19 stations) and moderately (6 stations) non-linear time series, at one or more components. Stations with gentle and moderately non-linear time series mainly concentrates offshore on northern Adriatic Sea and onshore, along the coastal belt; the stations with linear time series are mainly located along the coastal onshore belt.
As a last step, the time-series were compared between them after correcting the offsets (Table 3) and removing the linear trend as defined in Eq. (1). An example of such a comparison is reported in Fig. 4, while the remaining plots are reported in the online repository. This simple comparison highlights how the main features characterizing the time-series (e.g. seasonal signals, noisier time intervals, short-term transients, etc.) show a good agreement between the different solutions.
Velocity field analysis and comparison. In the following, we performed some simple comparisons between all the computed linear trend values; however, these analyses (and related results) must be considered with caution because of the non-linear motion detected at many stations.
As a first step, we calculated the residual values with respect to the mean ones for each solution and for each component (Table 6). Results are reported as frequency histograms in Fig. 5   As a second step, we computed for each solution the velocity field with respect to a Eurasian 35 reference frame; achieved results are reported in Fig. 6. Because of its simple definition, such a rotation affects only the horizontal velocities; therefore, the vertical ones remained constrained in their previous reference frame. Considering the horizontal velocity (Fig. 6), all the solutions show a general agreement in terms of azimuthal pattern and vector magnitudes, however the ONT_BERNESE solution show a systematic counterclockwise rotation with respect to the other ones, probably related to the Helmert constraints imposed for the reference frame transformation ( Table 2). The ONT_GIPSY solution, although referred to IGS14, is highly coherent with the other solutions. Most of the stations installed along the onshore coastal area show a prevailing NNW-oriented velocity pattern which is highly coherent with the regional deformation field 14,27,36 . Stations installed along the Adriatic offshore show a complex deformation field characterized by large variations both in the azimuthal pattern and the vector magnitudes (Fig. 6).
The vertical velocities are reported in Fig. 7 as averaged values (panel a) and as differences with respect to the average values (panel b). Figure 7a shows that most of the stations installed along the onshore coastal area are characterized by subsidence with rates up to 2 mm/yr in agreement with recent studies 8 , while all the stations installed on the offshore platforms exhibit a general subsidence with rates up to ~17 mm/yr (see also Table 6). The site-by-site comparisons reported in Fig. 7b confirm the features previously recognized: the ONT_GIPSY and ONT_BERNESE solutions are generally smaller than the mean values, while the other solutions are larger than www.nature.com/scientificdata www.nature.com/scientificdata/ the mean values. These differences seem related to the adopted reference (the ONT_GIPSY solution is referred to IGS14) and/or to the constraints imposed for the reference frame transformation (the ONT_BERNESE solution adopts a Helmert transformation based only on 4-parameters). We detected some offsets corresponding to equipment changes and to co-seismic deformation. Moreover, some detected offsets, classified as "unknown" in Table 3, are characterized by motion probably related to service operations on platform. All detected offsets have been identified in all solutions and are characterized by very similar values (see related log files on the online repository).
Examined position time series contain significant seasonal signals: annual amplitudes range between 0.5-5.3 mm horizontally and 3-8 mm vertically while semi-annual amplitudes are usually < 0.5 mm horizontally and ~1 mm vertically. Higher annual amplitude values were observed only at sites installed on the offshore platforms, with maximum amplitudes on January, September and August for North, East and Up components, respectively. The horizontal oscillations occur with a prevailing NNW-SSE attitude while the vertical oscillation are associated with uplift during summer and subsidence during winter. Conversely, sites installed along the coastal area are characterized by maximum amplitudes largely scattered and no clear oscillation patterns can be recognized. All these observations suggest that offshore stations would be affected by variations in hydrostatic pressure and  www.nature.com/scientificdata www.nature.com/scientificdata/ buoyancy on tubular members of the platforms caused by changes in tide levels. Moreover, since the Adriatic Sea is an almost land-locked basin where atmospheric conditions and vary considerably with the seasons 37,38 , environmental loads caused by wind (e.g. periodic seiche-like effects), current flows and waves would substantially contribute to the observed seasonal oscillations. Due to the main metallic-fabric of platform, seasonal thermal expansion-contraction cycles of the tiny structures could also contribute to the seasonal oscillations.
A large number of stations is characterized also by position time series with non-linear behavior at one or more components, as quantified by the computed I NL index (Table 5) (Table 5).
www.nature.com/scientificdata www.nature.com/scientificdata/ are generally characterized also by high subsidence rates (Table 6 and Fig. 7), clearly suggesting that, the tectonic deformation is superseded by the local sources ones (e.g. reservoir depletion, sediment compaction, etc.). Moreover, these local sources affected also the horizontal deformation rates leading to a complex pattern, strictly depending by the relative position of the station (i.e. the platform) with respect to the sources of deformation. Similar considerations can be done for the stations with 0.5 ≤ I NL ≤ 1, where however, the effects of local sources have had a minor influence on station motion. Conversely, stations with I NL < 0.5 show a velocity pattern which is highly coherent with the regional deformation field of the investigated area.
Based on the analysis and the comparisons performed in this study, the different time series solutions are highly consistent between them despite the use of different software, models, strategy processing and frame www.nature.com/scientificdata www.nature.com/scientificdata/ realizations. The analyzed dataset represents an invaluable dataset since it allows to improve the current GNSS stations density along the Adriatic coastal belt. In addition, this dataset allowed to discover a complex interplay between regional and local sources of deformation, whose relationships would be addressed in future analysis. Indeed, because of their uniqueness, these data could improve and promote further studies in offshore exploiting contexts including: (i) measurement and modelling of induced subsidence patterns and their spatial and temporal correlation with production data, (ii) estimation of natural and anthropic contributions to the overall ground subsidence, (iii) improvement on the long-term regional crustal deformation on offshore areas, (iv) coast-line dynamic and impact on human activities and natural ecosystems (e.g. coast line setback).