## 2 Data Processing
### GPS Data Processing  
To investigate the impact of MOT on the position of global IGS stations, we implement the above recommended conventional OTL correction method in GAMIT 10.4 by substituting the embedded OTL module with hardisp.f. Then a contrast experiment has been designed by reprocessing the GPS data of 109 evenly distributed global IGS stations using the original (Experiment A) and improved GAMIT (Experiment B) [21]. The difference between experiment A and B only comes from the OTL modeling. In experiment A, the OTL correction considered the effects of 11 main ocean tides, while in experiment B, we account for the MOT induced surface displacement by expanding the 11 main tides into 342 ocean tidal constituents. In this way, the contribution of MOT to daily GPS solution andcontinuous GPS time series of the global IGS stations could be determined by differencing between experiments A and B. Spatial distribution of the selected IGS stations is shown in Figure 1. To reduce time, only data on Wednesdays during the period of June, 1998 to December, 2010 are processed. The global grid of the amplitudes and phases for the 11 main tides with spatial resolution \(0.125^{\circ}\times 0.125^{\circ}\) is obtained from the FES2004 model [35].  
During the data reprocessing, satellite orbits, earth orientation parameters (EOPs), site coordinates together with the tropospheric delay and the horizontal gradient parameters are resolved simultaneously. Loose constraints are implemented on the stations, among which the constraints of IGS core stations are set as 5 cm, and non-core stations are 1 dm. Where possible, ambiguities are fixed to integers [36]. Satellite cut-off elevation angle is chosen as 10 degree, and site-specific, elevation-dependent weighting of the observations based on an assessment of the post phase residuals has been applied [27]. Corrections of solid Earth tide, ocean tide, pole tide as well as atmospheric tide have been implemented [6,37]. Vienna Mapping Function 1 (VMF1) troposphere mapping functions are used to calculate tropospheric delay [38,39]. Where possible, receiver-independent exchange (RINEX) meteorological files (.m) are used to provide pressure and temperature for the a priori zenith hydrostatic delay; otherwise, values from VMF1 are used [40]. Absolute antenna phase center offsets and variations are used (igs08_1636.atx) [41]. Impacts of second and third-order ionospheric delay are considered, during which the International Geomagnetic Reference Field 11 (IGRF 11) is selected to calculate the second-order ionospheric delay [36]. Non-tidal loading corrections are not implemented according to the IERS Conventions 2010 [4]. Finally, gross error rejection and datum transformation are implemented to the daily baseline resolutions using GLOBK to obtain stations' coordinate time series and velocities under the ITRF08 [42]. During the datum transformation, only six parameters, that is three rotation and three translation parameters are estimated to reduce the aliasing effects of the unmodeled surface mass loadings, e.g., the non-tidal atmospheric loading [43].  
Figure 1: Selected Global IGS reference station distribution. Black texts indicate the name of stations that mentioned in this paper.