## 1 Introduction  
The vertical movements of the Earth's crust are caused by environmental, tectonic, geological, and climatic changes. Cartographic models of the vertical movements of the Earth's crust are used in economics, Earth sciences, and reference systems. Data for cartographic models are acquired through the use of direct (geodetic measurements) and indirect (geological, geophysical, geomorphological, and hydrological data) techniques [1; 2; 3; 4; 5]. Geodetic datasets differ in their sources, measurement methods, spatial and temporal resolutions, levels of reference, data processing methods, and data trends. Cartographic models are also influenced by the applied interpolation method [6]. The choice of interpolation method is particularly important when different datasets are used [7]. The main aim of this study was to evaluate the applicability of the co-kriging method for modeling the vertical movements of the Earth's crust based on data acquired with the use of precision leveling techniques and measurements conducted using permanent GNSS stations. The advancements in measurement techniques, the development of measurement systems with a longer life cycle, and higher spatial and temporal resolution of measurements support the determination of the vertical movements of the Earth's crust with unprecedented accuracy. The oldest datasets (Figure 1) originated from precision leveling campaigns that were carried out every 20 years. These data were used to develop vertical movement networks. The node points in vertical movement networks are evenly distributed and mutually interconnected. Vertical movements are observed along leveling segments (lines) with a length of several dozen to several hundred meters. The vertical movements of the Earth's crust are usually observed relative to a point on the Earth's surface that is fixed or has a \"known\" velocity. The advantage of this method is that it minimizes the influence of local terrain deformations. Models of the vertical movements of the Earth's crust based on daily data from permanent Global Navigation Satellite System (GNSS) stations are the most widely used (Figure 1). These movements represent the absolute motion relative to an ellipsoid of revolution. The quality of GNSS data is determined mainly by the station's operating time (minimum of 3 years) and how data are processed by the system operator. The calculated velocities can represent local movements when GNSS stations are separated by a considerable distance (several dozen kilometers) and when they are not interconnected (do not form a network). The vertical movements of the Earth's crust are also analyzed based on satellite images, laser measurements, radar interferometry, altimetry measurements, gravimetry missions, and tide gauge data [8; 9; 10; 11; 12; 13; 14]. The combination of leveling data and GNSS data from the permanent stations is difficult to achieve due to various systems and working methods used [15] and the application of various solutions in the process of creating vertical crustal movement models, ranging from mathematical models to geophysical models [16]. The largest possible amount of available data (leveling and GNSS data) should be used to increase the spatial resolution of the network points. Statistical and mathematical methods for adjusting several networks, as well as interpolation techniques, are used for this purpose [17]. All datasets must have a network structure to be processed with the use of statistical methods. As mentioned earlier, vertical movements based on GNSS data do not have a network structure. In our previous study [18], the vertical movements of the Earth's crust were determined based on GNSS data as relative movements rather than absolute movements. A network structure based on Delaunay triangulation was proposed, where the observations were based on the vertical movements between GNSS stations that constitute triangle vertices [19]. This approach was adopted to adjust the network.