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

Abstract

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.

1. 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 clearly revolutionized research in many disciplines of Earth sciences. Besides GRACE-type satellite missions, the Global Navigation Satellite System (GNSS) has proven itself as one of the most powerful space geodesy tools for the determination of temporal mass variations in the Earth’s system and related geodynamic processes. It is successfully used to study the elastic ground loading deformation in response to the Earth’s hydrological mass variations. Since the beginning of this century, numerous studies on the relation between the mass loading and the Earth’s surface deformations obtained from GNSS data have been conducted on a global scale (e.g., [3,4,5,6,7]), on a continental scale, e.g., South America [8]; and Europe [9], as well as on a local/regional scale, e.g., Japan [10]; the Amazon basin [11]; Bangladesh [12]; Greenland [13]; the West Africa [14]; North America [15,16,17,18]; China [19,20,21,22]; Tibet [23,24]; Poland [25,26]; and the East Africa [27]. Except van Dam et al. [9] that indicated poor correlation between the annual vertical crustal deformation from GRACE and GPS at 36 GPS sites over Europe, all those studies revealed very good agreement (correlation coefficient >0.6) between the Earth’s surface deformations obtained from GNSS data and the corresponding ones determined from GRACE data. In addition to GNSS and GRACE data, other space geodetic techniques, e.g., the very long baseline interferometry (VLBI), indicated that the vertical crustal deformations due to continental water loading reach up to 15 mm [28].

The latest official GRACE/GRACE-FO level-2 products, i.e., RL06 (release 6) GRACE/GRACE-FO-based Global Geopotential Models (GGMs), developed by GFZ (GeoForschungsZentrum), CSR (Center for Space Research) and JPL (Jet Propulsion Laboratory) centers (cf. [29,30,31,32]) are currently available for scientific users. In addition to GRACE/GRACE-FO-based GGMs, global mass concentration (mascon) products of NASA (National Aeronautics and Space Administration), GSFC (Goddard Space Flight Center), JPL and CSR have also been released to the scientific community [33,34,35]. On the other hand, national GNSS CORS (Continuously Operating Reference Stations) networks are rapidly growing and being continuously expanded around the world. Currently, national GNSS CORS networks are in operation in Australia, North America and almost all European countries (e.g., networks of the EUPOS multifunctional reference stations) as well as in some countries in Africa, Asia and South America. National GNSS CORS networks are regarded as a multifunctional precise positioning system that provides very valuable information for: (1) static and Real-Time Kinematic (RTK) positioning; (2) monitoring crustal deformations; (3) definition and realization of global, regional and national geodetic reference frames; (4) atmosphere related studies; etc. Besides national GNSS CORS networks, a number of dense regional/local GNSS CORS networks have been established for a specific scientific research, e.g., the PBO H₂O (Plate Boundary Observatory Studies of the Water Cycle) established to study the 3D strain field resulting from deformation across the active boundary zone between the Pacific and North American plates in the western United States as well as remotely measured changes in soil moisture content, snow depth, water content of the troposphere, and ground motions related to changes in the water table (e.g., [36]).

The coarse spatial resolution of GRACE/GRACE-FO products is one of the main weaknesses of GRACE-type missions. Moreover, GRACE/GRACE-FO level-2 products (e.g., GRACE/GRACE-FO-based GGMs) as well as GRACE-based mascon products being the results of a complex processing of GRACE mission data, are released to scientific users with a delay of a few months. On the other hand, temporal mass variations can be obtained using GNSS data in real or near real time, but they represent the local effect for a very limited area within the radius of a few kilometers from a GNSS station. Furthermore, coordinate variations obtained from GNSS data are affected by many errors, e.g., the draconitic error that systematically affects the seasonal periods in GNSS series (e.g., [37]). However, investigation concerning the use of national GNSS CORS networks to study temporal mass variations within the Earth’s system is essentially needed. The main objective of this study is to investigate the use of national GNSS CORS networks for monitoring temporal mass variations within the Earth’s system. It is also aimed at the use of national GNSS CORS networks data to improve temporal mass variations within the Earth’s system on a regional/local scale determined from GRACE/GRACE-FO satellite missions’ data. The study area and data used are described in 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.

2. 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.

3. 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 Section 3.1 and Section 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.

3.1. 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]:

$$\mathsf{\Delta }{h}_{\mathrm{m}-\mathrm{GRACE}}=a{\displaystyle \sum _{n=0}^{N_{\max }}{\displaystyle \sum _{m=0}^n{P}_{nm}(\sin \phi )}\frac{h_n}{1+k_n}}(\mathsf{\Delta }{C}_{nm}\cos m\lambda +\mathsf{\Delta }{S}_{nm}\sin m\lambda )$$

(1)

$$\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}=\frac{a\times {\rho }_{avg}}{3{\rho }_w}{\displaystyle \sum _{n=0}^{N_{\max }}{\displaystyle \sum _{m=0}^n{P}_{nm}(\sin \phi )}}\frac{2n+1}{1+k_n}(\mathsf{\Delta }{C}_{nm}\cos m\lambda +\mathsf{\Delta }{S}_{nm}\sin m\lambda )$$

(2)

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], $\mathsf{\Delta }{C}_{nm}^{}$ and $\mathsf{\Delta }{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]):

$$\begin{array}{l}{P}_{00}(\sin \phi )=1.0\\ P_{11}(\sin \phi )=\sqrt{3}\cos \phi \\ P_{nn}(\sin \phi )=\sqrt{\frac{(2n+1)}{n}}\cos \phi P_{(n-1)(n-1)}(\sin \phi )\\ P_{nm}(\sin \phi )=\sqrt{\frac{(2n-1)(2n+1)}{(n-m)(n+m)}}\cos \phi P_{(n-1)m}(\sin \phi )\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\sqrt{\frac{(2n+1)(n+m-1)(n-m-1)}{(2n+3)(n+m)(n-m)}}\cos \phi P_{(n-2)m}(\sin \phi )\\ \end{array}\}$$

(3)

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.

3.2. 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]).

$$\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}=\mathsf{\Delta }{h}_{\mathrm{m}-\mathrm{GNSS}}\cdot {\left(2{\rho }_w\cdot R_0\cdot g\cdot \left(\frac{1-{\upsilon }^2}{E}\right)\right)}^{-1}$$

(4)

where R₀ 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.

3.3. 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:

$$\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{CombSol}}^i=\frac{w_{{}_{\mathrm{GNSS}}}^i\cdot \mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}^i+w_{{}_{\mathrm{GRACE}}}^i\cdot \mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}^i}{w_{{}_{\mathrm{GNSS}}}^i+w_{{}_{\mathrm{GRACE}}}^i}$$

(5)

where q denotes the total number of GNSS stations of the ASG-EUPOS CORS network investigated and $w_{\mathrm{GNSS}}$ and $w_{\mathrm{GRACE}}$ are weights for ΔEWT determined from GNSS and GRACE/GRACE-FO satellite missions’ data, respectively. These weights were estimated as follows:

$$\begin{array}{l}{w}_{{}_{\mathrm{GNSS}}}^i={\left(\sqrt{\frac{{\displaystyle \sum _{t=1}^k{\left(\left(\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}^i(t)-\mathsf{\Delta }EW{T}_{\mathrm{WGHM}}^i(t)\right)-\overline{\left(\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}^i(t)-\mathsf{\Delta }EW{T}_{\mathrm{WGHM}}^i(t)\right)}\right)}^2}}{k-1}}\right)}^{-1}\\ w_{{}_{\mathrm{GRACE}}}^i={\left(\sqrt{\frac{{\displaystyle \sum _{t=1}^k{\left(\left(\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}^i(t)-\mathsf{\Delta }EW{T}_{\mathrm{WGHM}}^i(t)\right)-\overline{\left(\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}^i(t)-\mathsf{\Delta }EW{T}_{\mathrm{WGHM}}^i(t)\right)}\right)}^2}}{k-1}}\right)}^{-1}\end{array}\}$$

(6)

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.

3.4. 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

$$\begin{array}{l}\mathsf{\delta }\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}^i=\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{WGHM}}^i-\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GRACE}}^i\\ \mathsf{\delta }\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}^i=\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{WGHM}}^i-\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{GNSS}}^i\\ \mathsf{\delta }\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{CombSol}}^i=\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{WGHM}}^i-\mathsf{\Delta }EW{T}_{\mathrm{m}-\mathrm{ComSol}}^i\end{array}\}$$

(7)

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.

4. Results

4.1. 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 Δ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. Statistics of peak-to-peak variations of Δh_m-GNSS and Δh_m-GRACE [mm].

Δh	Min	Max	Median	Std
Δhm-GNSS	10.2	38.2	21.7	4.0
Δhm-GRACE	13.8	21.4	20.0	1.6

.

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. Statistics of linear trends of vertical deformation of the Earth’s surface estimated at ASG-EUPOS sites (except sites SKSL, SKSK and SKSV) using GRACE/GRACE-FO satellite missions and GNSS data [mm/year].

Δh	Min	Max	Median	Std
Δhm-GNSS	−1.47	2.16	0.07	0.58
Δhm-GRACE	0.67	1.86	0.83	0.25

. 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) Δ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. Statistics of correlation coefficient values presented in Figure 6.

Δh	Min	Max	Mean
GNSS–GRACE	0.13	0.91	0.64
GNSS–GRACE (detrended)	0.31	0.83	0.72

. 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. Figure 4 and Figure 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.

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. Statistics of standard deviation values presented in Figure 7 [mm].

Δh	Min	Max	Mean
GNSS–GRACE	2.6	7.3	4.2
GNSS–GRACE (detrended)	2.5	5.2	3.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.

4.2. 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_{\mathrm{GNSS}}$ and $w_{\mathrm{GRACE}}$ estimated by Equation (4) and ΔEWT determined from GRACE/GRACE-FO and GNSS data. The ratio between these weights (i.e., $w_{\mathrm{GNSS}}$:$w_{\mathrm{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_{\mathrm{GNSS}}$ and $w_{\mathrm{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. Statistics of peak-to-peak variations of ΔEWT_m-GNSS, ΔEWT_m-GRACE, ΔEWT_m-CombSol and ΔEWT_m-WGHM [cm].

ΔEWT	Min	Max	Median	Std
ΔEWTm-GNSS	11.1	56.9	38.5	8.5
ΔEWTm-GRACE	11.1	22.2	18.7	3.2
ΔEWTm-CombSol	6.7	28.4	22.2	4.3
ΔEWTm-WGHM	5.7	37.1	21.4	4.9

. 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.

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. Statistics of correlation coefficient values presented in Figure 10.

ΔEWT	Min	Max	Mean
GNSS	−0.41	0.71	0.51
GRACE	−0.04	0.79	0.58
Combined solution	0.27	0.83	0.67

). 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 δΔ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. Statistics of standard deviation values presented in Figure 11 [m].

ΔEWT	Min	Max	Mean
GNSS	0.03	0.12	0.07
GRACE	0.03	0.07	0.05
Combined solution	0.02	0.06	0.04

. The results presented in Figure 11 and Table 7 reveal that the mean values of standard deviations of differences δΔEWT_m-GNSS, δΔEWT_m-GRACE and δΔEWT_m-CombSol are 7 cm, 5 cm and 4 cm, respectively. The values of standard deviations of δΔEWT_m-CombSol decreased at 44 sites of the ASG-EUPOS CORS network in comparison to standard deviations of δΔEWT_m-GRACE. At 48 sites of the ASG-EUPOS CORS network, standard deviations of δΔEWT_m-CombSol remained the same as the values of standard deviations of δΔEWT_m-GRACE. There are only 4 sites of the ASG-EUPOS CORS network where standard deviations of δΔEWT_m-CombSol increased in comparison to that of the δΔEWT_m-GRACE. Although there are improvements in the median values of standard deviations of δΔEWT_m-CombSol in comparison to the corresponding ones of δΔEWT_m-GRACE, they are limited to only 46% of GNSS sites investigated. 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 ΔEWT_m-GRACE and ΔEWT_m-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 ΔEWT_m-GNSS to ΔEWT_m-GRACE may not reduce standard deviations of δΔEWT_m-CombSol in comparison to that of δΔEWT_m-GRACE. Since, there are only 4 GNSS stations where the addition of ΔEWT_m-GNSS to ΔEWT_m-GRACE increase the standard deviation values of δΔEWT_m-CombSol, addition of ΔEWT_m-GNSS to ΔEWT_m-GRACE, in general, either reduced or did not change the values of standard deviations in comparison to the corresponding ones of δΔEWT_m-GRACE.

5. 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.