The Use of National CORS Networks for Determining Temporal Mass Variations within the Earth’s System and for Improving GRACE/GRACE-FO Solutions

Temporal mass variations within the Earth’s system can be detected on a regional/global scale using GRACE (Gravity Recovery and Climate Experiment) and GRACE Follow-On (GRACE-FO) satellite missions’ data, while GNSS (Global Navigation Satellite System) data can be used to detect those variations on a local scale. The aim of this study is to investigate the usefulness of national GNSS CORS (Continuously Operating Reference Stations) networks for the determination of those temporal mass variations and for improving GRACE/GRACE-FO solutions. The area of Poland was chosen as a study area. Temporal variations of equivalent water thickness ∆EWT and vertical deformations of the Earth’s surface ∆h were determined at the sites of the ASG-EUPOS (Active Geodetic Network of the European Position Determination System) CORS network using GRACE/GRACE-FO-based GGMs and GNSS data. Moreover, combined solutions of ∆EWT were developed by combining ∆EWT obtained from GNSS data with the corresponding ones determined from GRACE satellite mission data. Strong correlations (correlation coefficients ranging from 0.6 to 0.9) between detrended ∆h determined from GRACE/GRACE-FO satellite mission data and the corresponding ones from GNSS data were observed at 93% of the GNSS stations investigated. Furthermore, for the determination of temporal mass variations, GNSS data from CORS network stations provide valuable information complementary to GRACE satellite mission data.


Introduction
The determination of temporal mass variations in the Earth's system with high accuracy as well as high spatial and temporal resolutions using space geodetic data is one of the main scientific problems in the Earth science-related disciplines. Since the last decade, the Gravity Recovery and Climate Experiment (GRACE; [1]) satellite mission operated between March 2002 and October 2017, brought a unique opportunity for the determination of temporal mass variations within the Earth's system. Its concept, including theoretical background as well as impressive scientific results, has widely been demonstrated by many authors (e.g., see the review given by [2]). The tremendous success achieved from the mission emphasized the need of launching a GRACE-type mission for a sustainable long-term monitoring of temporal mass variations within the Earth's system. The GRACE Follow-On (GRACE-FO; https://gracefo.jpl.nasa.gov/) satellite mission of a designed life of five years has been launched in May 2018. The GRACE and GRACE-FO missions, without any doubt, are a state-of-the-art space geodetic technique for monitoring temporal mass variations within the Earth's system. They have Section 2. In Section 3, the methods implemented for the determination of vertical deformations of the Earth's surface ∆h and temporal variations of equivalent water thickness using GRACE/GRACE-FO satellite mission and GNSS data as well as the determination of ∆EWT combined solutions using these data are specified. The results obtained are presented and analyzed in Section 4. Moreover, in Section 4, the improvement of ∆EWT determined from GRACE satellite mission data when adding ∆EWT determined from GNSS data is illustrated. Finally, in Section 5, discussions and conclusions concerning the usefulness of GNSS data from CORS network stations for the determination of ∆h and ∆EWT as well as for improving GRACE/GRACE-FO solutions are given.

Study Area and Data Used
The area of Poland has been chosen as a study area. A national GNSS CORS network in this area, called the Active Geodetic Network of the European Position Determination System (ASG-EUPOS), consists recently of 102 GNSS stations apart from each other by about 70 km; http://www.asgeupos.pl, are operating since 2008. Daily GNSS data from 96 ASG-EUPOS stations covering the period 2008-2018 were used within the course of this study. Figure 1 illustrates the study area as well as the ASG-EUPOS sites used.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 22 mission and GNSS data as well as the determination of ΔEWT combined solutions using these data are specified. The results obtained are presented and analyzed in Section 4. Moreover, in Section 4, the improvement of ΔEWT determined from GRACE satellite mission data when adding ΔEWT determined from GNSS data is illustrated. Finally, in Section 5, discussions and conclusions concerning the usefulness of GNSS data from CORS network stations for the determination of Δh and ΔEWT as well as for improving GRACE/GRACE-FO solutions are given.

Study Area and Data Used
The area of Poland has been chosen as a study area. A national GNSS CORS network in this area, called the Active Geodetic Network of the European Position Determination System (ASG-EUPOS), consists recently of 102 GNSS stations apart from each other by about 70 km; http://www.asgeupos.pl, are operating since 2008. Daily GNSS data from 96 ASG-EUPOS stations covering the period 2008-2018 were used within the course of this study. Figure 1 illustrates the study area as well as the ASG-EUPOS sites used. Beside the GNSS data from ASG-EUPOS stations, CSR RL06 GRACE/GRACE-FO-based GGMs were utilized in this study. In order to reduce the noise included in RL06 GRACE/GRACE-FO-based GGMs, the DDK3 (decorrelation 3) filter (cf. [38]) that indicates a good performance over the area investigated [39], was applied. Moreover, CSR RL06 GRACE/GRACE-FO-based GGMs developed as monthly solutions obtained from GRACE/GRACE-FO satellite missions were truncated at d/o 60 that corresponds to the spatial resolution, i.e., 3° × 3° at the equator.
The WaterGAP Global Hydrology Model (WGHM) Version 2.2d (e.g., [40,41]) was used for the evaluation purposes. This hydrological model represents the terrestrial hydrological cycle with a spatial resolution of 0.5 degree. Within the context of the investigation conducted in this study, monthly variations of equivalent water thickness for the period 2008-2016 were obtained at the ASG-EUPOS sites from the WGHM. Beside the GNSS data from ASG-EUPOS stations, CSR RL06 GRACE/GRACE-FO-based GGMs were utilized in this study. In order to reduce the noise included in RL06 GRACE/GRACE-FO-based GGMs, the DDK3 (decorrelation 3) filter (cf. [38]) that indicates a good performance over the area investigated [39], was applied. Moreover, CSR RL06 GRACE/GRACE-FO-based GGMs developed as monthly solutions obtained from GRACE/GRACE-FO satellite missions were truncated at d/o 60 that corresponds to the spatial resolution, i.e., 3 • × 3 • at the equator.
The WaterGAP Global Hydrology Model (WGHM) Version 2.2d (e.g., [40,41]) was used for the evaluation purposes. This hydrological model represents the terrestrial hydrological cycle with a spatial resolution of 0.5 degree. Within the context of the investigation conducted in this study, monthly variations of equivalent water thickness for the period 2008-2016 were obtained at the ASG-EUPOS sites from the WGHM.

Methods
The methods implemented for the determination of monthly vertical deformations of the Earth's surface and monthly variations of equivalent water thickness from GRACE/GRACE-FO-based GGMs and GNSS data are described in Sections 3.1 and 3.2, respectively. In Section 3.3, the method applied to determine combined solutions of monthly variations of equivalent water thickness ∆EWT m-CombSol is presented. In Section 3.4, uncertainties of monthly variations of equivalent water thickness determined from GNSS and GRACE/GRACE-FO satellite missions' data as well as their combined solutions are estimated. The general steps of the method implemented within this investigation are summarized in Figure 2.

Methods
The methods implemented for the determination of monthly vertical deformations of the Earth's surface and monthly variations of equivalent water thickness from GRACE/GRACE-FO-based GGMs and GNSS data are described in Sections 3.1 and 3.2, respectively. In Section 3.3, the method applied to determine combined solutions of monthly variations of equivalent water thickness ∆EWTm-CombSol is presented. In Section 3.4, uncertainties of monthly variations of equivalent water thickness determined from GNSS and GRACE/GRACE-FO satellite missions' data as well as their combined solutions are estimated. The general steps of the method implemented within this investigation are summarized in Figure 2.

The Determination of ∆h m-GRACE and ∆EWT m-GRACE
Monthly vertical deformations of the Earth's surface ∆h m-GRACE and monthly variations of equivalent water thickness ∆EWT m-GRACE were determined at the sites of the ASG-EUPOS CORS network using GRACE/GRACE-FO-based GGMs as follows [9,20,42,43]: where ϕ, λ are spherical geocentric coordinates of the computation point, a is the radius of the Earth (semi-major axis), ρ w is the water density, ρ avg is the Earth's average density, k and h are load Love numbers of degree n based on the Preliminary Reference Earth Model (PREM; [44]) obtained from [45], ∆C nm and ∆S nm are the differences between fully normalized spherical harmonic coefficients from monthly RL06 GRACE/GRACE-FO-based GGMs and the corresponding ones obtained from a reference GGM, P nm are the fully normalized Legendre functions of degree n and order m, and N max is the maximum degree applied. The P nm (sinϕ) are determined using a recursive algorithm as follows (e.g., [46]): In this study, ∆h m-GRACE and ∆EWT m-GRACE were determined using CSR RL06 GRACE/GRACE-FO-based GGMs specified in Section 2 and the IGiK-TVGMF software [43]. These GGMs were obtained from the International Centre for Global Earth Models (ICGEM; cf. http://icgem.gfz-potsdam.de/series). The degree-1 and degree-2 spherical harmonic coefficients of these GGMs were replaced by the corresponding ones obtained from the solution described in [47] and from Satellite Laser Ranging (SLR) observations [48], respectively. The monthly estimates of these degree-1 and degree-2 spherical harmonic coefficients are available via GRACE TN-13 (GRACE Technical Note 13; https://podaac-tools.jpl.nasa.gov/drive/files/allData/ grace/docs/TN-13_GEOC_CSR_RL06.txt) and GRACE TN-11 (GRACE Technical Note 11; cf. https://podaac-tools.jpl.nasa.gov/drive/files/allData/ grace/docs/TN-11_C20_SLR.txt), respectively. The EGM2008 (Earth Gravitational Model 2008; [49]) was used as a reference GGM and the WGS84 (World Geodetic System 1984) was chosen as a geodetic reference system.

The Determination of ∆h m-GNSS and ∆EWT m-GNSS from GNSS Data
Daily vertical deformations of the Earth's surface ∆h d-GNSS for the period 2008-2018 were determined first at the location of 96 sites of the ASG-EUOPS CORS network using daily GNSS observations data processed by the GAMIT/GLOBK software version 10.7 [50,51]. The TEQC (Translation, Editing, and Quality Checking) software developed by UNAVCO [52] was used to check the quality of GNSS data from 96 stations of the ASG-EUPOS network. The linear combination of L1 and L2 was implemented to eliminate the signal delay caused by first-order ionospheric refraction. The delay induced by second-and third-order ionospheric refraction of the carrier wave was also removed using daily IONEX (IONosphere EXchange) files provided by the Centre for Orbit Determination in Europe (CODE; [53]). These IONEX files represent the ionosphere in the form of global maps of the Vertical Total Electron Content (VTEC). The troposphere delay was estimated using the Vienna Mapping Function (VMF; [54]) with a priori hydrostatic delay estimates based on the global pressure and temperature (GPT2) model [55]. The troposphere delay was estimated for every 2 h from GNSS data acquired at ASG-EUPOS sites investigated. The Atmospheric Loading Model (ATML; [56]) at the observation level provided by the MIT (Massachusetts Institute of Technology) was applied to remove deformation caused by atmospheric loading. The Finite Element Solutions 2004 tidal model (FES2004; [57]) was used to remove deformations induced from oceanic tidal loading.
From the GAMIT software, the daily loosely constrained solutions were obtained. These solutions were combined with global solutions provided by the MIT using the GLOBK software. In this step of processing, the GLORG was used to define the reference frame. Time series of daily vertical deformations of the Earth's surface were tight to the ITRF2014 (International Terrestrial Reference Frame 2014; [58]). Daily vertical deformations of the Earth's surface that exceed the daily repeatability within the same month were considered as outliers. These outliers were detected using GAMIT/GLOBK MATLAB TOOL [50] and subsequently removed. Finally, all daily vertical deformations of each month were averaged to estimate monthly vertical deformations of the Earth's surface ∆h m-GNSS . Time series of daily and monthly vertical deformations for an exemplary ASG-EUPOS station KLCE are shown in Figure 3. Then, ∆h m-GNSS were inverted into monthly variations of equivalent water thickness ∆EWT m-GNSS using the Green's function elastic half-space model (cf. [59]).
where R 0 presents the radius that has been estimated on the basis of the spatial resolution of GRACE/GRACE-FO satellite missions' data (i.e., 167 km), g denotes the gravity value (estimated to 9.81 Gal), υ is the Poisson's ratio, E is the Young's modulus. It should be noted that on the basis of the review of the literature (e.g., [59]), the Poisson's ratio and the Young's modulus used in this study were estimated to 0.25 and 70 Gpa, respectively.  [57]) was used to remove deformations induced from oceanic tidal loading. From the GAMIT software, the daily loosely constrained solutions were obtained. These solutions were combined with global solutions provided by the MIT using the GLOBK software. In this step of processing, the GLORG was used to define the reference frame. Time series of daily vertical deformations of the Earth's surface were tight to the ITRF2014 (International Terrestrial Reference Frame 2014; [58]). Daily vertical deformations of the Earth's surface that exceed the daily repeatability within the same month were considered as outliers. These outliers were detected using GAMIT/GLOBK MATLAB TOOL [50] and subsequently removed. Finally, all daily vertical deformations of each month were averaged to estimate monthly vertical deformations of the Earth's surface ∆hm-GNSS. Time series of daily and monthly vertical deformations for an exemplary ASG-EUPOS station KLCE are shown in Figure 3. Then, ∆hm-GNSS were inverted into monthly variations of equivalent water thickness ∆EWTm-GNSS using the Green's function elastic half-space model (cf. [59]).
where R0 presents the radius that has been estimated on the basis of the spatial resolution of GRACE/GRACE-FO satellite missions' data (i.e., 167 km), g denotes the gravity value (estimated to 9.81 Gal), υ is the Poisson's ratio, E is the Young's modulus. It should be noted that on the basis of the review of the literature (e.g., [59]), the Poisson's ratio and the Young's modulus used in this study were estimated to 0.25 and 70 Gpa, respectively.

The Detrmination of ∆EWT m-CombSol
Combined solutions of monthly variations of equivalent water thickness ∆EWT m-CombSol were developed for each (i = 1, 2, 3, . . . , q) GNSS station of the ASG-EUPOS CORS network: where q denotes the total number of GNSS stations of the ASG-EUPOS CORS network investigated and w GNSS and w GRACE are weights for ∆EWT determined from GNSS and GRACE/GRACE-FO satellite missions' data, respectively. These weights were estimated as follows: where (t = 1, 2, 3, . . . , k) denotes the month, k is the total number of months investigated, ∆EWT m-WGHM are time series of monthly variations of the equivalent water thickness obtained from the WGHM and the overbar operator presents the average over all k months.

Estimation of ∆EWT m-GRACE , ∆EWT m-GNSS , and ∆EWT m-CombSol Uncertainties
Monthly variations of equivalent water thickness determined from GNSS data and GRACE/GRACE-FO-based GGMs as well as their combined solutions determined with the use of methods specified in Section 3.1, Section 3.2, Section 3.3, were evaluated using ∆EWT m-WGHM . The differences were determined and standard deviations of these differences were computed. Moreover, correlations between equivalent water thickness obtained from the WGHM and the corresponding ones determined from GNSS and GRACE/GRACE-FO satellite missions' data as well as the combined solutions ∆EWT m-CombSol developed were estimated.

Monthly Vertical Deformations of the Earth's Surface
Time series of monthly vertical deformations of the Earth's surface ∆h at the stations of the ASG-EUPOS network obtained from GNSS data ∆h m-GNSS and ∆h m-GRACE determined from GRACE/GRACE-FO satellite missions' data are illustrated in Figure 4. The results presented in Figure 4 exhibit a seasonal pattern over the area investigated. It is observed that maximum values of ∆hm-GNSS and ∆hm-GRACE occur during the summer-fall seasons and minimum values during the winter-spring seasons [26]. This pattern can be attributed to the change in the mass loading caused by seasonal variations in the water mass with maximum values in March and minimum values in July-September [26,60]. For all GNSS stations investigated, the phase difference between ∆hm-GNSS and ∆hm-GRACE is less than one month, which can be considered negligible. The results presented in Figure 4 exhibit a seasonal pattern over the area investigated. It is observed that maximum values of ∆h m-GNSS and ∆h m-GRACE occur during the summer-fall seasons and minimum values during the winter-spring seasons [26]. This pattern can be attributed to the change in the mass loading caused by seasonal variations in the water mass with maximum values in March and minimum values in July-September [26,60]. For all GNSS stations investigated, the phase difference between ∆h m-GNSS and ∆h m-GRACE is less than one month, which can be considered negligible. Statistics of peak-to-peak variations of the seasonal pattern of the monthly vertical deformations are given in Table 1. The results presented in Table 1 indicate that although the medians of peak-to-peak variations of the seasonal pattern in ∆h m-GNSS and ∆h m-GRACE are in a good agreement, the dispersion and the standard deviation of those variations for ∆h m-GNSS with respect to the corresponding ones obtained from ∆h m-GRACE are higher by a factor of approx. 3.5, and approx. 2.5, respectively. The higher dispersion and standard deviation of peak-to-peak variations in ∆h m-GNSS may reveal that monthly vertical deformations obtained from GNSS data include the local deformation signal which cannot be detected from GRACE/GRACE-FO satellite missions' data. The sources of such local deformation signal may be (i) thermoelastic deformation of the monuments caused by seasonal variation of the Earth's temperature [61] and (ii) poroelastic deformations caused by large variation of the water table [62]. Besides these geophysical sources, systematic errors in GNSS observations such as mismodeling of daily and subdaily tidal signals [63], presence of GPS draconitic signal [37,64] as well as the variation in the phase center and the local multipath, can also induce the specific ∆h m-GNSS pattern for a station [65].
Furthermore, along with the seasonal pattern of ∆h obtained from GNSS and GRACE/GRACE-FO satellite missions data illustrated in Figure 4, there are also linear trends in the monthly vertical deformation of the Earth's surface ( Figure 5). Statistics of these linear trends are given in Table 2. It should be mentioned that due to limited GNSS data at SKSL, SKSK and SKSV stations of the ASG-EUPOS network, linear trends of ∆h have not been estimated at these sites. The results presented in Figure 5 and Table 2 indicate positive linear trends for ∆h m-GRACE at the ASG-EUPOS sites, ranging from 0.67 mm/year to 1.86 mm/year. These trends can be ascribed to the annual water depletion over the territory of Poland, observed beyond the extreme land hydrology events, i.e., flooding and increased precipitation, occurred in 2010-2011 (e.g., [66]). It should be noted that magnitudes of linear trends estimated from ∆h m-GRACE depend on both the covering period and number of months used. The ∆h m-GNSS reflects variable linear trends ranging from −1.47 mm/year to 2.16 mm/year. The reasons for these different linear trends can be attributed to the unstable monument [67], anthropogenic factors such as subsidence of mining areas and over exploitation of groundwater [68] and the presence of diverse geological structures such as Eastern European Precambrian platform, young western and middle Paleozoic platform as well as Carpathian region [69]. As the factors causing linear trends in ∆h m-GNSS and ∆h m-GRACE are different, the distinctive mismatch between linear trends of ∆h m-GNSS and ∆h m-GRACE for some ASG-EUPOS GNSS stations such as PRZM, SWKI occurs. The linear trends in ∆h m-GRACE in the case of ASG-EUPOS GNSS stations PRZM, SWKI can primarily be attributed to the manifestation of the elastic response of the water depletion over the area. On the other hand, the corresponding linear trends in ∆h m-GNSS can be ascribed to the monumental instability of the ASG-EUPOS GNSS station, anthropogenic activity and water mass change occurring in limited spatial scales, e.g., few kilometers, as well as the presence of geological features. The uncommon linear trends in ∆h m-GNSS and ∆h m-GRACE may result in disagreements between vertical deformations of the Earth's surface determined from GNSS and GRACE/GRACE-FO satellite missions' data. Time series of ∆h m-GNSS and ∆h m-GRACE were thus detrended for further analysis.  In order to compare monthly vertical deformations of the Earth's surface ∆h determined from GNSS and GRACE/GRACE-FO satellite missions' data, correlations and standard deviations of the differences between (1) ∆hm-GNSS and ∆hm-GRACE, and (2) detrended ∆hm-GNSS and detrended ∆hm-GRACE were computed. Figure 6 depicts coefficients of correlations between ∆h obtained from GNSS data and the corresponding ones obtained from GRACE/GRACE-FO satellite missions' data at the locations of the ASG-EUPOS stations investigated. Statistics of these correlation coefficients are given in Table 3. The results presented in Figure 6 and Table 3 reveal that at 68% of the ASG-EUPOS GNSS stations investigated, strong correlations (correlation coefficients ranging from 0.60 to 0.90) between ∆hm-GNSS and ∆hm-GRACE were obtained. Moderate/weak correlations (correlation coefficients ranging from 0.30 to 0.59) between ∆hm-GNSS and ∆hm-GRACE were observed at 30% of the stations investigated. At two ASG-EUPOS GNSS stations (PRZM and NWT1) correlations coefficients between ∆hm-GNSS and ∆hm-GRACE were found to be 0.13 and 0.29, respectively. After detrending ∆hm-GNSS and ∆hm-GRACE, strong correlations between ∆hm-GNSS and ∆hm-GRACE were obtained at 93% of the GNSS stations investigated. At 7% of the stations, moderate/weak correlations between ∆hm-GNSS and ∆hm-GRACE were obtained. Thus, values of correlations coefficients have clearly improved after detrending. For the ASG-EUPOS stations PRZM the value of the correlation coefficient increased to 0.71 after the removal of linear trends from ∆hm-GNSS and ∆hm-GRACE (cf. Figures 4 and 5). On the contrary the value of the correlation coefficient at the station NWT1 does not change much despite removing the linear trend. However, it should be noted that the value of the coefficient of correlation between ∆hm-GNSS and ∆hm-  In order to compare monthly vertical deformations of the Earth's surface ∆h determined from GNSS and GRACE/GRACE-FO satellite missions' data, correlations and standard deviations of the differences between (1) ∆h m-GNSS and ∆h m-GRACE , and (2) detrended ∆h m-GNSS and detrended ∆h m-GRACE were computed. Figure 6 depicts coefficients of correlations between ∆h obtained from GNSS data and the corresponding ones obtained from GRACE/GRACE-FO satellite missions' data at the locations of the ASG-EUPOS stations investigated. Statistics of these correlation coefficients are given in Table 3. The results presented in Figure 6 and Table 3 reveal that at 68% of the ASG-EUPOS GNSS stations investigated, strong correlations (correlation coefficients ranging from 0.60 to 0.90) between ∆h m-GNSS and ∆h m-GRACE were obtained. Moderate/weak correlations (correlation coefficients ranging from 0.30 to 0.59) between ∆h m-GNSS and ∆h m-GRACE were observed at 30% of the stations investigated. At two ASG-EUPOS GNSS stations (PRZM and NWT1) correlations coefficients between ∆h m-GNSS and ∆h m-GRACE were found to be 0.13 and 0.29, respectively. After detrending ∆h m-GNSS and ∆h m-GRACE , strong correlations between ∆h m-GNSS and ∆h m-GRACE were obtained at 93% of the GNSS stations investigated. At 7% of the stations, moderate/weak correlations between ∆h m-GNSS and ∆h m-GRACE were obtained. Thus, values of correlations coefficients have clearly improved after detrending. For the ASG-EUPOS stations PRZM the value of the correlation coefficient increased to 0.71 after the removal of linear trends from ∆h m-GNSS and ∆h m-GRACE (cf. Figures 4 and 5). On the contrary the value of the correlation coefficient at the station NWT1 does not change much despite removing the linear trend. However, it should be noted that the value of the coefficient of correlation between ∆h m-GNSS and ∆h m-GRACE in the case of NWT1 can be less reliable because of the shorter time span and gaps in ∆h time series. Overall, the results presented in Figure 6 show strong correlations between ∆h obtained from GNSS data and corresponding ones determined using GRACE/GRACE-FO data. Strong correlations between ∆h suggest that surface mass loading effects are the dominant contributor to the vertical deformation of the Earth's surface over the area investigated and also indicate consistency of the GNSS and GRACE/GRACE-FO satellite missions' data in monitoring surface mass loading effects. The cases of weak correlations between ∆h obtained from GNSS data and the corresponding ones determined using GRACE/GRACE-FO satellite missions' data may be attributed to variety of reasons such as local temporal variations of water masses and GNSS station dependent errors. However, in this study the values of correlation coefficients got improved after removing the linear trends and there are only two ASG-EUPOS stations (NWT1 and WIEL) for which values remain under 0.5. The case of NWT1 has already been stated above, and regarding WIEL the authors believe that the value of the correlation coefficient may get improved with the use of the longer data span.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 22 GRACE in the case of NWT1 can be less reliable because of the shorter time span and gaps in ∆h time series. Overall, the results presented in Figure 6 show strong correlations between ∆h obtained from GNSS data and corresponding ones determined using GRACE/GRACE-FO data. Strong correlations between ∆h suggest that surface mass loading effects are the dominant contributor to the vertical deformation of the Earth's surface over the area investigated and also indicate consistency of the GNSS and GRACE/GRACE-FO satellite missions' data in monitoring surface mass loading effects. The cases of weak correlations between ∆h obtained from GNSS data and the corresponding ones determined using GRACE/GRACE-FO satellite missions' data may be attributed to variety of reasons such as local temporal variations of water masses and GNSS station dependent errors. However, in this study the values of correlation coefficients got improved after removing the linear trends and there are only two ASG-EUPOS stations (NWT1 and WIEL) for which values remain under 0.5. The case of NWT1 has already been stated above, and regarding WIEL the authors believe that the value of the correlation coefficient may get improved with the use of the longer data span.  The standard deviations of the differences between ∆h obtained from GNSS data and the corresponding ones obtained from GRACE/GRACE-FO satellite missions' data were illustrated in Figure 7. Statistics of these standard deviations are given in Table 4. They reveal that the mean value of standard deviations of the differences between ∆hm-GNSS and ∆hm-GRACE decreased by ca. 18% (i.e., from 4.2 to 3.4 mm) after removing the linear trend of ∆h specified in Figure 5. The reduced value of standard deviations of ∆h differences after removing linear trend in ∆h signifies the reduction of amplitude discrepancy between detrended ∆hm-GNSS and detrended ∆hm-GRACE in relation to that of ∆hm-  The standard deviations of the differences between ∆h obtained from GNSS data and the corresponding ones obtained from GRACE/GRACE-FO satellite missions' data were illustrated in Figure 7. Statistics of these standard deviations are given in Table 4. They reveal that the mean value of standard deviations of the differences between ∆h m-GNSS and ∆h m-GRACE decreased by ca. 18% (i.e., from 4.2 to 3.4 mm) after removing the linear trend of ∆h specified in Figure 5. The reduced value of standard deviations of ∆h differences after removing linear trend in ∆h signifies the reduction of amplitude discrepancy between detrended ∆h m-GNSS and detrended ∆h m-GRACE in relation to that of ∆h m-GNSS and ∆h m-GRACE . Overall, the improvement, in terms of correlation coefficients ( Figure 6) and standard deviations of the ∆h differences (Figure 7), obtained after removing linear trends in ∆h is observed. The results also indicate that although ∆h obtained from GNSS and GRACE/GRACE-FO data agree well, linear trends in ∆h can vary significantly. Therefore, detrended ∆h m-GNSS are further used for estimating ∆EWT.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 22 GNSS and ∆hm-GRACE. Overall, the improvement, in terms of correlation coefficients ( Figure 6) and standard deviations of the ∆h differences (Figure 7), obtained after removing linear trends in ∆h is observed. The results also indicate that although ∆h obtained from GNSS and GRACE/GRACE-FO data agree well, linear trends in ∆h can vary significantly. Therefore, detrended ∆hm-GNSS are further used for estimating ∆EWT.

Monthly Variations of Equivalent Water Thickness
Monthly variations of equivalent water thickness ∆EWTm-GRACE and ∆EWTm-GNSS were determined with the use of Equations (2) and (4) as well as data from GRACE/GRACE-FO satellite missions and GNSS data from the ASG-EUPOS CORS network, respectively. Moreover, monthly variations of equivalent water thickness ∆EWTm-WGHM were obtained from the WGHM. In order to be consistent with the detrended ∆EWTm-GNSS, linear trends in ∆EWTm-GRACE and ∆EWTm-WGHM were also removed. Then, combined solutions ∆EWTm-CombSol were determined using Equation (5) with the weights GNSS w and GRACE w estimated by Equation (4) and ∆EWT determined from GRACE/GRACE-FO and GNSS data. The ratio between these weights (i.e., GNSS w : GRACE w ) is illustrated in Figure 8. The maximum (i.e., 0.73) and minimum (i.e., 0.14) values of this ratio were observed at the GNSS stations CHOJ and SIDZ, respectively. For ca. 72% of GNSS stations investigated (i.e., 69 GNSS stations of ASG-EUPOS network), the ratio between GNSS w and GRACE w is at the level of 0.46 ± 0.13. Time series of ∆EWTm-GRACE,

Monthly Variations of Equivalent Water Thickness
Monthly variations of equivalent water thickness ∆EWT m-GRACE and ∆EWT m-GNSS were determined with the use of Equations (2) and (4) as well as data from GRACE/GRACE-FO satellite missions and GNSS data from the ASG-EUPOS CORS network, respectively. Moreover, monthly variations of equivalent water thickness ∆EWT m-WGHM were obtained from the WGHM. In order to be consistent with the detrended ∆EWT m-GNSS , linear trends in ∆EWT m-GRACE and ∆EWT m-WGHM were also removed. Then, combined solutions ∆EWT m-CombSol were determined using Equation (5) with the weights w GNSS and w GRACE estimated by Equation (4) and ∆EWT determined from GRACE/GRACE-FO and GNSS data. The ratio between these weights (i.e., w GNSS :w GRACE ) is illustrated in Figure 8.
The maximum (i.e., 0.73) and minimum (i.e., 0.14) values of this ratio were observed at the GNSS stations CHOJ and SIDZ, respectively. For ca. 72% of GNSS stations investigated (i.e., 69 GNSS stations of ASG-EUPOS network), the ratio between w GNSS and w GRACE is at the level of 0.46 ± 0.13. Time series of ∆EWT m-GRACE , ∆EWT m-GNSS and ∆EWT m-CombSol as well as the corresponding ∆EWT m-WGHM obtained from the WGHM are illustrated in Figure 9. It is observed from the results presented in Figure 9 that likewise ∆h, there is also a seasonal pattern in ∆EWT. This seasonal pattern can be attributed to seasonal mass variations related with the hydrological cycle over the territory of Poland. Statistics of peak-to-peak variations of the seasonal pattern of the monthly equivalent water thickness variations are given in Table 5. From the values of medians and standard deviations, it is clear that overall ∆EWT m-CombSol are in better agreement with ∆EWT m-WGHM than those of ∆EWT m-GNSS and ∆EWT m-GRACE .
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 22 are illustrated in Figure 9. It is observed from the results presented in Figure 9 that likewise ∆h, there is also a seasonal pattern in ∆EWT. This seasonal pattern can be attributed to seasonal mass variations related with the hydrological cycle over the territory of Poland. Statistics of peak-to-peak variations of the seasonal pattern of the monthly equivalent water thickness variations are given in Table 5.
From the values of medians and standard deviations, it is clear that overall ∆EWTm-CombSol are in better agreement with ∆EWTm-WGHM than those of ∆EWTm-GNSS and ∆EWTm-GRACE.     In order to assess the performance of GNSS data from ASG-EUPOS CORS network for the determination of ∆EWT as well as for improving ∆EWT determined from GRACE/GRACE-FO satellite missions' data, correlations between (1) ∆EWT m-GNSS and ∆EWT WGHM , (2) ∆EWT m-GRACE and ∆EWT m-WGHM , and (3) ∆EWT m-CombSol and ∆EWT m-WGHM have been estimated ( Figure 10 and Table 6). Strong correlations (correlation coefficients ranging from 0.60 to 0.70), between ∆EWT m-GNSS and ∆EWT m-WGHM were obtained at 32% of GNSS sites investigated; moderate/weak correlations (correlation coefficients ranging from 0.30 to 0.59), were obtained at 56% of the GNSS sites investigated; at about 12% of GNSS sites investigated, no correlations (correlation coefficient ranging from −0.41 to 0.29) have been obtained. In terms of correlations between ∆EWT m-GRACE and ∆EWT m-WGHM , strong correlations (correlation coefficients ranging from 0.60 to 0.79) have been obtained at 57% of the stations investigated; moderate/weak correlations (correlation coefficients ranging from 0.30 to 0.59) were obtained in the case of 38% of the stations investigated; at the remaining 5% of the GNSS sites investigated correlation coefficients range from −0.04 to 0.29. It is clear that the obtained correlations between ∆EWT m-GRACE and ∆EWT m-WGHM are much stronger than the corresponding ones between ∆EWT m-GNSS and ∆EWT m-WGHM . The reason for this may be ascribed to the fact that ∆EWT m-GRACE is obtained directly from the monitoring of the gravitational effect of mass variations within the Earth's system, while ∆EWT m-GNSS is determined by inverting the response of the Earth's surface to the mass variation at a point. The response of the Earth's surface to the mass variation depends upon the underlying crustal properties. Moreover, ∆EWT m-GNSS is sensitive to the local signal and the presence of any such signal could affect the ∆EWT m-GNSS . At 83% of the GNSS sites investigated, strong correlation (correlation coefficients from 0.60 to 0.80) between ∆EWT m-CombSol and ∆EWT m-WGHM were obtained. At 16% of the GNSS stations investigated, the correlations are moderate/weak. For the remaining 1% of the GNSS stations investigated (one station) the correlation coefficient is 0.27. The values of coefficients of correlation between ∆EWT m-CombSol and ∆EWT m-WGHM have evidently been improved at 85% of the GNSS stations in comparison to values of coefficients of correlation between ∆EWT m-GRACE and ∆EWT m-WGHM . Furthermore, combining ∆EWT obtained from GNSS data with the corresponding ones determined from GRACE/GRACE-FO satellite missions' data does not change, however, the correlation coefficients' values at 3% of the GNSS stations investigated. At 12% of the GNSS stations, the addition of GNSS data worsen the ∆EWT in terms of correlations with the ∆EWT m-WGHM . The good agreement, in terms of correlations, between ∆EWT m-CombSol and ∆EWT m-WGHM , indicates that GNSS data from CORS network stations provides in general valuable information to improve GRACE/GRACE-FO solutions for monitoring mass variation within the Earth's system. With the use of Equation (7), the differences δ∆EWTm-GNSS, δ∆EWTm-GRACE and δ∆EWTm-CombSol were determined. Standard deviations of these differences were depicted in Figure 11. Statistics of these standard deviations are given in Table 7. The results presented in Figure 11 and Table 7  This can be attributed to the amplitude differences, due the fact that GNSS data present ∆EWT at points, while GRACE/GRACE-FO satellite mission data and the WGHM present ∆EWT over an area. Furthermore, amplitude differences between ∆EWTm-GRACE and ∆EWTm-WGHM may be due to (i) the coarse spatial resolution of the GRACE/GRACE-FO satellite missions data (ca. 3° × 3° at the equator), and (ii) the spatial resolution of the WGHM (ca. 0.5° × 0.5°), as well as uncertainties of the WGHM related to modeling of storage capacity, fluxes, and human intervention [70]. As the data from GNSS, GRACE/GRACE-FO satellite missions and the WGHM characterize with differences between each other, so incorporation of ∆EWTm-GNSS to ∆EWTm-GRACE may not reduce standard deviations of δ∆EWTm-CombSol in comparison to that of δ∆EWTm-GRACE.  Table 6. Statistics of correlation coefficient values presented in Figure 10. With the use of Equation (7), the differences δ∆EWT m-GNSS , δ∆EWT m-GRACE and δ∆EWT m-CombSol were determined. Standard deviations of these differences were depicted in Figure 11. Statistics of these standard deviations are given in Table 7. The results presented in Figure 11 and Table 7 Figure 11 [m].

Conclusions and Discussions
In this study, monthly vertical deformations of the Earth's surface ∆h as well as temporal variations of equivalent water thickness ∆EWT were determined using GNSS data from the ASG-EUPOS CORS network as well as the CSR RL06 GRACE/GRACE-FO-based GGMs over the area of Poland. Combined solutions of ∆EWT were developed by combining ∆EWT obtained from GNSS data and the corresponding ones determined from GRACE satellite mission data. The ∆EWT obtained from GNSS data, GRACE/GRACE-FO satellite missions' data and from the determined combined solutions of ∆EWT were evaluated using the corresponding values obtained from the WGHM hydrological model. The main findings reveal: Over the area investigated, the seasonal patterns in ∆h determined from GNSS data of the ASG-EUPOS CORS network ∆hm-GNSS are in a good agreement with the corresponding ones determined from GRACE/GRACE-FO satellite missions' data ∆hm-GRACE. The maximum and minimum values of ∆h occur during the summer-fall and winter-spring seasons, respectively. Linear trends in the time series of ∆hm-GNSS and ∆hm-GRACE are, however, not fully consistent. It is found that ∆hm-GRACE show   Figure 11 [m].

Conclusions and Discussions
In this study, monthly vertical deformations of the Earth's surface ∆h as well as temporal variations of equivalent water thickness ∆EWT were determined using GNSS data from the ASG-EUPOS CORS network as well as the CSR RL06 GRACE/GRACE-FO-based GGMs over the area of Poland. Combined solutions of ∆EWT were developed by combining ∆EWT obtained from GNSS data and the corresponding ones determined from GRACE satellite mission data. The ∆EWT obtained from GNSS data, GRACE/GRACE-FO satellite missions' data and from the determined combined solutions of ∆EWT were evaluated using the corresponding values obtained from the WGHM hydrological model. The main findings reveal: Over the area investigated, the seasonal patterns in ∆h determined from GNSS data of the ASG-EUPOS CORS network ∆h m-GNSS are in a good agreement with the corresponding ones determined from GRACE/GRACE-FO satellite missions' data ∆h m-GRACE . The maximum and minimum values of ∆h occur during the summer-fall and winter-spring seasons, respectively. Linear trends in the time series of ∆h m-GNSS and ∆h m-GRACE are, however, not fully consistent. It is found that ∆h m-GRACE show positive linear trends which can be related with water depletion over the territory of Poland. In the case of ∆h m-GNSS , the nature of linear trends significantly differs at the sites of the ASG-EUPOS CORS network. These differences in linear trends of ∆h m-GNSS can be due to monumental instability, anthropogenic activity, and geological features. Strong correlations between ∆h m-GNSS and ∆h m-GRACE were obtained at 68% of GNSS sites investigated. When removing the linear trends in ∆h m-GNSS and ∆h m-GRACE , the number of stations with strong correlations between detrended ∆h m-GNSS and detrended ∆h m-GRACE increased to 93% of GNSS sites investigated. Median values of standard deviations of differences between ∆h m-GNSS and ∆h m-GRACE as well as between detrended ∆h m-GNSS and detrended ∆h m-GRACE are 4.1 mm and 3.4 mm, respectively. Overall, removing the linear trends in ∆h m-GNSS and ∆h m-GRACE results in clear improvements of the fit, in terms of both correlations and standard deviations of the difference between ∆h determined from GNSS data of the ASG-EUPOS CORS network to the corresponding ones determined from GRACE/GRACE-FO satellite mission data.
Furthermore, ∆EWT determined from GRACE/GRACE-FO satellite missions' data (∆EWT m-GRACE ) and from GNSS data of the ASG-EUPOS CORS network (∆EWT m-GNSS ) as well as their combined solutions ∆EWT m-CombSol are also in a good agreement with the ∆EWT m-WGHM obtained from the hydrological model WGHM. Strong correlations between (1) ∆EWT m-GNSS and ∆EWT m-WGHM , (2) ∆EWT m-GRACE and ∆EWT m-WGHM , as well as (3) ∆EWT m-CombSol and ∆EWT m-WGHM were obtained at 32%, 57% and 83% of the GNSS sites investigated, respectively. It has been found that the combination of ∆EWT m-GNSS and ∆EWT m-GRACE improves the determination of ∆EWT, in terms of correlations coefficients, at 85% of GNSS sites investigated in comparison to the one obtained with the use of GRACE/GRACE-FO satellite missions' data only. The median values of standard deviations of differences ∆EWT m-GNSS , ∆EWT m-GRACE , and ∆EWT m-CombSol with respect to ∆EWT m-WGHM are 7 cm, 5 cm and 4 cm, respectively. The improvements, in terms of standard deviations of the differences, obtained for ∆EWT m-CombSol in comparison to ∆EWT m-GRACE are limited to 46% of the GNSS stations investigated, however, there are only 4% of the GNSS sites where standard deviations of the difference between ∆EWT m-CombSol and ∆EWT m-WGHM (i.e., δ∆EWT m-CombSol ) are larger than the corresponding ones obtained from the differences between ∆EWT m-GRACE and ∆EWT WGHM (i.e., δ∆EWT m-GRACE ). In general, the combination of ∆EWT m-GNSS and ∆EWT m-GRACE either reduce or do not change standard deviations of δ∆EWT m-CombSol in comparison to the corresponding ones of δ∆EWT m-GRACE .
Although linear trends in ∆h estimated from GNSS data of the ASG-EUPOS CORS network and the corresponding ones determined from GRACE/GRACE-FO satellite missions data are not fully consistent, the agreement in the seasonal pattern of ∆h and ∆EWT obtained from these both data types demonstrates the capability of using stations of the ASG-EUPOS CORS network as sensors for the determination of temporal mass variations within the Earth's system. Moreover, incorporation of ∆EWT m-GNSS to ∆EWT m-GRACE provides improved solutions for sensing mass variation of the Earth's surface. The ∆EWT m-CombSol developed from this incorporation can be thought as optimization of the determination of mass variations using GNSS data from the ASG-EUPOS CORS network and GRACE/GRACE-FO satellite mission data. This is due to the coarse spatial resolution of GRACE/GRACE-FO satellite mission data while GNSS provides local information about mass variations. Thus, ∆EWT m-GNSS and ∆EWT m-GRACE are complementary to each other.
Overall, the results obtained indicate that although national GNSS CORS networks are mainly developed for precise positioning related applications, they evidently provide a valuable information complementary to the one obtained from satellite data for the determination of temporal mass variations within the Earth's system. Data from the ASG-EUPOS CORS network in Poland improve models of temporal mass variations within the Earth's system obtained from GRACE/GRACE-FO satellite missions' data. However, the use of other national GNSS CORS networks operated worldwide for determining temporal mass variations within the Earth's system as well as for improving GRACE/GRACE-FO satellite missions' solutions is a subject of future investigation.