The renewed interest in lunar exploration and the development of future lunar communication and navigation services highlight the need for a precise, stable, and interoperable geodetic and timing infrastructure on the Moon. NovaMoon, proposed as a scientific and navigation payload for ESA's Argonaut lander, is designed as a lunar-based local differential, geodetic, and timing station supporting both operational needs in the Moon's south polar region and a broad range of scientific investigations. The payload integrates a lunar laser retroreflector, a Very Long Baseline Interferometry transmitter, a receiver for navigation signals compatible with LunaNet standards, high-stability atomic clocks, and direct-to-Earth radio links – making it the first lunar station to co-locate multiple ranging, tracking, and timing techniques. NovaMoon will enable sub-metre to decimetre positioning, provide local differential corrections for lunar users, and ensure an accurate and stable realisation of position and time. Preliminary simulation studies show that this multi-technique dataset improves the lunar reference frame, orientation and ephemerides, and estimates of interior parameters like tidal response and core properties. NovaMoon will also provide the first long-duration physical realisation of a lunar time reference. Beyond its primary goals, it supports improved cartography, precise surface geolocation, and higher-resolution topography, contributing to safer landings and operations. It also enables new tests of fundamental physics, including constraints on relativity and possible deviations from classical gravity.

The zero-difference (ZD) and double difference (DD) approach for processing GNSS observations is mathematically equivalent – meaning that in DD case the huge number of clock parameters are just pre-eliminated. The results for all remaining parameters are identical, if no numerical shortcut is done (e.g., ignoring parts of the correlations introduced by the DD approach). The advantage of processing DD observations is the reduced number of parameters, easier detection and potential correction of cycle slips, and direct access to the integer ambiguities without any phase bias parameter. On the other hand, ZD processing offers greater flexibility in network configuration and parameter handling. This is particularly advantageous when modifying the list of stations in the processing, including Low Earth Orbiting satellites (LEOs), or implementing advanced clock models, e.g., for Galileo satellites.Based on an experimentally developed ZD-based ambiguity resolution method that introduces ambiguity clusters and satellite-wise consistency corrections, the original research prototype was translated into a robust and automated routine processing chain which is suitable for operational use.The new procedure follows the structure of the established CODE DD strategy but adapts the individual processing steps, e.g. pre-processing, estimation of global parameters, handling of receiver-dependent parameters, and ambiguity resolution. Special emphasis is placed on numerical stability, the reliable handling of real-valued ambiguities, and the introduction of quality-control mechanisms designed for long-term autonomous operation. The resulting procedure enables efficient parallelization and delivers consistent orbit, clock, and ambiguity products. We have investigated the requirements on the station density for all these steps in order to optimize also the processing time. Stations, that are not needed in this context can get pre-processed based on the PPP approach and get added to the final solution only.Initial results show that the operational ZD processing chain reaches the accuracy and stability of CODE DD-based products while offering greater flexibility for future extensions. The results demonstrate that the ZD-based GNSS processing is sufficiently mature to generate stable global products on a daily basis and therefore represents a promising foundation for next-generation GNSS solutions computed at CODE.

In order to meet the demanding accuracy and availability requirements, GPS has introduced L5 signals that are compatible with Galileo E5a signals. These signals are designed to mitigate multipath issues and poor performance in challenging environments such forests and areas affected by jamming. As they are also intended to replace L2 signals in the future, steps must be taken to exploit the modern signal type that is currently broadcast by 20 GPS satellites of blocks IIF and III. This is considered particularly important for some LEO satellites, which rely exclusively on L1/L5 observations.In addition to standard analysis products based on L1/L2, the CODE (Center for Orbit Determination in Europe) IGS Analysis Center is in the process of testing a prototype processing chain to generate L1/L5-based products, paving the way for GNSS processing to be fully based on L1/L5 signals. The presentation addresses the application of these products to LEO orbit determination and PPP processing of the ground stations, considering different antenna calibrations for IIF satellites. A quantitative validation and comparison with L1/L2-based solutions is also discussed.

Abstract At the Astronomical Institute of the University of Bern (AIUB) the Celestial Mechanics Approach for the analysis of satellite orbital motion is currently employed for two space geodetic techniques, namely GNSS and Satellite Laser Ranging (SLR). Both are not only used for station network analyses but also for precise orbit determination (GNSS) and validation (SLR) of Low Earth Orbiting (LEO) satellites. In addition, GNSS-based LEO orbits and SLR observations to the spherical geodetic satellites are used in Earth gravity field recovery. Taking GRACE Follow-On as an example, we present temporal gravity field results, where we derive dynamic and reduced-dynamic orbits together with gravity field parameters from GPS carrier phase observations and combine them with SLR observations to the satellites. Hence, using the satellites as a space tie with respect to the gravity field. We validate our results by comparing them with solutions of superior quality derived from the operational GRACE Follow-On data processing at the AIUB, which include the K-band range-rate data, and we investigate the quality of the co-estimated orbits.

As low-Earth orbit (LEO) satellites are increasingly equipped with highly precise Global Navigation Satellite Systems (GNSS) receivers, new opportunities for GNSS-based precise orbit determination (POD) arise. The availability of GNSS observations from satellites that are geometrically well distributed over the entire Earth may enable the formation of an independent GNSS network in LEO, which can be processed solely on the basis of double-difference (DD) observations. This study presents a first attempt to compute a network solution of GNSS DD observations in space that combines data from eight satellites, including GRACE-FO C/D, Sentinel-3A/B, Sentinel-6A and Swarm A/B/C. Using satellite laser ranging (SLR) validation, we show that the absolute position accuracy of an ambiguity-float network solution increases from about 4 cm 1D RMS in a constellation of three satellites to about 2 cm–2.5 cm in a constellation of eight satellites. Similar results are obtained by comparing the baseline lengths from our DD network solution with those derived from external ambiguity-fixed zero-difference (ZD) orbits, with differences at a level of 2 cm to 3 cm RMS for a network of eight satellites. For both the SLR residuals and the baseline consistency, an improvement of up to 0.5 cm can be achieved when resolving 60

The European Space Agency’s BepiColombo mission continues its pioneering voyage to Mercury, the innermost planet of the Solar System. Among the advanced instruments onboard the Mercury Planetary Orbiter (MPO), one of the two spacecraft that comprise the BepiColombo mission, is the Italian Spring Accelerometer. The instrument’s primary scientific goals are closely linked to the Mercury Orbiter Radio-Science Experiment. Together, these instruments aim to provide valuable data on the spacecraft’s orbit, as well as Mercury’s gravity field and internal structure. This simulation study examines how accelerometer measurements affect orbit and gravity field recovery, and it explores strategies to overcome the challenges associated with the co-estimation of all parameters. In it, we evaluate the accuracy of the integrated retrieval of spacecraft orbit, gravity field, and accelerometer parameters under different noise levels and varying observation geometries during the mission. To achieve these goals, we propagate the orbit of MPO and simulate Doppler observations and accelerometer measurements based on the available noise models. We consider two scenarios by including either a realistic or a more conservative accelerometer noise level. Our results indicate that postponing the estimation of accelerometer biases until a preliminary gravity field is established helps prevent gravity field mismodelings from being absorbed into the accelerometer parameters. The daily estimation of bias was found to be essential. Using one year of Doppler tracking data and applying Kaula regularization, we carry out orbit determination and gravity field recovery under varying assumptions regarding the accelerometer noise and the observation geometry. Improved recovery of low-degree coefficients was observed under realistic noise assumptions. Moreover, the strong influence of observation geometry in shaping the error behavior was observed. namely, errors in cross-track bias estimation increased sharply when the β-Earth angle dropped below 45°, corresponding to the degraded Doppler observability periods. Orbit determination achieved decimeters-level accuracy in the radial direction and metre-level accuracy in the along- and cross-track directions, with a substantial rise in cross-track orbit errors observed during periods of reduced Doppler observability. This study also demonstrates the applicability of the planetary extension of the Bernese GNSS Software as an independent tool for orbit determination and gravity recovery for the BepiColombo mission.

In recent years, the use of non-scientific low Earth orbiting (LEO) satellites for gravity field determination has been increasingly explored. In principle, the same methods as for non-dedicated scientific missions are applied by analysing the satellites’ orbital perturbations to gain information about the Earth’s gravity field. With currently over 100 CubeSats in orbit, the Spire commercial constellation provides a huge amount of GNSS tracking data in the same timeframe as scientific missions. This can be exploited to increase the spatio-temporal resolution of estimated gravity field solutions. Although the data quality is limited, previous analyses using data from a 2020 ESA Announcement of Opportunity project have shown that a combined processing of Spire CubeSats can achieve monthly gravity field solutions of similar quality as non-dedicated scientific missions. In our work, we make use of both the Spire data from 2020 and a new dataset from 2023 provided by EUMETSAT. At the Astronomical Institute of the University of Bern, non-gravitational forces acting on the satellite have so far not been explicitly modelled when determining the gravity field, but instead absorbed by pseudo-stochastic parameters such as piecewise constant accelerations. As these forces are prominent for CubeSats due to their low altitude and their large area-to-mass ratio, this study explores how the estimation of the orbits and gravity field solutions improves when explicitly modelling these forces. To achieve this, we determine kinematic orbit positions from Spire CubeSat GNSS phase and code data and introduce them as pseudo-observations in an orbit and gravity field recovery step, where dedicated non-gravitational force modelling is applied using macro-model information provided by Spire. The combination of the different Spire satellites is performed at the normal equation level using a least-squares approach.

MESSENGER’s Mercury Laser Altimeter (MLA) provided high-resolution topographic models and geodetic insights into Mercury’s internal structure, though its data was limited to the northern hemisphere. Due to its less elliptical orbit, the BepiColombo Laser Altimeter (BELA) will soon be delivering global altimetry coverage. We present expected improvements on our knowledge of Mercury’s geodetic parameters, such as its obliquity, short- and long-period librations, rotation rate, and the amplitude of tidal deformations, by comparing crossover analyses of the MLA and BELA altimetry datasets within consistent closed-loop simulations. We discuss the potential of cross-dataset crossovers (consisting of a BELA track crossing with an MLA track) to refine the determination of key geodetic parameters by providing extended coverage at lower latitudes. Our simulations incorporate realistic altimetric ranges, orbital dynamics, and perturbations due to the imperfect knowledge of, e.g., MESSENGER orbit and attitude, and of Mercury’s orientation. We show the improvement of crossover-derived uncertainties of the geodetic parameters, calibrated using variance component estimation. In particular, rotation parameters, including long-period forced librations, significantly benefit from combining inter-instrument crossovers. These uncertainties are finally mapped through Markov Chain Monte Carlo to study BELA’s impact on constraining plausible models for Mercury’s internal structure.

Coronal mass ejections (CMEs) from the Sun can cause geomagnetic storms which cause the thermosphere to expand. This leads to enhanced air drag for satellites in low Earth orbit (LEO). This work focuses on the evaluation of orbital decay with a focus on selected geomagnetic storm events. Using the Bernese GNSS Software (BSW), reduced-dynamic orbits of different LEO satellites are computed from GNSS data of on-board receivers, where non-gravitational accelerations are modelled by means of estimated empirical piecewise-constant accelerations (PCAs). The orbital decay is then calculated by using the PCAs, or, in case of the GRACE Follow-On satellite, calibrated accelerometer data, to solve Gauss’s perturbation equation for the satellite’s semi-major axis. This method is compared to an approach where a fit model is applied to the osculating semi-major axis derived from the reduced-dynamic orbits computed by BSW. The fit model consists of a piece-wise linear model of the time-varying mean orbital decay and the time-varying amplitudes of the most dominant periodic oscillations. The results of both methods are compared and found to be similar for large orbital decays induced by CMEs. But the fit model struggles with low orbital decay. The Gaussian perturbation equation approach is far more precise than the fit model and can react, e.g., instantaneously to satellite maneuvers which change the semi-major axis. For satellites without an on-board accelerometer, PCAs can be estimated and used for the numerical integration.

An accurate knowledge of the orientation, the tidal deformability, and the gravity field of a celestial body is fundamental to provide constraints on its internal structure. These quantities may be retrieved by processing radiometric tracking and altimetry data from a probe in orbit around such body. This paper presents a method to combine altimetry crossovers with two-way Doppler tracking observations at normal equation level, using the Bernese GNSS Software and the pyXover software library. This method was applied to a proposed 200km altitude orbiter around Callisto, a privileged destination for the upcoming phase of Solar System exploration. Enhancing “standard” Doppler tracking with altimetry generally benefited both orbit determination and a joint estimation of the orientation of the north pole and of planetary librations. The retrieval of low-degree gravity field parameters was also improved by the addition of altimetry data. However, the improvements on the estimated parameters were highly dependent on the characteristics of the simulation, e.g., the underlying topography roughness. Overall, combining radioscience with altimetry data accounted for a visible reduction of correlations among estimated parameters, while also allowing for a consistent estimation of the “vertical” Love number h2 along with gravity.

The ESA GENESIS mission, which obtained green light at ESA's Council Meeting at Ministerial Level in November 2022 and which is expected to be launched in 2027, aims to significantly enhance the accuracy and stability of the Terrestrial Reference Frame (TRF). This shall be achieved by equipping one satellite at approximately 6000 km altitude with well-calibrated instruments for all four space-geodetic techniques contributing to TRF realizations, i.e., Global Navigation Satellite Systems (GNSS), Satellite Laser Ranging (SLR), Very Long Baseline Interferometry (VLBI) and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), and by exploiting the such realized very precise space collocations. The GENESIS satellite is foreseen to carry a zenith- and nadir-pointing GNSS antenna to track (at least) GPS and Galileo signals. Because of its very high altitude, GENESIS will cover much larger nadir angles as seen from the GNSS transmitting antennas, compared to receivers at ground or in low Earth orbit. This will partly result in GNSS observations with lower signal-to-noise ratios. Furthermore, to date, only little reliable information is available on the GNSS transmit antenna gain and carrier phase patterns at very large nadir angles. These problems lead to questions with respect to the best possible exploitation of GNSS data by GENESIS, e.g., whether phase pattern calibrations will need to be performed by means of tracked GNSS data (which would weaken the GNSS contribution to GENESIS). The aim of this study is to assess the impact of GNSS transmit antenna phase pattern errors on the GNSS-based POD of GENESIS, as well as global GNSS network solutions for GNSS orbits and clock corrections, Earth rotation and geocenter parameters and station coordinates based on GNSS observations from GENESIS and terrestrial stations. To accomplish this, we employ simulated GNSS pseudo-range and carrier phase data for GENESIS and ground stations, which have been generated based on detailed link-budget computations and a comprehensive set of transmit antenna gain patterns. The data are used for closed-loop simulation investigations, which allow to compare the reconstructed orbit and geodetic parameter solutions to the simulation truth and offer a quantification of the impact of transmit antenna phase pattern uncertainties on the estimated parameters.

To meet the demanding requirements in terms of accuracy and availability, GPS has introduced the signals on L5 that are compatible to Galileo E5a signals. The L5 signal was designed to mitigate the multipath and poor performance in harsh environments such as indoor, forests, and areas affected by jamming. As the L2 signal will become obsolete in the future, action must be taken to take advantage of the modern signal type which is currently broadcast by 19 out of GPS satellites. This is particularly important for some future LEO satellites (e.g. EPS-SG) which will rely exclusively on L1/L5. CODE (Center for Orbit Determination in Europe) is building up a prototype processing chain to generate L1/L5-based clock and bias products in addition to the classic L1/L2-based processing chain. First results regarding analysis product consistency are presented.

Satellite Laser Ranging (SLR) observations are provided by a global station network of the International Laser Ranging Service (ILRS). Since almost all SLR stations are unique in terms of utilized equipment, e.g., laser system, photo-detector or timing devices, their measurement performances may slightly differ. Furthermore, the tracked satellites have various properties, e.g., material composition (number and type of retro-reflectors or mantel material), diameter, area-to-mass ratio, altitude or inclination, which have an impact on the measurement precision and the requirements on the background force modeling. A reliable estimation of geodetic parameters solely based on SLR can only be performed by combining SLR observations to several spherical satellites provided by the entire ILRS network. To take the quality of each individual SLR measurement into account, several stochastic models, e.g., static or time-variable station- and/or satellitespecific weights, are introduced and their impact on the parameter estimation is studied.

The European Space Agency's BepiColombo mission continues its pioneering voyage to Mercury, the innermost planet of the Solar System. Among the advanced instruments onboard the Mercury Planetary Orbiter (MPO) is the Italian Spring Accelerometer (ISA), whose scientific objectives are closely linked to the Mercury Orbiter Radio-Science Experiment (MORE). Together, these instruments aim to provide high-precision data on the spacecraft's orbit, as well as Mercury's gravity field and internal structure. This simulation study investigates how the modeling and parametrization of accelerometer measurements influence orbit and gravity field recovery and explores strategies to overcome the challenges associated with the co-estimation of all parameters. We assess the integrated retrieval accuracy of the spacecraft orbit, gravity field, and accelerometer parameters under different noise levels and varying observation geometries during the mission. The MPO orbit is propagated, and Doppler and accelerometer measurements are simulated based on available noise models. Our results indicate that postponing the estimation of accelerometer biases until a preliminary gravity field is established helps prevent gravity field mismodelings from being absorbed into the accelerometer parameters. The daily estimation of bias was found to be essential. Using one year of Doppler tracking data and applying Kaula regularization, we demonstrate the potential to recover Mercury's gravity field up to degree and order 40. Low-degree coefficients improve under more optimistic noise conditions. Errors in cross-track bias estimation rise sharply when the beta-Earth angle falls below 45°, corresponding to the degraded Doppler observability. Orbit determination achieved centimetre-level radial accuracy and metre-level along- and cross-track accuracy, with increased cross-track errors during low-observability periods.

A growing number of Low Earth Orbiting (LEO) satellites are collecting GNSS tracking data that allows to recover the long-wavelength part of the Earth’s time-variable gravity field. Besides scientific LEO missions, commercial satellite constellations consisting of a huge number of nano-satellites are moving into focus. Due to an improved ground track coverage, such constellations offer the opportunity to increase the spatio-temporal resolution of derived gravity field models and can contribute to reduce temporal aliasing errors of dedicated gravity field missions. The Spire constellation is of particular interest as it consists of more than 100 nano-satellites (standardized CubeSats), all equipped with high-quality GNSS receivers. Furthermore, the Spire constellation offers a variety of orbital characteristics with different inclinations at altitudes of about 400–650 km. In this study, we use GNSS data from nine Spire CubeSats to derive monthly gravity field solutions covering a six-month period. The orbit and gravity field recovery is performed with the Bernese GNSS Software, which applies the Celestial Mechanics Approach. We demonstrate that the 2–3 times larger noise level of the Spire GNSS observations affects the quality of the retrieved gravity field solutions in the same order of magnitude. Therefore, a single Spire CubeSat solution cannot compete with those obtained from scientific LEO missions. However, with an increasing number of CubeSats, the performance improves so that a combination based on data from all nine Spire CubeSats can achieve a quality level comparable to a solution derived from ESA’s Swarm-B satellite.

We study temporal gravity field determination from GRACE Follow-On satellite-to-satellite tracking data using the inter-satellite link of the K-Band Ranging System (KBR) and kinematic positions of the satellites as observations and pseudo-observations, respectively. In addition, we introduce the range-rate observations collected by the Laser Ranging System (LRI) as further observation group. We compute our solutions with the Celestial Mechanics Approach using next to the kinematic positions either only LRI range-rate observations or combining them with the KBR derived range-rate data by the means of Variance Component Estimation (VCE).The stochastic noise of the KBR and LRI data is modelled with an empirical description of the noise based on the post-fit residuals between the final GRACE Follow-On orbits, that are co-estimated together with the gravity field, and the observations, expressed in position residuals to the kinematic positions and in KBR and/or LRI range-rate residuals. We validate the LRI-only and KBR+LRI monthly solutions by comparing them with the KBR-only derived models from the operational GRACE Follow-On processing at the AIUB, by examining the stochastic behaviour of respective post-fit residuals and by inspecting areas where a low noise is expected. Last but not least, we investigate the influence of LRI-only and KBR+LRI monthly gravity field solutions in a combination of monthly gravity fields based on other approaches as it is done by the Combination Service for Time-variable Gravity fields (COST-G) and make use of noise and signal assessment applying the quality control tools routinely used in the frame of COST-G.

Since the launch of the ESA Swarm satellites in 2014, GPS data has been used to monitor the monthly changes in Earth’s gravity field at a spatial resolution larger than 1,500 km (half-wavelength). We have chosen processing strategies that do not rely on assumptions of temporal and spatial correlations, producing models that are independent of GRACE and GRACE-FO data. We are a team composed of the Astronomical Institute of the University of Bern, the Astronomical Institute of the Czech Academy of Sciences, the Delft University of Technology, the Institute of Geodesy of the Graz University of Technology, and the School of Earth Sciences of the Ohio State University. We are supported by the European Space Agency and the Swarm Data, Innovation, and Science Cluster. Given the good health of the Swarm satellites, we will continue to produce high-accuracy hl-SST gravity field solutions in the foreseeable future.We produce individual gravity field models following independent gravity field inversion strategies. The International Combination Service for Time-variable Gravity Fields (COST-G) combines these individual models using weights derived with Variance Component Estimation, with the main objective of keeping the combined model unbiased to any individual solution. The models are published quarterly, pending the successful processing of kinematic orbits in the increasingly challenging environment resulting from the increased solar activity. The models are accessible at ESA’s Swarm Data Access server (https://swarm-diss.eo.esa.int) as well as at the International Centre for Global Earth Models (https://icgem.gfz.de/sp/02_COST-G_/Swarm).These data monitor geophysical processes in parallel to the ll-SST technique, offering the opportunity to validate the ll-SST models, act as an alternative in the event of gaps in their data records or to bridge the potential gap between the GRACE-FO and GRACE-C/MAGIC mission period.We show that the signal variability over the oceans of our models is a reliable measure of their accuracy by comparing it with the differences over land between Swarm and the GRACE and GRACE-FO solutions. Despite GRACE/GRACE-FO’s higher spatial resolution, we demonstrate that our Swarm models are able to resolve large-scale hydrological signals. Finally, we show that our effort has mitigated the challenges of estimating gravity field models derived from GPS data during the occurrence of high solar activities, which degraded the GPS signals.

Within the IGS, it was agreed that Precise Point Positioning (PPP) based on satellite orbit and clock corrections of the IGS analysis centers allow a direct access to the IGS realization of the International Terrestrial Reference Frame (ITRF). This convention is considered convenient for all PPP users and should not be changed in future. On the other hand, the groups determining the GNSS satellite orbits do prefer an origin of the frame that is related to the Earth instantaneous center of mass since this is the reference of the gravitational orbit force model. With this background, some discussions take place to change the convention related to the origin of the terrestrial frame because it is convenient for the orbit determination. This convention is, however, in contradiction to the expected needs of the PPP users that do prefer a stable coordinate origin in time. We will introduce a strategy to serve both needs by applying the center of mass corrections for the orbit determination only.

Satellite Laser Ranging (SLR) is essential for the geodetic parameter determination, e.g., geocenter and station coordinates, and, therefore, for long-term stable reference frame realizations. However, the orbit modeling and the quality of the parameter estimation partially depend on the background models. This study analyses the impact of static and time-variable a priori gravity field models provided by the Center for Space Research (CSR), the International Laser Ranging Service (ILRS), and the Combination Service for Time-variable Gravity Fields (COST-G) on the SLR data processing of spherical geodetic SLR satellites (LAGEOS-1/2 and LARES) at different orbital altitudes, by comparing the estimates of Earth rotation parameters, station coordinates and observation residuals. The COST-G model is further used to examine the impact of the mean pole model and the replacement of the spherical harmonic coefficients C_21/S_21 according to convention provided by the International Earth Rotation Service (IERS). While for LAGEOS-1/2 SLR data processing the a priori gravity field model has only a minor impact, the lower flying LARES satellite is more sensitive to the Earth’s gravity field and requires a more sophisticated gravity field modeling, e.g., COST-G Fitted Signal Model (FSM). In order to achieve higher consistency and thus improved solutions, the same mean pole model should be used in the SLR data processing as for the generation of the used a priori gravity field model. This study confirms the high quality of the COST-G FSM and demonstrates its suitability for potential use in the ILRS operational SLR processing.

Commonly, Global Positioning System (GPS)-based precise orbit determination (POD) of a low Earth orbiting (LEO) satellite is conducted by introducing fixed GPS orbits and clock corrections that were derived in a previous independent network solution using ground-based GPS receivers only. There have been a number of studies showing that the integration of LEO GPS observations can be advantageous in global network solutions, particularly when it comes to estimating geodetic parameters, such as the Earth’s center-of-mass. There has been an increase in GPS data available from LEO CubeSats over the past few years, such as from the Spire satellites, which are equipped with dual-frequency GPS receivers. The goal of the present study is to determine how GPS observations collected by specific Spire satellites affect global network solutions when combined with GPS observations from ground stations of the International GNSS Service. To obtain a combined GPS-LEO solution, GPS observations collected from receivers of selected Spire satellites and additional scientific LEOs are used to obtain combined GPS-LEO solutions by performing one joint least-squares adjustment process. Using this method, LEO satellite orbit parameters are determined along with GPS orbit parameters and geodetic parameters, namely GNSS station coordinates, Earth rotation parameters (ERPs), and the Earth’s center-of-mass. The influence of the Spire LEOs is assessed by analyzing different solutions that use GPS observations from different Spire and LEO satellites. Besides comparing the different global network solutions with one another, a comparison to a solution using only terrestrial GPS data is performed as well. As part of the quality assessment of the GPS-LEO solution, quality characteristics are analyzed, which are the geodetic parameters, orbit overlaps at arc boundaries, and formal errors determined for the estimated parameters. The results show that the estimation of the Earth’s center-of-mass coordinates benefits from including GPS-LEO observations. The solution is more improved by integrating GPS-LEO observations from scientific LEOs than from GPS-Spire-included solutions. The analysis reveals that GPS orbits are improved using the LEO-integrated approach, while orbit misclosures of LEO and Spire trajectories currently indicate a decrease in LEO orbit precision due to LEO orbit modeling deficiencies.