Study on the inter-annual hydrology-induced deformations in Europe using GRACE and hydrological models

Earth’s crust deforms in various time and spatial resolutions. To estimate them, geodetic observations are widely employed and compared to geophysical models. In this research, we focus on the Earth’s crust deformations resulting from hydrology mass changes, as observed by GRACE (Gravity Recovery and Climate Experiment) gravity mission and modeled using WGHM (WaterGAP Global Hydrological Model) and GLDAS (Global LandData Assimilation System), hydrological models. We use the newest release of GRACE Level-2 products, i. e. RL06, provided by the CSR (Center for Space Research, Austin) analysis center in the form of a mascon solution. The analysis is performed for the European area, divided into 29 river basins. For each basin, the average signal is estimated. Then, annual amplitudes and trends are calculated. We found that the eastern part of Europe is characterized by the largest annual amplitudes of hydrology-induced Earth’s crust deformations, which decrease with decreasing distance to the Atlantic coast. GLDAS largely overestimates annual amplitudes in comparison to GRACE and WGHM. Hydrology models underestimate trends, which are observed by GRACE. For the basin-related average signals, we also estimate the non-linear variations over time using the Singular Spectrum Analysis (SSA). For the river basins situated on the southern borderline of Europe and Asia, large interannual deformations between 2004 and 2009 reaching a fewmillimeters are found; they are related to high precipitation and unexpectedly large drying. They were observed byGRACEbutmismodelled in theGLDASandWGHMmod*Corresponding author: Artur Lenczuk, Faculty of Civil Engineering and Geodesy, Military University of Technology, Warsaw, Poland, e-mail: artur.lenczuk@wat.edu.pl, ORCID: https://orcid.org/0000-0002-9573-1302 Grzegorz Leszczuk, Anna Klos, Wieslaw Kosek, Janusz Bogusz, Faculty of Civil Engineering and Geodesy, Military University of Technology, Warsaw, Poland, e-mails: grzegorz.leszczuk@wat.edu.pl, anna.klos@wat.edu.pl, wieslaw.kosek@wat.edu.pl, janusz.bogusz@wat.edu.pl, ORCID: https://orcid.org/0000-0001-6920-3429 (G. Leszczuk), https://orcid.org/0000-0001-6742-8077 (A. Klos), https://orcid.org/0000-0002-1798-7097 (W. Kosek), https://orcid.org/0000-0002-0424-7022 (J. Bogusz) els. Few smaller inter-annual deformations were also observed by GRACE between 2002-2017 for central and eastern European river basins, but these have been also wellcovered by the WGHM and GLDAS hydrological models.


Introduction
The GRACE (Gravity Recovery and Climate Experiment) twin satellites delivered changes in the Earth's gravity field from 2002 until 2017. 15 years of global data coverage enabled many geodetic interpretations in terms of hydrology mass redistribution (e. g. [1,8,34,36,41,50]). Scientists compared GRACE-derived Total Water Storage (TWS) with values obtained from hydrological models to decide on the sensitivity of GRACE to changes of various TWS components or assimilated them into hydrological models to improve the TWS estimates (e. g. [9,39,51]). To compare GRACE and geophysical models' outputs with other geodetic techniques, as e. g. GPS (Global Positioning System), redistribution of masses has to be transformed into the Earth's crust deformations using load Love numbers. GRACE/hydrological models versus GPS comparison has been already widely performed (e. g. [4,6,11,21,25,27,47]). Authors focused mainly on seasonal signals being similar to each other or reduced once GRACE-observed or hydrologically-modeled values are removed from GPS displacements [12,15,23,24]. Klos et al. [19] showed that hydrologically-induced deformations are also observed in GPS displacements in short time periods. Han [16], Pan et al. [26], Tiwari et al. [42], or Wu et al. [45] pointed that hydrology loading affects Earth's crust deformations not only in the seasonal scale but also in the inter-annual scale. These inter-annual deformations driven by hydrology changes are also present in the GPS displacements. If this is the case, they may affect GPS station velocities and their uncertainties if not accounted for.
Examining various phenomena that contribute to the Earth's crust deformations as well as estimating their magnitudes is a challenging task. Once all contributors are properly recognized, the sensitivity of various geodetic techniques can be assessed. As a consequence, systematic errors of those techniques can be also examined. In this research, we decided to focus on the Earth's crust deformations, induced by hydrology loading. Our motivation is based on the fact that providing the magnitudes of inter-annual deformations would be of a benefit for numerous interpretations, where vertical displacements of the Earth's crust are taken into consideration. This may apply for future GRACE versus GPS comparisons, or deciding on the phenomena which induce long-term signals present in the GPS displacements. This would also definitely help in interpretations of GPS station velocities and understanding a broad range of phenomena which evoke their uncertainties.
We focus on the inter-annual signals present in hydrologically-induced vertical displacements provided for the European area. Vertical displacements are estimated using GRACE gravity observations and two hydrological models, namely: Water GAP Global Hydrological Model (WGHM) and Global Land Data Assimilation System (GLDAS) 2.1, Noah version. We estimate mean signals within the major European river basins and subject them to the Singular Spectrum Analysis (SSA) algorithm to estimate non-linear long-term changes. We show that large non-linear signals are present in vertical hydrologicallyinduced displacements for eastern Europe.

Datasets and methodology
We use vertical displacements induced by hydrological changes estimated from three different datasets. Firstly, the GRACE mascon solution is utilized. Then, two hydrological models are employed, namely WGHM and GLDAS. Both differ by a number and type of TWS components they include.

GRACE mascon solution
Mass concentration blocks or so-called mascons represent changes in the Earth's gravity field. They can be used alternatively to spherical harmonic functions, having basically a few advantages over them. First, geophysical constraints applied within the mascon solution support noise filtering. This is more rigorous than empirical methods used to remove north-south stripes present in spherical harmonic solutions. Besides better noise filtering, mascons also prevent information leakage between land and ocean grids. This leads to coastal areas being better represented [33]. Here, we use GRACE Level-2 RL06 (Release-06, RL06) mascon solution provided in the form of monthly gridded TWS values by the Center for Space Research (CSR) in Austin, U. S. [33]. The solution spans from April 2002 until June 2017. The CSR mascon solution is chosen, as it contains more signal variance than other GRACE solutions for corresponding degrees and orders [22]. In the CSR solution, degree 1 corrections were applied using GRACE TN13 values. The C20 coefficients were replaced using estimates from Satellite Laser Ranging (SLR). Glacial Isostatic Adjustment (GIA) effect has been removed using the ICE6G_C (VM5a) model [28].
TWS changes observed by GRACE represent a sum of all water storage components: where components are abbreviated by, respectively, CAN: plant canopy water storage, SN: snow water storage, SM: soil moisture, GW: groundwater, SW: surface water and ICE: water contained in glaciers.
We transform monthly gridded TWS values, provided originally in 0.25°per 0.25°grid nodes, into spherical harmonic coefficients C nm and S nm up to degree and order 120, using [43]: whereĈ nm andŜ nm are dimensionless components explained by Wahr et al. [43] within Eq. (11). The dσ term equals to sin θ dθ dφ and means the basic area element. Constants ρ w = 1000 kg/m 3 and ρ e = 5496 kg/m 3 stand for, a water density and a mean Earth's density, respectively. We use the methodology described by Farrell [10] to estimate Earth's crust deformations induced by changes in mass load. We utilize here only vertical displacements, which are computed by [44]: where R stands for Earth's radius, θ is colatitude, λ is a longitude,P nm are fully normalized associated Legendre functions for n degree and m order,C nm andS nm are spherical harmonic coefficients of the time-variable Earth's gravity field relative to the average gravity field model GGM05S [29], h ὔ n and k ὔ n are the load Love numbers provided by Farrell [10]. Using Eq. (4), we estimate vertical displacements in a 1°per 1°grid for European area, which we define by longitudes between 10°W and 60°E and latitudes between 34°N and 80°N (Fig. 1).
To show the differences between mascon solutions provided for the two consecutive releases: RL05 and RL06 (used here), and, also, between mascon solutions and spherical harmonic coefficients, we use the CSR RL05 mascon solution and the CSR RL06 spherical harmonic coefficients. For the CSR RL05 mascon solution, degree 1 corrections were applied using GRACE TN07 values. The C20 coefficients were replaced using SLR estimates. The GIA effect was removed using Geruo A et al. [13] model. TWS changes derived for the CSR RL05 mascon solution are then transformed into vertical displacements using Eqs. (2)(3)(4). For the CSR RL06 spherical harmonic coefficients, we applied degree 1 corrections using TN13 [40] and we replaced the C20 coefficients using SLR estimates. We also removed the average gravity field using the GGM05S model, to stay consistent with the mascon solution. The DDK3 filter was used to reduce the northsouth stripes [20]. Finally, the GIA effect was removed using Geruo A et al. [13] model. Spherical harmonic coefficients were then transformed into vertical displacements using Eq. (4).
We computed correlation coefficients between vertical displacements estimated in every 1°per 1°grid node for all types of observations we used. Correlation coefficients estimated between vertical displacements obtained from spherical harmonic coefficients and mascon solution in grid nodes are the smallest for coastal areas and they approach 1 inland (Fig. 1). There is a clear ocean-land leakage in spherical harmonic solution. This can be also observed for the Baltic coast, where the smallest correlation coefficients are noticed for Scandinavia. Comparing two mascon solutions, i. e. RL05 versus RL06, we notice the largest differences between them along the Atlantic coast. The difference between both solutions is mainly seen for Portugal and Spain, for which correlation coefficients decrease to 0.75. Differences between RL05 and RL06 CSR mascon solutions arise probably from the definition of a mascon. Mascons from coastal areas that covered both ocean and land in the RL05 version are divided in the RL06 version into two separate areas. Owing to that, better signal representation has been achieved. Moreover, for the RL06 solution, the recent background models for atmospheric and oceanic impact have been employed. Also, degree 1 coefficients and C20 parameter have been re-defined [32,33].

WGHM hydrological model
We make use of the WGHM global hydrological model developed at the University of Frankfurt. Version 2.2 of the WaterGAP model, which we use, contains information on CAN, SN, SM, GW and SW parameters (Tab. 1). Additionally, it is also supplemented with data about irrigation, animal breeding, household water usage (households and small businesses), industry, and power plant layers from five different global models. The model has also been implemented to provide information on ground and surface waters and water exchange between them (GWSWUSE, GroundWater Surface Water USE) [7]. Five different uncertainty sources have been identified in WaterGAP: cli- They were transformed into vertical displacements using Eqs. (2)(3)(4) in a 1°per 1°grid.

GLDAS Noah hydrological model
and are then transformed into vertical displacements using Eqs. (2)(3)(4) in a 1°per 1°grid. These have been estimated for the GRACE period, including also data gaps. In comparison to WGHM, the GLDAS model contains less hydrological layers. It does not include surface water layer, groundwater, and any information on irrigation.

European river basins
We distinguished the main river basins in Europe using the United Nations (UN) Global Compact initiative datasets (Fig. 2). We make use of catchments of areas larger than 50 000 km 2 . In this way, we include 29 basins, from which Rhine, Danube, Dniepr, Vuoksi-Neva, Northern Dvina, Volga, Don, Kura-Ozero Sevan, and Tigris Euphrates rivers are the largest; their areas are larger than 150 000 km 2 .
For each catchment, we estimate a mean vertical displacement based on GRACE RL06 mascon solution, WGHM, and GLDAS datasets (Fig. 3). This is performed by stacking 1°p er 1°grid nodes included in a certain basin and estimating a mean signal plus its standard deviation, showing a spatial variability within the basin. Those mean signals are then analyzed, as described in a Methodology section.

Methodology
To characterize changes of hydrology-induced deformations within individual river basins, we employed the Least-Squares Estimation (LSE) and Singular Spectrum Analysis (SSA) algorithms. The LSE is used to provide the estimates of trend and annual amplitudes within each catchment, basing on the average deformations. For this purpose, we employ the following time series model: where x 0 is the initial value, b is the mean trend, S i and C i are the sine and cosine terms of annual and semi-annual signals.
Then, average deformations for each basin were subjected to the SSA algorithm. The application of SSA to the analysis of geodetic time series has been already extensively described in many papers (e. g. [3,5,14,48] or [52]). The SSA allows estimating of signals which parameters may vary over time. It means that trends estimated with the SSA will no longer be linear, which is not provided with the LSE method. Also, amplitudes and phases of seasonal curves may change from year to year.
The SSA algorithm consists of two basic steps: decomposition and reconstruction. In a decomposition step, the vector of observations {X(t) : t = 1, . . . , N} is decomposed into a few sub-vectors of the length M. M stands here for the length of the moving window which slides along the whole time span of data, i. e. month-by-month. In this way, the trajectory vector is obtained:   Then, the lag-covariance matrix C x is computed using the Broomhead and King [3] algorithm: where D is the product of the trajectory matrix. Further, the lag-covariance matrix C x is decomposed into a set of eigenvectors and eigenvalues. Eigenvectors are also known as the Empirical Orthogonal Functions (EOFs). A sum of all eigenvalues allows reconstruction of the total variance of the original observations. Eigenvalues are sorted in a decreasing order. Then, eigenvectors which correspond to those eigenvalues are estimated and Principal Components (PCs) are obtained as: Since eigenvalues were sorted in a decreasing order, PCs we estimate are sorted according to the amount of total variance of the original observations they explain, e. g. the first PC explains the maximum amount of variance. Based on the PCs and EOFs, Reconstructed Components (RCs) are estimated using: We tested various lengths of time moving windows and we finally adopted a 2-year moving window. This length of time window enables that inter-annual signals of periods longer than 2 years are reliably estimated. For longer time windows, inter-annual signals are not well pronounced, i. e. signals are too smoothed.
Principal components PCs are sorted according to the decreasing amount of variance they explain, i. e. the first PC explains the maximum variance of the original time series. We applied the Welch periodogram to derive periods of individual PCs. For the European river catchments, inter-annual and annual signals are always contained in the first four PCs. It means that these signals are the most powerful signals contained in the vertical deformations induced by hydrology changes. The first four PCs are then combined, according to the period they represent. In this way, annual and inter-annual reconstructed components RCs are obtained.

Results and discussion
Vertical displacements induced by hydrological loading are not constant over time, which can be noticed from the average signals estimated for individual catchments. The ones we present in Figure 3 were chosen randomly. The Loire river basin is situated in western Europe. For this basin, annual amplitudes estimated for GRACE, WGHM and GLDAS Noah are similar to each other with no clear long-term trend being noticed. Danube and Volga rivers are situated in central and eastern parts of Europe. For these, the differences in annual amplitudes from year to year can be noticed. There also exist changes in a longterm trend, with evident jumps occurring in 2006 and 2010 for, respectively, Danube and Volga rivers. From this simple plot, we can conclude that the time-variability of hydrology-induced deformations is changing depending on the area. Central and eastern parts of Europe are known to be affected by larger signals related to continental hydrosphere than western coastal parts [39]. Hydrological signals are more prominent for eastern areas in terms of either annual amplitudes or time-variability of the long-term trend, which can be noticed from Figure 3.
The above statement is also supported by the rootmean-square (RMS) values estimated for average signals provided for individual basins (Fig. 4). The RMS values are quite similar for GRACE-and WGHM-derived deformations, while they are larger for those provided for the GLDAS hydrology model. The RMS values are significantly larger for eastern Europe than they are for remaining areas. It was also presented by Scanlon et al. [35] or Springer et al. [39]. The RMS values are mainly influenced by the annual signal which overwhelms and covers other signals within displacement time series; note the similarity between results presented in Figure 4 and Figure 5. It can also be observed between results presented in Figure 4  and Figure 6, especially for southeastern Europe, where the trend estimated from GRACE-derived displacements is well-pronounced.
Annual amplitudes and their distribution over the European area are similar to the results presented before by Jin [18] and Scanlon et al. [37]. If the GRACE solution is treated as a reference, we notice that the WGHM model overestimates annual amplitudes for central Europe (Fig. 5), noticed before by Scanlon et al. [36]. It also underrates changes observed by GRACE in basins for southern European areas, as Tigris Euphrates, Kura-Ozero Sevan, Sakarya, and Kizilirmak river basins. Land surface models (LSMs), i. e. GLDAS, overestimate seasonal amplitudes in high latitudes [35,37], compared to the WGHM Table 2: Trend values (in mm/yr) estimated for European river basins with the LSE approach. Grey shaded rows correspond to the "red" and "yellow" rivers from Figure 7. estimates. This is visible in Figure 5, where the GLDAS Noah model significantly revalues annual amplitudes by even 3 mm for east European river basins. This revaluation of amplitude values in the GLDAS model is attributed to the combined effect of mismodelling of SN in winter and SM in summer. This might be also related to the overestimates of evapotranspiration in GLDAS Noah compared to other GLDAS versions, as CLM, VIC, and Mosaic [49]. Scanlon et al. [37] showed that previous versions of the GLDAS model agree better with GRACE observations. Linear trends estimated for GRACE-derived displacements with the LSE approach vary from −0.06 to +0.91 mm/yr (Fig. 6, Tab. 2). Positive values relate to lesser water stored in the catchment, while negative values correspond to more water being accumulated. Almost all river basins we analyze are characterized by positive trends. The largest values (>0.2 mm/yr) are observed for Danube, Dniestr, Northern Dvina, Volga, Dniepr, Kuban, Kizilirmak, Sakarya, Don, Kura-Ozero Sevan, and Tigris Euphrates; all these river basins are situated in central and eastern Europe. These large positive trends indicate that coastal areas of the Black Sea and Arabian Peninsula are drying, which is also showed by Scanlon et al. [37]. A positive trend values are the result of abnormal groundwater depletion and droughts [15,31] and by precipitation changes in eastern Europe [17,31]. These values are similar to the ones presented by Humprey et al. [17], Jin [18], Rodell et al. [31] and Scanlon et al. [35,36].

RIVER BASIN
Trends found by the LSE for GRACE displacements are underestimated by both hydrological models we employed. WGHM trends are slightly negative for almost all river basins we considered. Significantly negative trends derived from the WGHM model for Finland and Russia are caused by the overestimation of precipitation [37]. Only small positive trends are retrieved for the east coast of the Black Sea and the Arabian Peninsula. For the GLDAS model, the European area is divided into halves. Central and eastern parts are characterized by positive trends, while western areas are characterized by negative trend values. From this perspective, we can say that geophysical models underestimate trends derived from the GRACE observations. This may limit the comparison between them, and also a combined comparison to other geodetic techniques, as GPS.
Trends presented so far were estimated using the LSE approach assuming their time-constancy. To get an idea of how these may differ and change over the years, we use the SSA algorithm. The inter-annual signals we obtain are a combination of corresponding reconstructed components retrieved from EOFs. These are sorted according to the variance of the original time series they explain. Basing on that variance, we can state that the percentage of time series variance explained by inter-annual signals for the GRACE-derived vertical displacements exceeds the variance explained by the annual signal for basins: Dniepr, Dniestr, Don, Guadiana, Kuban, Kura-Ozero Sevan, Tigris Euphrates, and Volga. In other words, for these basins inter-annual changes are more powerful than the annual signal itself, i. e. they explain the maximum amount of variance of the original time series. This is however no longer observed for WGHM and GLDAS, for which the largest part of time series variance is explained by the annual signal. There are two exceptions for the WGHM displacements; Kura-Ozero Sevan and Tigris Euphrates basins have prominent inter-annual changes. Figure 7 presents a comparison of inter-annual signals for European river basins we analyze. Plotted are GRACE, WGHM, and GLDAS estimates. Comparing all three plots, it is obvious that GRACE-derived displacements are characterized by the largest inter-annual signals and, also, the most apparent trends. There is a distinct drop in the vertical displacements noticed between 2004 and 2006 (Fig. 7). Then, Earth's crust is recovering and strongly uplifting up to 2009. This is observed for all eastern European rivers, i. e. Tigris Euphrates, Kuban, Kura-Ozero Sevan, Don, Kizilirmak, Sakarya, Dniestr, Danube, and Dniepr. Central European rivers, i. e. Volga, Vuoksi-Neva, Narva, Western Dvina, Neman, Wisla, and Oder, have similar down-and uplifting, but this recovers faster of three years; it lasts from 2004 until 2006 (Fig. 7). To prove that long-term trends are mismodelled by both WGHM and GLDAS, we use the same colors and mark those rivers on WGHM and GLDAS plots. From these, we can notice that the large down-and uplifting of the Earth's crust, which has been retrieved from GRACE is not covered by any of the two models we analyzed. Also, trends that dominate GRACE displacements are mismodelled. Few inter-annual peaks observed by GRACE in 2002, 2006, 2012, and 2014 are covered well by both WGHM and GLDAS models. As they are retrieved by GLDAS, it may be concluded that these variations are mostly generated by changes in SM. Other compartments included in the GLDAS should not be so prominent. CAN is constant and including this compartment brings no difference to the TWS values and it is also hardly detectable by GRACE. The SN may be considered only for northeastern areas of Europe, as Onega or Northern Dvina river basins.

Conclusions
In this study, we showed the analysis of hydrologyinduced deformations over the European area. Both GRACE-observed and WGHM-and GLDAS-modelled deformations were analyzed. We found that GLDAS overestimates the seasonal changes, in comparison to GRACE observations and WGHM model. This has been also already noticed by Blitzkow et al. [2] and Sliwinska et al. [38]. Both GLDAS and WGHM models underestimate trends observed by GRACE, shown also by Scanlon et al. [36], but for a global scale. We used the SSA algorithm to estimate inter-annual signals in the displacement time series from GRACE, WGHM, and GLDAS. For that, we observed that the eastern part of Europe is characterized by large long-term non-linear variations present in the displacement time series, captured by all three datasets. We discovered, however, that GRACE observations contain also large subsidence and recovering between 2004 and 2009, which is unmodelled in both WGHM and GLDAS. This is probably related to the high precipitation period, happening in Europe in 2004. Once this precipitation period stopped in 2005, a constant uplift is observed, up to present.
Analysis of such long-term non-linearities is possible thanks to a long and reliable time series provided by the GRACE mission. These time series are compared with other techniques, as GPS, to decide whether these are sensitive to hydrological mass changes or not. Inter-annual signals will limit the commonly performed GRACE versus GPS comparison, when not accounted for. It also needs to be investigated whether GPS European stations are sensitive to those inter-annual deformations.