Contribution of GRACE Satellite Mission to the Determination of Orthometric/Normal Heights Corrected for Their Dynamics—A Case Study of Poland

: Physical heights were traditionally determined without considering the dynamic processes of the Earth induced from temporal mass variations. The Gravity Recovery and Climate Experiment (GRACE) mission provided valuable data that allow the estimation of geoid/quasigeoid height changes and vertical deformations of the Earth’s surface induced from temporal mass loading, and thereby temporal variations of physical heights. The objective of this investigation is to discuss the determination of orthometric/normal heights considering mass transports within the Earth’s system. An approach to determine such heights was proposed. First, temporal variations of or-thometric/normal heights ( ∆ H / ∆ H * ) were determined using the release 6 GRACE-based Global Geopotential Models together with load Love numbers obtained from the preliminary reference Earth model. Then, those variations were modelled and predicted using the seasonal decomposition (SD) method. The proposed approach was tested on the territory of Poland. The main results obtained reveal that ∆ H / ∆ H * over the area investigated are at the level of a couple of centimetres and that they can be modelled and predicted with a millimetre accuracy using the SD method. Orthomet-ric/normal heights corrected for their dynamics can be determined by combining modelled ∆ H / ∆ H * with orthometric/normal heights referred to a speciﬁc reference epoch.


Introduction
The determination of heights is one of the major tasks in geodesy. Knowledge of heights is needed in various engineering applications as well as in scientific research associated with the Earth's shape and its dynamics. Physical heights, e.g., orthometric/normal heights (H/H * ), were traditionally obtained using the expensive and timeconsuming precise spirit levelling. Nowadays, they are more and more frequently determined as a combination of ellipsoidal heights from precise satellite navigation techniques with geoid/quasigeoid heights derived from gravity data. Extensive investigations have thus been conducted to determine geoid/quasigeoid heights, ellipsoidal heights, and the resulting orthometric/normal heights of high accuracy. The issue of modelling gravimetric geoid/quasigeoid of a sub-cm accuracy is currently discussed by many authors, e.g., [1][2][3]. It is also considered as one of the main activities of the International Association of Geodesy [4]. The ellipsoidal heights can be determined with a couple of millimetres accuracy using long GNSS (global navigation satellite system) observation sessions and a high-precision scientific GNSS data processing software [5]. The H/H * of the mean square error of the adjusted levelling network less than 1 mm/km can be obtained from the first-order levelling network established using spirit levelling.

The Determination of Orthometric/Normal Heights Corrected for Their Dyn
The orthometric/normal heights corrected for their dynamics H(t)/H * (t) can fined at arbitrary time t as a sum of the orthometric/normal height H(t0)/H * (t0) at erence epoch t0, and temporal variations of orthometric/normal heights ∆H(t)/∆H The mathematical formulation of the proposed approach for the determin H(t)/H * (t) at epoch t is expressed as follows: where ∆H M (t)/∆H *M (t) are values of ∆H/∆H * at time t calculated from the ∆H/∆H This approach is described as follows (cf.

The Determination of Orthometric/Normal Heights Corrected for Their Dynamics
The orthometric/normal heights corrected for their dynamics H(t)/H * (t) can be defined at arbitrary time t as a sum of the orthometric/normal height H(t 0 )/H * (t 0 ) at the reference epoch t 0 , and temporal variations of orthometric/normal heights ∆H(t)/∆H * (t) at t.
The mathematical formulation of the proposed approach for the determination of H(t)/H * (t) at epoch t is expressed as follows: H * (t) = H * (t 0 ) + ∆H *M (t) (2) where ∆H M (t)/∆H *M (t) are values of ∆H/∆H * at time t calculated from the ∆H/∆H * model. This approach is described as follows (cf.  Precise levelling networks, i.e., first order levelling networks, on a national as well as a regional scale, were commonly established by levelling measurement campaigns conducted over a long-time span in different epochs. Thus, a proper reduction of those levelling measurements to a common epoch is applied to determine H(t0)/H * (t0) at the reference epoch t0. Moreover, H(t0)/H * (t0) are, nowadays, determined by combining the ellipsoidal height of high accuracy obtained from GNSS measurements with geoid/quasigeoid height obtained from a precise geoid/quasigeoid model developed with the use of high-quality gravity data. It should be noted that establishing/re-measuring a national precise levelling network usually takes many years. Thus, secular (i.e., linear trend) variations of orthometric/normal heights that occurred during these years must be considered for the establishment of the 1st order levelling network referenced to a specific epoch such as t0. Although ∆H(t)/∆H * (t) can theoretically be obtained with the use of geometric data, e.g., spirit levelling, in practice, it is not possible to repeat levelling campaigns sufficiently frequently to provide as frequent as monthly solutions for heights [9]. Thus, monthly variations of orthometric/normal heights can practically be estimated using data from dedicated gravity satellite missions, in particular the GRACE/GRACE-FO missions, combined with longer time series data such as repeated levelling, glacial isostatic adjustment (GIA), GNSS time series, and tide gauge records, that also contain variations of physical heights signal. In this study, ∆H(t)/∆H * (t) are estimated as follows Precise levelling networks, i.e., first order levelling networks, on a national as well as a regional scale, were commonly established by levelling measurement campaigns conducted over a long-time span in different epochs. Thus, a proper reduction of those levelling measurements to a common epoch is applied to determine H(t 0 )/H * (t 0 ) at the reference epoch t 0 . Moreover, H(t 0 )/H * (t 0 ) are, nowadays, determined by combining the ellipsoidal height of high accuracy obtained from GNSS measurements with geoid/quasigeoid height obtained from a precise geoid/quasigeoid model developed with the use of high-quality gravity data. It should be noted that establishing/re-measuring a national precise levelling network usually takes many years. Thus, secular (i.e., linear trend) variations of orthometric/normal heights that occurred during these years must be considered for the establishment of the 1st order levelling network referenced to a specific epoch such as t 0 . Although ∆H(t)/∆H * (t) can theoretically be obtained with the use of geometric data, e.g., spirit levelling, in practice, it is not possible to repeat levelling campaigns sufficiently frequently to provide as frequent as monthly solutions for heights [9]. Thus, monthly variations of orthometric/normal heights can practically be estimated using data from dedicated gravity satellite missions, in particular the GRACE/GRACE-FO missions, combined with longer time series data such as repeated levelling, glacial isostatic adjustment (GIA), GNSS time series, and tide gauge records, that also contain variations of physical heights signal. In this study, ∆H(t)/∆H * (t) are estimated as follows where t G i (i = 1, 2, . . . , n) is the epoch of the ith (of n considered) monthly GRACE/GRACE-FO global geopotential models (GGMs), T(t G i )/T * (t G i ) present the linear trend of ∆H(t)/∆H * (t) obtained from GRACE/GRACE-FO data, and T ∆H (t G i )/T * ∆H * (t G i ) denote the linear trend (i.e., secular variation) of ∆H(t)/∆H * (t), estimated with the use of data from repeated levelling, GIA, GNSS, or tide gauges. It should be noted that when the epoch t 0 , to which orthometric/normal heights of the levelling network are referenced does not correspond to one of the epochs t G i , a proper transformation of H(t 0 )/H * (t 0 ) to the closest t G i epoch would be required.
The ∆N(t G i ), ∆ζ(t G i ) and ∆h(t G i ) can be estimated from GRACE/GRACE-FO-based GGMs as follows [33,39]: and surface density coefficients where ϕ, λ, r are the latitude, longitude, and the geocentric radius of the computation point Q on the physical surface of the Earth, a is the semi-major axis of the reference ellipsoid, GM is the product of the Newtonian gravitational constant and the Earth's mass, C U lm and S U lm are spherical harmonic coefficients of the normal gravity field, C W lm and S W lm are fully normalized spherical harmonics from GRACE/GRACE-FO-based GGMs, γ denotes the normal gravity at Q, P lm is a fully normalized Legendre function of degree l and order m, l max is the maximum degree applied, ρ avg is the Earth's average density, ρ w is the water density, and h l and k l are the load Love numbers (LLN) of degree l.
Geoid heights N(t G i ) and N(t 0 ) are obtained using Equations (5) and (6), respectively, considering the reduction of the potential induced by topographic masses, i.e., the topographic bias [40], and assuming r = R and γ = γ 0 , where R is the mean radius of the Earth and γ 0 is the normal gravity at Q 0 on the ellipsoid, i.e., the corresponding point of Q projected along the plumb line onto the ellipsoid. The ∆N(t G i ) can be considered equal to ∆ζ(t G i ) [35]. The ∆H/∆H * can be modelled using a suitable method such as the seasonal decomposition (SD) method [41] and the principle component analysis/empirical orthogonal function (PCA/EOF) method [42]. In this research, the SD method, following the algorithm and computation steps described in detail in [29], was applied to model ∆H/∆H * . The ∆H c (t)/∆H *c (t) time series is obtained by centering the ∆H(t)/∆H * (t) time series to zero. Then, the ∆H c (t)/∆H *c (t) is decomposed into a seasonal component S(t)/S * (t), a long-term component LT(t)/LT * (t), and an unmodelled residual component ε(t)/ε * (t). The model of temporal variations of orthometric/normal heights ∆H SD (t)/∆H *SD (t) is obtained as a sum of long-term and seasonal components. The ∆H M (t)/∆H *M (t) were determined by referring the ∆H SD (t)/∆H *SD (t) to the reference epoch t 0 . Then, ∆H M (t)/∆H *M (t) values were interpolated/extrapolated at the height datum benchmarks investigated. The H(t)/H * (t) were calculated by combining these interpolated/extrapolated values with H(t 0 )/H * (t 0 ) at those sites.
In order to mitigate the effect of the gaps in data from the GRACE/GRACE-FO satellite missions and the latency of GRACE/GRACE-FO products, prediction of the ∆H M /∆H *M values for months where these products are unavailable, would be required. Within this study, the prediction of ∆H M /∆H *M values were investigated using the empirical approach proposed in Godah et al. [29] based on the results obtained from the SD method. The predicted values of ∆H M /∆H *M were calculated as a sum of the S(t)/S * (t) component and the ∆H/∆H * values obtained from an appropriate mathematical model fitted to LT(t)/LT * (t) component.

A Case Study
The area of Poland, as a unique one covered with high-quality terrestrial gravity datasets, GNSS/levelling data, and gravimetric geoid/quasigeoid models, has been chosen as a study area. The accuracy of quasigeoid heights obtained from the recent gravimetric quasigeoid model and GNSS/levelling data as well as combined GGMs, e.g., the EGM2008 (Earth Gravitational Model 2008; [43]), over this area is estimated at the level of 1 cm [44]. Over this area, ∆N/∆ζ were reliably determined, analyzed, and modelled using GRACE satellite mission data [29], and ∆h were detected at the ASG-EUPOS (Polish Active Geodetic European Position Determination System; http://www.asgeupos.pl/, accessed on 1 August 2022) network sites using GNSS data [21,22].

Orthometric/Normal Heights at Reference Epoch
In this study, H(t 0 )/H * (t 0 ) from the ASG-EUPOS network sites were utilized. The values of normal heights referred to the Kronstadt86 vertical datum for the territory of Poland (PL-KRON86-NH) are accessible via the ASG-EUPOS website (http://www.asge upos.pl/index.php accessed on 1 August 2022). These normal heights were transformed into the PL-EVRF2007-NH (European Vertical Reference Frame 2007 recently applied for the area of Poland). The reference epoch t 0 of PL-EVRF2007-NH is 1 January 2008, which corresponds to its realization epoch [45]. The PL-EVRF2007-NH is based on results of the adjustment of the 4th levelling campaign in Poland (1998-2012), EUVN (European Vertical Reference Network), and EUVN-DA (European Vertical Reference Network-Densification Action) solutions as well as solutions for the benchmarks of eccentric points of permanent GNSS stations of the ASG-EUPOS network. Heights are given in the zero tidal system. The accuracy of levelling measurements was estimated as 0.74 mm/km. Based on the adjustment of the levelling network, the estimated accuracy of a single benchmark is 3.5 mm [46].
The orthometric heights at the ASG-EUPOS sites at the reference epoch t 0 were determined using the well-known basic relation for the geoid-quasigeoid separation [47]: where ∆g B is the Bouguer gravity anomaly and γ is the mean normal gravity along the normal plumb line between the ellipsoid and the telluroid.
The H and H * of ASG-EUPOS sites referred to the PL-EVRF2007-NH at the reference epoch t 0 are shown in Figure 3. It should be noted that the differences between orthometric and normal heights over the area of Poland do not exceed a couple of decimeters [44].
where ∆g B is the Bouguer gravity anomaly and  is the mean normal gravity along th normal plumb line between the ellipsoid and the telluroid.
The H and H * of ASG-EUPOS sites referred to the PL-EVRF2007-NH at the referenc epoch t0 are shown in Figure 3. It should be noted that the differences between ortho metric and normal heights over the area of Poland do not exceed a couple of decimeter [44].

Temporal Variations of Orthometric/Normal Heights
Temporal variations of orthometric/normal heights ∆H(t)/∆H * (t) over Poland in duced by TMVES were considered. In Section 4.2.1, temporal variations of an equivalen water thickness ∆EWT(t) and ∆H(t)/∆H * (t) over the area investigated were estimated, an the relation between ∆EWT(t) and ∆H(t)/∆H * (t) was discussed. The territory of Poland was divided into equal area blocks of a size equivalent to th spatial resolution of GRACE data, i.e., ~334 × 334 km 2 or 3° × 3° on the equator. Thos equal area blocks are localized in the same configuration as the Jet Propulsion Laborator (JPL) mascons (mass concentrations), global 3° × 3° equal area spherical cap mascon so lutions. Figure 4 shows the area investigated, including the location of four mascon used. The main TMVES over this area can be ascribed to temporal hydrological mas variations [28]. Temporal variations of equivalent water thickness at the location of those mascons obtained from monthly release 06 (RL06) JPL GRAC mascon solutions using the mascon visualization tool (c

Temporal Variations of Orthometric/Normal Heights
Temporal variations of orthometric/normal heights ∆H(t)/∆H * (t) over Poland induced by TMVES were considered. In Section 4.2.1, temporal variations of an equivalent water thickness ∆EWT(t) and ∆H(t)/∆H * (t) over the area investigated were estimated, and the relation between ∆EWT(t) and ∆H(t)/∆H * (t) was discussed.  Figure 4 shows the area investigated, including the location of four mascons used. The main TMVES over this area can be ascribed to temporal hydrological mass variations [28]. Temporal variations of equivalent water thickness ∆EWT(t G i ) over t G i -t 0 at the location of those mascons obtained from monthly release 06 (RL06) JPL GRACE mascon solutions using the mascon visualization tool (cf. https://ccar.colorado.edu/grace/index. html, accessed on 1 August 2022) are depicted in Figure 5.
The results presented in Figure 5 indicate a distinctive pattern of seasonal water mass variations with minimum and maximum values in July-September and March, respectively. This seasonal variations pattern agrees well with the corresponding one documented in Krynski et al. [28]. It reveals that the decrease/increase in water masses over Poland results from intensive water evaporation during dry months in the summer season and the accumulation of snow/water in the winter season. https://ccar.colorado.edu/grace/index.html, accessed on 1 August 2022) are depicted in Figure 5.  The results presented in Figure 5 indicate a distinctive pattern of seasonal water mass variations with minimum and maximum values in July-September and March, respectively. This seasonal variations pattern agrees well with the corresponding one documented in Krynski et al. [28]. It reveals that the decrease/increase in water masses over Poland results from intensive water evaporation during dry months in the summer season and the accumulation of snow/water in the winter season. The , were estimated at the location of mascons (cf. Figure 4) using monthly CSR (Computation Centre of the University of Texas Center of Space Research) RL06 GRACE-based GGMs [48], LLN calculated using the preliminary reference Earth model (PREM; [49]), and the IGiK-TVGMF (Instytut Geodezji i Kartografii-Temporal Variations of Gravity/Mass Functionals; [50]) software. The degree-1 and degree-2 harmonic coefficients from the solution described in Swenson et al. [51] and the corresponding ones obtained from satellite laser ranging (SLR) [52], respectively, replaced those of CSR RL06 GRACE-based GGMs. In order to reduce the noise https://ccar.colorado.edu/grace/index.html, accessed on 1 August 2022) are depicted in Figure 5.  The results presented in Figure 5 indicate a distinctive pattern of seasonal water mass variations with minimum and maximum values in July-September and March, respectively. This seasonal variations pattern agrees well with the corresponding one documented in Krynski et al. [28]. It reveals that the decrease/increase in water masses over Poland results from intensive water evaporation during dry months in the summer season and the accumulation of snow/water in the winter season. The , were estimated at the location of mascons (cf. Figure 4) using monthly CSR (Computation Centre of the University of Texas Center of Space Research) RL06 GRACE-based GGMs [48], LLN calculated using the preliminary reference Earth model (PREM; [49]), and the IGiK-TVGMF (Instytut Geodezji i Kartografii-Temporal Variations of Gravity/Mass Functionals; [50]) software. The degree-1 and degree-2 harmonic coefficients from the solution described in Swenson et al. [51] and the corresponding ones obtained from satellite laser ranging (SLR) [52], respectively, replaced those of CSR RL06 GRACE-based GGMs. In order to reduce the noise The ∆N(t G i )/∆ζ(t G i ), as well as ∆h(t G i ), were estimated at the location of mascons (cf. Figure 4) using monthly CSR (Computation Centre of the University of Texas Center of Space Research) RL06 GRACE-based GGMs [48], LLN calculated using the preliminary reference Earth model (PREM; [49]), and the IGiK-TVGMF (Instytut Geodezji i Kartografii-Temporal Variations of Gravity/Mass Functionals; [50]) software. The degree-1 and degree-2 harmonic coefficients from the solution described in Swenson et al. [51] and the corresponding ones obtained from satellite laser ranging (SLR) [52], respectively, replaced those of CSR RL06 GRACE-based GGMs. In order to reduce the noise included in CSR RL06 GRACE-based GGMs, the decorrelation filter DDK3 [53] was utilized. This filter was chosen as it compromises between reducing the noise and keeping the signal over the area investigated [29]. Moreover, those GGMs were truncated at d/o 60 that corresponds to the spatial resolution of the mascon grid ( Figure 4) and the DDK3 filter. Then, ∆H(t G i )/∆H * (t G i ) were estimated as a sum of detrended ∆N(t G i )/∆ζ(t G i ), detrended ∆h(t G i ), and the linear trend (i.e., secular variations) of H/H * obtained from the NKG2016LU model [54]. The resulting ∆N(t G i )/∆ζ(t G i ), ∆h(t G i ), and ∆H(t G i )/∆H * (t G i ) for the period April 2002-August 2016 are depicted in Figure 6. lized. This filter was chosen as it compromises between reducing the noise and keeping the signal over the area investigated [29]. Moreover, those GGMs were truncated at d/o 60 that corresponds to the spatial resolution of the mascon grid ( Figure 4) and the DDK3 filter. Then, were estimated as a sum of detrended , and the linear trend (i.e., secular variations) of H/H * obtained from the NKG2016LU model [54]. The resulting ( ) / ( ) for the period April 2002-August 2016 are depicted in Figure 6.
were obtained at the beginning of spring (in February-April) and minimum values at the end of sum- The maximum values of ∆N(t G i )/∆ζ(t G i ) were obtained at the beginning of spring (in February-April) and minimum values at the end of summer/beginning of autumn (in August-October); they are consistent with the corresponding ones estimated by Krynski et al. [28] and Godah et al. [29,30,32]. For ∆h(t G i ), maximum values were obtained in August-October and minimum values in February-April, which correspond to minima and maxima of ∆N(t G i )/∆ζ(t G i ), respectively. For a subarea, the dispersions of ∆N(t G i )/∆ζ(t G i ) and ∆h(t G i ) reach 8 mm and 17 mm, respectively. Minimum and maximum values of ∆H(t G i )/∆H * (t G i ), obtained as the combination of ∆N(t G i )/∆ζ(t G i ) and ∆h(t G i ), were observed in February-April and in August-October, respectively. The ∆H(t G i )/∆H * (t G i ) for the same subarea reach up to 23 mm, and 4 mm from one subarea to another at the same epoch. The results presented in Figure 6 also indicate that approximately 66% of the ∆H(t G i )/∆H * (t G i ) signal is due to ∆h, while the remaining part of this signal is induced from ∆N/∆ζ. The coefficients of correlation between ∆EWT(t G i ) (cf. Figure 5) and ∆N(t G i )/∆ζ(t G i ), ∆h(t G i ), and ∆H(t G i )/∆H * (t G i ) are presented in Figures 7-9, respectively, showing strong correlations between water mass changes and ∆N/∆ζ as well as ∆h. The correlations between ∆EWT(t G i ) and ∆N(t G i )/∆ζ(t G i ) are slightly stronger (correlation coefficients of 0.89 ± 0.02) compared to those between ∆EWT(t G i ) and ∆h(t G i ), characterized by correlation coefficients of −0.71 ± 0.04. This might be due to the fact that N/ζ is basically associated with the disturbing potential generated from the Earth's mass distribution, while the vertical deformations of the Earth's surface may depend on other factors, e.g., the elasticity of the Earth's crust. Overall, the results obtained reveal that the resulting ∆H/∆H * are strongly correlated with ∆EWT (correlation coeffiecients of −0.80 ± 0.04).
of ( )/ ( ) were observed in February-April and in August-October, respectively. The for the same subarea reach up to 23 mm, and 4 mm from one subarea to another at the same epoch. The results presented in Figure 6 also indicate that approximately 66% of the H t H t signal is due to Δh, while the remaining part of this signal is induced from ΔN/Δζ. The coefficients of correlation between Figure 5) and       Figure 11 shows models developed using the SD method (see Equation (15)). It also illustrates the differences δ∆H/δ∆H * between and their models The seasonal decomposition (SD) method was used for the analysis and modelling of the ∆H/∆H * determined. The ∆H c (t G i )/∆H * c (t G i ) time series, obtained by centering the ∆H(t G i )/∆H * (t G i ) time series to zero for the period between January 2004 and December 2010 exhibiting no gaps in monthly CSR RL06 GRACE-based GGMs, was applied. The long-term and seasonal components of ∆H SD (t G i )/∆H * SD (t G i ), depicted in Figure 10, were used in Equation (15) for modelling ∆H/∆H * . Figure 11 shows ∆H c (t G i )/∆H * c (t G i ) and ∆H SD (t)/∆H * SD (t) models developed using the SD method (see Equation (15)). It also illustrates the differences δ∆H/δ∆H * between ∆H c (t G i )/∆H * c (t G i ) and their models ∆H SD (t)/∆H * SD (t).   The ∆H SD (t)/∆H *SD (t) models are centered to zero (see Figure 11). Thus, in order to determine ∆H M (t)/∆H *M (t), the ∆H SD (t)/∆H *SD (t) had to be referred to the reference epoch t0, i.e., 1 January 2008, that has been used as the realization epoch of orthometric/normal heights. In this study, ∆H M (t)/∆H *M (t) were determined by shifting those ∆H SD (t)/∆H *SD (t) Figure 11. Time series of ∆H c (t G i )/∆H * c (t G i ) (black line) and ∆H SD (t)/∆H * SD (t) models obtained using the SD method (red line), and differences between ∆H c (t G i )/∆H * c (t G i ) and their corresponding ∆H SD (t)/∆H * SD (t) models (green line). Figure 10 indicates that the dominant signal of temporal variations of orthometric/ normal heights comes from seasonal components of ∆H c (t G i )/∆H * c (t G i ). The amplitudes of these seasonal components reach the level of 6.5 ± 0.5 mm. The long-term components of ∆H c (t G i )/∆H * c (t G i ) for the area and time period investigated exhibit a non-linear changing pattern. They fluctuate within the range of ±2 mm. The results presented in Figure 11 indicate a good agreement between ∆H SD (t)/∆H * SD (t) models and ∆H c (t G i )/∆H * c (t G i ) data. The coefficients of correlation between ∆H SD (t)/∆H * SD (t) and the respective ∆H c (t G i )/∆H * c (t G i ) are at the level of 0.97. The standard deviations of δ∆H/δ∆H * are in the range of 1.3 ± 0.1 mm.
The ∆H SD (t)/∆H *SD (t) models are centered to zero (see Figure 11). Thus, in order to determine ∆H M (t)/∆H *M (t), the ∆H SD (t)/∆H *SD (t) had to be referred to the reference epoch t 0 , i.e., 1 January 2008, that has been used as the realization epoch of orthometric/normal heights. In this study, ∆H M (t)/∆H *M (t) were determined by shifting those ∆H SD (t)/∆H *SD (t) considering the offset (τ) that is equal to: 1.9 mm, 1.9 mm, 2.6 mm, and 2.0 mm for the mascons #400, #401, #470, and #471, respectively (see Figure 12). considering the offset (τ) that is equal to: 1.9 mm, 1.9 mm, 2.6 mm, and 2.0 mm for the mascons #400, #401, #470, and #471, respectively (see Figure 12). The GRACE/GRACE-FO level 2 products, i.e., GGMs or mascon solutions, are generally available with a latency of a couple of months from taking measurements by the GRACE/GRACE-FO satellites. Thus, ∆H M /∆H *M obtained from these products for the present time should be predicted. With the use of ΔH c /ΔH *c data from 6 years preceding the prediction period, S/S * and LT/LT * components of ∆H/∆H * were determined with the The GRACE/GRACE-FO level 2 products, i.e., GGMs or mascon solutions, are generally available with a latency of a couple of months from taking measurements by the GRACE/GRACE-FO satellites. Thus, ∆H M /∆H *M obtained from these products for the present time should be predicted. With the use of ∆H c /∆H *c data from 6 years preceding the prediction period, S/S * and LT/LT * components of ∆H/∆H * were determined with the SD method. Several mathematical models, e.g., exponential models, Fourier series, Gaussian models, polynomial models, and sum of sines models, were investigated for their fit to the LT/LT * component of ∆H/∆H * using Matlab [55]. The 7th degree polynomial model was chosen as the best fitted one. The predicted values of ∆H M /∆H *M were determined by combing (1) S/S * component of ∆H/∆H * , (2) values of ∆H/∆H * resulted from the 7th degree polynomial model fitted to LT/LT * component of ∆H/∆H * , and (3) the offset (τ). Figure 13 depicts the example of ∆H M /∆H *M together with their predicted values ∆H P /∆H *P as well as S/S * and LT/LT * components of ∆H/∆H * . It also illustrates the LT/LT * component of ∆H/∆H * obtained using the 7th degree polynomial model. The ∆H P /∆H *P variations presented in Figure 13 were computed for the period of January 2010-June 2010 using ∆H c /∆H *c from January 2004 to December 2009. The ∆H c /∆H *c shifted by τ and ∆H P /∆H *P , as well as the predicted values of S/S * and LT/LT * for this prediction period, were also shown in Figure 13. To obtain ∆H P /∆H *P in different periods of the year, this procedure was repeated seven times, shifting the beginning of the time series by 1 month. Each time, ∆H P /∆H *P were predicted for six epochs coinciding with the respective 6 months ( Figure 14). The differences as well as the differences were obtained. The statistic of the differences δP 1 and δP 2 for the predicted six months are given in Table 1.     The results given in Table 1 and presented in Figures 13 and 14 reveal uneven differences between the predicted ∆H M /∆H *M and their corresponding ∆H c /∆H *c data shifted by the offset τ for all investigated cases. Differences δP 2 for the first 3 months range from −3.1 to 1.4 mm, when predicting the values of ∆H M /∆H *M for 6 months from January 2010 to June 2010 using ∆H c /∆H *c data from January 2004 to December 2009; they range from −2.9 to 0.0 mm, when predicting the values of ∆H M /∆H *M for the period from March 2010 to August 2010 using ∆H c /∆H *c data from March 2004 to February 2010. This may indicate that the fit of the predicted values of ∆H M /∆H *M to ∆H c /∆H *c data shifted by the offset τ strongly depends on the magnitude and the character of ∆H/∆H * within the predicted period. The statistics given in Table 1 indicate that the prediction accuracy in terms of the standard deviations of δP 1 and δP 2 ranges from 0.2 to 1.2 mm, and from 1.2 to 2.3 mm, respectively. It shows that the prediction accuracy of ∆H/∆H * in the case of δP 1 is higher compared to the case of δP 2 . This might be ascribed to the fact that differences δP 2 include not only the error of the predicting procedure, such as δP 1 , but also the error resulting from the modelling procedure.

The Determination of H(t)/H * (t)
The inverse distance to a power interpolation method was used to interpolate/extrapolate the ∆H M (t)/∆H *M (t) referred to the epoch of 1 January 2008 at the sites of the ASG-EUPOS network (see Figure 3).

The Determination of H(t)/H (t)
The inverse distance to a power interpolation method was used to interpolate /extrapolate the ∆H M (t)/∆H *M (t) referred to the epoch of 1 January 2008 at the sites of the ASG-EUPOS network (see Figure 3). Figure 15, showing these changes, exhibits clear variability of H/H * at the ASG-EUPOS sites in both time and space domains. This variability can reach 20 mm from March to September. Figure 15

Conclusions
The paper discusses the determination of orthometric/normal heights H/H * considering their dynamics of both secular and seasonal character. An approach based on the GRACE satellite mission data are proposed to determine such heights. The results obtained from the implementation of this approach over the area of Poland that was chosen as a case study indicated that:

1.
Temporal variations of orthometric/normal heights ∆H/ ∆H * for the period from April 2002 to August 2016, obtained as a combination of temporal variations of geoid/quasigeoid heights ∆N/∆ζ and vertical deformations of the Earth's surface ∆h reach up to 23 mm.

2.
The major part of the signal, i.e., ca. 66%, of ∆H/ ∆H * results from ∆h, while its remaining part is due to ∆N/∆ζ.

4.
The use of the seasonal decomposition method makes possible modelling ∆H/ ∆H * with one millimetre accuracy at the confidence level of 97%; it also makes possible predicting them for the next six months with the accuracy of ca. 1-2 mm.
Overall, the research conducted within the course of the paper emphasizes the need for the H/H * to be corrected for their dynamics. Such heights will be required for fulfilling the contemporary geodetic scientific purposes and high-precision applications associated with the physical heights. For example, ∆H/ ∆H * can particularly be significant to mitigate artifacts and aliasing of repeated levelling measurements. Let us assume that two repeated levelling campaigns were conducted in different epochs over the area of Poland, the first campaign in winter/spring seasons (February-April) and the second campaign in summer/autumn seasons (i.e., August-October). In order to merge and integrate H/H * determined from these levelling campaigns, ∆H/ ∆H * should be considered, otherwise, all levelling measurements associated with the first campaign will be biased by ca. 2 cm with respect to those from the second campaign. Furthermore, H/H * are nowadays, more and more frequently determined by combining ellipsoidal heights from GNSS data with geoid/quasigeoid heights obtained from highly accurate gravimetric geoid/quasigeoid model. Taking into the consideration ∆N/∆ζ and ∆h in Poland, the dynamics of orthometric/normal heights will essentially be needed for the determination of H/H * in Poland when using GNSS solutions combined with a precise geoid/quasigeoid model.
The gravity-dedicated satellite missions, e.g., GRACE, provide valuable information for the determination of H/H * corrected for their dynamics. In particular, they provide unique information concerning long-wavelength components of physical height changes. Since the spatial resolution of ∆H/ ∆H * obtained from GRACE satellite mission data is limited to 3 • × 3 • at the equator, complementary data of high spatial resolution from other sources, e.g., high-resolution hydrological models, would be required to gain more information about short-medium wavelength components, i.e., beyond d/o 60, of ∆H/ ∆H * .