## 2 Materials and Methods  
Four precision leveling campaigns have been carried out in Poland (#1: 1926-1937, #2: 1953-1955, #3: 1974-1982, #4: 1997-2003). Most of the collected data covers the territory of Poland after World War II. The first campaign was weakly documented and the relevant data were not used in this study: The remaining data were used to develop maps of the vertical crustal movements (Table 2, datasets A and B). The vertical movements were determined at points in the double-leveling network (vertical movement network). Successive movements were determined at network points based on data from the second, third, and fourth campaign (Table 2, dataset C), which determined the relative vertical movements of the Earth's crust based on data from permanent GNSS stations (with a life cycle of 5 years) (Table 2, dataset D). For this article, the first joint alignment of vertical motion networks based on leveling and GNSS data was carried out. (Table 2, dataset E). All values in the analyzed datasets relate to the same point with a constant velocity. The same adjustment method was used to facilitate the analysis and compare the generated models. Datasets differed in their spatial and temporal resolutions (Figure 2). \begin{table}
\begin{tabular}{c c c c c} \hline \hline \multirow{2}{*}{**Data Set**} & \multirow{2}{*}{**Leveling Data**} & \multicolumn{2}{c}{**GNSS Stations**} & \multicolumn{1}{c}{**Number of Nodal**} & \multicolumn{1}{c}{**Number of Leveling**} \\  & & **Data** & **Points** & **Campaign** \\ \hline A & + & – & 98 & \#2, \#3 \\ \hline B & + & – & 222 & \#3, \#4 \\ \hline C & + & – & 228 & \#2, \#3, \#4 \\ \hline D & – & + & 123 & – \\ \hline E & + & + & 345 & \#3, \#4 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Characteristics of datasets used in the analysis. GNSS: Global Navigation Satellite System. Figure 2: Graphical presentation of the point distributions of the datasets. The spatial distribution of points in each dataset was determined by the density of the double-leveling networks (A, B, C) and the location of permanent GNSS stations (D). Due to the shape of the leveling network, dataset A was composed of sparsely distributed points on leveling lines. The remaining datasets had a higher spatial resolution. Datasets C and E were used as the reference. A fully functional model can be derived from a limited set of representative data with the use of kriging interpolation [30]. The interpolated values are determined based on the correlations between the value and the location of the analyzed points. Kriging differs from other interpolation methods in that it relies on interpolation algorithms to identify, describe, and eliminate trends from representative data [23]. Kriging is a non-uniform method. The same data can be processed in many ways to produce radically different results that highlight and expose trends and information according to the preset parameters. Kriging methods can be classified based on the sampling interval for interpolation (local kriging, global kriging) and their application (ordinary kriging, block kriging, indicator kriging, and co-kriging, i.e., kriging with an additional variable). Simple kriging, ordinary kriging, and universal kriging methods are the most popular. Similar to least-squares collocation, kriging relies on the assumptions of the least-squares adjustment method. In kriging, the estimated value is expressed as the weighted average of the observation vector, and weights are assigned based on the distance between observations and the distance between observations and the location of the interpolated points [23]. This solution is used to identify significant discontinuities and their parameters, to determine the predicted value, and to evaluate the estimated value (variance estimation). In co-kriging, various datasets can be merged and analyzed in terms of their characteristic features. The relationships between random variables are typically determined with the use of variant-covariant matrices.