Combination of Soviet-Era Surface Gravity and Modern Satellite Data for Geoid Model Computation: A Case Study for Kazakhstan

1. Introduction

The precise determination of the geoid, which is defined as the equipotential surface that closely approximates the mean sea level, represents a fundamental challenge in geodetic science. It is of prime importance for various scientific, engineering, and geophysical applications. Geoid modeling involves the integration of global geopotential models with terrestrial gravity anomalies and local topographic data to achieve high accuracy [1]. Traditionally, gravimetric geoid determination has relied on extensive ground-based gravity surveys, which can be resource-intensive and geographically sparse in certain regions [2]. Consequently, the integration of Global Navigation Satellite System (GNSS) and leveling data has emerged as a crucial approach, offering a direct means to evaluate and refine geoid models by providing geometric geoid undulations at specific points [3].

A high-accuracy geoid model offers a cost-effective alternative, enabling direct conversion of ellipsoidal heights obtained from GNSS derived ellipsoidal heights into orthometric heights, without the need for extensive leveling campaigns [4]. Such an approach is essential for supporting a wide range of geoscientific and engineering applications, including large-scale infrastructure projects, cadastral surveying, environmental monitoring, and national mapping programs [5,6]. Despite significant advancements in satellite gravimetry and computational techniques, achieving greater geoid accuracy remains a persistent challenge, particularly in regions with complex topographical and gravimetric characteristics. Expensive field surveys in remote and rugged terrain, along with a lack of GNSS stations, often result in inadequate data density for robust geoid determination [7,8]. Moreover, temporal variations of the Earth’s gravity field, caused by mass redistribution processes such as ice sheet melting or hydrological changes, further necessitate the requirement for continuous monitoring and periodic updates of geoid models to maintain their precision over time [9,10,11].

Regional geoid models are particularly critical for practical applications within specific geographical areas, where local geodynamic processes and geological structures significantly influence the gravity field [12,13]. They offer enhanced resolution and accuracy compared to global models, addressing the limitations imposed by sparse data in challenging terrains like mountainous regions [14]. Global geopotential models (GGMs) like EGM2008 [15], EIGEN-6C4 [16], and XGM2019e [17] provide adequate accuracy for many regional and global-scale applications, they cannot fully capture the short-wavelength gravity field features associated with geographies having complex terrain. Consequently, their direct use for high-precision engineering and surveying in the region can lead to systematic biases and reduced accuracy, particularly in mountainous areas [18]. Developing a dedicated local geoid model mitigates these shortcomings by incorporating detailed local gravity observations, refined topographic modeling, and advanced computation methods tailored to the country’s geophysical characteristics.

In addition to the geodetic benefits, establishing a local geoid model carries strategic significance. It ensures independence from foreign datasets, supports national sovereignty in spatial data management, and underpins all engineering and infrastructure projects with a unified and precise vertical reference system [19]. For large-scale projects such as the construction of transportation corridors, hydrotechnical structures, and energy infrastructure high-accuracy orthometric heights derived from a local geoid model improve design reliability and reduce construction risks [20].

Modern geoid modeling techniques, including least-squares collocation and the least-squares modified Stokes’ formula (LSMS), provide the necessary tools to achieve the required accuracy under complex terrain conditions [21,22]. The LSMS approach, in particular, offers flexibility in incorporating heterogeneous gravity datasets, topographic corrections, and stochastic modeling of residual signals, making it well-suited for mountainous and seismically active regions [23].

Kazakhstan is a country characterized by its vast and diverse topography, encompassing flat plains, deserts, and the high mountains of the Tian Shan range, which presents unique challenges for precise geoid determination. Furthermore, the country's limited coverage of uniformly distributed gravimetric data, coupled with the sparse availability of GNSS/leveling benchmarks, exacerbates the complexity of establishing a high-precision vertical datum [24]. The complex terrain and substantial variations in elevation necessitate a robust regional geoid model that can accurately reflect the intricate gravity field behavior across the region. Although global models provide a starting point, their accuracy remains insufficient for engineering and surveying in such a geodynamically complex region. Despite the importance of a unified vertical datum, Kazakhstan currently lacks a dedicated, validated national geoid model. This deficiency hinders various geodetic and practical applications, making the development of such a model a pressing necessity for the country's infrastructural development and scientific advancement [25].

To address this gap, the present study develops the first high-accuracy gravimetric geoid model for Kazakhstan, termed QazGM2025. The model is computed using the Least-Squares Modification of Stokes’ formula (LSMSA/KTH method), which allows the optimal integration of heterogeneous datasets. Specifically, we combined digitized Soviet-era gravimetric maps, recently acquired absolute gravity measurements, GNSS/levelling benchmarks from the national geodetic network, and high-resolution topographic data (FABDEM). The model is evaluated against independent GNSS/levelling data and benchmarked against leading global models (EGM2008, EIGEN-6C4, and XGM2019e). This study not only delivers a precise reference surface with ~6 cm accuracy but also establishes the foundation for Kazakhstan’s modern vertical datum, enabling robust support for surveying, engineering, and geoscientific applications.

2. Study Area

The focus area is the Republic of Kazakhstan (Figure 1), which, despite encompassing vast expanses of lowland terrain, is characterized by highly diverse topography ranging from extensive steppe and desert plains to some of the highest mountain ranges in Central Asia. The western and northern parts of the country are dominated by the Caspian Depression and the West Siberian Plain, while the central Kazakh Uplands feature hilly and plateau landscapes of complex geological structure. The southeastern and southern margins of Kazakhstan include the towering mountain systems of the Tien Shan and Zhetysu Alatau, with elevations exceeding 7000 m, presenting significant challenges for precise geoid modeling due to steep relief gradients and active tectonics [26]. Pression and the West Siberian Plain, while the central Kazakh Uplands feature hilly and plateau landscapes of complex geological structure. The southeastern and southern margins of Kazakhstan include the towering mountain systems of the Tien Shan and Zhetysu Alatau, with elevations exceeding 7000 m, presenting significant challenges for precise geoid modeling due to steep relief gradients and active tectonics [26].

The geographic boundaries of Kazakhstan extend from a minimum latitude of approximately 40° N to a maximum latitude of about 56° N, and from a minimum longitude of roughly 46° E to a maximum longitude of about 88° E. However, to minimize edge effects during geoid modeling for the main focus area, the computation region was extended beyond the national boundaries, covering the approximate range from 38° N to 58° N in latitude and from 44° E to 90° E in longitude. This expanded area ensures that boundary effects are properly accounted for, allowing for a more accurate representation of the geoid across Kazakhstan’s complex and varied terrain [27].

3. Methodology

The geoid modeling methodology adopted in this study integrates multi-source gravimetric data and advanced computational techniques to achieve high precision, particularly in areas with complex topography. Figure 2 illustrates the structured workflow for geoid determination and evaluation using a hybrid approach that combines global gravity models with regional terrestrial data. We used the Least -Squares Modification of Stokes' formula to optimally integrate diverse datasets, ensuring a robust and accurate geoid determination. The process begins with the digitization of historical gravimetric maps, including datum transformation from Pulkovo to WGS84, application of systematic corrections, validation against terrestrial gravity observations, and integration with global datasets. GNSS/levelling data at first- to third-order levelling benchmarks are employed for input data evaluation, optimisation of Stokes’ kernel modification parameters, and final model validation. Normal heights of control points are converted to orthometric heights. The optimal GGM and digital terrain model are selected based on statistical consistency with GNSS/levelling data.

The geoid is computed with additive corrections - topographic, atmospheric, ellipsoidal, and downward continuation. The final stage comprises rigorous accuracy assessment using independent GNSS/levelling benchmarks and trend removal of systematic effects.

3.1. Data Collection and Database Compilation

3.1.1. Gravity Data

The regional gravity database compiled for QazGM2025 incorporates a heterogeneous collection of historical gravimetric data for Kazakhstan and neighboring countries, as well as contemporary terrestrial observations from the Modern Gravimetric Framework of Kazakhstan, where absolute gravity observations were conducted in 2023–2024 by the National Centre for Geodesy and Spatial Information [28].

In this study, the core dataset comprises digitized Soviet-era gravimetric maps at a scale of 1:200,000, providing dense coverage across nearly 90% of Kazakhstan (Figure 3). Priority areas, particularly those of strategic or scientific importance, were surveyed at finer resolutions (1:50,000 and larger), ensuring enhanced detail in complex terrains. The maps were originally compiled from ground-based gravity surveys conducted primarily with SN-3, GAK, and GR-2K gravimeters, which were widely used during that period for regional-scale geodetic campaigns. For standardization, the gravity observations were reduced to Bouguer anomalies using intermediate crustal densities of 2.30 g/cm³ or 2.67 g/cm³, accompanied by records of free-air anomalies, terrain corrections, and observation heights [29]. This systematic framework provided a consistent baseline for the construction of Kazakhstan’s first digital gravity database.

The maps were digitized through a multi-stage workflow to ensure positional and numerical accuracy. Original sheets were scanned at 600 dpi, georeferenced to the Pulkovo 1942/Gauss–Krüger system, and subsequently vectorized to extract gravity anomaly values, observation heights, and terrain corrections. A preliminary anomaly surface was generated via Kriging interpolation (250 m grid) to identify and correct discrepancies, after which the dataset was transformed into WGS84 coordinates using standard transformation parameters (EPSG:15865). Preference was given to Bouguer reductions based on 2.67 g/cm³ density, while 2.30 g/cm³ data were used where necessary.

Although the digitized terrestrial data provide dense coverage across most of Kazakhstan, the mountainous and border regions remained undersampled. To fill these gaps and maintain spatial continuity, global gravity anomalies from the WGM2012 model (5′ × 5′ resolution) were incorporated (Figure 4). As a result, the final gravity database consists of 651,862 terrestrial measurement points (0.24 points/km², ~13.8 per 5′ grid cell), complemented by 103,854 global points from WGM2012, yielding a total of 755,716 data points. This integration ensured both dense coverage and improved homogeneity across the entire study area, forming the basis for subsequent geoid computation.

To assess the reliability of the digitized historical dataset, we compared it with recent absolute gravity measurements obtained from the Modern Gravimetric Framework of Kazakhstan. This network includes 76 absolute stations and 142 nearby relative stations (within ~50 km), observed using a Micro-g LaCoste A10 absolute gravimeter and Scintrex CG-6 Autograv relative gravimeters. The expanded uncertainty of these measurements was below 10 μGal at all 218 stations. Statistical comparison revealed a mean difference of only 0.66 mGal between the reconstructed values from the Soviet-era maps and the modern absolute observations, confirming the overall consistency of the historical dataset for use in regional geoid modeling.

3.1.2. GNSS/Levelling Points

A dataset of 3,653 GNSS/levelling points, evenly distributed across Kazakhstan (Figure 5), was selected as the primary reference for validation of all datasets used in this study. These control points were established under the National Spatial Data Infrastructure (NSDI) project, which led to the creation of the official national reference frame QazTRF-23, replacing the outdated SK-42 system with an ITRF-consistent framework [28].

From approximately 45,600 GNSS stations surveyed nationwide, over 10,000 were associated with geometric levelling benchmarks. However, because the majority of these belonged to the 4th accuracy class, unsuitable for high-precision applications, they were excluded. After filtering, only benchmarks along first- and second-order levelling lines, supplemented by a subset of reliable third-order stations, were retained. This selection resulted in 3,653 points, comprising 243 first-order, 974 second-order, and 2,436 third-order stations.

The geodetic framework follows the “general-to-specific” principle, transferring high-accuracy coordinates from 11 Fundamental Astronomical and Geodetic Network (FAGN) stations to the Astronomical and Geodetic Network (AGN), Geodetic Densification Network (GDN), and the State Levelling Network (SLN). At FAGN stations, continuous 24-hour GNSS sessions with dual-frequency receivers were processed using Bernese GNSS Software, incorporating precise satellite ephemerides, Earth orientation parameters, tidal models, and IGS reference stations. The adjustment achieved absolute accuracies of 0.10–0.15 m, with relative accuracies better than 0.02 m horizontally and 0.03 m vertically. Validation against AUSPOS online GNSS processing services confirmed differences within 0.015 m horizontally and 0.03 m vertically.

By establishing this unified and high-accuracy coordinate framework, Kazakhstan’s geodetic network provides a robust basis for geoid model evaluation. Importantly, the GNSS/levelling benchmarks were also used to validate the DEM and to assess the performance of global geopotential models, ensuring consistency in the validation framework across all datasets employed in this study.

3.1.3. Digital Elevation Model

For the topographic component of the geoid modeling, the FABDEM dataset was used. FABDEM is a global digital elevation model derived from the Copernicus GLO-30 DEM [30,31,32], corrected to represent the bare-earth surface by removing vegetation and building artifacts [33]. It provides a spatial resolution of 1 arc-second (~30 m), enabling accurate representation of terrain variations across Kazakhstan, including the steep gradients of the Tien Shan and Zhetysu Alatau ranges [34].

The quality of FABDEM over the study area was assessed using GNSS-derived elevations from control points and EGM2008 geoid height. The comparison showed a mean bias of 0.397 m, a standard deviation of 1.365 m, and an RMSE of 1.422 m. These results confirm that FABDEM provides sufficient vertical accuracy for use in computations of the geoid model.

3.1.4. Global Geopotential Models

To establish the optimal long-wavelength gravity field for geoid modeling in Kazakhstan, we evaluated several recent high-resolution global geopotential models (GGMs) developed by the International Centre for Global Earth Models [35]. These included XGM2019e_2159, an enhanced combination model incorporating GRACE, GOCE, and terrestrial gravity data up to degree/order 2190 [17]; SGG-UGM-2, derived from GOCE satellite gravity gradients combined with surface gravity [36]; GECO, a hybrid model integrating satellite-only solutions with EGM2008 long-wavelength information [18] EIGEN-6C4, which combines GRACE and GOCE data with extensive surface gravity observations up to degree/order 2190 [16]; and the widely used EGM2008, which incorporates both satellite and terrestrial data to the same spectral resolution [15].

These models were tested against GNSS/levelling benchmarks distributed across the study area to identify the model that shows the highest consistency with local observations and could therefore serve as the reference for subsequent geoid computations. For each benchmark, height anomalies derived from GNSS/leveling were compared with the corresponding model-predicted values.

3.2. Geoid Model Development & Validation

3.2.1. Least – Squares Modification of Stokes’ (LSMSA) Approach

The gravimetric geoid is commonly determined using Stokes’ formula, which relates gravity anomalies on the Earth’s surface to geoid undulations [37]. However, direct application of Stokes’ integral requires globally complete and error-free data, a condition that is rarely met in practice. In addition, truncating the integral to regional domains introduces errors due to unmodelled long-wavelength signals [2,6,38] .

To overcome these limitations, the Least-Squares Modification of Stokes’ Formula (LSMSA or KTH method) was developed [39]. This approach modifies the classical Stokes kernel by introducing least-squares optimal coefficients that minimize the expected error in the estimated geoid height [40,37,41,42]. The method achieves this by optimally balancing contributions from three main sources: the long-wavelength gravity field provided by a Global Geopotential Model (GGM), the medium- and short-wavelength gravity anomalies derived from terrestrial and airborne observations, and the residual terrain effects obtained from high-resolution digital elevation models (DEMs).

In this framework, the geoid undulation $\widehat{N}$ at a given location is expressed as:

$$\widehat{N}= \tilde{N}+ {\delta N}_{comb}^{Top}+ {\delta N}_{DWC}+ {\delta N}_{comb}^{Atm}+{\delta N}_{ell}$$

(1)

where, $\tilde{N}$ is a fundamental term representing the primary approximation of the geoid height, obtained from Stokes’ formula using gravity anomalies referenced to a mean spherical surface before applying specific correction terms. ${\delta N}_{comb}^{Top}$ accounts for topographic correction, ${\delta N}_{DWC}$ corrects for the difference between gravity anomalies observed at the Earth’s surface and those referenced to the geoid level, mitigating distortions in mountainous regions. ${\delta N}_{comb}^{Atm}$ and ${\delta N}_{ell}$ are parameters accounting for atmospheric and ellipsoidal corrections.

The key principle of LSMSA is the modification of Stokes’ kernel to minimize the propagation of data errors into the computed geoid heights [43].This is achieved by estimating optimal modification parameters from the empirical covariance functions of residual gravity anomalies, fitted to analytical models. The modification degree, spectral weighting factors, and combination coefficients are then determined to balance the spectral content of terrestrial anomalies and the long-wavelength information from the GGM.

Once the modification parameters are established, the gravimetric geoid is computed directly by applying the least squares optimized Stokes’ integral to the prepared set of terrestrial gravity anomalies. The formulation inherently allows for the rigorous incorporation of additive corrections, which compensate for effects not modeled in the main integral [44,45].

These features make LSMSA particularly well-suited for Kazakhstan’s geodynamically complex terrain, where data coverage is irregular and large elevation variations necessitate careful treatment of residual terrain effects. All stages of geoid computation described herein were implemented within the LSMSsoft software environment [1]. This sophisticated approach ensures that the resulting geoid model is not only accurate but also robust against inherent data limitations and topographic complexities [46,47].

3.2.2. Parameter Testing Approach

A critical aspect of the LSMSA method is the selection of spectral and stochastic parameters that govern the balance between contributions from the global geopotential model and terrestrial gravity anomalies. Two parameters are particularly important: the maximum spherical harmonic degree/order (L, M) used to truncate the Stokes’ integral, and the a priori error variance of the terrestrial gravity anomalies (C(0)).

The truncation degree controls the spectral resolution of the residual gravity field included in the modified Stokes’ integral. Lower degrees emphasize long-wavelength signals but risk omitting shorter-wavelength information, while higher degrees capture finer detail but may amplify observational noise. To identify a suitable compromise, we tested a range of truncation degrees from 180 to 760, in line with previous studies of regional geoid modelling [3,41,48].

The parameter C(0) represents the stochastic characterization of gravity anomaly errors, influencing the relative weight assigned to terrestrial gravity data in the least-squares solution [49,50]. To account for different assumptions of observational uncertainty, values of 1, 3, and 9 mGal² were considered, consistent with the accuracy level of the digitized Soviet-era gravimetric dataset and its validation against modern absolute gravity measurements.

For each combination of parameters, the LSMSsoft software was used to compute preliminary geoid solutions. Their performance was subsequently evaluated against GNSS/levelling benchmarks, and the optimal configuration was selected based on statistical measures of residual differences.

3.2.3. Geoid Computation

The development of the QazGM2025 geoid model followed a structured workflow designed to ensure consistent integration of heterogeneous datasets within the LSMSA framework. All computations were carried out using the LSMSsoft software package [1], which implements the least-squares modification of Stokes’ formula together with additive corrections and covariance modeling.

The computation began with the assembly of input datasets, including terrestrial gravity anomalies from the digitized Soviet-era maps (supplemented by WGM2012 in undersampled regions), the selected global geopotential model (XGM2019e) for the long-wavelength gravity field, the FABDEM digital elevation model for terrain effects, and GNSS/levelling benchmarks for external evaluation.

The terrestrial gravity anomalies were first pre-processed to remove the GGM contribution, producing residual anomalies suitable for modified Stokes integration. Additive corrections were then applied to account for topographic, atmospheric, and ellipsoidal effects, as well as for downward continuation of anomalies from the Earth’s surface to the geoid. The corrected anomalies were integrated using the LSMSA kernel, truncated at the tested degrees and weighted according to the assumed gravity error variances, to obtain preliminary geoid surfaces.

Finally, the resulting geoid models were compared with GNSS/levelling benchmarks to assess consistency and to detect possible systematic trends. The best-performing parameter configuration was retained for the computation of the final QazGM2025 model.

3.2.4. Validation of Geoid Model

The development of the QazGM2025 geoid model followed a structured workflow designed to ensure consistent integration of heterogeneous datasets within the LSMSA framework. All computations were carried out using the LSMSsoft software package [1], which implements the least-squares modification of Stokes’ formula together with additive corrections and covariance modeling.

This rigorous validation process is crucial for verifying the accuracy and reliability of the QazGM2025 model for practical geodetic applications within Kazakhstan, especially given the historical nature of some input data [51]. To assess the accuracy of the computed geoid model, comparisons were made against the GNSS/levelling benchmarks described in Section 3.1.2. Since the official vertical reference system of Kazakhstan is realized in terms of normal heights, the comparison required a consistent transformation. Following [52], the geoid heights were derived from GNSS/levelling anomalies using:

$$\mathrm{N}= \mathsf{\zeta }-{\displaystyle \frac{\mathsf{\Delta }{\mathrm{g}}_{\mathrm{B}}}{{\mathsf{\gamma }}_0}}\mathrm{H},$$

(2)

where $\underset{\_}{\gamma }$ in the denominator is replaced by the normal acceleration of gravity for an arbitrary standard latitude (${\gamma }_0$). usually 45. and H is the topographic height of the computation point.

The accuracy assessment of the derived geoid model was carried out using the classical approach, by comparing the modelled geoid heights with their observed counterparts, obtained as the difference between ellipsoidal and orthometric heights at benchmarks of the national levelling network of Kazakhstan. As control points, benchmarks from first- and second-order precise levelling, as well as lower-accuracy third-order points, were employed, totaling 3,653 stations. Residuals were calculated as:

$$\mathsf{\Delta }\mathrm{N}=\left(\mathrm{h}-\mathrm{H}\right)-{\mathrm{N}}^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}},$$

(3)

The scatter of the points selected in this way revealed that the variations in the differences ΔN exhibit a pronounced shift and a slope from northwest to southeast, indicating the accumulation of systematic errors in the data. To approximate this constant trend, we employed a linear model dependent on latitude and longitude coordinates. The model coefficients were determined using linear regression, which, for a spherical surface, takes the form:

$$\mathsf{\Delta }\mathrm{N}(\phi ,\lambda =\mathrm{a} \cos \phi \cos \lambda +\mathrm{b} \cos \phi \sin \lambda +\mathrm{c} \sin \phi +\mathrm{d},$$

(4)

where φ and λ – latitude and longitude, a, b, c and d – model coefficients.

The coefficients of this model were estimated using the least-squares method, providing an empirical correction surface that links the gravimetric geoid to the official height reference system.

4. Results

4.1. Evaluation of Global Geopotential Models

To identify the most suitable long-wavelength gravity field for geoid modeling in Kazakhstan, five recent high-resolution global geopotential models (GGMs) were evaluated against the GNSS/levelling benchmarks described in Section 3.1.2. The models tested included XGM2019e_2159, EIGEN-6C4, GECO, SGG-UGM-2, and EGM2008. For each benchmark, the observed geometric height anomaly derived from GNSS/levelling was compared with the corresponding model-predicted value, and the differences were analysed statistically. The resulting parameters are presented in Table 1.

The comparative analysis demonstrates that XGM2019e_2159 and EIGEN-6C4 yield the lowest RMSE values (0.140 m), with identical standard deviations and mean biases. GECO and SGG-UGM-2 perform only marginally worse (RMSE = 0.146 m and 0.143 m, respectively), while EGM2008 shows the largest dispersion (RMSE = 0.148 m). Although EIGEN-6C4 exhibits a slightly narrower residual range than XGM2019e, the differences are not statistically significant.

Considering its broader spectral content, inclusion of the most recent GRACE and GOCE satellite data, and assimilation of dense terrestrial gravity information, XGM2019e_2159 was selected as the long-wavelength reference model for subsequent LSMSA computations in this study. This selection aligns with practices in other regions where XGM2019e has demonstrated robust performance, even in areas with limited terrestrial data availability [17]. Therefore, the choice of XGM2019e_2159 is expected to provide a highly accurate and consistent long-wavelength gravity field for the QazGM2025 model.

4.2. LSMSA Parameter Optimization

The selection of optimal parameters for the Least-Squares Modification of Stokes' Formula with Additive Corrections is critical for achieving a precise gravimetric geoid model [37]. The LSMSA method requires the specification of two key parameters: the maximum spherical harmonic degree/order (L, M) used to truncate the Stokes’ integral and the error variance of terrestrial gravity anomalies, denoted as C(0). To evaluate their influence on the computed geoid, a series of test solutions were generated by systematically varying these parameters. Truncation degrees from 180 to 760 were tested in combination with gravity anomaly error variances of 1, 3, and 9 mGal². The performance of each configuration was assessed against GNSS/levelling benchmarks using the statistical measures of mean difference, standard deviation, and RMSE.

The results (Table 2) indicate that geoid accuracy is sensitive to the choice of truncation degree. At very low truncation (L = 180), the omission of short-wavelength components results in larger discrepancies, while at higher degrees (L ≥ 500) the amplification of noise leads to increased scatter. A balanced performance was obtained for L = M = 300, which consistently produced the lowest RMSE values across the tested error variances.

Similarly, the assumed error variance of terrestrial gravity anomalies (C(0)) influenced the weighting of input data in the least-squares solution. The best performance was achieved for C(0) = 1 mGal², which reflects the reliability of the digitized Soviet-era gravimetric dataset after validation with modern absolute measurements. Larger values of C(0) (3–9 mGal²) reduced the influence of terrestrial gravity and degraded the fit to GNSS/levelling benchmarks.

Based on the evaluations described above, the final configuration for computing the QazGM2025 model was defined as follows. The terrestrial gravity anomalies, digitized from historical maps and supplemented by WGM2012 data in undersampled regions, were combined with the long-wavelength spherical harmonic coefficients from XGM2019e. Topographic information was provided by the FABDEM DEM, while the LSMSA was implemented with a modification degree of L = M = 300, an assumed terrestrial gravity data error variance of C(0) = 1 mGal², and an integration cap radius of Ψ = 1°. This configuration was used to generate the final gravimetric geoid surface for Kazakhstan.

4.3. Computed Geoid Surface

Using the optimized LSMSA configuration, the QazGM2025 gravimetric geoid model was computed for the territory of Kazakhstan. The resulting surface reflects the combined influence of the long-wavelength field from XGM2019e, the medium- and short-wavelength signal from terrestrial gravity anomalies, and the residual terrain effects derived from the FABDEM digital elevation model.

The applied corrections (Figure 6) — topographic, analytical downward continuation, ellipsoidal, and atmospheric — improved the initial geoid estimates by accounting for terrain mass effects, gravity field continuation, reference surface geometry, and atmospheric mass distribution. The largest adjustments were due to topographic effects (up to −3.86 m), while ellipsoidal and atmospheric terms contributed at the millimeter level. The resulting geoid model computed using the LSMS method (Figure 7) captures the major gravity field patterns across Kazakhstan, with heights ranging from −67.75 m in the southeast to 20.70 m in the northwest, clearly reflecting the influence of complex topography and crustal density variations. The lowest values occur in the western Caspian Depression, where negative geoid undulations are characteristic of thick sedimentary basins and low crustal density. The highest values are observed in the Tien Shan and Altai mountain ranges, where steep relief and crustal thickening dominate the gravity field. Intermediate values are found across the central steppes and plateaus, producing a smooth transition between extremes.

Furthermore, the computed geoid surface of Kazakhstan has been presented at a resolution of 5′ × 5′, highlighting the major spatial patterns (Figure 7). The model captures the expected correlation between geoid undulations and large-scale tectonic and topographic structures, providing a realistic representation of the regional gravity field. The final geoid surface serves as the foundation for subsequent accuracy assessment and validation against GNSS/levelling benchmarks.

4.4. Validation of the Developed Geoid Model

The computed QazGM2025 model was validated against 3,653 GNSS/levelling benchmarks distributed across Kazakhstan. The comparison of observed and modelled geoid heights revealed systematic differences, with a pronounced northwest–southeast tilt in the residuals. This large-scale distortion is consistent with cumulative errors in normal heights introduced by the realization of the Baltic 1977 vertical reference system. To account for this effect, the residuals were modelled using a linear regression function of latitude and longitude. The estimated coefficients were:

a=0.442795±0.070324,
b=1.681212±0.153129,
c=0.504583±0.185160,
d=-1.661341±0.246897,

(5)

Application of the regression correction significantly improved the agreement between QazGM2025 and GNSS/levelling observations. The RMSE decreased from 0.185 m to 0.066 m, while the mean bias was effectively eliminated (Table 3).

The spatial trend in residuals is shown in Figure 8, while the corrected residual distribution is presented in Figure 9. The histograms in Figure 10 illustrate that, after correction, the residuals closely follow a normal distribution centered near zero, confirming the absence of systematic bias.

Importantly, the analysis showed no clear accuracy advantage for first-order benchmarks compared to second- and third-order sites. The dispersion of residuals within each order was similar, indicating that the QazGM2025 model achieves a consistent level of accuracy across the heterogeneous national levelling network. This robust validation process, leveraging a comprehensive dataset of GNSS/levelling benchmarks, underscores the high precision and reliability of the QazGM2025 model in representing the regional geoid for Kazakhstan.

4.5. Comparison with Existing Global Geopotential Models

To evaluate the performance of QazGM2025 relative to existing global solutions, its accuracy was compared with that of five widely used GGMs: XGM2019e, EIGEN-6C4, GECO, SGG-UGM-2, and EGM2008. The validation was performed using the same set of GNSS/levelling benchmarks, ensuring consistency across datasets.

The comparison (Table 4) demonstrates that while the best-performing GGMs (XGM2019e and EIGEN-6C4) achieve RMSE values of approximately 0.14 m, the corrected QazGM2025 model reduces this error to 0.066 m, representing an improvement of more than 50%. The reduction in bias and the near-normal distribution of residuals further confirm the superiority of the regional model over global solutions.

These results confirm that incorporating dense terrestrial gravity anomalies, refined DEM-based terrain corrections, and optimal LSMSA parameterization significantly enhances the accuracy of geoid determination in Kazakhstan compared to global models alone.

5. Discussion

The QazGM2025 model demonstrates a high-accuracy fit to ground truth, with GNSS/levelling comparisons yielding a root-mean-square error of 0.066 m after applying a regression correction. This post-correction performance indicates that most systematic differences (e.g. bias and tilt between the gravimetric geoid and the leveling datum) have been removed, leaving only small random residuals.

The LSMSA approach itself is a proven technique in modern geoid modeling – it facilitates the optimal combination of terrestrial gravity measurements with a geopotential model to mitigate truncation [46,53]. The success of QazGM2025 (reflected in its ~0.07 m agreement with GNSS/levelling) attests to the effectiveness of this method and the improved data inputs. After removing datum biases, the remaining geoid discrepancies are on the order of a few centimeters, which is a strong validation of the model’s internal consistency. Minor residuals can be attributed to localized gravity field complexities (e.g. rugged topography, geologic density variations) and the uneven distribution of input data, but overall, the results confirm that the model achieves the expected precision for a regional gravimetric geoid.

The performance and methodology of QazGM2025 can be benchmarked against global and regional geoid models from other countries with comparable geographical challenges. While global geopotential models such as EGM2008 and EIGEN-6C4 provide advancements in geodetic applications, their accuracy can vary significantly depending on data availability for specific regions, with Central Asia exhibiting errors up to 35 cm [8]. The previous Kazakh geoid model (KazGM2015) showed GPS/levelling discrepancies ranging from –0.295 m to +0.327 m (std. dev. ≈0.17 m) [54]. Similarly, Turkey’s recent gravimetric geoid (Turkey-ITUGRG2021) was computed with a very similar approach (LSMSA with remove-restore using XGM2019 gravity model), and the model achieved an accuracy of about 10.1 cm before any fitting, and ~8.6 cm after a 7-parameter corrector surface was applied [55]. The IRG2018 model of Iran, incorporating a remove-restore technique, though using least-squares collocation and a hybrid solution attained roughly 20 cm RMS agreement with GNSS/leveling after incorporating over 1,200 leveling points [56]. The current QazGM2025 model's enhanced performance underscores the criticality of integrating local gravimetric data with improved global models for regional geoid refinement, surpassing the capabilities of purely global solutions [57].

A reliable high-resolution geoid model such as QazGM2025 provides immediate and concrete benefits to a wide range of stakeholders in Kazakhstan. Geodetic surveyors and civil engineers are among the primary beneficiaries: with this model, they can directly obtain orthometric heights (physical elevations) from GNSS readings, streamlining surveying workflows for infrastructure projects. Urban planning and cadastral authorities benefit by easily tying development projects to the official height system, while hydrologists and disaster managers can integrate LiDAR, river gauge, and coastal data within a unified datum for accurate flood mapping and reservoir management. The model also supports environmental monitoring and geoscience, enabling consistent tracking of Caspian Sea level fluctuations, glacier retreat, snowpack dynamics, and land subsidence. In effect, QazGM2025 ensures that accurate, consistent height information is readily available to a broad range of stakeholders, strengthening infrastructure, resource management, and environmental resilience.

Despite these advances, several limitations remain. The geoid’s absolute accuracy is constrained by data quality, particularly the absence of absolute gravity control points and sparse coverage in mountainous and border regions, where errors may exceed the national RMSE. Future campaigns, such as airborne gravimetry and new GNSS benchmarks, are needed to densify observations. The model is also static, while long-term processes (tectonic uplift, hydrological mass change, post-glacial rebound) can cause centimeter-level geoid variations over decades. Periodic updates, potentially incorporating GRACE/GRACE-FO time-variable gravity data, will be essential to maintain precision. Finally, global integration is required: QazGM2025 should be tied to the International Height Reference Frame (IHRF) to ensure cross-border consistency. Methodologically, refinements to LSMSA assumptions and tests with alternative approaches (e.g. least-squares collocation, multi-resolution models) could further improve local performance.

In conclusion, while QazGM2025 is a significant leap forward, and provides a solid foundation in regional geoid building for Kazakhstan. However, ongoing efforts in data collection, methodological refinement, and global integration will be essential to further enhance the geoid model. Such future work will ensure that Kazakhstan’s vertical reference frame remains accurate, robust, and interoperable, supporting the country’s needs well into the future.

6. Conclusions

This study presents the first dedicated solution to the problem of geoid determination for Kazakhstan, addressing long-standing challenges in regional gravity field modelling and vertical datum unification. A comprehensive geodetic information base was established, combining digitized Soviet-era gravimetric maps, modern absolute gravity control, and GNSS/levelling benchmarks. The validation of the digitized gravity dataset confirmed its suitability for regional geoid modeling, while a comparative evaluation of global geopotential models identified XGM2019e_2159 as the optimal long-wavelength reference field. For the topographic component, the FABDEM DEM was verified as sufficiently accurate for geoid computation.

The geoid was computed using the Least-Squares Modification of Stokes’ formula (LSMSA) with optimized parameters (L = M = 300; C(0) = 1 mGal²), integrating terrestrial gravity anomalies, global model coefficients, and DEM-based corrections. Validation against 3,653 GNSS/levelling benchmarks demonstrated that systematic distortions in the national height system could be corrected through regression adjustment, reducing the RMSE from 0.185 m to 0.066 m. The resulting QazGM2025 geoid model outperforms global models, achieves stable centimeter-level accuracy across diverse terrain, and provides a reliable and consistent vertical reference surface.

These results establish QazGM2025 as a robust foundation for modernizing Kazakhstan’s vertical datum, enabling GNSS-based height determination, unification of the national height system, and future alignment with the International Height Reference Frame (IHRF). Ongoing improvements in data coverage, absolute gravity control, and dynamic model updates will further enhance its accuracy and ensure its utility for surveying, infrastructure development, and environmental monitoring.