Traditional Very Long Baseline Interferometry (VLBI) faces limitations in providing comprehensive geodetic products due to its reliance solely on ground-based observations of extragalactic radio sources. This study investigates the feasibility and benefits of a novel three-level integrated VLBI observation mode that combines quasar, satellite, and space VLBI observations to enhance geodetic parameter estimation. We conducted Monte Carlo simulation using GALILEO, BDS, GPS satellites and ETALON-1/-2 satellites as the Medium Earth Orbit space segment platform, integrating space and satellites VLBI observations with existing continuous geodetic VLBI campaign schedules (CONT17). Our analysis focused on assessing the impact of space VLBI observation quality and quantity on derived geodetic parameters, particularly station positions and Earth Rotation Parameters (ERPs). Results demonstrate that space VLBI observation could perform comparably to satellite VLBI for geodetic references due to providing superior geometric configuration. The integration of space VLBI into combined quasar-satellite schedules, yielded a 10% improvement in station position quality. Assuming 300-picosecond precision for space VLBI observations and GALILEO satellite scans replacing every fourth quasar scan, the derived ERP precision achieved 23.2 las (xp), 20.4 las (yp), 1.5 ls (dUT1). Space VLBI observations may offer the potential to improve the positioning of satellite orbital planes within the celestial reference frame by providing additional spatiotemporal information and broader sky coverage. These observations can also help compensate for the limited arc length of satellite VLBI tracking, particularly when the visible segment of the satellite orbit is short, thereby contributing to amore robust reconstruction of the orbital geometry. This integrated approach demonstrates substantial potential for advancing space geodetic reference frame precision and establishing more consistent space ties between geodetic techniques. (c) 2026 COSPAR. Published by Elsevier B.V. All rights are reserved, including those for text and data mining, AI training, and similar technologies.

Abstract. River discharge provides one of the clearest measures of how Arctic freshwater systems respond to climate change, yet long-term and continuous monitoring remains challenging. Here we present a gauge-based monthly natural discharge reconstruction dataset for 94 gauging stations in the Eurasian Arctic from 1952 to 2025. The dataset was generated using a basin-scale framework driven by precipitation and air temperature and aided by snow information. Terrestrial water storage (TWS) serves as the key link in the framework. Historical TWS was first reconstructed using an improved state-update empirical model constrained by TWS observations from the Gravity Recovery and Climate Experiment (GRACE), and monthly natural discharge was then estimated from reconstructed TWS using a refined runoff–TWS relationship calibrated against gauge discharge observations. The reconstructed monthly natural discharge agreed well with observations at minimally regulated gauges and during pre-regulation periods at regulated gauges, with median monthly Kling–Gupta efficiency (KGE) values exceeding 0.88. When monthly values were aggregated to annual discharge, the reconstruction still captured observed variability at minimally regulated gauges, with a median annual KGE of 0.73. The dataset generally agreed better with gauge observations than discharge estimates derived from existing GRACE-like TWS products and an existing gridded runoff reconstruction product. Example applications illustrate spatially heterogeneous trends in annual discharge and seasonal allocation over the past seven decades and show how the dataset can be used to assess the effects of human regulation on river discharge. This dataset provides a long-term natural-discharge baseline for characterizing historical discharge variability and assessing the roles of climate variability and human regulation in poorly monitored Arctic basins. The reconstructed monthly natural discharge dataset is available at https://doi.org/10.5281/zenodo.21157816 (Liu et al., 2026).

Potential theory provides rigorous solutions to geodetic boundary value problems involving Laplace's equation, often expressed through global integral transformations. However, the practical application of these solutions is hindered by the incomplete availability of highresolution gravity data. This limitation necessitates the separate evaluation of far-zone contributions-commonly referred to as truncation errors. These errors are typically characterized using truncation coefficients, whose computation via recurrence formulas becomes numerically unstable at high altitudes and degrees when using double precision. This study presents a robust and comprehensive method for accurately computing truncation coefficients for the Poisson, Hotine, and Stokes integrals across all degrees and altitudes. The proposed approach integrates three specialized strategies: (1) an asymptotic expansion for high-degree terms to ensure both accuracy and efficiency; (2) a two-regime altitude model for low-degree terms, with binomial expansions applied at low altitudes and stabilized recurrence series used at high altitudes; and (3) double-precision numerical experiments that validate the method, showing relative errors below 10-5 for the Poisson integral and below 10- 8 for the Hotine and Stokes integrals. These results demonstrate that the method is both accurate and reliable for geodetic applications at all altitudes-from airborne surveys to satellite missions-and across the full range of spherical harmonic degrees.

The Volga River Basin is one of the most important water systems on the Eurasian continent. Variations in its terrestrial water storage (TWS) not only exert a profound influence on regional hydrological processes but also dominate the fluctuation of the Caspian Sea level (Correlation coefficient: 0.90). However, due to the limited temporal coverage of Gravity Recovery and Climate Experiment (GRACE) observations, understanding of the long-term evolution of TWS in this region remains insufficient. This study introduces a multi-source data fusion approach that does not rely on any a priori information from GRACE. By integrating ERA5-Land, FLDAS, and GLDAS products with the terrestrial water balance equation and the Extended Triple Collocation (ETC) algorithm, we reconstructed a long-term TWS time series for this region spanning more than four decades (1982-2023). Benefiting from the higher spatial resolution (0.1 degrees & times; 0.1 degrees) of fusion model, a pronounced declines trend (approximately-15 mm/yr) in TWS within certain rivers and lakes was successfully identified. The decomposition results of terrestrial water storage anomalies reveal a progressively increasing amplitude of its seasonal component, while the long-term trend is primarily modulated by climatic wetness and dryness conditions. Snow water equivalent (SWE) and Soil moisture (SM) are the dominant contributors to TWS change, accounting for 48.3% and 36.1%, respectively. The contribution of SWE is largely concentrated in high-latitude regions with an increasing trend. This study makes up for the limitation of the insufficient temporal span in GRACE, providing an objective basis for understanding long-term TWS variations under climate change.

The Newtonian volume integral is the fundamental mathematical framework for modeling the gravitational effects of mass bodies with arbitrary geometry. Tesseroids provide a natural discretization of spherical shells and are therefore well suited for global-scale applications. However, the curved geometry of tesseroids makes it difficult to derive analytical solutions for the Newtonian integral. Numerical approaches, such as Taylor expansions and Gauss–Legendre quadrature, suffer from two major limitations: near-zone errors and polar-region effects, the latter being closely related to the choice of integration elements. In this study, we derive linearly analytical solutions (LASs) for the gravitational effect of tesseroids by applying linear approximations to both the integrands and variables, based on multiple integration elements with and without the cosφ^' term. The main results are as follows: (1) With respect to the integration elements, the cosφ^' term preserves the meridian convergence effect, and both numerical and analytical methods constructed with integration elements containing this term are inherently free from polar-region errors. (2) For Newtonian integration with arbitrary-order polynomial density, numerical integration does not introduce additional errors; in particular, the error does not increase with the polynomial degree. (3) We characterize, for the first time, the radial extent of a tesseroid—defined by its latitudinal and longitudinal coordinates—as a composite function, and derive LAS for the Newtonian integral over a non-rectangular domain. This formulation enables accurate modeling of tesseroids with inclined top or bottom surfaces. Building upon previous developments of LAS in spherical polar coordinates, this study establishes a solid mathematical foundation for modeling tesseroid gravitational effects using linearly analytical expressions. It also advances the understanding of Newtonian integration under different integration elements and represents the first attempt to formulate quadrature rules for numerical integration when the integration elements include the cosφ^' term.

The combination of satellite gravimetry measurements with other techniques can promote the integration of the "Three Pillars" of geodesy, namely, the Earth's shape, gravity field, and rotation. Combined processing of Gravity Recovery and Climate Experiment (GRACE) gravimetry and ground-based Global Positioning System (GPS) measurements at the observation level can theoretically improve parameter accuracy. However, because of the implementation challenges, existing combination experiments are limited to low-degree (up to degree and order 20) gravity field determination or daily and weekly solutions. Here we present a one-step method for determining higher monthly gravity field solutions and the resulting models are up to degree and order 60. A comprehensive analysis of theoretical models and experimental results demonstrates the advantages of the one-step approach over the conventional two-step method, which processes ground and GRACE observations separately. Compared to two-step results, the one-step approach yields improvements in low-degree gravity field coefficients, with a 60 % improvement in the signal of the C20 coefficient. The polar motion x-component shows an average improvement of 58 %, while GPS satellite orbits achieved approximately 35 % reduction of Satellite Laser Ranging (SLR) residuals RMS. The GRACE satellite orbit exhibits a 22 % reduction in SLR residuals RMS for GRACE-A and 17 % for GRACE-B. In summary, the one-step method produces more consistent gravity and geometric products, improving the accuracy of shared parameters between ground-based and GRACE observations. (c) 2025 The Authors. Published by Elsevier B.V. on behalf of COSPAR. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/).

Gravity satellites play an essential role in the recovery of the Earth’s gravitational field. The satellites are equipped with kinds of instruments, each installed at predetermined locations within the satellite’s framework, designed to meet various mission requirements. It is necessary to correct the observations from different instruments location to the satellite’s center of mass for the precise orbit determination. This implies that the accurate determination of various baseline lengths within the satellite is crucial, such as the correction from the antenna phase center and the accelerometer’s center to the satellite’s center of mass.This study focuses on the GOCE satellite and discusses the satellite gravimetry measurement mode when there is an offset between the accelerometer and the satellite’s center of mass. In this study, the satellite’s orbit, accelerometer calibration coefficients, and gravitational field potential coefficients are simultaneously estimated using the dynamic orbit determination methods. Preliminary results indicate the influence of the GOCE satellite’s HL-SST observations on the gravitational field, and accurate center offset parameters play a positive role in the recovery of gravitational field model. Combining the HL-SST observations and LL-SST KBR observations from the GRACE satellites, the accuracy of the upgraded gravity field model is consistent with the state-of-the-art international gravity field model.

The gravitational effect of tesseroid, expressed as a Newtonian volume integral, is generally not analytically solvable across most of the globe, in spherical coordinates. A closed-form solution exists only when the computation point is on the polar axis. In this case, spherical and spherical polar coordinates coincide, meaning that a tesseroid regularly defined in spherical coordinates remains regular in spherical polar coordinates. Thus, the analytical solvability of the tesseroid’s gravitational effect depends on whether its domain remains regular in spherical polar coordinates. This study decomposes a tesseroid with an irregular domain in spherical polar coordinates into a set of regular triangular prism elements and derives an analytical solution for the gravitational effect of a triangular prism. This allows us to obtain an analytical solution for a tesseroid at any computation point. To validate this decomposition, we examine the linear relationships along each tesseroid edge and demonstrate that it can be well approximated by four regular triangular prism elements. To assess accuracy, the relative errors of the analytical solution are compared globally against two reference methods: the 3D Gauss–Legendre quadrature (3D-GLQ) and the closed-form solution at the polar axis. Results show that, except near the computation point and its symmetric point located on the opposite side of the earth, our analytical solution surpasses the 3D-GLQ in accuracy and achieves 26.2

The quality of low Earth orbit (LEO) satellite-borne global positioning system (GPS) observation data has an important influence on the precise orbit determination of LEO satellite and the study of the ionosphere. In this contribution, we develop a geometry-free orbit-dynamic-constraint approach to analyze the satellite-borne GPS observations errors, including code errors, carrier phase errors, receiver clock errors, ionospheric delay, to detect and repair cycle slips for satellite-borne dual-frequency carrier phase observations. Different from the traditional methods suitable to satellite-borne GPS observations, we take full advantage of the correlations between the dual-frequency observations, between the satellite orbit variations and the orbit dynamic background models, without a prior precise orbit. The proposed approach can not only remove the geometry variations between GPS satellites and LEO satellite, but also have partial characteristics of other signals separated from GPS observations, such as the receiver clock offset. Then, the performance of the proposed approach is verified by GPS measurement data of GOCE/GRACE satellite with different sampling interval. The results indicate that the proposed approach preforms well on observations error analysis, can effectively detect and fix cycle slips with high success rates, even in real-time applications.

Cycle-slip detection is crucial in achieving high-accuracy Global Navigation Satellite System (GNSS) data processing for Low Earth Orbit (LEO) satellites. The detector based on the dynamic force model has emerged as a promising approach to cycle-slip detection, as it is insensitive to the number of visible tracked satellites and insufficient accuracy of a prior orbit. However, the ionospheric delay limits the applications of the method for high-speed LEO satellites under different scenarios. To complete the method, we propose a new com-bination of two test parameters: the phase ionosphere-free (IF) combinations in the second-order time-difference model and the Mel-bourne-Wu & BULL;bbena (MW) combinations. Meanwhile, we highlight the identification of whether the fake-positive cycle slips are derived from the flicker noises and jumps in the receiver quartz-based clock. Considering highly, moderately active and quiet ionospheric activ-ities, comparisons are performed among our proposed two test parameters and the phase second-order time-difference wide-lane (WL) combinations. The results indicate that the detector based on the IF combinations performs optimally under the three levels of iono-spheric activities. The time-difference ionospheric delays vary dramatically even under moderately active and quiet ionospheric activities, which leads to the inferiority of the detector based on the WL combinations to one based on the MW combinations. The performances of detectors based on the IF and MW combinations are further evaluated through real-time kinematic orbital determinations comparisons with the Jet Propulsion Laboratory (JPL) precise science orbits (PSO). After using the IF detector, the mean 3-Dimensional (3D) root-mean-squares (RMS) of the orbit differences have slightly decreased by 12% and 6% for the Gravity Recovery and Climate Experiment (GRACE) A and B satellites, respectively. In the case of the GRACE follow-on (GRACE-FO) satellites, the IF detector has achieved a more than 40% improvement over the MW detector. Further, the use of both detectors simultaneously results in a slight improvement in orbit accuracy.& COPY; 2023 COSPAR. Published by Elsevier B.V. All rights reserved.

Ambiguity resolution (AR) is critical for enhancing the orbit accuracy in precise orbit determination (POD) for low earth orbit (LEO) satellites. While orbit estimation using single-difference (SD) AR has been widely researched, the investigation of orbit estimation using undifferenced (UD) AR for LEO satellites is limited due to time-varying hardware biases at the LEO receiver end. To address this deficiency, we propose an improved UD AR method for LEO satellite orbit determination. The method employs the optimal integer datum ambiguity for a 1-day observation arc, and the random-walk clock model is utilized to transfer the integer ambiguity datum to the other epochs, which enables the arc-wise ambiguities to regain the integer property within the time frame. To validate the effectiveness of our improved method, numerical experiments are conducted. Moreover, we assess the performance of the random-walk clock model for a spaceborne ultra-stable oscillator (USO). The high frequency stability of the USO establishes the satisfactory requirements for this study. Both in the kinematic and dynamic modes, the AR success rates of our improved method are - 2% higher than that of the current SD AR method, and the orbit results indicate a slight improvement. The 3-Dimensional (3D) root-mean-squares (RMS) of the orbit differences between JPL precise science orbits (PSO) and our orbits are reduced by up to 4% and 5% for kinematic and dynamic orbits, respectively. The subtle benefit is also proven by K-band ranging (KBR) system validations. In the case of the Gravity Recovery and Climate Experiment Follow-on (GFO) satellites, the results indicate that the effect of our improved UD AR can be equivalent to that of SD AR. Our improved UD AR method can be a good alternative for its superior ability to generate hardware delay at the LEO receiver end. & COPY; 2023 COSPAR. Published by Elsevier B.V. All rights reserved.

Precise calibration of the GOCE gravity gradiometer is one of the premises to determine a state-of-the-art global gravity field model. Provided the center of the gradiometer coincides to the GOCE center of mass, the common-mode accelerometer observations can generally represent the satellite’s non-gravitational accelerations. However, during the satellite flight, it is necessary to consider the possible offset of the gradiometer center with the center of the satellite towards the end of the mission due to orbital maneuvers and other reasons. In this research, the gradiometer calibration model is established by a priori gravity field model with six accelerometer observations and three star trackers angular rates. Besides the inverse calibration matrices, it includes the mass center offset as new calibration parameters and the impact of the parameter drifts is discussed. The feasibility of this method to solve this parameter is verified using the nominal mode measured L1b data, and the impact of this offset on GOCE gravity gradient data quality are reduced to a large extent. In addition, the analysis of time series of the calibration parameters and the estimated PSD of the gravity gradient trace demonstrate that the accuracy of the calibrated gravity gradients is improved.

During the processing of GRACE observation data, the GNSS satellite ground station observations are added, which can directly transfer the datum from the ground stations to achieve a more theoretically rigorous parameter estimate. The software platform independently developed by our research group is used to realize the integrated solution at the observation level of GRACE high-low and low-low tracking data and ground station GNSS observations. AOD1B product is a prior correction model for ocean and atmosphere de-aliasing used in GRACE data processing. AOD1B product is an atmospheric and ocean mass changes high-frequency signal and belongs to non-tidal time-varying gravity field signals. The effect of the application of the AOD1B product in the orbital dynamics model of GNSS satellites on the solution of the satellite precision orbit determination and monthly earth gravity field time-varying models are considered in the integrated solution of multi-source data combination. The differences between the estimated satellite orbit parameters and time-varying gravity field models' time series with and without the AOD1B product in the GNSS satellite orbit dynamics model are displayed, along with the time series comparison of estimated tropospheric parameters for the several GNSS tracking stations.

Abstract The cycle-slip detection is the prerequisite for achieving the utmost precise carrier-phase measurements for the low earth orbit (LEO) spaceborne Global Navigation Satellite System (GNSS) receivers. However, the issue has not yet been solved entirely due to the limitations of traditional cycle-slip detection methods. In the paper, we proposed a novel cycle-slip detection method based on the second-order time-difference phase ionosphere-free (IF) or wide-lane (WL) observations with dynamic orbital constraints. Nominal orbits at the dm-level accuracy and the earth gravity field of at least 30 degrees and orders (d/o) within 30 s or higher sampling rate data could guarantee the proposed method reliable for the flying altitude of near 500 km of LEO satellites. To validate the proposed method effectiveness, numeric experiments with comparisons to the traditional cycle-slip detector based on the Melbourne–Wübbena (MW) combination are carried out with well-chosen measurements for the LEO satellite under ionospheric activities of different levels. The results, in terms of the correctness and misjudgment of the cycle slips detected and the precision of estimated orbit, indicate that, on the one hand, the detector based on the proposed phase IF test parameters are superior to one based on the proposed phase WL test parameters or the MW combinations under any level of ionospheric activities; on the other hand, under highly active ionospheric activities, the proposed phase WL cycle-slip detector discovers a fewer number of real cycle slips, which are flagged by the reference method introduced in the paper, but more misjudgments than the MW cycle-slip detector does, which leads to the inferior of the WL cycle-slip detector. However, the phase WL cycle-slip detector performs slightly better than the MW cycle-slip detector under moderately active and very quiet ionospheric activities. In view of the developing multi-frequency multi-GNSS uncombined observation data scanning, the proposed phase IF test parameter accompanied with the MW combinations is recommended to be applied to the cycle-slip detection if the predictable ionospheric refraction lies at the strongly intensive level of the ionospheric activity. Otherwise, the proposed phase WL test parameter would be involved.

The accurate calibration of the GOCE gravity gradiometer is one of the premises for determining the earth's gravity field model with high precision. L1b dataset of GOCE gradiometer and star sensors are used in the internal calibration. The combination of different star sensors by least-squares adjustment can prevent the propagation of the less accurate component due to the reference frame transformation, thus the accuracy of angular rates used by internal calibration can be improved. In this paper, GOCE data in November 2009 are used to verify the ESA's calibration method. On this foundation, considering the rotation matrices between star sensors and gradiometer are time-varying, improved calibration model is presented by using parameters of the rotation matrices. The analysis shows that the parameters are about 100 arcseconds with a liner drift of 3~30 arcseconds in this month. Based on ESA's internal calibration model considering parameters of three accelerometer pairs, a calibration by using transformation matrix from star sensors to gradiometer and the parameters of three accelerometer pairs is presented in this paper. The accuracy of gravity gradients after calibration shows the effectiveness of this method below frequency of 0.005 Hz. Possible developments of GOCE gradiometer calibration based on this method are discussed in this paper which provides foundation for processing of GOCE and other gravity satellites.

High mountain glaciers (HMGs), called the water towers of the world, are vulnerable to the effects of climate change and thus are rapidly shrinking. Monitoring and evaluating large-scale glacier and snow (GS) mass changes are critical for humans and ecosystems. Although modern gravity satellites monitor GS on a global scale, the contemporary Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow−On (GRACE−FO) missions are only able to do so at coarse spatial resolutions and are restricted by their limited observation periods. Moreover, because of the complex mechanism behind a glacier's mass balance, existing reconstruction methods for GRACE-like mass anomalies, including machine learning approaches and statistical models, cannot be applied to HMGs. Here, we propose a precipitation and temperature data-driven statistical model combining hydrology and GS processes to reconstruct GRACE-like mass anomalies of HMGs, from which the hydrological and GS signals can be further separated. The prediction and reconstruction performance of this method was consistent with the GRACE/GRACE−FO observations (the median correlation coefficient/Nash-Sutcliffe efficiency: ~0.92/0.85). Additionally, the separated GS signals agreed with the independent external data used for comparison. Compared with machine learning methods, this method better reconstructed the long-term trend of GRACE-like mass anomalies. Our study provides an acceptable method for expanding mass anomaly time series on HMGs, thereby assisting in the sustainable management and protection of water resources in their downstream areas.

基于DE405、DE421、DE430、DE440星历,计算各大行星在地球质心及太阳系质心惯性系中的位置,比较其他星历相对于DE440星历的行星位置精度.分析讨论各历表下月球的地心位置和速度精度,以及历表对于月球探测器的位置和速度从月心惯性系转为月固系的影响,并给出使用建议.结果表明,各大行星由于受观测数据等因素影响,位置精度差异较大,跨度从m级到106 m量级.对于行星的位置精度,DE421和DE430相对于DE405有1～2个量级的提高,DE430相对于DE421提高50％.DE405历表的月球地心位置、速度精度分别为7 m和0.02 mm/s,DE421历表分别为1.5 m和0.004 mm/s,DE430历表分别为1.3 m和0.0035 mm/s.对于星历用于月球探测器从月惯系转为月固系产生的坐标和速度误差,DE405分别为30 m和3 cm/s,DE421分别为1.3 m和1.2 mm/s,DE430分别为1 m和0.9 mm/s.历表对于月球坐标系转换的影响为m级.对于月球探测器导航及相关任务,推荐使用DE430或DE440行星历表.

In this study, we composed the velocity and acceleration of GRACE Follow-On satellites directly estimated from the onboard carrier phase observations based on the FIR (finite impulse response) differentiator instead of the orbit differentiation method. Using the level 1B data of GRACE Follow-On from DOY 305 to 314 in 2018, CODE(Center for Orbit Determination in Europe)precise ephemeris and 5 s clock corrections, the results showed that based on the ninth-point differentiator, the 3D root mean square (RMS) of satellite C and D velocity can achieve the accuracy at 0.227 6 mm/s and 0.238 4 mm/s (differentiator interval=60 s); Meanwhile, the 3D RMS values of acceleration can achieve the accuracy at 4.1 μm/s2 and 4.5 μm/s2 for satellite C and D respectively (differentiator interval=90 s). The GNSS phase direct differentiation method does not require a fixed ambiguity, which weakens the influence of the correlation between orbit epochs, respect to the kinematic orbit difference method and provides high-precision velocity and acceleration information for gravity field model determination.

Based on a satellite constellation composed of two GRACE-type satellite formations with different inclinations (near polar orbit + low inclination) and the theory of repeat orbit cycle, we discuss the methods for selecting medium-low inclinations for global and local gravity fields. The effects of this constellation configuration on gravity field inversion are comparatively analyzed using a whole-course dynamics simulation. The results show that compared with the single GRACE-type satellite formation, the use of satellite constellations with different inclination configurations improves the gravity solution precision by 34%. The inclusion of multi-directional observations can improve the spatio-temporal resolution of the satellite missions, and yield gravity field solutions with higher isotropic sensitivity. Furthermore, it is necessary to select the optimal low inclination according to the study area, which will have a significant influence on the gravity field solution.