Assessment of Temporal Variations of Orthometric/ Normal Heights Induced by Hydrological Mass Variations over Large River Basins Using GRACE Mission Data

Almost half of the Earth’s land is covered by large river basins. Temporal variations of hydrological masses induce time-varying gravitational potential and temporal mass loading that deforms the Earth’s surface. These phenomena cause temporal variations of geoid/quasigeoid and ellipsoidal heights that result in temporal variations of orthometric/normal heights ∆H/∆H*. The aim of this research is to assess ∆H/∆H* induced by hydrological masses over large river basins using the Gravity Recovery and Climate Experiment (GRACE) satellite mission data. The results obtained reveal that for the river basin of a strong hydrological signal, ∆H/∆H* reach 8 cm. These ∆H/∆H* would be needed to reliably determine accurate orthometric/normal heights. The ∆H/∆H* do not exceed ±1 cm in the case of the river basin of the weak hydrological signal. The relation between hydrological mass changes and ∆H/∆H* was investigated. Correlations between ∆H/∆H* and temporal variations of equivalent water thickness were observed in 87% of river basins subareas out of which 45% exhibit strong correlations. The ∆H/∆H* determined over two river basins that characterize with the strongest and weakest temporal variations were analysed using the Principal Component Analysis method. The results obtained reveal that ∆H/∆H* in subareas of the same river basin can significantly differ (e.g., ±2 cm in the Amazon basin) from each other, and are strongly associated with different spatio-temporal patterns of the entire river basin.


Introduction
The orthometric/normal heights, defined as a distance between the geoid/quasigeoid surface and a point on the Earth's surface measured along the real/normal plumb line, are substantially needed for scientific research, e.g., for studying crustal deformation, as well as for a wide range of engineering applications such as construction and infrastructure projects (i.e., roads, bridges, railways, dams, etc.). These heights are traditionally obtained using levelling measurements, and nowadays, more often by the combination of a high accuracy precise gravimetric geoid/quasigeoid height with an ellipsoidal height. In the last years, research of numerous teams, e.g., [1][2][3], conducted by the Commission 2 "Gravity Field" of the International Association of Geodesy (IAG), the Intercommission Committee on Theory (ICCT), Joint Study Group 0.15 (JSG 0.15) [4] was focused on determining the gravimetric geoid/quasigeoid height of a sub-centimetre accuracy. The current multi-constellation Global Navigation Satellite System (GNSS) allows determining ellipsoidal heights at a few millimetres accuracy level, e.g., [5]. Taking into geoid/quasigeoid and ellipsoidal heights are changing due to temporal mass variations within the Earth system, temporal variations of orthometric/normal heights should be considered to determine orthometric/normal heights of high accuracy.
Temporal variations of hydrological masses are considered one of the main sources for temporal variations of orthometric/normal heights. They generate the time-varying gravitational potential and thereby temporal variations of the geoid/quasigeoid surface. They also induce mass loadings that deform the Earth's surface in both horizontal and vertical directions. Temporal variations of hydrological masses lead to different scenarios of temporal variations of orthometric/normal heights (see Figure 1). For example, when hydrological masses increase on the Earth's surface (Figure 1a), the gravitational potential on the geoid/quasigeoid surfaces also increases, which leads to ascending the geoid/quasigeoid surface. This is due to the fact that according to the Newton's low of gravitation the relation between the mass and the gravitational potential is the directly proportional (cf. [6] (Ch. 1, p. 1)). On the other hand, the increase of the hydrological masses on the Earth's surface induces mass load that results in downward elastic deformation. The opposite scenario will occur when hydrological masses decrease ( Figure 1b). Therefore, temporal variations of hydrological masses indicate that temporal variations of orthometric/normal heights might increase/decrease by two factors: variations of geoid/quasigeoid heights and vertical deformations of the Earth's surface. Moreover, negative coefficients of correlation between temporal variations of hydrological masses and temporal variations of orthometric/normal heights are expected. This is because an increase in hydrological masses would result in a decrease in orthometric/normal heights (e.g., [7][8][9][10]). The GRACE (Gravity Recovery and Climate Experiment [11]) satellite mission, operated in the years 2002-2017, provided unique data for the determination of temporal mass variations within the Earth system including temporal variations of hydrological masses on a global/regional scale. The success of the GRACE satellite mission has emphasized the need for such satellite missions. On 22 May 2018, the GRACE Follow-On (GRACE-FO) satellite mission, designed for a nominal lifetime of five years, was launched for further monitoring temporal mass variations within the Earth's system [12]. Several authors demonstrated the usefulness of GRACE satellite mission data for investigating temporal variations of geoid/quasigeoid heights, e.g., [13][14][15][16][17][18][19]. Moreover, many authors have proved that the estimated vertical deformations of the Earth's surface using GRACE satellite mission data are in a good agreement with the corresponding ones determined from other space geodetic techniques such as GNSS on a global scale, e.g., [20][21][22], on a continental scale, e.g., [23,24] in Europe, as well as on a local/regional scale, e.g., [25,26] in the Amazon basin, [27] in the West Africa, [28] in the East Africa, [29] in Japan, [30] in Bangladesh, [31][32][33] in China, [34,35] in Tibet, [36,37] in Poland, [38] in Greenland, and [39][40][41][42] in the North America, and Very Long Baseline Interferometry (VLBI) The GRACE (Gravity Recovery and Climate Experiment [11]) satellite mission, operated in the years 2002-2017, provided unique data for the determination of temporal mass variations within the Earth system including temporal variations of hydrological masses on a global/regional scale. The success of the GRACE satellite mission has emphasized the need for such satellite missions. On 22 May 2018, the GRACE Follow-On (GRACE-FO) satellite mission, designed for a nominal lifetime of five years, was launched for further monitoring temporal mass variations within the Earth's system [12]. Several authors demonstrated the usefulness of GRACE satellite mission data for investigating temporal variations of geoid/quasigeoid heights, e.g., [13][14][15][16][17][18][19]. Moreover, many authors have proved that the estimated vertical deformations of the Earth's surface using GRACE satellite mission data are in a good agreement with the corresponding ones determined from other space geodetic techniques such as GNSS on a global scale, e.g., [20][21][22], on a continental scale, e.g., [23,24] in Europe, as well as on a local/regional scale, e.g., [25,26] in the Amazon basin, [27] in the West Africa, [28] in the East Africa, [29] in Japan, [30] in Bangladesh, [31][32][33] in China, [34,35] in Tibet, [36,37] in Poland, [38] in Greenland, and [39][40][41][42] in the North America, and Very Long Baseline Interferometry (VLBI) (e.g., [43]). The use of GRACE satellite mission data for estimating physical height changes, i.e., temporal variations of orthometric/normal heights, was investigated on local and regional scales for the area of Turkey, Poland, Central Europe and Greenland [7][8][9][10]44].
In the majority of land areas over the world, temporal variations of orthometric/normal heights induced by the changes of hydrological masses have not been investigated. In other words, orthometric/normal heights are, so far, practically considered static. The only exception are areas of an evident land subsidence and post-glacial land uplift. For such areas, the secular change of orthometric/normal heights was considered (e.g., [45][46][47]). The main aim of this research is to assess temporal variations of orthometric/normal heights induced by temporal variations of hydrological masses over almost a half of Earth's land areas, i.e., over twenty-four worldwide large river basins (Figure 2), using GRACE satellite mission data. Temporal variations of hydrological masses deform both the geoid/quasigeoid surfaces and the surface of the Earth's which result in orthometric/normal height changes (see Figure 1). Such assessment is vital for obtaining reliable, accurate orthometric/normal heights to monitor, interpret, analyse and model changes in orthometric/normal heights used in scientific and engineering works. The river basins chosen represent main sources of temporal variations of hydrological masses over the continents. They cover land areas that range from 5.89 × 10 6 km 2 for the Amazon basin to 0.44 × 10 6 km 2 for the Don basin. Considering the spatial resolution of GRACE satellite mission data, these area sizes are suitable to reliably investigate the differences between orthometric/normal height changes over the entire river basin. Furthermore, the large river basins selected characterize with huge diversities in terms of climate, environment, geological structure, topography, tectonics and geodynamics. Such diversities allow the assessment of temporal variations of orthometric/normal heights on broad ranges. The data used and the methods implemented to determine and analyse orthometric/normal height changes are described in Section 2. In Section 3, results concerning orthometric/normal height changes in the time domain, as well as time series of temporal variations of geoid/quasigeoid heights and vertical deformations of the Earth's surface over the selected large river basins were presented and analysed. Moreover, the relation between temporal variations of orthometric/normal heights obtained from GRACE data and temporal variations of equivalent water thickness obtained from an independent hydrological model over large basins are discussed. Furthermore, in Section 3, orthometric/normal height changes in the space-time domain for large river basins of the strongest and weakest water change signal are analysed. In Section 4, discussions and conclusions concerning orthometric/normal heights changes assessed over large basins using GRACE satellite mission data are given.
In the majority of land areas over the world, temporal variations of orthometric/normal heights induced by the changes of hydrological masses have not been investigated. In other words, orthometric/normal heights are, so far, practically considered static. The only exception are areas of an evident land subsidence and post-glacial land uplift. For such areas, the secular change of orthometric/normal heights was considered (e.g., [45][46][47]). The main aim of this research is to assess temporal variations of orthometric/normal heights induced by temporal variations of hydrological masses over almost a half of Earth's land areas, i.e., over twenty-four worldwide large river basins (Figure 2), using GRACE satellite mission data. Temporal variations of hydrological masses deform both the geoid/quasigeoid surfaces and the surface of the Earth's which result in orthometric/normal height changes (see Figure 1). Such assessment is vital for obtaining reliable, accurate orthometric/normal heights to monitor, interpret, analyse and model changes in orthometric/normal heights used in scientific and engineering works. The river basins chosen represent main sources of temporal variations of hydrological masses over the continents. They cover land areas that range from 5.89 × 10 6 km 2 for the Amazon basin to 0.44 × 10 6 km 2 for the Don basin. Considering the spatial resolution of GRACE satellite mission data, these area sizes are suitable to reliably investigate the differences between orthometric/normal height changes over the entire river basin. Furthermore, the large river basins selected characterize with huge diversities in terms of climate, environment, geological structure, topography, tectonics and geodynamics. Such diversities allow the assessment of temporal variations of orthometric/normal heights on broad ranges. The data used and the methods implemented to determine and analyse orthometric/normal height changes are described in Section 2. In Section 3, results concerning orthometric/normal height changes in the time domain, as well as time series of temporal variations of geoid/quasigeoid heights and vertical deformations of the Earth's surface over the selected large river basins were presented and analysed. Moreover, the relation between temporal variations of orthometric/normal heights obtained from GRACE data and temporal variations of equivalent water thickness obtained from an independent hydrological model over large basins are discussed. Furthermore, in Section 3, orthometric/normal height changes in the space-time domain for large river basins of the strongest and weakest water change signal are analysed. In Section 4, discussions and conclusions concerning orthometric/normal heights changes assessed over large basins using GRACE satellite mission data are given.

Data and Methods Used
In this study, temporal variations of orthometric ∆H and normal heights ∆H* were estimated at subareas coincided with the Jet Propulsion Laboratory (JPL) mascons solution grid using the release 6 (RL06) of monthly GRACE-based GGMs from the University of Texas at Austin, Center of Space Research (CSR) computation centre [48], load Love numbers which are calculated using the Preliminary Reference Earth Model (PREM) [49], and the IGiK-TVGMF (Instytut Geodezji i Kartografii-Temporal Variations of Gravity/Mass Functionals) package [9]. It should be noted the JPL mascons solution grid was chosen due to the fact that this grid was defined by one of the official GRACE Science Data System centers, i.e., the JPL center, as well as the JPL 3 • equal-area mascon grid corresponds to the spatial resolution of GRACE data. The degree-1 and degree-2 harmonic coefficients of CSR RL06 GRACE-based Global Geopotential Models (GGMs) were replaced by degree-1 from the solution described in [50] and the degree-2 obtained from Satellite Laser Ranging (SLR) [51]. The DDK3 decorrelation filter [52] that compromises between reducing noise and keeping the signal [18], was utilized to reduce the noise included in CSR RL06 GRACE-based GGMs. Moreover, these GGMs were truncated at d/o 60 that corresponds to the spatial resolution of the DDK3 filter. GRACE-based GGMs RL06, developed by the CSR, GFZ (German Research Centre for Geosciences) and JPL were intercompared by several researchers. For example, Göttl et al. [53] illustrated consistency between these GRACE-based GGMs RL06 developed by CSR, JPL and GFZ. Kvas et al. [54] concluded that GRACE-based GGMs RL06 developed by the CSR, GFZ and JPL exhibit the same signal content. Adhikari et al. [55] revealed that differences, in terms of sea-level fingerprints, between these three solutions are within 1-sigma uncertainties.

Data and Methods Used
In this study, temporal variations of orthometric ∆H and normal heights ∆H* were estimated at subareas coincided with the Jet Propulsion Laboratory (JPL) mascons solution grid using the release 6 (RL06) of monthly GRACE-based GGMs from the University of Texas at Austin, Center of Space Research (CSR) computation centre [48], load Love numbers which are calculated using the Preliminary Reference Earth Model (PREM) [49], and the IGiK-TVGMF (Instytut Geodezji i Kartografii-Temporal Variations of Gravity/Mass Functionals) package [9]. It should be noted the JPL mascons solution grid was chosen due to the fact that this grid was defined by one of the official GRACE Science Data System centers, i.e., the JPL center, as well as the JPL 3° equal-area mascon grid corresponds to the spatial resolution of GRACE data. The degree-1 and degree-2 harmonic coefficients of CSR RL06 GRACE-based Global Geopotential Models (GGMs) were replaced by degree-1 from the solution described in [50] and the degree-2 obtained from Satellite Laser Ranging (SLR) [51]. The DDK3 decorrelation filter [52] that compromises between reducing noise and keeping the signal [18], was utilized to reduce the noise included in CSR RL06 GRACE-based GGMs. Moreover, these GGMs were truncated at d/o 60 that corresponds to the spatial resolution of the DDK3 filter. GRACE-based GGMs RL06, developed by the CSR, GFZ (German Research Centre for Geosciences) and JPL were intercompared by several researchers. For example, Göttl et al. [53] illustrated consistency between these GRACE-based GGMs RL06 developed by CSR, JPL and GFZ. Kvas et al. [54] concluded that GRACE-based GGMs RL06 developed by the CSR, GFZ and JPL exhibit the same signal content. Adhikari et al. [55] revealed that differences, in terms of sea-level fingerprints, between these three solutions are within 1-sigma uncertainties.
Temporal variations of orthometric/normal heights ∆H/∆H* were obtained from GRACE-based GGMs as follows [8,9], * for orthometric heights for normal heights where ∆h denotes the vertical deformation of the Earth's surface, ∆N presents the temporal variation of a geoid height and ∆ζ is the temporal variation of a quasigeoid height. The approximate relation between geoid and quasigeoid heights can be expressed as follows [6] (Ch. 8, p. 328, Equation (8-104)), where H is the orthometric height (in kilometres) of the computation point and H av is the average height (in kilometres) over the computation area. On the basis of Equation (2), the difference δ between ∆ζ and ∆N is: For the highest point on the Earth (H ≈ 8.848 km-Mount Everest), and H av ≈ 6.100 km (the average elevation of Great Himalayas), assuming hydrological mass changes inducing maximum 10 centimetres variation in both the orthometric height of the computation point and the average height of the area, the δ will reach the level of 15 × 10 −5 m (0.15 mm), which is negligible. Therefore, under the assumption that orthometric height changes from one epoch to another do not exceed 10 cm, ∆N can merely be considered equal to ∆ζ.
The ∆N, ∆ζ and ∆h were estimated using GRACE-based GGMs as follows [9, [56][57]: where ∆h denotes the vertical deformation of the Earth's surface, ∆N presents the temporal variation of a geoid height and ∆ζ is the temporal variation of a quasigeoid height.
The approximate relation between geoid and quasigeoid heights can be expressed as follows [6] (Ch. 8, p. 328, Equation (8)-(104)), where H is the orthometric height (in kilometres) of the computation point and H av is the average height (in kilometres) over the computation area. On the basis of Equation (2), the difference δ between ∆ζ and ∆N is: For the highest point on the Earth (H ≈ 8.848 km-Mount Everest), and H av ≈ 6.100 km (the average elevation of Great Himalayas), assuming hydrological mass changes inducing maximum 10 centimetres variation in both the orthometric height of the computation point and the average height of the area, the δ will reach the level of 15 × 10 −5 m (0.15 mm), which is negligible. Therefore, under the assumption that orthometric height changes from one epoch to another do not exceed 10 cm, ∆N can merely be considered equal to ∆ζ. The ∆N, ∆ζ and ∆h were estimated using GRACE-based GGMs as follows [9, 56,57]: and surface density coefficients ∆C σ lm and ∆S σ lm defined as, with where ϕ, λ are spherical geocentric coordinates of the computation point, a is the semi-major axis of reference ellipsoid, GM is the product of the Newtonian gravitational constant G and the Earth's mass M, r is the geocentric radius of the computation point on the ellipsoid, C W lm , S W lm are fully normalised Stokes coefficients from GRACE-based GGMs, C U lm , S U lm denote spherical harmonic coefficients of the normal gravity field, ∆C lm , ∆S lm are the differences between C lm , S lm obtained from monthly RL06 GRACE-based GGMs and the corresponding ones obtained from a suitable reference model, P lm are the fully normalized Legendre functions of degree l and order m, l max is the applied maximum degree, ρ w is the water density, ρ ave is the Earth's average density, k l and h l are the load Love numbers of degree l, and γ denotes the normal gravity at the reference ellipsoid for ∆N and at the physical surface of the Earth's for ∆ζ.
The hydrological mass changes are the main signal of mass transport that causes ∆H/∆H* over large basins investigated. So that in order to verify H/∆H* determined from GRACE satellite mission data (cf. Equation (1)), the relation between those ∆H/∆H* and temporal variations of equivalent water thickness ∆EWT obtained from an independent source is investigated. As a source of ∆EWT values the WGHM (WaterGAP Global Hydrology Model) was used within this study. These ∆EWT consist of temporal variations of fast-surface and subsurface runoff, groundwater recharge and river discharge as well as storage variations of water in canopy, snow, soil, groundwater, lakes, wetlands and rivers. The WGHM delivers water storages and fluxes for land areas of the Earth, except Antarctica, at a spatial resolution of 0.5 • × 0.5 • . The basic temporal resolution of the model is one day, whereas predominantly the model outputs are used at monthly time steps [58]. Since 2003, the WGHM has been developed integrating GRACE total water storage for calibration and data assimilation [58]. Many investigations, e.g., [59][60][61][62][63], regarding the comparison of temporal mass variations within the Earth's system, obtained from GRACE data with those from WGHM, have been conducted. In general, these investigations reveal a good agreement, in terms of annual periods and phases, between ∆EWT obtained from WGHM and GRACE satellite mission data, however, the solutions, based on WGHM data, underestimate large decadal declining and rising water storage trends with respect to those obtained from GRACE satellite data (e.g., [64]). Therefore, in this study, linear trends in ∆EWT, and ∆H/∆H* were removed to avoid discrepancies between WGHM and GRACE data in terms of water storage trends. In this study, ∆EWT values from the recent version of WGHM (i.e., version 2.2d), developed in the Goethe University in Frankfurt [58], were used.
The Principal Component Analysis (PCA) method [65] was utilized to analyse the spatio-temporal variations of ∆H/∆H*. This method allows the analysis of ∆H/∆H* in the space-time domain. According to Rangelova and Sideris [16] and Li et al. [66], the PCA method is one of the methods recommended for the analysis and modelling of temporal mass variations within the Earth's system and temporal variations of gravity functionals. The fundamentals of the PCA method have widely been presented in the Earth science related textbooks (e.g., [65,67]). The general aim of PCA method is to reduce the dimensionality within a dataset. This is achieved by decomposing original data matrix X into empirical orthogonal functions (EOFs) and the corresponding principal components (PCs) time series obtained using the singular value decomposition as follows, where UD are PCs time series and V are EOFs which are the eigenvectors of the data autocovariance matrix. The detailed algorithm and computation steps used within this study to obtain PCs time series and EOFs can be found in Godah et al. [8]. In this investigation, ∆H/∆H* normalized with their standard deviations for subareas of the river basins were used to determine the matrix X. The analysis of the spatio-temporal variations of ∆H/∆H* was conducted using the IGiK-TVGMF package [9].
It should be noted that in the PCA method, the first PCs time series approximate the data. Therefore, in this study, only the 1st to 4th PCs time series were considered.

Orthometric/Normal Height Changes in the Time Domain
On the basis of GRACE satellite mission data, temporal variations of geoid/quasigeoid heights ( Figure 3) and vertical deformations of the Earth's surface ( Figure 4) were estimated over large river basins specified in Figure 2 using Equations (4) and (5), respectively. These height changes were estimated with a high accuracy, as GRACE satellite mission data are sufficient to determine geoid height changes for a short period (i.e., less than 3 months) with an accuracy of 0.25 mm (see [68]). The phase shift between mean temporal variations of geoid/quasigeoid heights and mean vertical deformations of the Earth's surface over large river basins investigated are depicted in Figure 5. The results presented in Figure 3 indicate that the amplitude of temporal variations of geoid/quasigeoid heights reach 8 mm over the Amazon basin, which agrees with the results presented in [26], and 2 mm over the Orange basin. The results presented in Figure 4 reveal that the amplitude of vertical deformations of the Earth's surface reach the level of 13 mm over the Amazon basin, which agrees with the results given in [13], and 3 mm over the Orange basin. The amplitudes of temporal variations of geoid/quasigeoid heights and vertical deformations of the Earth's surface for the remaining large river basins investigated are at the level of 5 mm, and 8 mm, respectively. Moreover, the results presented in Figures 3 and 4 reveal that, for most of river basins investigated, seasonal patterns of temporal variations of geoid/quasigeoid heights are shifted in phase with respect to vertical deformations of the Earth's surface by approximately 6 months (cf. Figure 5). This illustrates the relation between temporal variations of hydrological masses, temporal variations of geoid/quasigeoid heights and the crustal deformations in the vertical component. The ∆H/∆H* over large river basins investigated ( Figure 6) were determined by combining temporal variations of geoid/quasigeoid heights (cf. Figure 3) and vertical deformations of the Earth's surface (cf. Figure 4). The dispersions (Max-Min) of ∆h, ∆N/∆ζand ∆H/∆H* for the river basins investigated, taking into consideration all subareas within the river basin, are illustrated in Figure 7. The extreme values of ∆h, ∆N/∆ζand ∆H/∆H* are almost symmetric in relation to their mean values that are approximately equal to zero. Figures 6 and 7 show that pick-to-pick variations of ∆H/∆H* differ significantly (e.g., ca. 8 cm for the Amazon basin and ca. 2 cm for the Orange basin) for large river basins investigated within this study. This is due to the significant difference between temporal variations of hydrological masses for these large river basins. At different epochs, the amplitudes of ∆H/∆H* can reach up to ca. 40 mm, and at most 5 mm over the Amazon, and Orange river basins, respectively. For the remaining large river basins, these amplitudes do not exceed 13 mm. The results presented in Figure 6 also reveal that for the period between 2005 and 2016, secular patterns of ∆H/∆H* appear over some large river basins, e.g., Danube, Dnieper, Don and Tigris and Euphrates river basins. These secular patterns are estimated to 1.5 ± 0.5 mm/year, and they are ascribed to the drought pattern over these river basins. Furthermore, Figure 6 exhibits, for some river basins, a good agreement between the average of ∆H/∆H* over the whole area of the river basin and ∆H/∆H* from individual subareas (i.e., individual cell of the JPL 3 • equal-area mascon grid) of the river basin. For example, for the Danube, Dnieper, Don, Volga and Murray Darling river basins, differences between ∆H/∆H* at subareas of these basins with respect to their average of ∆H/∆H* values are at the level of ±3 mm. For the remaining river basins, these differences become larger reaching almost ±10 mm for the Nile river basin, ±15 mm for the Congo and La Plata river basins, and ±20 mm for the Amazon river basin. This indicates that ∆H/∆H* over the entire river basins investigated (except for Danube, Dnieper, Don, Volga and Murray Darling river basins) significantly differ from one subarea to another. Additionally, for each subarea of large river basins investigated root mean squares (RMS) of ∆h time series, RMS of ∆N/∆ζtime series and RMS of ∆H/∆H* time series were estimated. The maximum and minimum of the RMS obtained for each large river basin investigated are provided in Table 1. For Danube, Dnieper, Don, Volga and Murray Darling river basins the range of RMS does not exceed 0.7 mm for ∆N/∆ζ, 1.0 mm for ∆h and 1.6 mm for ∆H/∆H*. For the remaining river basins, the differences in RMS between subareas are from 1.0 mm for ∆N/∆ζ, 1.1 mm for ∆h and 2.3 mm for ∆H/∆H*for Orange river basins to 6.5 mm for ∆N/∆ζ, 10.7 mm for ∆h and 17.3 mm for ∆H/∆H*for the Amazon. This reveals that amplitudes of ∆h, ∆N/∆ζand ∆H/∆H* at subareas over entire each of large river basins investigated (except the Danube, Dnieper, Don, Volga and Murray Darling river basins) are inconsistent. mm for ∆N/∆ζ, 1.0 mm for ∆h and 1.6 mm for ∆H/ΔH*. For the remaining river basins, the differences in RMS between subareas are from 1.0 mm for ∆N/∆ζ, 1.1 mm for ∆h and 2.3 mm for ∆H/ΔH*for Orange river basins to 6.5 mm for ∆N/∆ζ, 10.7 mm for ∆h and 17.3 mm for ∆H/ΔH*for the Amazon. This reveals that amplitudes of ∆h, ∆N/∆ζ and ∆H/∆H* at subareas over entire each of large river basins investigated (except the Danube, Dnieper, Don, Volga and Murray Darling river basins) are inconsistent.      In order to investigate the impact of temporal variations of water mass, represented by ∆EWT, on ∆H/∆H*, coefficients of correlations between ∆H/∆H* and ∆EWT were computed. The WGHM was used as an independent source of data to estimate ∆EWT that are linearly related to the surface density coefficients used in Equation (5) (cf. [69]). In order to avoid the discrepancies in linear trends of mass changes estimated from WGHM and GRACE data, ∆EWT and ∆H/∆H* were detrended. Figure 8 shows coefficients of correlations between detrended ∆H/∆H* and detrended ∆EWT. It illustrates that at 48% of subareas (i.e., mascon locations) investigated, strong negative correlations between ∆H/∆H* and ∆EWT (correlation coefficients ranging from −0.97 to −0.70) were obtained. Weak/moderate negative correlations between ∆H/∆H* and ∆EWT (coefficients of correlation ranging from −0.69 to −0.30) were observed at 40% of subareas investigated. At ca. 9% of subareas investigated the coefficients of correlation obtained range from −0.29 to zero, while at 3% of subareas investigated positive values of correlation coefficients (from zero to +0.37) are observed. The subareas exhibiting no correlations between ∆H/∆H* and ∆EWT (ca. 12% of subareas investigated) are mainly clustered in areas of week hydrological signals such as in southern Sahara. For a few subareas in Lena, Mackenzie, Ganges-Brahmaputra and Yangtze (Chang Jiang) river basins the reason of the observed lack of correlation between ∆H/∆H* from GRACE satellite mission data and ∆EWT from the WGHM is not evident.  In order to investigate the impact of temporal variations of water mass, represented by ∆EWT, on ∆H/∆H*, coefficients of correlations between ∆H/ΔH* and ∆EWT were computed. The WGHM was used as an independent source of data to estimate ∆EWT that are linearly related to the surface density coefficients used in Equation (5) (cf. [69]). In order to avoid the discrepancies in linear trends of mass changes estimated from WGHM and GRACE data, ∆EWT and ∆H/ΔH* were detrended.

Orthometric/Normal Height Changes in the Space-Time Domain
As examples, spatio-temporal variations of orthometric/normal heights for large river basins of the strongest and weakest ∆EWT signal, i.e., for Amazon and Orange basins, normalized with their

Orthometric/Normal Height Changes in the Space-Time Domain
As examples, spatio-temporal variations of orthometric/normal heights for large river basins of the strongest and weakest ∆EWT signal, i.e., for Amazon and Orange basins, normalized with their standard deviations were investigated. These spatio-temporal variations were obtained using the PCA method. In order to avoid gaps in ∆H/∆H* time series, GRACE data from the period between January 2004 and December 2010 were used in the investigation. Figures 9 and 10 show the results of the analysis of temporal variations of orthometric/normal heights ∆H/∆H* using the PCA method, i.e., the percentage of a total variance of ∆H/∆H* and the first, second, third and fourth PC time series, for the Amazon, and Orange basins, respectively.   The Orange basin, in contrary to the Amazon basin, is characterized by the weakest signal of ΔH/ΔH* among the river basins investigated (cf. Figure 7). This is because ΔEWT obtained from GRACE satellite mission data over this basin are very weak [75]. This might be ascribed to the fact that the majority of the Orange basin is covered by desert areas, i.e., Kalahari Desert and Namib Desert. The weakest ΔH/ΔH* signals over the Orange basin are observed. The results presented in Figure 10 indicate that in terms of the percentage of a total variance of ∆H/∆H* reflected by PCAs, the PC time series 1 reflects ca. 85% of ΔH/ΔH* signal and the second one reflects ca.  The results presented in Figure 9 reveal that, in terms of the percentage of a total variance of temporal variations of orthometric/normal heights reflected by PCA/EOFs, the first two PCs time series reflect ca. 97% of ∆H/∆H* signal for the Amazon basin. The PC time series 1 accounts the largest variance, i.e., 77%, of the total signal. This is due to the fact that it recovers the dominant spatio-temporal component of the signal. It clearly illustrates seasonal variations of ∆H/∆H* with minimum values in March-May and maximum values in September-November. This is because over the Amazon basin, the minimum and maximum values of ∆EWT were observed in March-May and September-November, respectively (e.g., [70][71][72]). It also shows that the largest ∆H/∆H* can be seen in 2005 resulting from the extreme drought in 2005 [73], followed by minimum ∆H/∆H* values due to the larger than usual flood that occurred in 2006 and 2009 [74]. Temporal variations of the PC time series 2 precede those of PC time series 1 in phase by about two months. According to [70][71][72], this phase delay is due to the fact that for the Amazon basin, the average rainfall maximum is in February while the maximum of water storage is delayed by ca. 2 months, peaking in March/April. The EOF1 reveals that the maximum EOF loading pattern is within the downstream area of the Amazon basin, i.e., the eastern and south-eastern parts of the Amazon basin that combines the effect of rainfall-induced soil moisture and groundwater variations and seasonally varying storage in surface water for the river and inundation areas [61]. In terms of the EOF2 the Amazon basin is divided into two distinguished parts, northern and southern one. The positive values of the EOF2 loading pattern are obtained for the northern part, while negative values of the EOF2 loading pattern for the southern part. This north-south division can be ascribed to the fact that the rainfall seasonality induced by the seasonal north-south migration of the Intertropical Convergence Zone (ITC), and the fact that the minimum water storage in the northern part of the Amazon basin in observed from December to February while in the southern part-from July to September. Overall, the first and second EOFs obtained illustrate clear differences between EOFs loading patterns for subareas of the entire Amazon basin.
The Orange basin, in contrary to the Amazon basin, is characterized by the weakest signal of ∆H/∆H* among the river basins investigated (cf. Figure 7). This is because ∆EWT obtained from GRACE satellite mission data over this basin are very weak [75]. This might be ascribed to the fact that the majority of the Orange basin is covered by desert areas, i.e., Kalahari Desert and Namib Desert. The weakest ∆H/∆H* signals over the Orange basin are observed. The results presented in Figure 10 indicate that in terms of the percentage of a total variance of ∆H/∆H* reflected by PCAs, the PC time series 1 reflects ca. 85% of ∆H/∆H* signal and the second one reflects ca. 9% of ∆H/∆H* signal for the Orange basin. They recover the dominant spatio-temporal component of the ∆H/∆H* signal over the Orange basin. For the PC time series 1, a downward trend is observed and also a sharp fall at the end of 2005 and the beginning of 2006 is clearly noticeable. At that time there were heavy rains that caused flooding in the area of Southern Africa [76]. In the case of the PC time series 2, for the period 2005-2011, a distinct seasonal pattern is observed. This seasonal pattern coincides with the wet, and dry seasons, respectively, over the Orange basin, and it might be attributed to the seasonal runoff linked with climatic conditions within the Orange basin [77]. For 2004, this seasonal pattern is unclear. This might be ascribed to the fact that practically no rainfall was recorded in the fall season of 2004 over Orange basin [76]. The EOF1 reveals that the maximum EOF loading values occur in the north-eastern part of the Orange basin, whereas the minimum values were obtained for the south-western part. There is a slope from north-eastern to south-western in the EOF1 pattern. It means that subareas, located in the north-eastern part of the basin, have the major impact on PC time series 1. This is because in the north-eastern part of the Orange basin, the annual rainfall exceeds the annual evaporation, and thereby, induces significant run off [77]. On the contrary, in the south-western part of the Orange basin the annual evaporation exceeds the annual rainfall [77]. Moreover, the northern catchment areas of the Orange basin are characterized with a relatively high rainfall, while the central and southern catchments areas of the basin are characterized with lower rainfall and high evaporation. Therefore, relatively strong water mass variations signal induced from the northern catchment areas are not explained by PC time series 1 in the central and southern catchment areas of the basin [78]. The opposite situation takes place in the case of the EOF2. The descending slope of loading pattern from south-western to north-eastern is observed, and thus subareas located in the south-west have the biggest impact on PC time series 2. This is due to the fact that the topography and climate of the Orange basin varies significantly, i.e., from mountains to flat areas and from moist sub-humid to hyper-arid (deserts) climate zones, from northern to southern catchment areas [67]. Furthermore, 75% of the annual runoff occurs during summer months and 25% of the runoff occurs during winter [77]. The lower part of the Orange basin is also characterized by a complex geological structure [77], and therefore, the EOF2 pattern may be attributed to the geological structure and climatic conditions prevailing in southern Orange basin.
It should be noted that over the Amazon basin, the third and fourth PCs time series accounts for 3%, and 1% of total variance of ∆H/∆H* signal, respectively. For the Orange basin, total variances of ∆H/∆H* signal were at the level of 3% for the third PC time series and the level of 2% for the fourth PC time series. In terms of total variance of ∆H/∆H* signal, the third and fourth PCA/EOFs can be considered insignificant in comparison to the first and second PCs time series. Furthermore, the summation of total variances of ∆H/∆H* signal for the remaining 5th to 75th PCs time series for the Amazon basin and the summation of total variances of ∆H/∆H* signal for the remaining 5th to 15th PCs time series over the Orange basin are at the level of 1%. Therefore, the analysis of these remaining PCs time series were neglected in this investigation.

Discussion and Conclusions
Although, orthometric/normal heights are considered static over the majority of the Earth land areas, except in areas of land uplift, temporal mass variations within the Earth's system result in both temporal variations of geoid/quasigeoid heights and vertical deformations of the Earth's surface, and thereby, temporal variations of orthometric/normal heights. Considering the fact that approximately half of the Earth's land areas are covered by large river basins, in this study, temporal variations of orthometric/normal heights ∆H/∆H* for 24 large river basins over the world were determined using GRACE-based GGMs. These ∆H/∆H* are needed to correct orthometric/normal heights. Correlations between these ∆H/∆H* and temporal variations of equivalent water thickness ∆EWT obtained from the WGHM hydrological model were investigated. The spatio-temporal analysis of ∆H/∆H* for two river basinsn characterized by the strongest and the weakest ∆H/∆H* signaln was conducted using the PC method.
The main findings of the research are as follows: 1.
Amplitudes of ∆H/∆H* significantly differ from each other between large river basins investigated. The largest amplitudes were obtained for river basins located over tropical rainforests, whilst the lowest ones were observed over desert areas. This is due to the fact that hydrological mass changes within large river basins located in the tropical rain forest are substantially larger than the ones over desert areas. For example, the range of ∆H/∆H* for the same subarea at different epochs reaches 8 cm over the Amazon basin and only ca. 2 cm for the Orange basin.

2.
For some large river basins ∆H/∆H* time series obtained for particular subareas of the river basin are similar and close to each other, e.g., Danube, Dnieper and Don basins. This is because the hydrological mass changes patterns over the entire large river basins are consistent. However, in many cases differences between ∆H/∆H* obtained at subareas within the same river basin are significant, sometimes exceeding three times the amplitude of the average of ∆H/∆H* over the whole river basin, e.g., for the Congo basin. This is due to the fact that hydrological mass changes patterns substantially differ among subareas located in the same large river basin.

3.
For 88% of river basins subareas negative correlations between detrended ∆H/∆H* and detrended ∆EWT were observed, whereas they are strong for 48% of those subareas. This is because the increase of hydrological masses results in decrease of ∆H/∆H*, and vice versa, the decrease of hydrological masses results in increase of ∆H/∆H*. For the remaining 12% of river basins subareas there are no correlations between detrended ∆H/∆H* and detrended ∆EWT. The main reason for observing no correlations can be ascribed to the fact that the majority of those subareas are located in regions of a very weak hydrological signal (e.g., southern Sahara). However, further investigations concerning the correlation between ∆H/∆H* from GRACE satellite mission data and ∆EWT from the WGHM for some subareas such as the ones located in Lena, Mackenzie and Ganges-Brahmaputra river basins are recommended.

4.
For Amazon and Orange basins, the 1st and 2 nd PCs time series reflect together 95%, and 94% of a total variance of ∆H/∆H* signal, respectively. Although, these percentages are close to each other, the contribution of the first and second PCs time series is different in the case of both river basins, i.e., the first and second PCs time series are 77%, and 20%, respectively, for the Amazon basin, and the first and second PCs time series are 85% and 9%, respectively, for the Orange basin. The first and second PCs time series exhibit that ∆H/∆H* are not an artefact, but are a consequence of the processes inducing hydrological mass transport. These temporal variations of orthometric/normal heights are strongly associated to different spatio-temporal patterns of the entire river basins. They can be associated with the extreme drought, the unusual flood, the location of the upstream and downstream areas of the river basin, the rainfall seasonality induced by the seasonal migration of the Intertropical Convergence Zone (ITC) and other climatic conditions, the spatial distribution of the water storage of the entire river basin and the geological structure of the river basin.
Overall, the results obtained in this investigation clearly prove that GRACE satellite mission data are significant for determining ∆H/∆H* induced by temporal variations of hydrological masses over large river basins. These temporal variations are needed for precise determination of orthometric/normal heights to monitor, interpret, analyse and model changes in those heights used in scientific (e.g., crustal motion, subsidence and isostatic readjustment) and engineering (e.g., the deformation of bridges, dams and large constructions) works, especially for areas characterized with significant temporal mass variations within the Earth's system such as hydrological masses in large basins.