Layout 6 1 ANNALS OF GEOPHYSICS, 63, 3, GD327, 2020; doi:10.4401/ag-8180 M2 constituent of ocean tide loading displacements from VLBI CONT14 hourly sessions Kamil Teke1 (1) Hacettepe University, Department of Geomatics Engineering, 06800, Ankara, Turkey Article history: received April 16, 2019; accepted October 14, 2019 Abstract Several studies prove that ocean tide loading (OTL) displacements can be observed with space geodetic techniques. In this study, the amplitudes and Greenwich phase lags for each coordinate component, i.e., radial, west, and south of the principal lunar semidiurnal tide, M2 of OTL displacements were estimated at the very long baseline interferometry (VLBI) sites of the 15 days long continuous VLBI campaign, CONT14, carried out by the International VLBI Service for Geodesy and Astrometry (IVS). In the estimation of the amplitudes and Greenwich phase lags of the M2 tidal constituent, hourly VLBI station coordinate time series were used as observations derived through analyzing 1 hour VLBI sessions of the CONT14 campaign. In the analysis of hourly sessions of the CONT14 campaign, to derive accurate hourly station coordinates, troposphere delays estimated from daily sessions were reduced from the observations a priori to the analysis. The estimated amplitudes and Greenwich phase lags of the M2 constituent of OTL displacements were compared with the predictions the state-of-the-art ocean tide models, among others, FES2012 [Lyard et al., 2006; Carrère et al., 2012], FES2014 [Carrère et al., 2016], and TPXO8 [Egbert and Erofeeva, 2002; Egbert et al., 2010]. Both the amplitudes and the phases between CONT14 estimates and ocean tide models agree well for the M2 tide at all the sites and in most of the coordinate components. The RMS misfits of the M2 tide of OTL displacements in all coordinate components between CONT14 and ocean tide models over coastal sites are found about two times larger than those of inland sites. This result confirms the modeling insufficiencies in shallow waters of ocean tide models which cause an accuracy restriction of OTL displacement predictions around coastal regions. Keywords: VLBI; troposphere zenith delays; CONT14 hourly sessions; ocean tide loading displacements; M2 constituent. 1. Introduction The seafloor pressure variations due to the ocean tide loading (OTL) cause position and tidal frequency- dependent harmonic displacements on the Earth crust, the so-called OTL displacements. OTL displacements can be predicted by convolution software, e.g., OLFG/OLMPP [Scherneck, 1991], SPOTL [Agnew, 1996], NLOADF [Agnew, 1997], GOTIC2 [Matsumoto et al., 2001], CARGA [Bos and Baker, 2005], LoadDef [Martens et al., 2016] which convolve the elastic load Green’s function over the gridded values of global ocean tide models such as SCHW81 [Schwiderski, 1980], TPXO8 [Egbert and Erofeeva, 2002; Egbert et al., 2010], FES2012 [Lyard et al., 2006; Carrère et al., 2012], or FES2014 [Carrère et al., 2016] for the whole oceans. Green’s function formulates the elastic response of the Earth’s crust to the surface load that depends on the angular distance to the load as well as Earth structure [Farrell, 1972]. The sparse spatial resolution of global tide models and modeling insufficiencies in shallow waters have been proved to considerably restrict the accuracy of OTL displacements around coastal regions (up to ~150 km distance from the coast) [e.g. Agnew, 1997; Khan and Scherneck, 2003; Bos and Baker, 2005; Penna et al., 2008; Yuan and Chao, 2012] whereas tidal heights derived from global ocean models do agree within 2-3  cm in deep oceans (Shum et al., 1997). Further details and discussions on the calculation of OTL displacements from ocean tide models and errors in the OTL displacements propagated from ocean tide models are provided, e.g., in Farrell [1972, 1973], Pagiatakis [1990, 1992], Scherneck [1991, 1993], Agnew [1996, 1997], Scherneck and Bos [2002], and Baker and Bos [2003]. Several studies prove that OTL displacements can be observed with space geodetic techniques. In literature, two main approaches are proposed for resolving the tidal harmonic constituents from the observations of space geodetic techniques, static and kinematic. In the static approach - also called harmonic parameter estimation approach [e.g., Penna et al., 2015] real (in-phase) and imaginary (out-of-phase) components of OTL displacements tidal constituents are included in the observation model along with daily station coordinates as part of daily solutions. Then, the daily estimates of tidal constituents and their covariance information are stacked in a combined solution, e.g. using a Kalman Filter as often implemented in Global Navigation Satellite System (GNSS) analysis [e.g., Schenewerk et al., 2001; Dach and Dietrich, 2001; Allinson et al., 2004; King et al., 2005; Thomas et al., 2007; Yuan et al., 2013] and in several very long baseline interferometry (VLBI) analysis [e.g., Schuh and Moehlmann, 1989; Sovers, 1994; Haas and Schuh, 1998; Scherneck et al., 2000; Petrov and Ma, 2003]. Yuan et al., [2013] estimated OTL displacements 3D components of eight principal harmonics at semi-diurnal and diurnal tidal periods by means of analyzing the observations of globally distributed 456 continuous GPS stations since 1996 till 2011. They used the precise point positioning technique [PPP, Zumberge et al., 1997] on daily GPS batches where the tidal constituents were treated as additional parameters. Then, they combined the resultant daily estimates of tidal constituents with their variance-covariance matrices. Among other results, Yuan et al. [2013] found the accuracy of their lunar only constituent estimates, M2, N2, O1, and Q1 in horizontal and vertical components better than 0.12 mm and 0.24 mm, respectively. Analyzing more than 3 million VLBI observations of 3126 sessions observed from April 1980 to January 2002 and considering 40 stations, Petrov and Ma [2003] showed that the VLBI technique is capable of observing ocean tide loading displacements with the average accuracy of ~0.5 mm and ~1.7 mm in the horizontal and radial coordinate components, respectively. For 28 co- located sites Petrov and Ma (2003) found that VLBI has a better agreement in vertical displacement amplitudes at 8 principal tidal waves than GPS [Schenewerk et al., 2001] if compared to the GOT00.2 (Ray 1999) ocean tide model values. In the kinematic approach, station coordinate time series and their formal errors are estimated at sub-daily intervals (usually in 1-2 hour batches) without reducing the OTL displacements from the observations a priori to the parameter estimation. Then, the amplitudes and phases of OTL displacements tidal constituents are extracted through harmonic analysis of sub-daily station coordinate time series, thereby treating the station coordinates as independent tide gauge measurements. The kinematic approach is implemented in several GNSS analyses, [e.g, by Baker et al., 1995; Dragert et al., 2000; Khan and Tscherning, 2001; Khan and Scherneck, 2003; Vey et al., 2002; King, 2006; Yun et al., 2007; Melachroinos et al., 2008; Vergnolle et al., 2008; Penna et al., 2015; Martens et al., 2016]. Dragert et al., [2000] estimated the position time series of the coastal GPS site HOLB (Canada) using 4 weeks of continuous GPS data for 3 hour and daily solutions. According to Dragert et al. [2000], un-modeled parts of ocean loading displacements propagate to daily averages of station positions as minimal biases. Besides, they found that introducing ocean loading corrections in daily GPS solutions and constraining e.g. position to daily average, the vertical tidal motion propagates into the hourly troposphere delay estimates, thereby biasing the hourly zenith delay estimates by up to 1 cm. Khan and Tscherning [2001] estimated differential amplitudes of M2 and N2 between two GPS stations (Fair and Chi3) from hourly solutions of standard relative GPS positioning technique and found similar M2 and N2 differential amplitudes to those of the GOT99.2 [Ray, 1999] ocean tide model. Due to the fact that estimating troposphere zenith delay in hourly intervals absorbs and eliminates the loading signal [Dragert et al., 2000; Khan and Tscherning, 2001] fixed troposphere zenith delays to the Saastamoinen [1972] model in addition to daily L3 ambiguities in their hourly solutions. As a follow up study, Khan and Scherneck [2003] could also separate differential vertical and north-south amplitudes and phases of the Kamil Teke 2 M2 tide from zenith total delay (ZTD) estimates through hourly solutions of GPS observations at the sites, Fair and Chi3 in Alaska during a period of 49 days and validated their results using those of the GOT99.2 ocean tide model. Vey et al. [2002] investigated the impact of ocean tide loading corrections on GPS ZTD estimates by means of analyzing the observations of 3 hour sessions with the contribution of six GPS stations over a period of 3 days and estimated ZTD with and without introducing OTL displacements using the CSR4.0 ocean tide model [Eanes and Shuler, 1999]. By means of comparing the correlation between the OTL vertical position displacements with the ZTD differences of these two solutions as a rule of thumb, Vey et al. [2002] concluded that 4.4 cm station height error would cause an error of 1 cm ZTD. Similarly, based on the VLBI observations of the CONT11 campaign, Teke et al. [2013] found about 1 cm zenith delay differences between the analyses of 24 hours and 2-hour sessions as well as the radial components of station positions differ by about 2 to 4 cm depending on the number of observations and their sky coverage in 2-hour sessions. The main purpose of this study is to investigate the possibility of estimating M2 constituent of OTL displacements from the VLBI observations carried out during the 15 days long continuous VLBI campaign, CONT14 and investigate the level of agreement of the estimated M2 constituent of OTL displacements with the predictions of the ocean tide models FES2012 [Lyard et al., 2006; Carrère et al., 2012], FES2014 [Carrère et al., 2016] and TPXO8 [Egbert and Erofeeva 2002; Egbert et al., 2010], which were generated from global hydrodynamic simulations that assimilate tide gauge and satellite altimetry data. Troposphere delay is the largest portion of error in the reduced observations (observed minus computed vector) of radio-frequency space geodetic techniques. Relatively to the other geophysical effects, troposphere delays cannot be reduced from the observations at the desired level with the widely-used present troposphere delay models e.g. Saastamoinen [1973], Davis et al. [1993], Böhm et al. [2006], and Chen and Herring [1997]. The original value of this study is reducing the relatively accurate troposphere delays, estimated through the analysis of 24- hour sessions, from the observations of 1 hour sessions a priori to the parameter estimation. This contributes to the reliable estimation of sub-daily varying geodetic and geophysical parameters. For this study, among others the sub-daily varying geophysical parameter is selected as M2 tidal signal of OTL displacements and the investigation is focused on the possibility of estimating this sub-daily parameter with such a small data set (about 360 positions for a station evenly distributed at hourly epochs covering 15 days) i.e. CONT14 campaign. The reason of choosing the M2 tide as a metric for the success of observing the OTL displacements from the 15 days long continuous VLBI observations, CONT14 is that it has the largest amplitude among those of other principal semi-diurnal and diurnal tides of OTL displacements that leads to more distinct comparisons with those of ocean tide models. In this study, CONT14 hourly sessions were obtained by dividing each of the original 15 1-day sessions into 24 1-hour segments. An intrinsic parameter estimation strategy for the analysis of 1 hour VLBI sessions of the CONT14 campaign was used by means of reducing the troposphere delays estimated from 24-hour sessions and using a CONT14 specific terrestrial reference frame (TRF) and Earth orientation parameters (EOP) series as a priori values. It is worth to note that this paper does not aim to attain better accuracies of M2 OTL constituent from VLBI CONT14 campaign than those of ocean tide models. The experimental work presented in this paper provides investigations into how reliable and accurate M2 tidal signals of OTL displacements can be estimated from the sub-daily i.e. hourly VLBI station positions over 15 days. Only M2 tide is estimated in this study and the station displacements caused by the remaining; long period, semi-diurnal and diurnal tides are calculated from FES2014 model and reduced from the hourly station coordinates a priori to the estimation. Because the 15 days long sub-daily i.e. hourly station coordinate time series over CONT14 campaign does not satisfy the minimum period of Rayleigh criterion [Foreman, 1977] to distinguish between neighboring frequencies. As an example, according to the Rayleigh criterion to distinguish between M2 tide with a period of 12.42  hour and N2 tide (12.66 hours) at least 27.3 days of data is needed. In Section 2 the VLBI CONT14 campaign is briefly described, and then a detailed explanation of the analysis of VLBI hourly sessions is presented focusing on the reduction of the 24-hour session troposphere delays from the observations of 1 hour sessions a priori to the adjustment. In Section 3 the estimation procedure of the M2 tide amplitudes and phase lags of OTL displacements at VLBI sites of CONT14 from the hourly station coordinate time series is introduced, and the level of agreement between the estimates of this paper and those of the FES2012, FES2014, and TPXO8 ocean tide models is presented. In Section 4 the impact of a priori ocean tide model errors on the estimated M2 tide displacements is discussed in the scope of the analysis of 1 hour sessions over 15 days. 3 Ocean loading M2 tide displacements from VLBI Kamil Teke 4 2. CONT14 campaign and analysis of VLBI hourly sessions The International VLBI Service for Geodesy and Astrometry [IVS, Nothnagel  et al.,  2017; Behrend,  2013; Schuh and Behrend, 2012] carries out continuous VLBI campaigns over two weeks with the state-of-the-art VLBI technology every three years. One of the scientific opportunities enabled by these campaigns is to address the discrepancies between ocean tide models and VLBI observations at principal diurnal, semidiurnal, and even ter-diurnal periods. In this study, the CONT14 campaign is analyzed, which was continuously observed by the IVS from 6 (0  UT) to 20 (24  UT) May 2014 to provide the highest accuracy of geodetic parameters that VLBI can currently provide. The CONT14 observation network consists of 17 VLBI stations located at 16 sites. Homogeneous distribution of the network sites is provided as far as possible through selecting 6 sites in the southern hemisphere (Figure 1). For further information on the CONT14 campaign, e.g. scheduling, performing correlation, technical and scientific prospects, readers are referred to the web site https://ivscc.gsfc.nasa.gov. The geodetic VLBI observatories (stations) contributing to the CONT14 campaign are listed in Table 1. The 17 VLBI stations are separated into inland (7) stations and coastal (10) stations, depending on whether the station is located within 150 km from the coastline or not (Table 1). Although the VLBI stations at Zelenchukskaya (Russia) and Matera (Italy) are located within 150 km from the coastline, they are considered as inland stations in this paper due to the very low tidal amplitudes at the Black Sea and the Mediterranean Sea. In contrast to 24-hour session analysis, the main disadvantage of analyzing the observations of 1 hour VLBI sessions is that reliable station positions cannot be estimated due to the high correlation between troposphere delays and station positions [e.g. Rothacher and Beutler, 1998; Teke et al., 2013]. To show the amount of shared variances (degree of linear relationship) between the ZWD and radial coordinates estimated simultaneously from 1 hour sessions of the CONT14 campaign, ZWD from 24 hour sessions were subtracted from those of 1 hour sessions, ΔZWD1H-24H=ZWD1H-ZWD24H. Then, Pearson correlation coefficients were calculated between ΔZWD1H-24H and the radial coordinates of the stations for each 1 hour session. All negative correlations between ZWD1H-ZWD24H and radial estimates from 1 hour session analysis at the VLBI stations of CONT14 are found between -0.74 and -0.93. The p-values of these correlations are below 0.05 indicating that they are all statistically significant. It is inferred Figure 1. The VLBI stations contributing to the CONT14 campaign. https://ivscc.gsfc.nasa.gov that, as a rule of thumb, 1 cm ZWD1H variation propagates into the radial positions approximately 1.5 to 2.5 cm when 1 hour session is analyzed as well as these two parameters are simultaneously estimated (see e.g. Figure 2 for Kokee VLBI station and supplementary material for the other stations). According to Teke et al. [2013] depending on the number and sky distribution of the observations over the hourly sessions, correlations between antenna radial positions and troposphere delays are varying. They also state that these correlations between these two parameters can be mitigated if homogeneously distributed adequate number of observations are carried out at each antenna and at each hourly session. Consequently, the classical Gauss-Markov least squares adjustment cannot decorrelate these two parameters when they are estimated simultaneously in sub-daily e.g. 1 hour intervals. Thus, troposphere delays and antenna coordinates propagate into each other and resulting in unreliable estimates. For example, the cyan and grey lines in Figure 3 show the very noisy ZWD and radial coordinate components estimated simultaneously from the analysis of 1 hour sessions at the VLBI station Kokee and for other stations readers are referred to the supplementary material of this paper. To overcome this restriction, external troposphere slant delays Δ𝐿 estimated from daily (24 hour) sessions were reduced from the observations of hourly sessions a priori to the adjustment and residual troposphere delays were not estimated. The troposphere delay model from Davis et al., (1993) Δ𝐿(𝛼,𝜀) = ZHD 𝑚ℎ (𝜀)+ZWD 𝑚𝑤 (𝜀)+ 𝑚𝑔 (𝜀)[𝐺𝑛 cos(𝛼)+𝐺𝑒 sin(𝛼)] (1) was used to compute the slant delays Δ𝐿 so as to reduce them from the VLBI observations of hourly sessions. In Equation (1) 𝛼 denotes the azimuth of the observation, 𝜀 the outgoing elevation angle from the local horizon, ZHD the troposphere zenith hydrostatic delay, ZWD the zenith wet delay, 𝑚ℎ the troposphere hydrostatic mapping 5 Ocean loading M2 tide displacements from VLBI Table 1. Geodetic VLBI stations participating in the CONT14 campaign with their geographical coordinates and tide-free ellipsoidal heights of VieTRF13b [Krásná et al., 2014]. The VLBI stations which are closer than 150 km distance from the nearest coastline are assumed as coastal stations. The VLBI sites are ordered according to the latitudes of the sites from north to south. Observatory Name Country VLBI Acronym Latitude (°) Longitude (°) Ellipsoidal height (m) Region (distance to the nearest coastline in km) Ny Ålesund Norway NYALES20 78.93 11.87 87.79 coast (2) Onsala Sweden ONSALA60 57.39 11.93 59.73 coast (1) Badary Russia BADARY 51.77 102.23 822.02 inland (1849) Wettzell Germany WETTZELL 49.14 12.88 669.56 inland (374) Zelenchukskaya Russia ZELENCHK 43.79 41.57 1175.43 inland (98) Westford USA WESTFORD 42.61 288.51 87.19 coast (44) Matera Italy MATERA 40.65 16.70 543.80 inland (29) Yebes Spain YEBES40M 40.52 356.91 989.17 inland (261) Tsukuba Japan TSUKUB32 36.10 140.09 85.14 coast (44) Kokee Park USA KOKEE 22.13 200.33 1168.72 coast (5) Fortaleza Brazil FORTLEZA -3.88 321.57 23.48 coast (5) Katherine Australia KATH12M -14.38 132.15 189.55 inland (227) Hartebeesthoek South Africa HART15M -25.89 27.68 1409.82 inland (482) Yarragadee Australia YARRA12M -29.05 115.35 248.47 coast (47) Warkworth New Zealand WARK12M -36.43 174.66 128.09 coast (3) Hobart Tasmania HOBART26 -42.80 147.44 65.52 coast (8) Hobart Tasmania HOBART12 -42.80 147.44 41.39 coast (8) function and 𝑚𝑤 the wet mapping function, 𝐺𝑛 and 𝐺𝑒 are the north and east total horizontal gradients, respectively. In Equation (1), the formulation by Chen and Herring (1997) with C=0.0032 as a gradient mapping function 𝑚𝑔, was used. For all solutions of hourly sessions, the ZHD were calculated at observation epochs with total surface pressure values measured at the VLBI sites [Saastamoinen, 1972; Saastamoinen, 1973; Davis et al., 1985]. ZWD were estimated at 20-minute intervals with relative loose constraints as 1.5 cm after 20 minutes from 24-hour sessions. Note that in the analysis of these 24-hour sessions FES2004 ocean tide model [Lyard et al., 2006] displacements are introduced to the a priori station coordinates (solution: Case 1, see Table 3 in Section 4). Then, these ZWD estimated from 24 hour VLBI sessions were linearly interpolated to observation epochs of hourly sessions. ZHD and ZWD were mapped with the Vienna Mapping Functions  1 [VMF1, Böhm  et al.,  2006] to get the slant hydrostatic and slant wet delays. Kamil Teke 6 Figure 2. The circles show the ZWD1H - ZWD24H versus hourly radial positions estimated simultaneously once for each 1 hour session of the CONT14 campaign at Kokee VLBI station. The correlation between ZWD1H - ZWD24H and radial positions is found as -0.80 at this station. Figure 3. Cyan and grey lines show the hourly ZWD and the radial coordinate components estimated simultaneously from the analysis of 1 hour sessions at the VLBI station Kokee, respectively. Blue line illustrates the ZWD estimated from 24 hour sessions. Blackline shows the radial coordinate components estimated from 1 hour sessions when ZWD from 24 hour sessions (blue line) are reduced from the observations a priori to the parameter estimation. Similar to ZWD, troposphere total north and east gradients at 1 hour intervals with relative loose constraints as 1 mm after 1 hour estimated from 24 hour VLBI sessions were linearly interpolated to the observation epochs of hourly sessions. Then, azimuthally asymmetric troposphere delays were calculated through mapping these horizontal total north and east gradients to slant direction using the third term of Equation (1) where the gradient mapping function by Chen  and  Herring  [1997] was used. Finally, troposphere slant delays Δ𝐿 from each VLBI observation (delay) of the hourly sessions were subtracted a priori to the parameter estimation. As previously mentioned, due to the small number and inhomogeneous sky coverage of observations as well as the weak global coverage of the stations, each hourly session was analyzed with a specific parameterization. The terrestrial reference frame (TRF) catalogue, series of daily nutation offsets and hourly Earth rotation parameters intrinsic to the CONT14 campaign (covering only the time span and observations of CONT14 campaign) were estimated from a series of global solutions, used as a priori values and held as fixed parameters in the analysis of hourly sessions. Specific to each global solution the following parameterization was introduced: • CONT14-TRF catalog: A CONT14 specific TRF catalog from a global solution with the observations of CONT14 (referred in this paper as CONT14-TRF) was estimated. In this global TRF solution no-net-translation (NNT) and no-net rotation (NNR), TRF datum conditions with respect to VieTRF13b catalog [Krásná et al., 2014] were introduced to the accumulated datum free normal equation system and station TRF velocities were fixed to those of VieTRF13b. The datum conditions were not imposed on the estimated coordinates of the VLBI stations TSUKUB32 and WARK12M since VieTRF13b coordinates of these stations are not available for the CONT14 period. The reason for creating CONT14 specific TRF, CONT14-TRF is to eliminate the biases on estimated coordinates resulting from e.g. the inadequacies in the a priori TRF coordinates and velocities as well as unmodeled episodic displacements. • CONT14-EOP-D series: Daily Earth orientation parameters for CONT14 were estimated in a global solution (referred in this paper as CONT14-EOP-D) of which a priori nutation values were taken from the IERS C04 08 corrections [Bizouard  and  Gambis,  2009] in addition to the IAU2006 precession-nutation model [Petit and Luzum, 2010]. Thus, VLBI only derived daily EOP corrections intrinsic to CONT14 were introduced to the a priori IERS C04 08 series. The reason for a global solution of daily EOP series, CONT14-EOP-D is to derive more accurate daily nutation parameters specific to the CONT14 campaign. • CONT14-ERP-H series: The hourly Earth rotation parameters (ERP: polar motion coordinates and Earth’s daily rotation phase angle, UT1-UTC) for CONT14 were estimated (referred in this paper as CONT14-ERP-H) where a priori daily nutation offsets are fixed to CONT14-EOP-D and a priori ERP were taken from CONT14-EOP-D plus high-frequency corrections. In all hourly session analyses, interpolated values of the CONT14-ERP-H series to observation epochs were fixed. The parameterization for the analysis of the hourly VLBI sessions is as follows: The CONT14 VLBI observations were analyzed using Vienna VLBI and Satellite Software [VieVS, Böhm et al., 2018]. The parameters were estimated using the classical Gauss Markov least-squares adjustment method. The observations were not removed below a certain elevation angle and not down-weighted. Source coordinates were fixed to ICRF2 (International Celestial Reference Frame 2, Fey et al., 2009). The high-frequency variations of Earth rotation parameters were modeled as recommended by the IERS Conventions 2010 [Petit and Luzum, 2010]. Tidal and nontidal atmospheric loading [Petrov and Boy, 2004], as well as ocean tide loading based on the ocean model FES2014 [Carrère et al., 2016], were introduced to each observation a priori to the adjustment. A priori station coordinates were not corrected for hydrological loading and non-tidal ocean loading. Troposphere delays from 24-hour sessions were reduced from the observations a priori to the adjustment and residual delays upon the a priori ZWD, and corrections to the a priori north and east gradients were not estimated. Daily nutation offsets and hourly ERP were fixed to CONT14-EOP-D and CONT14-ERP-H series, respectively. For each VLBI hourly session, clock synchronization errors as one offset and a rate of a polynomial for each clock with respect to a fixed clock and antenna coordinates with NNT/NNR conditions with respect to CONT14-TRF were estimated. Analyzing CONT14 hourly sessions were repeated 17 times, i.e., the same number as stations in CONT14, using the same parameterization at each run with the difference that the FES2014 ocean tide model displacements were not introduced to the a priori coordinates of one of the VLBI stations and NNT/NNR datum conditions were not imposed on those station coordinates. At least by doing so the transformation defect of NNT/NNR TRF datum constraints on the estimated hourly coordinates of the concerning station was recovered. To prevent singularity of the normal equation system solution of hourly sessions, as an 7 Ocean loading M2 tide displacements from VLBI experience, the antennas of which the number of observations is less than 20, were excluded from the analysis. A posteriori variances (Chi2) of 9 sessions out of 360 were found larger than 2. Thus, these sessions were excluded from the hourly session set so the remaining 351 hourly sessions were considered for each run. Through introducing the solid Earth tide displacements, which can be modeled with an often-stated accuracy of around 1%, as well as pole tide, tidal and non-tidal atmosphere loading corrections to the a priori station coordinates, at the final point, OTL displacements were assumed to be unveiled by the estimated hourly station positions (e.g. see Figure 4 for the VLBI site, Onsala in Sweden). 3. M2 constituent of OTL displacements from CONT14 hourly sessions and assessment on the level of agreement with those of the selected ocean tide models The tidal frequency and position-dependent harmonic displacements on the Earth crust, caused by the seafloor pressure variations due to the ocean tide loading (OTL), Δ𝑛,𝑘 for the 𝑘’th coordinate component (i.e. radial, west, or south) at a site, 𝑛 and at a particular time 𝑡 can be formulated with the harmonic function Δ𝑛,𝑘 =�𝐴𝑛,𝑗,𝑘 cos(𝜒𝑗(𝑡)‒𝜑𝑛,𝑗,𝑘) (2) � where 𝜒𝑗 denotes the astronomical phase of the 𝑗’th tidal constituent which is computed from fundamental astronomical arguments [Petit and Luzum, 2010]. In Equation (2), 𝐴𝑛,𝑗,𝑘 and 𝜑𝑛,𝑗,𝑘 are the amplitude and the phase lag with respect to the 𝑗’th tidal potential at Greenwich meridian, respectively. The Darwinian designations of the principal tides at semi-diurnal periods M2, S2, N2, K2, at diurnal periods K1, O1, P1, Q1, and at longer periods Mf, Mm, and Ssa describe more than 95% of the tidal signal [Lambeck, 1988]. Principal lunar semi-diurnal tide, M2 is the strongest tidal constituent that has a frequency of 2 cycles per lunar day (period: 12.42 hour) equals to twice the first Doodson argument and zero times all of the others [Doodson, 1921]. Kamil Teke 8 Figure 4. Hourly coordinate estimates of ONSALA60 (Sweden) station when the ocean tide model corrections are not introduced to the coordinates a priori to the adjustment are depicted in black (hourly estimates). The differences in hourly estimates with respect to FES2014 ocean tide model displacements are shown in orange. The median formal errors of the hourly coordinate estimates of ONSALA60 station are 5.1, 3.2 and 4.5 mm in radial, east and north components, respectively. Readers are referred to the supplementary material of this paper for the plots of the hourly OTL displacement estimates of all CONT14 stations. For the calculations of OTL displacements from ocean tide models according to Equation (2), the amplitudes and phase lags of principal tides were derived from “Ocean Tide Loading Provider”, a web-based facility provided by M.S. Bos and H.G. Scherneck (http://holt.oso.chalmers.se/loading, accessed on July 2017 for FES2012, FES2014 and TPXO8) that uses OLFG/OLMPP algorithm [Scherneck, 1991] and the Green’s functions using Gutenberg-Bullen standard Earth model to determine the deformation due to point loads [Farrell,  1972]. Tidal coefficients were selected as not to correct for geocenter motion due to ocean tides. For all calculations of OTL displacements from the tidal constituents, the sidelobes of the principal tides, caused by the nodal modulation of these tides with the 18.6-year period of the lunar node [McCarthy, 1996; Scherneck, 1999; Tamura 1987] were not introduced. Instead, fully interpolated loading tide spectrum with 342 tidal constituents was computed by spline interpolation of the real (in-phase) and imaginary (quadrature) components of principal tides (Hartmann and Wenzel 1995) as recommended by Petit  and  Luzum  [2010]. Matlab functions of VieVS based on a Fortran program, “hardisp.f” provided by Petit and Luzum [2010] were used. The CONT14 hourly VLBI station coordinate time series Δ𝑛,𝑘 were handled as observations to estimate the amplitudes (𝐴𝑛,𝑗,𝑘) and phase lags (𝜑𝑛,𝑗,𝑘) of the M2 constituent of OTL displacements using the linearized form of Equation (2) Δ𝑛,𝑘 =�𝑅𝑛,𝑗,𝑘 cos 𝜒𝑗(𝑡) + 𝐼𝑛,𝑗,𝑘 sin 𝜒𝑗(𝑡) (3) � where 𝑛 is the station index, 𝑅𝑛,𝑗,𝑘 = 𝐴𝑛,𝑗,𝑘 cos 𝜑𝑛,𝑗,𝑘 and 𝐼𝑛,𝑗,𝑘 = 𝐴𝑛,𝑗,𝑘 sin 𝜑𝑛,𝑗,𝑘 are the real (in-phase) and imaginary (quadrature) parts of the 𝑗’th tide (M2) and the 𝑘’th coordinate component. Note that no constraint was introduced on the estimated real and imaginary parts of the OTL displacements M2 constituent. The real and imaginary parts were then converted to an amplitude and Greenwich phase lag with the convention of phase lags taking negative values. The amplitudes, 𝐴𝑛,𝑗,𝑘 = �𝑅²𝑛 ,𝑗,𝑘 + 𝐼²𝑛 ,𝑗,𝑘 and the Greenwich phase lags, 𝜑𝑛,𝑗,𝑘 = arctan (𝐼𝑛,𝑗,𝑘 / 𝑅𝑛,𝑗,𝑘) of the M2 tide, were computed from the direct estimates of 𝑅𝑛,𝑗,𝑘 and 𝐼𝑛,𝑗,𝑘, respectively. The errors of the amplitudes and phase lags were calculated from the one sigma formal errors of the estimated real and imaginary parts of the M2 constituent using the law of propagation of variances. The displacements due to long-period tides i.e. Mf, Mm and Ssa as well as the semi-diurnal and diurnal tides other than M2 tide provided by the FES2014 model were reduced from the hourly station coordinates a priori to the estimation. The estimated amplitudes and phases of M2 tide, as well as their formal errors at the CONT14 VLBI sites, are presented in Table 2. Depending on the number and quality of the observations per station, the formal errors of the estimated amplitudes of M2 tidal constituent were found similar, ranging from 0.04 to 0.30 mm in all coordinate components except for the VLBI sites Hartebeesthoek and Fortaleza with formal errors of the amplitudes as large as about 0.40- 0.48 mm in radial and 0.20-0.35 mm in horizontal components. The formal errors of the phase angles of M2 tide are found mostly on the order of 0.4-9.6 degrees except for the VLBI sites at Wettzell and Hartebeesthoek that reach up to 12.7 degrees (south component Wettzell) and 14.6 degrees (west component Hartebeesthoek). The reason for the larger formal errors of the south component phase of Wettzell and the west component phase of Hartebeesthoek is that the amplitudes of these components are small. The best agreement of the estimated radial amplitudes of M2 tide with respect to those of ocean tide models were found at the sites Onsala, Badary, Wettzell, Zelenchukskaya, Matera, Tsukuba, and Yarragade varying between 0.01 and 0.40  mm whereas the worst agreement is seen at Ny-Ålesund, Yebes, and Fortaleza sites with the radial amplitude differences of 1.48  mm with respect to FES2014, 0.93  mm and 1.50  mm with respect to TPXO8, respectively. At most of the stations, the estimated tangential amplitudes of M2 tide with respect to those of ocean tide models vary from 0.03 mm (at Matera) to 0.75 mm (at Fortaleza). It is worth to note that there is a large M2 tide radial amplitude difference of about 0.7 mm between the TPXO8 model and both FES2014 and FES2012 models at Ny-Ålesund site. When the Greenwich phase lags of M2 tide are considered across the CONT14 sites, the agreement between the estimates and the ocean tide models do agree within 1 to 5 degrees for most of the stations in all coordinate components. However, the largest M2 tide phase difference is found between the CONT14 estimates and ocean tide models at Hartebeesthoek site in west component with the value of about 27 degrees. To see the bar plots of the amplitudes and the phases of M2 constituent OTL displacements derived from the ocean tide models with respect to the estimates of the CONT14 campaign at all VLBI sites and in all coordinate components readers are referred to the supplementary material of this paper. 9 Ocean loading M2 tide displacements from VLBI http://holt.oso.chalmers.se/loading From the phasor vector plots of the M2 tidal constituent of OTL displacements (see Figure 5 and the supplementary material of this paper) it can be seen that the agreement of the amplitudes and phases between CONT14 estimates and ocean tide models do agree well for most of the sites and coordinate components except Fortaleza. The amplitudes of the vector differences of the phasor vectors were used as a metric for assessing the level of agreement of the estimated amplitudes and phases of M2 tidal constituent derived from CONT14 observations to those of the ocean tide models i.e. FES2014, FES2012, and TPXO8. The amplitudes of the vector differences were calculated from ⃒ 𝑑𝑛,𝑗,𝑘 ⃒ = ⃒ 𝐴𝑉𝐿𝐵𝐼(cos𝜑𝑉𝐿𝐵𝐼+𝑖 sin𝜑𝑉𝐿𝐵𝐼)𝑛,𝑗,𝑘−𝐴𝑀𝑂𝐷𝐸𝐿(cos𝜑𝑀𝑂𝐷𝐸𝐿 +𝑖 sin𝜑𝑀𝑂𝐷𝐸𝐿 )𝑛,𝑗,𝑘 ⃒ (4) where 𝑗 denotes to the tidal constituent i.e. M2, 𝑘 the coordinate component (in radial, west, and south directions) at the station 𝑛. The amplitudes of M2 tide phasor vector differences between FES2014, FES2012, TPXO8 ocean tide models and those estimated from the CONT14 campaign are shown in Figure 6. For most of the sites, the amplitudes of the vector differences in radial components of the M2 tide are up to two times larger than those of the tangential components. The amplitudes of the M2 tide phasor vector differences in all coordinate components are below 0.5 mm at most of the inland sites i.e. Badary, Wettzell, Zelenchukskaya, Kamil Teke 10 Stations Radial West South A M2 [mm] φ M2 [°] A M2 [mm] φ M2 [°] A M2 [mm] φ M2 [°] NYALES20 9.26±0.13 [7.78] 175.7±0.8 [178.5] 2.84±0.05 [2.57] -19.5±1.1 [-12.5] 1.52±0.05 [1.34] -27.5±1.9 [-18.1] ONSALA60 3.28±0.06 [3.67] -63.1±1.1 [-63.9] 1.47±0.04 [1.49] 89.8±1.7 [85.0] 0.67±0.05 [0.58] 115.8±4.1 [111.3] BADARY 0.86±0.08 [0.85] 85.5±5.5 [78.0] 0.63±0.06 [0.49] 149.6±5.9 [156.2] 0.42±0.07 [0.27] 81.2±9.6 [87.6] WETTZELL 5.44±0.06 [5.14] -72.20±0.7 [-71.3] 2.06±0.04 [1.91] 81.30±1.2 [77.6] 0.23±0.05 [0.39] -55.30±12.7 [-49.8] ZELENCHK 2.81±0.10 [2.64] -68.8±1.9 [-64.9] 0.77±0.08 [0.76] -2.7±5.9 [-3.0] 0.70±0.08 [0.74] -136.2±6.2 [-122.9] WESTFORD 7.65±0.15 [7.24] -167.9±1.1 [-167.4] 3.59±0.08 [3.40] -130.3±1.3 [-129.4] 1.85±0.09 [1.73] -22.5±3.0 [-15.8] MATERA 5.76±0.08 [5.41] -79.8±0.8 [-78.5] 1.24±0.06 [1.24] 72.3±2.7 [62.8] 0.93±0.07 [0.96] -81.2±4.4 [-78.6] YEBES40M 14.19±0.10 [13.35] -87.3±0.4 [-87.4] 4.22±0.08 [3.91] 70.5±0.8 [69.0] 3.07±0.06 [2.88] -64.2±1.1 [-62.2] TSUKUB32 7.66±0.15 [7.57] 49.7±1.1 [50.2] 2.57±0.08 [2.35] -14.0±1.8 [-9.5] 1.90±0.06 [1.95] -80.4±1.9 [-76.1] KOKEE 11.97±0.29 [12.6] -120.6±1.4 [-121.3] 2.97±0.18 [2.49] 150.8±3.5 [159.0] 4.44±0.14 [4.09] 96.5±1.7 [97.3] FORTLEZA 36.88±0.48 [35.36] 32.0±0.8 [34.0] 5.20±0.35 [4.45] 27.9±3.9 [20.7] 5.77±0.33 [5.15] 49.4±3.3 [52.1] KATH12M 2.57±0.15 [3.02] 111.8±3.3 [112.6] 1.97±0.11 [1.60] 85.3±3.1 [87.5] 1.15±0.10 [1.19] 44.4±5.0 [47.3] HART15M 17.22±0.40 [16.66] -129.6±1.3 [-131.9] 0.78±0.20 [0.58] 69.2±14.6 [43.1] 1.84±0.26 [1.39] 79.7±8.1 [69.6] YARRA12M 3.61±0.19 [3.48] 143.9±3.0 [146.4] 2.28±0.13 [2.24] -129.8±3.3 [-128.7] 2.52±0.14 [2.12] -161.20±3.3 [-154.1] WARK12M 26.41±0.24 [25.63] 57.9±0.5 [56.6] 8.58±0.14 [8.27] -21.8±0.9 [-20.8] 5.70±0.20 [5.01] 14.6±2.0 [20.1] HOBART26 10.38±0.18 [9.57] 150.7±1.0 [152.7] 4.26±0.13 [3.97] 104.7±1.7 [103.4] 1.10±0.16 [1.21] 107.3±8.2 [103.5] HOBART12 10.20±0.18 [9.56] 151.30±1.0 [152.7] 4.36±0.14 [3.97] 101.3±1.8 [103.4] 1.06±0.17 [1.21] 109.8±9.3 [103.5] Table 2. The amplitudes and phase lags of OTL displacements M2 constituent with their formal errors estimated from CONT14 hourly sessions at all contributing stations in millimeter and arc degrees, respectively. The FES2014 ocean tide model values are written in brackets. Matera, and Katherine as well as at the coastal sites Onsala and Westford. At Ny-Ålesund, Yebes, Yarragadee, and Hobart (HOBART12 and HOBART26), the tangential amplitudes of vector differences are below 0.6 mm. However, at these sites, except Yarragadee, the radial amplitudes of vector differences are larger relative to tangential components and range from 0.7 to 1.6 mm. The worst agreement between CONT14 estimates and the ocean tide models in terms of M2 tide phasor vector differences is found at the site Fortaleza varying in all coordinate components between 0.7 mm to 1.6 mm (see Figure 6). 11 Ocean loading M2 tide displacements from VLBI Figure 5. Phasor vectors of OTL displacements M2 constituent with one sigma error ellipses at Wettzell and Matera VLBI sites. The formal errors of the estimated amplitudes and phase-lags are written to the top of each subplot, respectively. The horizontal axis represents the in-phase component, and the vertical axis out-of-phase component relative to the tidal potential at Greenwich. The Greenwich phase lag is zero along the positive direction in the horizontal axis and the phase angle increases counterclockwise. Figure 6. The amplitudes of the M2 tide phasor vector differences between FES2014, FES2012, TPXO8 ocean tide models and those estimated from the CONT14 campaign. 4. The impact of introducing different ocean tide models on the estimated OTL displacements and the M2 tide when analyzing CONT14 campaign The errors of the geophysical models, e.g. Earth tides, pole tide, ocean loading, atmosphere loading, and hydrological loading, introduced a priori to the station coordinates, propagate into the residual estimates of station positions in the estimation process. To quantify the effect of ocean tide model errors on the estimated hourly positions, besides FES2014 (see Case 1, in Table 3), the OTL displacements calculated from another ocean tide model, GOT00.2 [Ray, 1999] (Case 2, in Table 3) is introduced to the station positions a priori to the analysis of hourly sessions. When the a priori ocean tide model is switched between FES2014 and GOT00.2, the level of change of the estimated OTL displacements of the corresponding site is tested. A relatively older and less accurate model, GOT00.2 is selected to get more significant differences to FES2014. It is worth to note that for both of the hourly solutions (Case 1 and Case 2), the reduced external troposphere delays estimated from the 24 hour sessions are the same and FES2004 ocean tide model [Lyard et al., 2006] is used for predicting the OTL displacements in both of the analyses of 24 hour sessions. Using the same observations and parameterization in the analysis of 1 hour sessions and only changing the a priori ocean tide model results in estimating the same position errors in magnitudes and directions. Thus, as expected, all the errors in the estimated OTL displacements are removed after having the differences, except those effects caused by the discrepancies of the ocean tide models (see e.g. Figure 7). The differences of the OTL displacement estimates reveal a tide-like harmonic behavior (see Figure 7). The peak- to-peak amplitudes of the differences of OTL displacements range from 0.1 to 1 mm for all coordinate components of all stations while the medians of formal errors of the estimated displacements for all coordinate components are varying in 5 to 10 mm. This finding suggests that in the solutions of hourly VLBI sessions, the effects of changing a priori ocean tide model on the station coordinate estimates would not exceed 1  mm which is far below the magnitude formal errors of the hourly OTL displacement estimates. To quantify the ocean tide model errors on the estimated amplitudes and phases of M2 tide, amplitudes of phasor vector differences are used as a metric. Different ocean tide models introduced in various stages of the analysis, are shown in Table 3 and explained in the below items. • To quantify the effect of using different a priori ocean tide models i.e. FES2014 or GOT00.2 in the analysis of 1 hour sessions the amplitudes of phasor vector differences of M2 tide between Case 1 and Case 2 are calculated (see Table 3 and Figure 9). • To get only the M2 tide OTL displacements, in addition to the FES2014, the GOT00.2 displacements, except M2, are reduced from the estimated hourly position time series. If the FES2014 or the GOT00.2 models would contain errors, these tides would not be completely reduced at the end and would affect the estimated M2 tide amplitudes and phases especially due to the short time period of 15 days. The effect of reducing FES2014 or GOT00.2 tidal displacements (except M2) from the hourly position series, on the estimated M2 tide, is quantified through producing the amplitudes of the phasor vector differences between Case 1 and Case 3 (see Table 3 and Figure 9). Kamil Teke 12 Table 3. Introduction of different ocean tide models into various stages of the analysis. For correcting the a priori station positions in the analysis of 24 hour sessions For correcting the a priori station positions in the analysis of 1 hour sessions For reducing all the other tides, except M2, from the estimated OTL hourly displacements Case 1 FES2004 FES2014 FES2014 Case 2 FES2004 GOT00.2 FES2014 Case 3 FES2004 FES2014 GOT00.2 Case 4 FES2014 FES2014 FES2014 • To quantify the effect of using different a priori ocean tide models i.e. GOT00.2 or FES2014 both in the analysis of 1  hour sessions and in the reduction of all the other tides, except M2, from the estimated OTL hourly displacements, the amplitudes of phasor vector differences between Case 2 and Case 3 are calculated (see Table 3 and Figure 9). • When FES2014 or FES2004 ocean tide model displacements are introduced to the a priori station coordinates in the analysis of 24 hour sessions, two different sets of troposphere delays are estimated (e.g. see Figure 8 and readers are referred to the supplementary material for the plots of the ZWD differences at the other VLBI sites when different a priori ocean tide models are used in the analysis of 24  hour sessions). The level of influence of these two types of troposphere delays estimated from 24 hour sessions and reduced from the observations of 1 hour sessions, on the estimated M2 tide of 1 hour sessions is evaluated. Thus, the estimated M2 tides are compared in terms of the amplitudes phasor vector differences between Case  1 and Case  4 (see Table 3 and Figure 9). The smallest variations in the amplitudes of the phasor vector differences of M2 tide are found between the results of Case 1 and Case 3 solutions. For most of the sites and coordinate components the amplitudes of the phasor vector differences between Case 1 and Case 3 are not larger than 0.1 mm, except the radial components of 13 Ocean loading M2 tide displacements from VLBI Figure 7. The differences of OTL displacements in radial, east, and north coordinate components at Kokee estimated from 1 hour sessions of CONT14 when the a priori ocean tide model is switched between FES2014 and GOT00.2 (Case 1 – Case 2). Readers are referred to the supplementary material for the differences of OTL displacements at the other sites. Figure 8. The ZWD differences estimated from the analysis of 24  hour sessions at Ny  Alesund when a priori station coordinates are corrected with the ocean tide models; FES2014, FES2004 or GOT00.2. Ny  Alesund, Warkworth, and Hobart (Hobart26) VLBI sites (see the red bars in Figure  9). The phasor vector differences for most of the sites and coordinate components between Case 1 and Case 2 range within 0.1-0.2 mm. The phasor vector differences between Case 2 and Case 3 are a bit larger than those of Case 1 and Case 2. As an inference from these findings, when 15 days long, 1 hour VLBI sessions are considered, selection of ocean tide model used for correcting the station positions a priori to the analysis of hourly sessions as well as for reducing the other tides (except M2) from the position time series would make a difference of a few sub-millimeters, i.e. about 0.1 to 0.3 mm in the differences of M2 tide phasor vectors. The largest differences of M2 tide phasor vectors are found between the solutions of Case 1 and Case 4 in the radial components of the coastal VLBI sites. The amplitudes of the differences of M2 tide radial phasor vectors of these sites are found as 1.1 mm at Ny Alesund, 0.6 mm at Onsala, 0.5 mm at Kokee, 0.4 mm at Fortaleza, 0.7 mm at Katherine (assumed to be an inland site, not too far away from sea about 230 km), 1.4 mm at Warkworth and 0.3 mm at Hobart (Hobart26). In the analysis of 24 hour sessions, the station displacements due to the change of a priori ocean tide models, i.e. FES2004 or FES2014, cause the estimated troposphere delays i.e. zenith wet delays (ZWD) to vary in -2 and 2 mm (see e.g. Figure 8 for the station Ny Alesund and the supplementary material for the other stations). Root-mean-square of the estimated ZWD differences (purple horizontal lines in Figure 9) and the M2 tide phasor vector differences in radial components between Case 1 and Case 4 reveal a statistically significant strong positive correlation of 0.89. These findings suggest that the larger errors of ocean tide models (e.g. FES2014 and FES2004) at coastal sites relative to the inland sites, propagate into the troposphere delays estimated from 24 hour sessions, which in turn these errors in the troposphere delays cause larger differences in the M2 tide phasor vectors of OTL displacements estimated from the 1 hour sessions. Kamil Teke 14 Figure 9. The amplitudes of the estimated M2 tide phasor vector differences from the CONT14 campaign between four different cases of parameterization (see Table 3). The effect of using different a priori ocean tide models i.e. FES2014 or GOT00.2 in the analysis of 1 hour sessions are plotted as black bars (Case 1 - Case 2). The effect of reducing FES2014 or GOT00.2 tides, except M2, from the estimated OTL hourly displacements are shown in red bars (Case 1 - Case 3). The effect of using different a priori ocean tide models i.e. FES2014 or GOT00.2 both in the analysis of 1 hour sessions and in the reduction of all the other tides, except M2, from the estimated OTL hourly displacements are shown in grey bars (Case 2 -Case 3). The effect reducing the two types of troposphere delays, when FES2014 or FES2004 is introduced a priori in the analysis of 24h sessions, from the observations of 1  hour sessions is depicted as orange bars (Case  1 -  Case  4). Root-mean-square (RMS) of the differences of ZWD estimates from the analysis of 24 hour sessions when a priori ocean tide model is used as FES2014 or FES2004 are plotted as purple horizontal lines. To assess the overall agreement between CONT14 estimates and ocean tide models the root-mean-square (RMS) of the amplitudes of the vector differences between observed and modeled phasor vectors across the stations for the 𝑗’th tidal constituent (M2) and the 𝑘’th coordinate component were calculated using (5) and are given in Table 4. Hereupon, the term “RMS misfits” refers to the RMS of the amplitudes of phasor vector differences between the CONT14 estimates and those of ocean tide models i.e. FES2014, FES2012 and TPXO8 across all the stations. The RMS misfits of the M2 tide in the radial component between CONT14 and ocean tide models over coastal sites are found about 0.70-0.96 mm both from the solutions of Case 1 and Case 4 which is significantly larger than those of inland sites i.e. varying in 0.34-0.51 mm. The RMS misfits in both radial and tangential components over coastal sites are found about two times larger than those of inland sites which is valid for both solutions of Case 1 and Case 4 (see Table 4). When all sites are considered, the RMS misfits of the M2 tide between CONT14 and ocean tide models are within 0.66-0.81 mm in radial and 0.21-0.40 mm in tangential components. The best agreement between CONT14 and the selected ocean tide models for M2 tide is seen in all the components of inland sites estimated from the Case 4 solution (see Table 3) with the RMS misfits ranging in 0.10-0.41 mm. 5. Conclusions In this study, hourly VLBI station coordinate time series derived from the analysis of 1 hour VLBI sessions of 15 days long continuous VLBI campaign, CONT14 were used to estimate the principal lunar semidiurnal tide, M2 of OTL displacements. Only M2 tide is estimated in this study and the station displacements caused by the remaining; long period, semi-diurnal and diurnal tides were calculated from FES2014 model and reduced from the hourly station coordinates a priori to the estimation. Intrinsic to this study, the troposphere delays estimated from 24-hour sessions were reduced from the observations of hourly sessions a priori to the analysis. 𝑅𝑀𝑆𝑗,𝑘 =� � ⃒ 𝑑𝑛,𝑗,𝑘 ⃒ 21𝑁 𝑚𝑖𝑠𝑓𝑖𝑡𝑠 � �⁼¹ 15 Ocean loading M2 tide displacements from VLBI Table 4. RMS misfits of M2 tide calculated from Equation (5) in mm between CONT14 phasor vectors from the solutions of Case 1, Case 4 and those of ocean tide models (i.e. FES2014, FES2012, and TPXO8) over coastal (10 stations, see Table 1), inland (7 stations otherwise, see Table 1) and all stations. Results from the Case 4 solution is written in italics. Case 1 Case 4 Radial West South Radial West South coast CONT14-FES2014 0.96 0.45 0.42 0.84 0.37 0.27 CONT14-FES2012 0.92 0.45 0.41 0.70 0.32 0.22 CONT14-TPXO8 0.89 0.46 0.40 0.89 0.39 0.25 inland CONT14-FES2014 0.45 0.24 0.15 0.36 0.20 0.10 CONT14-FES2012 0.46 0.23 0.15 0.34 0.19 0.10 CONT14-TPXO8 0.51 0.25 0.16 0.41 0.21 0.12 all CONT14-FES2014 0.81 0.39 0.35 0.70 0.32 0.22 CONT14-FES2012 0.78 0.38 0.34 0.66 0.32 0.21 CONT14-TPXO8 0.77 0.40 0.33 0.75 0.33 0.21 The estimated amplitudes and phase lags of the M2 constituent of OTL displacements at each site are compared with the predictions of the considered ocean tide models, i.e., FES2012 [Lyard et al., 2006; Carrère et al., 2012], FES2014 [Carrère et al., 2016], and TPXO8 [Egbert and Erofeeva, 2002; Egbert et al., 2010]. The amplitudes of M2 tide phasor vector differences in all coordinate components between CONT14 estimates and those of ocean tide models are found below 0.5 mm at the inland sites i.e. Badary, Wettzell, Zelenchukskaya, Matera, and Katherine as well as at the coastal sites Onsala and Westford. The worst agreement between CONT14 estimates and ocean tide models for M2 tide in terms of phasor vector differences is found at the site Fortaleza where the values in all coordinate components are between 0.7 mm to 1.6 mm. This is due to the extremely humid weather at Fortaleza which results in the unmodelled troposphere delays propagate into coordinate components in the least-squares adjustment. Besides, the Fortaleza VLBI site is located at about 5 km far away from the Atlantic ocean and its M2 amplitude of OTL displacements is about 36 mm in the radial component, the largest radial amplitude relative to the other CONT14 VLBI sites. As an overall assessment, the RMS of the amplitudes of M2 tide phasor vector differences between the CONT14 estimates and those of ocean tide models i.e. FES2014, FES2012, and TPXO8 across all the stations, namely RMS misfits are considered. The best agreement between CONT14 and the selected ocean tide models for M2 tide is seen in all the components of inland sites with the RMS misfits ranging in 0.1-0.4 mm. The RMS misfits of the M2 tide of OTL displacements in all coordinate components between CONT14 and ocean tide models over coastal sites are found about two times larger than those of inland sites. This result confirms the modeling insufficiencies in shallow waters of ocean tide models which results in an accuracy restriction of OTL displacements around coastal regions. The larger errors of ocean tide models (e.g. FES2014 and FES2004) at coastal sites relative to the inland sites, propagate into the troposphere delays estimated from 24 hour sessions, which in turn these errors in the troposphere delays cause larger differences in the M2 tide phasor vectors of OTL displacements estimated from the 1 hour sessions. The larger discrepancies between the VLBI observed and model-predicted M2 tides at coastal sites may not only be resulted from the errors of ocean tide models at coastal regions but also additional errors can become from VLBI data and analysis. Besides, hydrological loading and non-tidal ocean loading, which are not introduced a priori to the station coordinates in the VLBI analysis of 1 hour sessions in this study, would contribute to the increase in the misfits of phasor vectors between ocean tide models and VLBI at both coastal and inland sites. Because, even these geophysical effects have periods much longer than one day, given the short time span of CONT14 their contributions may not be negligible. To conclude, the level of agreement between the VLBI observed and model-derived OTL displacements can be significantly increased through reducing the 24-hour session troposphere delays from the observations of hourly sessions. As an outlook to improve the estimated accuracies of semidiurnal and diurnal tidal constituents of OTL displacements from the analysis of the observations of VLBI, further developments should be introduced to the present-day troposphere delay models as well as to the geophysical models for the station position displacements. Besides, the number of VLBI observations should be increased. Their temporal distribution over a session and sky coverage should be homogenized. Acknowledgments. This work is supported by The Scientific and Technological Research Council of Turkey (TÜBİTAK), project number: 115Y244. The author acknowledges the International VLBI Service for Geodesy and Astrometry (IVS) for providing the observations of the CONT14 campaign. The author is grateful to the anonymous reviewers for their comments which were helpful for improving the paper. References Agnew, D. C. (1996). SPOTL: some programs for ocean-tide loading, SIO Reference Series 96-8, Scripps Institution of Oceanography, La Jolla, California, 35. Agnew, D. C. (1997). NLOADF: A program for computing ocean-tide loading, J. Geophys. Res., 102(B3), 5109-5110. Allinson, C. R., P. J. Clarke, S. J. Edwards, M. A. King, T. F. Baker and P.R. Cruddace (2004). Stability of direct GPS estimates of ocean tide loading, Geophys. Res. Lett., 31, L15603. Baker, T. F., D. J. Curtis and A. H. Dodson (1995). Ocean tide loading and GPS, GPS world, vol March, 54–59. Kamil Teke 16 Baker, T. F. and M. S. Bos (2003). Validating Earth and ocean tide models using tidal gravity measurements, Geophys. J. Int., 152(2), 468-485. Behrend, D. (2013). Data Handling within the International VLBI Service, Data Science Journal, Vol. 12, WDS81–WDS84. Bizouard, C. and D. Gambis (2009). The combined solution C04 for Earth orientation parameters consistent with International Terrestrial Reference Frame, in: Drewes, H., (Ed.), Geodetic Reference Frames, IAG Symp, 134, 265–270. Bos, M. S. and T. F. Baker (2005). An estimate of the errors in gravity ocean tide loading computations, J. Geod., 79(1), 50-63. Böhm, J., B. Werl and H. Schuh (2006). Troposphere mapping functions for GPS and Very Long Baseline Interferometry from European Centre for Medium-Range Weather Forecasts operational analysis data, J. Geophys. Res., 111, B02406. Böhm, J., S. Böhm, J. Boisits, A. Girdiuk, J. Gruber, A. Hellerschmied, H. Krasna, D. Landskron, M. Madzak, D. Mayer, J. McCallum, L. McCallum, M. Schartner and K. Teke (2018). Vienna VLBI and Satellite Software (VieVS) for Geodesy and Astrometry, Publications of the Astronomical Society of the Pacific, 130(986), 044503, 1 – 6. Carrère, L., F. Lyard, M. Cancet, A. Guillot and L. Roblou (2012). A new global tidal model taking advantage of nearly 20 years of altimetry, in: Proceedings of Meeting “20 Years of Altimetry”, Venice. Carrère, L., F. Lyard, M. Cancet, A. Guillot, N. Picot (2016). Finite Element Solution FES2014, a new tidal model - Validation results and perspectives for improvements, presentation to ESA Living Planet Conference, Prague. Chen, G. and T. A. Herring (1997). Effects of atmospheric azimuthal asymmetry on the analysis from space geodetic data, J. Geophys. Res., 102(B9), 20489-20502. Dach, R. and R. Dietrich (2001). The Ocean Loading Effect in the GPS Analysis: A Case Study in the Antarctic Peninsula Region, Marine Geodesy, 24:1, 13-25. Davis, J. L., T. A. Herring, I. I. Shapiro, A. E. E. Rogers and G. Elgered (1985). Geodesy by radio interferometry: effects of atmospheric modeling errors on estimates of baseline length, Radio Sci., 20 (6), 1593–1607. Davis, J. L., G. Elgered, A. E. Niell and C. E. Kuehn (1993). Ground-based measurements of gradients in the “wet” radio refractivity of air, Radio Sci., 28 (6), 1003–1018. Doodson, A. T. (1921). The Harmonic Development of the Tide-Generating Potential, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 100 (704): 305-329. Dragert, H., T. S. James and A. Lambert (2000). Ocean loading corrections for continuous GPS: A case study at the Canadian coastal site Holberg, Geophys. Res. Lett., 27(14), 2045– 2048. Eanes, R. J. and A. Shuler (1999). An improved global ocean tide model from TOPEX/Poseidon altimetry: CSR4.0, in EuropeanGeophysical Society, 24th General Assembly, The Hague. Egbert, G. D. and S. Y. Erofeeva (2002). Efficient inverse modeling of barotropic ocean tides, Journal of Atmospheric and Oceanic Technology, 19.2 (2002), 183-204. Egbert, G. D., S. Y. Erofeeva and R. D. Ray (2010). Assimilation of altimetry data for nonlinear shallow-water tides: Quarter-diurnal tides of the Northwest European Shelf, Cont. Shelf Res., 30(6), 668-679. Farrell, W. E. (1972). Deformation of the Earth by surface loads, Rev. Geophys. Space Phys., 10(3), 761-797. Farrell, W. E. (1973). Earth tides, ocean tides and tidal loading, Philosophical Transactions of the Royal Society of London, Series A. 274, 253-259. Fey, A., D. Gordon and C.S. Jacobs (2009). The second realization of the International Celestial Reference Frame by Very Long Baseline Interferometry, IERS Technical Note, vol. 35. Verlag des Bundesamt für Kartographie und Geodäsie, Frankfurt am Main. Foreman, M. (1977). Manual for tidal heights analysis and prediction, Pac. Mar. Sci. Rep., 77–10, 1–58. Haas, R. and H. Schuh (1998). Ocean loading observed by geodetic VLBI, in Proceedings of the 13th International Symposium on Earth Tides, edited by B. Ducarme and P. Paquet, 111-120. Hartmann, T. and H-G. Wenzel (1995). The HW95 tidal potential catalogue, Geophys. Res. Lett., 22(24), 3553-3556. Krásná, H., J. Böhm, L. Plank, T. Nilsson and H. Schuh (2014). Atmospheric effects on VLBI-derived terrestrial and celestial reference frames, Earth on the Edge: Science for a Sustainable Planet Proceedings of the IAG General Assembly, Melbourne, Australia, in Chris Rizos, Pascal Willis (Editors), Series: International Association of Geodesy Symposia, 139, 203–208. Khan, S. A. and C. C. Tscherning (2001). Determination of semi-diurnal ocean tide loading constituents using GPS in Alaska, Geophys. Res. Lett., 28(11), 2249-2252. 17 Ocean loading M2 tide displacements from VLBI http://publik.tuwien.ac.at/files/publik_268662.pdf http://publik.tuwien.ac.at/files/publik_268662.pdf http://publik.tuwien.ac.at/files/publik_268662.pdf Khan, S. A., H. G. Scherneck (2003). The M2 ocean tide loading wave in Alaska: Vertical and horizontal displacements, modelled and observed, J. Geod., 77(3-4), 117–127. King, M. A., N. T. Penna, P. J. Clarke and E. C. King (2005). Validation of ocean tide models around Antarctica using inshore GPS and gravity data, J. Geophys. Res., 110, B08401. King, M. A. (2006). Kinematic and static GPS techniques for estimating tidal displacements with application to Antarctica, J. Geodyn., 41(1-3), 77-86. Lambeck, K. (1988). Geophyscial geodesy-The slow deformations of the Earth, Oxford University Press, Oxford. Lyard. F., F. Lefevre, T. Letellier and O. Francis (2006). Modelling the global ocean tides: Modern insights from FES2004, Ocean. Dyn., 56(5-6), 394-415. Martens, H. R., M. Simons, S. Owen and L. Rivera (2016). Observations of ocean tidal response in South America from subdaily GPS positions, Geophys. J. Int., 205, 1637-1664. Matsumoto, K., T. Sato, T. Takanezawa and M. Ooe (2001). GOTIC2: A program for computation of oceanic tidal loading effect, J. Geod. Soc. Japan, 47, 243-248. McCarthy, D. D. (Editor) (1996). IERS Technical Note 21, IERS Conventions (1996), International Earth Rotation and Reference Systems Service (IERS), Verlag des Bundesamt für Kartographie und Geodäsie, Frankfurt am Main. Melachroinos, S. A., R. Biancale, M. Llubes, F. Perosanz, F. Lyard, M. Vergnolle, M. N. Bouin, F. Masson, J. Nicolas, L. Morel and S. Durand (2008). Ocean tide loading (OTL) displacements from global and local grids: comparisons to GPS estimates over the shelf Brittany, France, J. Geod., 82(6), 357-371. Nothnagel, A., T. Artz, D. Behrend, and Z. Malkin (2017). International VLBI Service for Geodesy and Astrometry - Delivering high-quality products and embarking on observations of the next generation, Journal of Geodesy, 91(7), 711–721. Pagiatakis, S. D. (1990). The response of a realistic earth to ocean tide loading, Geophys. J. Int., 103, 541-320. Pagiatakis, S. D. (1992). Program LOADSDP for the calculation of ocean load effects, Manusc. Geod., 17, 315-320. Penna, N. T., M. S. Bos, T. F. Baker and H-G. Scherneck (2008). Assessing the accuracy of predicted ocean tide loading displacement values, J. Geod., 82(12), 893-907. Penna, N. T., P. J. Clarke, M. S. Bos and T. F. Baker (2015). Ocean tide loading displacements in western Europe: 1. Validation of Kinematic GPS estimates, J. Geophys. Res. Solid Earth, 120, 6523-6539. Petit, G. and B. Luzum (2010). IERS Conventions 2010, IERS Technical Note; 36. Verlag des Bundesamt für Kartographie und Geodäsie, Frankfurt am Main. Petrov, L., C. P. Ma (2003). Study of harmonic site position variations determined by very long baseline interferometry, J. Geophys. Res., 108(B4), 2190. Petrov, L. and J. P. Boy (2004). Study of the atmospheric pressure loading signal in Very Long Baseline Interferometry observations, J. Geophys. Res., 109 (B3), B03405. Ray, R. D. (1999). A global ocean tide model from TOPEX/POSEIDON Altimetry: GOT99.2, NASA/TM-1999-209478, Greenbelt, 58. Rothacher, M. and G. Beutler (1998). The Role of GPS in the study of Global Change, Physics and Chemistry of the Earth, 23(9/10), 1029-1040. Saastamoinen, J. (1972). Atmospheric correction for the troposphere and stratosphere in radio ranging of satellites, in The use of artificial satellites for geodesy, Geophys. Monogr. Ser. 15, Amer. Geophys. Union, 274–251. Saastamoinen, J. (1973). Contribution to the theory of atmospheric refraction (in three parts). Bull. Geod., 105–107, 279–298. Schenewerk, M. S., J. Marshall and W. Dillinger (2001). Vertical ocean-loading deformations derived from a global GPS network, J. Geod. Soc. Japan, 47(1), 237–242. Scherneck, H-G. (1991). A parameterized solid earth tide model and ocean tide loading effects for global geodetic baseline measurements, Geophys. J. Int., 106(3), 677-694. Scherneck, H-G. (1993). Ocean tide loading: propagation of errors from the ocean tide into loading coefficients, Manuscripta Geodaetica, 18, 59-71. Scherneck, H-G. (1999). Explanatory supplement to the section “Local site displacement due to ocean loading” of the IERS Conventions 1996 Chapters 6 and 7. DGFI Report 71. Schuh H., ed. Verlag des Bundesamt für Kartographie und Geodäsie, Frankfurt am Main, 19-23. Scherneck, H-G., R. Haas and A. Laudati (2000). Ocean loading for, in and from VLBI, in IVS 2000 General Meeting Proceedings, Greenbelt, Vandenberg N. and K. Baver (Editors), 257-262. Kamil Teke 18 Scherneck, H-G and M. S. Bos (2002). Ocean Tide and Atmosphere Loading, IVS 2002 General Meeting Proceedings, 205-214. Schuh, H. and L. Möhlmann (1989). Ocean loading station displacements observed by VLBI, Geophys. Res. Lett., 16, 1105-1108. Schuh. H. and D. Behrend (2012). VLBI: A fascinating technique for geodesy and astrometry, Journal of Geodynamics, 61, 68–80. Schwiderski, E. W. (1980). On Charting Global Ocean Tides, Rev. Geophys. Space Phys., 18(1), 243-268. Shum, C. K., P. L. Woodworth, O. B. Andersen, G. D. Egbert, O. Francis, C. King, S.M. Klosko, C. Le Provost, X. Li, J- M. Molines, M. E. Parke, R. D. Ray, M. G. Schlax, D. Stammer, C. C. Tierney, P. Vincent, C. I. Wunsch (1997). Accuracy assessment of recent ocean tide models, J. Geophys. Res. 102(C11), 25173-25194. Sovers, O. J. (1994). Vertical ocean loading amplitudes from VLBI measurements, Geophys. Res. Lett., 21, 357-360. Tamura, Y. (1987). A harmonic development of the tide-generating potential, Bull. Inf. Mar. Terr., 99, 6813-6857. Teke, K., J. Böhm, T. Nilsson and H. Krasna (2013a). Sub-daily antenna position estimates from the CONT11 campaign, Proceedings of the 21st Meeting of the European VLBI Group for Geodesy and Astrometry (EVGA), Reports of the Finnish Geodetic Institute, N. Zubko and M. Poutanen (Editors), Helsinki, Finland, 131-134. Thomas, I. D., M. A. King M. A. and P. J. Clarke (2007). A comparison of GPS, VLBI and model estimates of ocean tide loading displacements, J. Geod., 81(5), 359-368. Vergnolle, M., M-N. Bouin, L. Morel, F. Masson, S. Durand, J. Nicolas and S. A. Melachroinos (2008). GPS estimates of ocean tide loading in NW-France: Determination of ocean tide loading constituents and comparison with a recent ocean tide model, Geophys. J. Int., 173(2), 444-458. Vey, S., E. Calais, M. Llubes, N. Florsch, G. Woppelmann, J. Hinderer, M. Amalvict, M. F. Lalancette, B. Simon, F. Duquenne and J. S. Haase (2002). GPS measurements of ocean loading and its impact on zenith tropospheric delay estimates: a case study in Brittany, France, J. Geod., 76, 419-427. Yuan, L. and B. F. Chao (2012). Analysis of tidal signals in surface displacement measured by a dense continuous GPS array, Earth planet. Sci. Lett., 355, 255-261. Yuan, L., B. F. Chao, X. Ding and P. Zhong (2013). The tidal displacement field at Earth’s surface determined using global GPS observations, J. Geophys. Res. Solid Earth, 118, 2618-2632. Yun, H-S., D. H. Lee and D. S. Song (2007). Determination of vertical displacements over the coastal area of Korea due to the ocean tide loading using GPS observations, J. Geodyn., 43(3-5), 528-541. Zumberge, J. F., M. B. Heflin, D. C. Jefferson, M. M. Watkins and F.H. Webb (1997). Precise point positioning for the efficient and robust analysis of GPS data from large networks, J. Geophys. Res., 102(B3), 5005-5017. *CORRESPONDING AUTHOR: Kamil T E K E, Hacettepe University, Department of Geomatics Engineering, Ankara, Turkey e-mail: kteke@hacettepe.edu.tr © 2020 the Istituto Nazionale di Geofisica e Vulcanologia. All rights reserved. 19 Ocean loading M2 tide displacements from VLBI