# Combination of altimetry crossovers and Doppler observables for orbit determination and geodetic parameter recovery: Application to Callisto

William Desprats, Stefano Bertone, Daniel Arnold, Martin Lasser, Adrian

Jäggi, Michel Blanc

► To cite this version:

William Desprats, Stefano Bertone, Daniel Arnold, Martin Lasser, Adrian Jäggi, et al.. Combination of altimetry crossovers and Doppler observables for orbit determination and geodetic parameter recovery: Application to Callisto. Acta Astronautica, 2025, 226, pp.585-600. 10.1016/j.actaastro.2024.10.045. hal-04938255

HAL Id: hal-04938255

https://hal.science/hal-04938255v1

Submitted on 12 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research paper

# Combination of altimetry crossovers and Doppler observables for orbit determination and geodetic parameter recovery: Application to Callisto

William Desprats<sup>a, *</sup>, Stefano Bertone<sup>b, c</sup>, Daniel Arnold<sup>a</sup>, Martin Lasser<sup>a</sup>, Adrian Jäggi<sup>a</sup>, Michel Blanc<sup>d</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

### Keywords:

Altimetry

Planetary

Geodesy

Orbit

## ABSTRACT

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 200 km 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 $h_2$ along with gravity.

## 1. Introduction

### 1.1. Characterization of Callisto's interior

Among the four Galilean moons, Callisto is probably the most recently formed, and the one with the smallest amount of accretion heating. This may have caused Callisto to be only partly differentiated [1], making it one of the most “pristine” objects in the Jupiter system. Moreover, its heavily cratered surface has kept records of the episodes of intense bombardment that characterized the first billion years of the Solar System. Callisto can then be considered a key “celestial record keeper” of the outer solar system. For such reasons, it has been identified in the “Planetary Exploration, Horizon 2061” foresight exercise as the logical preferred target for a future space mission dedicated to unraveling the origins of the Jupiter system [2]. ESA’s Jupiter Icy Moons Explorer (JUICE) [3] will be the next mission to visit Callisto with 21 flybys planned between 2031 and 2034, to characterize the outer shells of Callisto (including a putative sub-surface ocean [4]), to determine the composition of the non-ice material, and to study

past activity. However, an orbiter would enable global and uniform mapping of Callisto, as well as recovering its gravity field and other geodetic parameters to much higher resolution, thus achieving a more detailed characterization of the presence and extent of the ocean and its relation to the deeper interior (see, e.g., [5]). In this sense, several missions have been proposed to send an orbiter around Callisto, such as the Magnetics, Altimetry, Gravity, and Imaging of Callisto (MAGIC) mission [6] and the Jupiter component of the Tianwen-4 [7] mission led by the Chinese National Space Agency (CNSA) [8,9].

As it is usual for planetary missions, the telecommunications system can be used both for tracking and for gravity and radio science investigations. For instance, the Gravity and Geophysics of Jupiter and the Galilean Moons (3GM) experiment [10,11] onboard JUICE will contribute to the characterization of the interior state of Ganymede, informing on the presence of a deep ocean and on other gravity anomalies, and studying Ganymede and Callisto surface properties. On the other hand, an onboard laser altimeter may be used to characterize the topography and the tidal deformation of Callisto, as it will be

<footnote><sup>a</sup> Astronomical Institute, University of Bern, Switzerland</footnote>

<footnote><sup>b</sup> University of Maryland, College Park, MD, USA</footnote>

<footnote><sup>c</sup> INAF, Astrophysical Observatory of Torino, Italy</footnote>

<footnote><sup>d</sup> IRAP, CNRS-Université Paul Sabatier, Toulouse, France</footnote>

<footnote>* Corresponding author.
E-mail address: william.desprats@unibe.ch (W. Desprats).</footnote>
done for Ganymede with the GAnymede Laser Altimeter (GALA) [12] carried by JUICE. Synergies between these two instruments are expected. For instance, in the case of JUICE and Ganymede, the ice crust thickness will be constrained by determining the tidal response using altimetry [13] and gravity potential measurements [14], as well as by estimating the amplitudes of physical librations [15]. Furthermore, the combination of altimetry and gravity measurements would be used to assess current hypotheses of non-hydrostatic components in the gravity field of Callisto [5].

## 1.2. Scope of the study

This paper discusses the advantages of enhancing Doppler tracking with a laser altimeter on board of a low-altitude orbiter, in terms of orbit determination and of the recovery of geodetic parameters. Despite having considered an orbiter around Callisto, most of the results and conclusions also apply to other planetary bodies. The altimeter ranges were analyzed in the form of crossovers (see Section 2.3.2), which may be seen in simplified terms as an observation of the spacecraft's radial position in the generalized orbit determination process [16]. Altimetry crossovers have already been used in the framework of planetary orbiters such as Mars Global Surveyor (MGS) [17–19], Lunar Reconnaissance Orbiter (LRO) [20–22], and Kaguya [23]. The formulation used in the previously mentioned studies consists in projecting each track on a sloped terrain and fitting them using three-dimensional polynomials [17]. The computed distance between the polynomials corresponding to two intersecting ground tracks is thus the quantity to minimize, providing sensitivity to directions also other than the radial one.

The inclusion of altimetry crossovers in the orbit determination of MGS allowed to improve the spacecraft's trajectory and attitude [17, 18], but also the low-degree gravity field coefficients of Mars [19]. For LRO, the five distinct altimetry tracks from the multibeam altimeter allowed to significantly increase the number of crossovers [20], which resulted in an improved orbit during the first mission phases [21], and in improved solutions for the lunar tides and gravity field later on [22]. However, Mazarico et al. [24] did not use crossovers for the orbit determination of LRO, as the improved gravity field model and the accuracy of radio tracking reduced the presumed benefit of crossovers with respect to the high computational effort needed to include them. Finally, the higher-altitude Kaguya orbiter benefited from the inclusion of crossovers in the orbit determination procedure [23].

Previous simulation studies combining altimetry crossovers and Doppler observations for planetary orbiters have been carried out, e.g., for BepiColombo around Mercury [25,26], and for JUICE around Ganymede [27], which will both benefit from a highly accurate multi-link tracking. Marabucci [25], concluded that the inclusion of altimetry crossovers does not improve the polar orbit enough to justify the increase of computational complexity inherent to this observation type, partially because there are only very few crossovers at latitudes < 65°. On the other hand, Villamil et al. [27] showed significant orbit improvements thanks to the inclusion of altimetry crossovers, especially when the probe's line-of-sight is occulted from Earth (by Jupiter or by the Sun), or when the orbital plane is close to perpendicular to the Earth direction ($\beta_{Earth} \sim 90^{\circ}$). Moreover, they found that the tidal Love number $k_2$ and the low degrees of Ganymede's gravity field benefit from the inclusion of crossovers, and so does the estimation of biases of an onboard accelerometer instrument. Differences with respect to BepiColombo can be explained by a reduced Doppler tracking time (8 h/day for JUICE), and by a less polar orbit, resulting in more crossovers at low latitude. Nonetheless, their crossover formulation was valid only for nadir pointing instruments [27]. Additionally, none of these studies estimated the rotation state of the orbited body. In case of an orbiter around Callisto, comprehensive numerical simulations and covariance analyses have been performed for MAGIC mission scenarios [28], showing the expected estimates of Callisto's gravity

field, orientation and geophysical parameters, and how they would contribute to the determination of the moon's internal structure. An altimeter was also considered for MAGIC to measure tidal deformations, but no combined analysis of radio tracking and altimetry data has been presented.

In this paper, a 200 km altitude Callisto orbiter was modeled for 5 Callisto days (i.e., ~83 Earth days). Three different orbit configurations with respect to Earth's direction were analyzed following the procedure described in [29]. While most previous studies assessed the performance of different mission scenarios by covariance analyses, we performed full closed-loop simulations. A development version of the Bernese GNSS Software (BSW) [30] adapted for planetary probes [29, 31–33] was combined with pyXover, a Python suite of altimetry analysis tools for planetary geodesy [34], in order to precisely assess how well the orbit and geodetic parameters, i.e., rotation, orientation, and gravity field coefficients, as well as the $k_2$ and $h_2$ Love numbers, can be recovered.

# 2. Simulation setup and observable generation

## 2.1. Force model

The spacecraft's orbit is governed by the equations of motion which can be written, in a quasi inertial frame, as

$$ \ddot{\mathbf{r}} = -\mu_C \frac{\mathbf{r}}{r^3} + \mathbf{f}_p(t, \mathbf{r}, \dot{\mathbf{r}}, q_1, \dots, q_d), \quad (1) $$

where $\mathbf{r}$ is the position vector of the spacecraft center of mass with respect to the center of mass of Callisto, $r$ is its absolute value, a dot denotes derivative in time $t$, $\mu_C$ is the standard gravitational parameter of Callisto and $\mathbf{f}_p$ collects all perturbing accelerations beyond the central term (acceleration due to central point mass only). The $q_i$, $i = 1, \dots, d$, denote the parameters used to represent different kinds of orbit perturbations.

The accelerations induced by the gravity field of Callisto and by the Jupiter-induced tides are described in Sections 2.1.1 and 2.1.2, respectively. The perturbations of all Solar System planets and of the other Galilean moons are modeled by considering them as point masses, except for Jupiter for which zonal gravity field coefficients up to degree 6 have also been included [35]. We considered accelerations due to direct solar radiation pressure and indirect radiation pressure, such as planetary radiation pressure due to reflected and emitted radiation from the Callisto. The positions of Solar System planets were retrieved from the DE430 JPL ephemerides [36], and the position of the bodies within the Jovian system from the ephemerides JUP310 [35]. The rotation models of Jupiter and Callisto follow the recommendations of the International Astronomical Union (IAU) [37]. The orientation parameters of Callisto were estimated, namely the constant terms of right ascension and declination of its north pole $\alpha_0$ and $\delta_0$, respectively, its rotation rate $\omega_1$ and the amplitude of its main longitudinal libration $W_1$.

### 2.1.1. Synthetic gravity field

The gravity field of a celestial body is determined by its internal mass distribution. The acceleration acting on a spacecraft can be computed from the gradient of the gravitational potential $V$, which can be decomposed in spherical harmonics such that [38]

$$ V(r, \lambda, \phi) = \frac{\mu_C}{r} \cdot \sum_{n=2}^{n_{\max}} \sum_{m=0}^{n} \left( \frac{R_C}{r} \right)^n \bar{P}_{nm}(\sin \phi)(\bar{C}_{nm} \cos m\lambda + \bar{S}_{nm} \sin m\lambda), \quad (2) $$

where $\lambda$ and $\phi$ are longitude and latitude in a central body-fixed reference frame, $\bar{P}_{nm}$ are the fully normalized Legendre functions of degree $n$ and order $m$, and $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are the fully normalized spherical harmonic coefficients. $R_C = 2410.3$ km [37] is the equatorial radius of Callisto.

Since only the degree-2 gravity field coefficients have been so far estimated for Callisto [1], we generated a synthetic gravity field from
degree 3 to $n_{max}$ = 100. The gravity field coefficients were derived from the Earth Moon's gravity field model GRGM900C [39], appropriately scaled by the squared ratio of the surface gravity of Callisto and of the Moon. Due to the Moon's strong mass concentrations, the chosen gravity field creates significant orbit perturbations at low altitudes. It can thus be considered as a worst case scenario for the currently (mostly) unknown Callisto's gravity field. In the following, this synthetic reference gravity field model is labeled CALGLMo.

### 2.1.2. Solid tides

Jupiter's gravitational force produces a deformation of the Galilean moons. The amplitude of the radial deformation induced by tides is parametrized by the Love number $h_2$ (see Section 2.3.1). The tidal secondary potential induced by the mass redistribution is instead modeled via a change of the degree-2 gravity field coefficients [40], as

$$ \Delta \bar{C}_{2m} - i \Delta \bar{S}_{2m} = \frac{k_{2m}}{5} \frac{\mu_J}{\mu_C} \left( \frac{R_C}{r_J} \right)^3 \bar{P}_{2m}(\sin \phi_J) e^{-im\lambda_J}, \quad (3) $$

where $k_{2m} = k_2 = 0.3$ is the Love number for degree 2 and order $m = 0, 1, 2, r_J, \lambda_J$ and $\phi_J$ denote the spherical coordinates of Jupiter in a Callisto-fixed reference system, and $\mu_J$ denotes the gravity constant of Jupiter.

The mean corrections for the coefficients $\bar{C}_{20}$ and $\bar{C}_{22}$ are not zero. The time independent part of these corrections produces a permanent deformation. These permanent tidal corrections were separated from the main terms, so that

$$ \Delta \bar{C}_{20}^{\text{perm}} = -\frac{k_{20}}{2\sqrt{5}} \frac{\mu_J}{\mu_C} \left(\frac{R_C}{a_C}\right)^3, \quad (4) $$

$$ \Delta \bar{C}_{22}^{\text{perm}} = \frac{k_{22}}{2} \sqrt{\frac{3}{5}} \frac{\mu_J}{\mu_C} \left( \frac{R_C}{a_C} \right)^3, $$

with $a_C = 1883134$ km [41] the mean semi-major axis of the orbit of Callisto around Jupiter. Permanent tides can be included in the $\bar{C}_{20}$ and $\bar{C}_{22}$ coefficients of gravity field model of the Galilean moons. In this case, the gravity field is labeled “zero tide”, similarly to the Earth’s geopotential conventions [40]. Thus, before applying the corrections from Eq. (3), the $\Delta \bar{C}_{2m}^{\text{perm}}$ must be removed from $\bar{C}_{2m}$, with $m = 0, 2$, so that the permanent tide is counted only once.

### 2.2. Doppler observables

In this work, only one ground station tracked the probe, similarly to JUICE's baseline mission scenario. As for Europa Clipper [42], two-way X-band Doppler observables were modeled, but emitted and received from the Jiamusi ground station (46°29'36"N, 130°46'13"E) of the Chinese Deep Space Network (DSN) [43]. Radioscience observables were computed in the BSW following the formulation of Moyer [44] for un ramped two-way Doppler observables, including the Shapiro gravitational delay [45] on light propagation due to the gravity fields of the Sun, Earth and Jupiter. The delays due to Earth's troposphere, to Earth's ionosphere, and to the Io Plasma Torus (IPT), not considered in our simulations, were assumed to be fully calibrated. Current and future models of the IPT (e.g., [46]) would help to at least partially calibrate the signal delay [42,47]. The residual noise due to the media through which the signal travels is discussed in the following.

Various sources of noise affect Doppler observations [48,49], which may be characterized by their Allan deviation $\sigma_y(\tau)$ [50], where $\tau$ is the Doppler integration time (see Fig. 1). The propagation of the signal through media with a variable refractive index causes frequency fluctuations. The noise due to solar plasma highly depends on the angular distance of the line of sight between the Earth and the probe with respect to the Sun and in this case, i.e., for a Sun–Earth-probe angle $> 130^\circ$, $\sigma_y$ was found to be $< 8.5 \times 10^{-14}$ for two-way X-band Doppler observations at $\tau = 60$ s, using the model described in [51].

The refractive index is also fluctuating in the troposphere, mainly due to water vapor fluctuations [48]. Iess et al. [51] established a

noise model with a seasonal variation, adopted in this work. The noise coming from the wet troposphere is maximum in summertime and minimum in wintertime. Due to wind, gravitational loading, and thermal deformation, the antennas of the ground stations are also a source of noise, which can be evaluated at $\sigma_y^{GS} = 1.6 \times 10^{-14}$ for $\tau = 60$ s [49]. We assumed $\sigma_y^{S/C} = 1.8 \times 10^{-14}$ for all the spacecraft electronics, which are generally dominated by noise from the Deep Space Transponder (DST) [49]. The simultaneous use of X- and Ka-band Doppler tracking data, as expected for BepiColombo and JUICE, would allow for a nearly complete calibration of the dispersive noise, leading to a high tracking accuracy of $4 \times 10^{-14}$ for $\tau = 60$ s. However, this is not always possible (see, e.g., Europa Clipper [42]).

The simulated Doppler noise was thus added to the Doppler observations, which were generated under the following conditions: minimum spacecraft elevation of 10° above the horizon and no occultations of the spacecraft line of sight by the Sun, Jupiter, and Callisto.

### 2.3. Laser altimetry data

In addition to Doppler tracking data, realistic laser altimeter measurements were simulated from the spacecraft to the surface of Callisto by considering a synthetic topography model (see Section 2.3.1), using pyXover, a Python suite of altimetry analysis tools for planetary geodesy [34]. Although the formulation is valid regardless of the altimeter pointing, this work assumes only nadir pointing measurements to avoid the additional complexity of defining an attitude model for the spacecraft or the instrument. We followed the same formulation that has been previously used for crossover analyses of the MESSENGER probe at Mercury [52], to compute crossover discrepancies at Callisto (see Section 2.3.2). In order to combine the crossover analysis with the generalized Doppler orbit determination, the crossover discrepancies and their partial derivatives were extracted from pyXover (Section 2.3.2) and interfaced to the BSW routines.

#### 2.3.1. Altimetric range simulation

The laser altimeter range $LA(t)$ is modeled as the distance between the spacecraft on its orbit and the topography $T(\phi, \lambda)$ of the moon at the sub-satellite point,

$$ LA(t) = r(t) - T(\phi, \lambda) - u_r(\mathbf{r}, t) + \varepsilon_{\text{range}}, \quad (5) $$

where $\phi$ and $\lambda$ are the latitude and longitude of the probe at $t$, respectively. $u_r$ is the radial deformation due to Jupiter-induced tides, which can be expressed as (e.g., [53])

$$ u_r(\mathbf{r}, t) = \frac{h_2}{g} \frac{GM_J}{2} \frac{R_C^2}{d_J^2} (3 \cos^2 \psi - 1), \quad (6) $$

with $\psi$ being the angle between $\mathbf{r}$ and the direction from the center of the moon to Jupiter, $g$ the gravitational acceleration at the surface, $d_J$ the distance between the centers of mass of the Jupiter and Callisto, and $h_2$ the Love number of degree 2. $\varepsilon_{\text{range}}$ represents the white noise contributions added for a realistic simulation of the altimetric range. The main sources of errors for an altimeter are the instrument intrinsic errors $\varepsilon_{\text{range}}^{\text{instr}}$ and the errors in pointing and alignment $\varepsilon_{\text{range}}^{\text{pointing}}$. The standard deviation of $\varepsilon_{\text{range}}^{\text{instr}}$ is $< 0.2$ m in the case of the BepiColombo Laser Altimeter (BELA) [26,54] and GALA [12], while for $\varepsilon_{\text{range}}^{\text{pointing}}$, it depends on the altitude $h$, spacecraft pointing uncertainties $\Delta\Phi$, and surface slope $S$, so that

$$ \varepsilon_{\text{range}}^{\text{pointing}} = h \tan S \tan \Delta\Phi. \quad (7) $$

For a probe orbiting at $h = 200$ km, the standard deviation of $\varepsilon_{\text{range}}^{\text{pointing}}$ is $\sim 1$ m, assuming in average $S = 6^\circ$ and $\Delta\Phi = 10''$ (see, e.g., [27]). In light of the significantly larger errors originating from the interpolation of altimetric ranges in the presence of small-scale topography (see below), $\varepsilon_{\text{range}}^{\text{instr}}$ and $\varepsilon_{\text{range}}^{\text{pointing}}$ were neglected in this study.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Noise model for two-way X-band Doppler received at a station in the northern hemisphere at $\tau = 60$ s. The shaded area corresponds to the time span considered.</figcaption></figure>

The absence of highly accurate topographic models for the Galilean moons prevents using altimetric ranges as direct measurements [24, 55]. Nevertheless, the topography can be canceled at crossover locations by defining differential measurements between ranges to the same surface location at different epochs (i.e., crossover discrepancies, see Section 2.3.2). An important source of uncertainty in the processing of altimetry crossovers is related to interpolating between the bouncing points of the laser altimeter ranges and the crossover location on the surface. These uncertainties can be taken into account as part of the altimetric ranges noise budget (e.g., [27]) as they depend on topography, on the velocity of the probe with respect to the surface of Callisto, and on the sampling rate of the altimeter instrument. This approach is compared with the interpolation error induced by a detailed topography model in Section 5.

Since no Digital Elevation Model (DEM) for the Galilean moons is currently available for our simulations, we considered a global DEM of Mercury [56] for the similar-sized Callisto. However, the chosen global DEM presents only large-scale features with an effective resolution down to several km/pixel. To also take into account smaller-scale features, we simulated a small-scale topography for the surface (see, e.g., [52] Fig. 10, left, and description), with 10 m resolution and an elevation Root Mean Square (RMS) of 7.8 m. The total topography $T$ considered for the Galilean moon when simulating the altimetry data set is then the sum of the large-scale topography $T_{DEM}$ and of the small-scale topography $T_{local}$, so that

$$T(\lambda, \phi) = T_{DEM}(\lambda, \phi) + T_{local}(\lambda, \phi). \qquad (8)$$

It is important to note that $T_{local}$ was based on past analyses at Mercury, which might not represent the actual (currently unknown) small-scale topography of Callisto. The magnitude of small-scale topography has a direct impact on the interpolation error for the crossover analysis. For instance, a smoother small-scale topography would yield a lower interpolation error.

Altimetric ranges were simulated using the pyXover software package, based on spacecraft orbits integrated in the BSW (see Section 2.1), with a 10 Hz sampling, consistent with BELA [57]. We also performed a few tests at 30 Hz, as foreseen for GALA [12] on board of JUICE. However, such a higher sampling increases the amount of data, and thus the computational load and related challenges.

### 2.3.2. Crossover discrepancy

A crossover point is a location at which two distinct ground-tracks, i.e., the projections of the orbit on the ground surface along the line-of-sight of the altimeter, intersect. The search for the location of all

crossovers during the mission duration is central to this analysis, as all combinations between all ground-tracks have to be analyzed. This is efficiently performed in pyXover, as detailed in [52]. Differential measurements at crossover locations are particularly interesting [16, 17], as they are due to a combination of errors in the spacecraft orbit and attitude, interpolation errors of the surface topography between the altimeter footprints, and geophysical signals (e.g., mismodellings in orientation or vertical tides model).

The search and computation of the crossover observables start with geolocating the laser altimetry bouncing points on the surface of the moon. Each altimetry measurement gets assigned a set of latitude $\phi$, longitude $\lambda$ and elevation $\eta$, corresponding to the radial component of the geolocalized range measurements. Then, the partial derivatives with respect to the orbit parameters, the rotation parameters $\alpha_0$, $\delta_0$, $w_1$ and $W_1$, and the tidal parameter $h_2$ are computed along each track by finite differencing (except for $h_2$, whose derivatives are computed analytically, see [52]). In pyXover, the orbit parameters are orbit corrections modeled as 3 constant offsets estimated on each track in the radial, along-track and cross-track directions.

Finally, the intersections between the ground tracks are identified by interpolating between the bouncing points, and the crossover discrepancy $v$ between two tracks A and B is computed as the difference in elevation at the crossover point, such that

$$v = \eta_A - \eta_B. \qquad (9)$$

This discrepancy is modeled as a function of errors in the a priori orbit, deviations from the rotation model, and mismodelling of tidal deformations. In an ideal case, one would expect $v = 0$. The crossover observables and their partial derivatives computed in pyXover, can either be processed as in [52] within pyXover to find corrections to the considered parameters, or they can now be exported to the BSW for a joint analysis with Doppler observables, as presented in this study (see Section 3.2). Since the BSW uses a different orbit parametrization, i.e., initial osculating orbital elements for each arc, the partial derivatives had to be adapted accordingly.

For each parameter $p_j$, the total derivative of the crossover discrepancy $v_i$ can be written as

$$\frac{dv_i}{dp_j} = \frac{\partial v_i}{\partial \Delta_{i,1}} \frac{\partial \Delta_{i,1}}{\partial p_j} + \frac{\partial v_i}{\partial \Delta_{i,2}} \frac{\partial \Delta_{i,2}}{\partial p_j} + \frac{\partial v_i}{\partial p_j}, \qquad (10)$$

with $\Delta_{i,1}$ and $\Delta_{i,2}$ the orbit offsets of respectively the first and second track of the crossover $i$. Using the chain rule, we derive

$$\frac{dv_i}{dp_j} = \frac{\partial v_i}{\partial \Delta_{i,1}} \frac{\partial \Delta_{i,1}}{\partial \mathbf{r}_1} \frac{\partial \mathbf{r}_1}{\partial p_j} + \frac{\partial v_i}{\partial \Delta_{i,2}} \frac{\partial \Delta_{i,2}}{\partial \mathbf{r}_2} \frac{\partial \mathbf{r}_2}{\partial p_j} + \frac{\partial v_i}{\partial p_j}, \qquad (11)$$
<caption>**Table 1**</caption>

List of estimated parameters. The osculating orbital elements are estimated at the start of each of the 80 ~25 h-arcs.

<table><thead><tr><th>Parameter</th><th>Number</th><th>Type</th></tr></thead><tbody><tr><td>Osculating orbital elements</td><td>6</td><td>Local</td></tr><tr><td>Tidal Love number k<sub>2</sub></td><td>1</td><td>Global</td></tr><tr><td>Tidal Love number h<sub>2</sub></td><td>1</td><td>Global</td></tr><tr><td>Gravity field coefficients (d/o 100)</td><td>8278</td><td>Global</td></tr><tr><td>Pole orientation offsets a<sub>0</sub>, δ<sub>0</sub></td><td>2</td><td>Global</td></tr><tr><td>Rotation rate w<sub>1</sub></td><td>1</td><td>Global</td></tr><tr><td>Main libration amplitude W<sub>1</sub></td><td>1</td><td>Global</td></tr></tbody></table>

with $\mathbf{r}_k = \mathbf{r}(t_k)$. The position of the spacecraft is reconstructed from the orbit offset according to $\mathbf{r}_k = \mathbf{r}_{0,k} + \mathbf{R}_k \Delta_{i,k}$, with $\mathbf{R}_k$ denoting the rotation matrix from the local orbital frame to the inertial frame at $t_k$, and $\mathbf{r}_{0,k}$ being the a priori value of $\mathbf{r}_k$. We then have

$$ \frac{\partial \Delta_{i,k}}{\partial \mathbf{r}_k} = \frac{\partial \mathbf{R}_k^T}{\partial \mathbf{r}_k} (\mathbf{r}_k - \mathbf{r}_{0,k}) + \mathbf{R}_k^T = \frac{\partial \mathbf{R}_k^T}{\partial \mathbf{r}_k} \mathbf{R}_k \Delta_{i,k} + \mathbf{R}_k^T. \quad (12) $$

The total derivatives (Eq. (11)) are evaluated with all the parameters set to their a priori value, i.e., in particular by setting $\Delta_{i,k} = 0$ in Eq. (12), so that

$$ \frac{dv_i}{dp_j} = \frac{\partial v_i}{\partial \Delta_{i,1}} R_1^T \frac{\partial \mathbf{r}_1}{\partial p_j} + \frac{\partial v_i}{\partial \Delta_{i,2}} R_2^T \frac{\partial \mathbf{r}_2}{\partial p_j} + \frac{\partial v_i}{\partial p_j}. \quad (13) $$

The partial derivatives $\frac{\partial v_i}{\partial \Delta_{i,k}}$ are computed in pyXover (see Section 2.3.2) together with the partials $\frac{\partial v_i}{\partial p_j}$ for $p_j \in \{\alpha_0, \delta_0, w_1, h_2\}$

and $\frac{\partial v_i}{\partial p_j} = 0$ for other parameters. The epochs of the crossovers, the crossover discrepancies, and their partial derivatives were exported from pyXover to XOVNEQ, a program within the BSW developed in the frame of this work. XOVNEQ reads the observation files from pyXover and then uses the previously integrated orbit to compute the partials $\frac{\partial \mathbf{r}_i}{\partial p_j}$ as detailed in Section 3.1, and then the total derivatives according to Eq. (13). The resulting Normal Equation Systems (NEQs) can then be combined with other observables within the BSW.

### 3. Combining radioscience and altimetry data for generalized orbit determination on normal equation level

The estimation of the orbit and geodetic parameters considered in this work (summarized in Table 1) is based on two software packages: the BSW [30], and pyXover [34]. This section describes our approach to interface and combine altimetry crossovers and Doppler observables across them, following the pipeline depicted in Fig. 2.

#### 3.1. Celestial mechanics approach

Orbit and geodetic parameters were estimated following the Celestial Mechanics Approach (CMA) [58], which has been successfully applied to planetary probes [29,31–33]. This approach is divided into an initial arc-wise orbit determination followed by a generalized orbit improvement.

The initial arc-wise orbit determination was performed in GRAVDET2, a program from a development branch of the BSW. The length of an arc was not fixed and varied between 24 h and 28 h as in [29]. These arc lengths provide a good balance between cumulative propagation error and dynamical consistency of the orbit for gravity field recovery. The equations of motion (1) were numerically integrated along with the variational equations for the six initial osculating orbital elements (the semi-major axis *a*, the eccentricity *e*, the inclination *i*, the longitude of the ascending node *Ω*, the argument of perigee *ω*, and the argument of latitude *u*<sub>0</sub>), based on the fixed set of background force models detailed in Section 2.1. The partial derivatives $\frac{\partial \mathbf{r}}{\partial p}$ of the orbit

with respect to the other parameters listed in Table 1 were then computed by solving definite integrals as described in [59]. This allowed for a first Doppler-based least squares adjustment of local parameters (i.e., the six initial osculating orbital elements, see Table 1) to minimize the observation residuals. This process is generally repeated for at least 4 iterations, and until the relative RMS change of residuals between one iteration and the previous one is below 0.5%.

Based on the updated orbits, the equations of motion and the variational equations were once more integrated in GRAVDET2 to compute the partial derivatives with respect to both local and global parameters (Table 1), thus allowing us to set up the NEQs for all parameters and for all arcs separately. The arc-wise NEQs were stacked for the whole mission span (e.g., 5 Callisto days), using another BSW program named ADDNEQ2 (Fig. 2). The accumulated NEQ was inverted to retrieve corrections to the global parameters as well as to the local parameters, along with their formal uncertainties. The complete procedure was iterated until convergence was reached.

#### 3.2. Combination of altimetry crossover and Doppler observables

There are two interface points to connect the two software packages, as depicted in Fig. 2. First, the orbits generated in the BSW were exported to pyXover, to simulate altimetry data in the pyAltsim program (see Section 2.3), and to geolocate the altimetry bouncing point in the pyGeoloc program. Second, all the relevant information concerning the crossovers computed in the pyXover program (location, discrepancies, partial derivatives, etc.) were exported to XOVNEQ, a BSW routine.

The generalized orbit determination, based on a non-linear least-square adjustment, is an iterative process (see Section 3.1). At every iteration, the probe's orbit is integrated from updated initial conditions and force models. Unlike Doppler observations, the crossover discrepancies and their locations change with respect to the orbit and the rotation model. Additionally, ground tracks may no longer intersect from one iteration to the other, or vice versa, new intersections might be found. At every iteration, a computationally heavy search for crossovers is thus, in principle, required, which can be eased by leveraging crossover locations from the previous iteration. For these reasons, the first step was to fit orbits using Doppler observables only within the GRAVDET2 program until convergence, which constituted the baseline for the results presented in this paper. The number of iterations requiring the computation of crossovers was therefore reduced as suggested, e.g., in [27]. The improved orbit was finally exported to pyXover.

In pyGeoloc, the orbit is divided into several segments, called gtracks. Each gtrack corresponds to one revolution of the probe around Callisto. The gtracks are separated into two hemispheres, so that pyXover can later search for crossovers between each gtrack separately in the northern hemisphere and in the southern hemisphere. In the end, each gtrack is the projection of a part of the orbit, up to half a period long, on the Callisto reference ellipsoid. One of the advantages offered by this splitting is the possibility to extensively parallelize this geolocation step for each gtrack. In this work, altimetry data were grouped in batches of 7 days for effective parallelization. For each altimetry observation in the gtrack object, pyGeoloc computes the latitude *φ*, longitude *λ*, and elevation *η* of the laser spot bouncing point on the surface, based on the input orbit and on light-time modeling. Additionally, the partial derivatives of *φ*, *λ* and *η* with respect to the estimated parameters are computed by finite differencing (except for $h_2$, see Section 2.3.2).

Once the gtracks were processed, combinations between the weekly batches were formed. pyXover thereby searched for crossovers within each of these combinations, and for each hemisphere. The search for crossovers is divided into two steps. A first rough search employs sub-sampled altimetry data points to limit the computational load. In case of a 10 Hz altimeter sampling rate operating at 200 km altitude above Callisto, we considered 1 observation every 20. The distance between
<figure><img src="image_4.png" /><figcaption>Fig. 2. Flowchart for Doppler and altimetry crossovers combination.</figcaption></figure>

two sub-sampled points is then 0.076°, which corresponds to 3.2 km on the surface of Callisto. The altimetry bouncing points are projected onto the surface using a polar stereographic projection adapted for each hemisphere.

The second step is a fine search based on the results of the rough search. The neighboring fully sampled data points were projected around the “roughly guessed” crossover to more precisely locate the crossover by interpolating between these points using a cubic spline. The coordinates of the crossovers $\phi$ and $\lambda$ were obtained from the interpolation, while the crossover discrepancy $v$ was calculated from the difference between the elevation of the surface $\eta$ interpolated along each of the two tracks. The partial derivatives of $v$ with respect to all estimated parameters were computed similarly using finite differences. More details of the entire process can be found in [52].

The crossover discrepancies and their partial derivatives were subsequently exported to XOVNEQ, which used the orbits previously generated by GRAVDET2 (and used to find the crossovers), to compute crossover NEQs for each combination. In particular, XOVNEQ computed the partial derivatives according to the parametrization used by the BSW, based on the parametrization used in pyXover for orbit corrections (see Section 2.3.2). For this study, each crossover observation was assigned a unit weight. However, our setup would also accommodate a more advanced weighting scheme, such as assessing the interpolation error based on the distance of a given crossover with respect to the surrounding altimetry bouncing points (e.g., [52]).

The Doppler NEQs and the crossover NEQs were then stacked within the ADDNEQ2 program to form a global NEQ to be solved, as for the Doppler-only case. There are several ways to stack these multiple NEQs. In this work, NEQs were first stacked over time for each observation type separately, thus retrieving one “global” Doppler NEQ and one “global” crossover NEQ, which were then stacked together with a given weight. Variance Component Estimation (VCE) (see, e.g., [60]) was particularly useful to determine “optimal” weights, but we also experimented using empirically determined weights, which is discussed later.

## 4. Doppler-only solution

We assumed the 200 km-altitude probe to be tracked only by the Jiamusi ground station of the Chinese DSN, for a maximum tracking time of 6 h/day at X-band, for 5 Callisto days (i.e., ~83 days or 80 arcs of ~25 h). Three scenarios were chosen, corresponding to different angles between the orbital plane and the Earth direction, i.e., the $\beta_{Earth}$ angle, which would not vary significantly during the chosen 5 Callisto days, starting on 1-May-2031 and with an orbit inclination of 89°. The first orbit is nearly edge-on, with $\beta_{Earth} = 1.0° \pm 4.2°$. In such

<figure><img src="image_5.png" /><figcaption>Fig. 3. Ground coverage of Callisto from a polar orbiter with $h = 200$ km occulted by Callisto for 83 days, for $\beta_{Earth} = 1.0° \pm 4.2°$. The 60 s-sampled Doppler observations are grouped in bins of 1.8° in longitude and 2.8° in latitude.</figcaption></figure>

a configuration, the probe would pass regularly behind Callisto with respect to Earth, which would prevent the probe from being directly tracked by ground stations on Earth (see Fig. 3). In the other two scenarios, the probe would never be occulted by Callisto as seen from Earth with $\beta_{Earth} \ge 67°$. The second scenario’s orbit is at the limit of this threshold, with $\beta_{Earth} = 71.8° \pm 4.7°$, while the last scenario corresponds to a near face-on configuration with $\beta_{Earth} = 84.0° \pm 2.6°$. We simulated the acquisition of 18 413, 29 305 and 29 429 two-way Doppler X-band observations, for $\beta_{Earth} = 1°, 72°, \text{ and } 84°$, respectively, and we present in this section the results for a Doppler-only solution, using the CAL-GLMo model up to degree and order (d/o) 100 as a priori gravity field model and considering a perfect force modeling. This choice allowed us to avoid a lengthy iterative retrieval process, as shown in [29], where they started the procedure with a degree-2 gravity field as a priori information, while still being able to estimate a full gravity field model within a few iterations from Doppler measurements.

### 4.1. Orbit parameters

The determination of a probe’s orbit around a celestial body using radio tracking from Earth is known to have degeneracies [61,62], which can be explained by geometrical symmetries in the system [63]. The weakly determined components of the position $\mathbf{r}$ with respect to the central body are contained in a plane orthogonal to the line
<figure><img src="image_6.png" /><figcaption>Fig. 4. Daily RMS of the orbit differences for $\beta_{Earth} = 1°, 72°$ and $84°$ considering only Doppler observations.</figcaption></figure>

of sight from Earth to the central body. This results in the direction perpendicular to the line of sight and to the position vector being the worst determined direction. In case of an edge-on orbit, the latter corresponds to the cross-track direction. This was confirmed by the daily RMS of the orbit differences between the true and the fitted orbit shown in Fig. 4. In the other two cases, the error was spread between the along-track and cross-track directions.

Wood [61] showed that some orbital elements may be indeterminate. Some configurations proved to be worse than others, e.g., when $\beta_{Earth}$ is close to 90° as shown by the large global RMS of the orbit differences of 7.0 m for $\beta_{Earth} = 84°$. Additionally, Russell and Thurman [62] quantified the influence of $\beta_{Earth}$ and showed that orbit errors increase with decreasing $\beta_{Earth}$. However, this was not clear when looking only at the global 3D orbit RMS of orbit differences, which is 1.9 m, 4.6 m and 7.0 m for the experiments at $\beta_{Earth} = 1°, \beta_{Earth} = 72°$ and $\beta_{Earth} = 84°$, respectively. Nevertheless, correlations between the estimated initial orbital elements were insightful. Fig. 5 represents the absolute correlations between each orbit parameter within each arc averaged over the 80 arcs and shows significant correlations between the osculating orbital elements in all three cases. The correlations were distributed across all orbital elements for low $\beta_{Earth}$. The overall lowest correlations were for the case $\beta_{Earth} = 72°$, except for the correlations between inclination $i$ and the argument of latitude $u_0$. The near face-on orbit suffered from very high correlations between $i, u_0$ and the longitude of the ascending node $\Omega$. For $\beta_{Earth} = 1°$ the correlations between the inclination $i$ and both the longitude of the ascending node $\Omega$ and the argument of latitude $u_0$ were significantly lower, < 0.5, compared to values > 0.9 in the other two cases. Such high correlations between these parameters may explain why the orbit accuracy was worse for $\beta_{Earth} = 72°$ than for $\beta_{Earth} = 1°$.

## 4.2. Global solution

Gravity field parameters were estimated up to d/o 100, together with the Love number $k_2$, the right ascension $\alpha_0$ and declination $\delta_0$ of the north pole and the amplitude of the main longitudinal libration $W_l$. The formal errors of the estimated gravity field coefficients (see Fig. 6) are theoretically independent from the a priori gravity field. Thus, choosing a different field in this study would mainly change the maximum degree for which the degree-amplitude formal errors are lower than the signal (i.e., the resolution of the gravity field solution).

As discussed in [29], the observation geometry from Earth is critical for gravity field recovery when radio tracking observations are used. An edge-on orbit ($\beta_{Earth} = 1°$) provides a better sensitivity to the gravity field, since most of the gravity signal is aligned with the line of sight from Earth. However, this benefit needs to be balanced with occultations of the line-of sight by Callisto, as seen in Fig. 3. This

effect, combined with the limited number of Doppler observations (only one tracking station) resulted in a gravity field solution suffering from aliasing. When estimating gravity field coefficients up to d/o 100, the signal from the high degree coefficients, which cannot be resolved because of the depleted ground coverage, is absorbed by the lower degree coefficients (see blue curve in Fig. 6).

Constraining the high degree coefficients to zero using, e.g., a Kaula law, limited aliasing due to the ground coverage (Fig. 6). Each coefficient of degree $n$ was constrained by $\sigma_K(n) = \frac{K}{n^2}$, with the constant $K = 0.018$ derived using VCE as in [64]. On the other hand, a nearly face-on orbit (e.g., $\beta_{Earth} = 84°$) is detrimental in terms of gravity field recovery, as the gravity signal acting on the probe is mostly perpendicular to the line-of-sight, resulting in a lower Doppler sensitivity. However, despite using a Kaula-constraint in the case $\beta_{Earth} = 1°$, the gravity field solution for the near face-on case performed better than the edge-on case in the low degrees. This may be explained by the very low number of observations for $\beta_{Earth} = 1°$, and may change, e.g., in case the probe would be tracked by more than one tracking station [29]. The intermediate case $\beta_{Earth} = 72°$ benefited from the absence of occultations along the line-of-sight from Callisto, and from a better sensitivity to the gravity signal than a face-on orbit, thus resulting in a better unconstrained gravity field solution.

The solutions for the orientation parameters and Love number $k_2$ are presented in Section 6. In general, the rotation rate $w_1$ was not estimated, as the short mission time span considered did not allow to reduce the high correlations with the non-zonal gravity field coefficients without applying tight constraints on $w_1$ (see Fig. 7). Nevertheless, in Section 5 a few solutions including the estimation of $w_1$ are presented to investigate how the combination with altimetry helps reducing correlations.

## 5. Impact of topography roughness on noise and parameter estimation

This section presents the added value of altimetry data on the recovery of a Doppler-fitted near edge-on orbit ($\beta_{Earth} = 1.0° \pm 4.2°$). We considered three different topography models (no topography at all, only a DEM, and synthetic small-scale topography on top of a DEM), two different altimeter sampling rates (10 Hz or 30 Hz), and a “control case” where all topography was replaced by a simple 12 m-white noise added to the altimetry ranges to match the crossover discrepancy distribution of the 10 Hz case. In addition, the following shows how the procedure was validated and discusses the robustness of the estimation of Callisto’s rotation rate.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Absolute correlations between initial osculating orbital elements averaged over 80 arcs for different $\beta_{Earth}$ using only Doppler observables.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Difference (solid) and error (dashed) degree amplitudes of gravity field solutions for different $\beta_{Earth}$ using only Doppler observations. A Kaula constraint (noted w/ K.) was considered in case of the nearly edge-on orbit, to compensate for the limited ground coverage. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Distribution of the crossovers along latitudes $\phi$ for $\beta_{Earth} = 1^{\circ}$. 95% of the crossovers are located at latitudes $>70^{\circ}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Correlations between the low-degree gravity field coefficients and rotation parameters using only Doppler observations.</figcaption></figure>

### 5.1. Distribution of the crossover discrepancies

As described in Section 3.2, altimetry data were grouped in batches of 7 days, so that the 83 days of data were split in 12 weeks, resulting in 760 gtracks in each hemisphere (one gtrack may be split into two gtracks if the ground track spreads over two different weeks). The 12 weeks of data then resulted in a total of 78 combinations. As the orbit is quasi-polar, most of the ~522 000 crossovers were located at latitudes >70° (see Fig. 8).

Crossover discrepancies can be due to several factors: orbit modeling errors, a geophysical signal (e.g., tidal deformation), but also, as discussed in Section 2.3, intrinsic errors in the instruments and pointing errors. Moreover, a potentially significant error is introduced when interpolating the ground tracks between the altimetry bouncing points and the crossover coordinates. This interpolation error is strongly dependent on topography. In case of a 10 Hz sampling rate, the inter-spot distance, i.e., the distance between the bouncing points, is 160 m for the orbits considered in this paper. The distance between the closest bouncing point and the crossover cannot be larger than 80 m (half the inter-spot distance) and this distance was found to be in average 26 m. The distributions of the pre-fit crossover discrepancies are presented in the following, for different topography models, different altimeter frequencies, and with or without considering additive white noise. The standard deviations of the crossover discrepancies are summarized in Table 2.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Distribution of the crossover discrepancies *v* for *β*<sub>Earth</sub> = 1° without topography or with a DEM at 10 Hz (left), with small-scale topography or with 12 m white noise on the DEM at 10 Hz (center), with small-scale topography at 10 Hz or 30 Hz (right).</figcaption></figure>

The first case is without any topography, i.e., an ellipsoid, and without other sources of noise, in order to validate the combination procedure. As expected, the crossover discrepancies were very low in this ideal case, with a standard deviation of 1 cm (see Fig. 9 left). As soon as a DEM with maximum effective resolution of several hundred meters to a few km was introduced, still without noise, the discrepancies increased by one order of magnitude. This difference is due to the variation of altitude between two bouncing points (i.e., the slope).

For a more realistic case, a simulated small-scale topography was added on top of large-scale features modeled by the DEM. The elevation of a given point on the surface of Callisto is then the sum of both elevations. As a result, the interpolation errors increased up to 20 m, as shown in Fig. 9 (center), because of surface roughness [52,65]. An additional case was examined without any topography, but with a white noise equivalent to the small-scale topography overlaid to the ellipsoid, to evaluate the impact of a structured topography on the recovery of orbit and geodetic parameters. As shown in Fig. 9 (center), the two distributions of pre-fit crossover discrepancies were similar but not identical, and their standard deviation for the white noise case was slightly lower.

The interpolation error depends on the topography, on the velocity of the probe with respect to the surface of Callisto, and on the sampling rate of the altimeter. The higher the probe, or the smoother the topography, the smaller the interpolation error. To test its impact on interpolation error and on the solution, the sampling rate of the altimeter was increased with the same orbit and topography. Increasing the sampling rate of the altimeter from 10 Hz (as for BELA [57]) to 30 Hz (as for GALA [12]) reduced the distance between the bouncing points, and thus reduced the a priori crossover discrepancies, as expected (see Fig. 9 right). The standard deviation of the discrepancies was reduced from 8.3 m to 6.4 m. However, because of the larger amount of data and computational burden, only one test was performed at 30 Hz with small-scale topography, keeping a 10 Hz sampling for the tests with different topography models.

## 5.2. Reduction of the orbit differences

Once the crossovers were located, and once their discrepancies and partial derivatives were computed, NEQs were derived in XOVNEQ for each of the 78 week-combinations, as described in Section 3.2. The global Doppler NEQ and the global crossovers NEQ were subsequently combined, mainly following VCE-provided relative weights. However, in more ideal cases (i.e., without any noise on the altimetry data and either without topography or with only a DEM) increasing the weight given by VCE was found to favor the crossovers dataset does help

<caption>**Table 2**</caption>

Observation residuals and orbit differences for several orbit and topography scenarios. The RMS of Doppler residuals is 0.776 mHz in the three orbit cases with Doppler only, while the RMS of the differences are 1.95 m for $β_{Earth} = 1°$.

<table><thead><tr><th rowspan="2">Topography</th><th rowspan="2">Weight ratio</th><th colspan="2">Stdev v [m]</th><th rowspan="2">Dop. [mHz]</th><th rowspan="2">Orb. [m]</th></tr><tr><th>Pre-fit</th><th>Post-fit</th></tr></thead><tbody><tr><td>No topography</td><td>1 × 10<sup>3</sup></td><td>0.00980</td><td>0.000259</td><td>0.787</td><td>0.350</td></tr><tr><td>DEM</td><td>1</td><td>0.123</td><td>0.0890</td><td>0.784</td><td>0.846</td></tr><tr><td>Small-scale 30 Hz</td><td>1.2 × 10<sup>-3</sup></td><td>6.40</td><td>6.34</td><td>0.777</td><td>0.934</td></tr><tr><td>Small-scale 10 Hz</td><td>6.9 × 10<sup>-4</sup></td><td>8.33</td><td>8.32</td><td>0.777</td><td>1.39</td></tr><tr><td>White noise</td><td>7.7 × 10<sup>-4</sup></td><td>7.84</td><td>7.84</td><td>0.776</td><td>1.93</td></tr></tbody></table>

reducing post-fit orbit differences. The ratios of the weights given to altimetry crossover discrepancies with respect to the weights given to the Doppler observations are summarized in Table 2. They were empirically determined in the three ideal cases, and derived following VCE in the other cases.

Considering a large-scale DEM already limited the orbit improvement given by crossovers with respect to a Doppler-only fit to ~50%, from ~80% when not considering topography. It is important to note that these tests are not realistic, since no source of noise was considered, which would have further reduced the benefit of altimetry crossovers with respect to more realistic Doppler observations. However, these results constitute a first necessary validation of the simulation setup. Only in the case of no topography or when using only a DEM, the crossover discrepancies were visibly reduced after the combination, at the price of a slight increase of the Doppler post-fit residuals.

In the more realistic cases, i.e., with small-scale topography or with additive white noise, the post-fit Doppler residuals and crossover discrepancies were overall very close to the values after the Doppler-only fit. However, the orbits were visibly improved (see Fig. 10) with respect to the Doppler-only solution when using a simulated small-scale topography, with the RMS of post-fit differences being reduced by ~26% (see Table 2). Increasing the sampling rate by a factor of 3 improved the orbit by a significant 33% with respect to the 10 Hz case. There was no visible orbit improvement in the near edge-on case when using white noise instead of a fixed small-scale topography, as expected. Indeed, the discrepancy of a perfectly modeled crossover located precisely at the bouncing points of two perfect ground tracks would have been affected only by the range noise, as the topography would have canceled (see Eq. (5)), and interpolation would not have been needed. Thus, if the interpolation error were part of the range noise budget, and not a consequence of the considered small-scale topography, it would affect
<figure><img src="image_12.png" /><figcaption>Fig. 10. Daily RMS of the orbit differences for $\beta_{Earth} = 1^{\circ}$ considering only Doppler observations (D) and combining with crossovers (X) for different topography models.</figcaption></figure>

the error on the estimated parameters by minimizing the crossover discrepancies.

Because of the nature of the altimetry crossovers data set, one could expect the main improvement to happen in the radial direction. However, because a Doppler-fitted orbit is most weakly determined in the direction perpendicular to the radial direction (see Section 4), the orbit differences in the radial direction were already smaller (<5 cm) than in the along-track and cross-track directions. Moreover, even small improvements in the radial direction induce improvements in the other directions. In particular for an orbit close to an edge-on configuration, the direction weakly determined by Doppler observations was the cross-track direction, with orbit differences >5 m in case of a Doppler-only solution (see Fig. 10).

### 5.3. Estimation of the geodetic parameters

After looking at the benefit of a combination between altimetry crossovers and Doppler observables in terms of orbit parameters, we present in this section its impact on the estimation of rotation and orientation parameters, and on the estimation of gravity field parameters and tidal Love numbers. Only the near edge-on case ($\beta_{Earth} = 1.0^{\circ} \pm 4.2^{\circ}$) is presented here, for which the depleted ground coverage may harm a freely estimated gravity field Doppler-only solution, as discussed in Section 4, and for which the coverage is even more depleted when considering Doppler observations from only one tracking station (Jiamusi).

#### 5.3.1. Rotation and orientation parameters only

We first evaluated the free estimation of the global parameters as performed in `pyXover` with altimetry data only, namely the constant terms of the right ascension and declination of the north pole $\alpha_0$ and $\delta_0$, the rotation rate of Callisto $w_1$, the main longitudinal libration $W_1$, and the tidal Love number $h_2$. Next, crossovers were combined with Doppler observations and, similarly to Section 5.2, an optimal weight between the observation types was determined using VCE, except for the ideal cases without small-scale topography and where no additional noise is considered, for which the weight was fixed to the values given in Table 2.

For each of the cases described in Section 5.2, the formal errors of the previously mentioned parameters are given in Table 3 and the true errors (between the estimated and the true values) are given in Table 4. We compared the errors when using only Doppler observations (except for the Love number $h_2$), when using altimetry crossovers only, and when combining both observation types (X+D). Formal and true errors are additionally represented in Fig. 11. The combination of altimetry crossover discrepancies and Doppler observables resulted in a general improvement of the formal errors with respect to each observation type taken separately, except for $h_2$ in the ideal case of no topography.

<caption>Table 3</caption>

Formal errors of the estimated rotation and orientation parameters and $h_2$ Love number, considering different observation types and different topography models for $\beta_{Earth} = 1^{\circ}$.

<table><thead><tr><th>Topography</th><th>Obs. type</th><th>α<sub>0</sub> [°] ×10<sup>4</sup></th><th>δ<sub>0</sub> [°] ×10<sup>4</sup></th><th>w<sub>1</sub> [°/day] ×10<sup>9</sup></th><th>W<sub>l</sub> [°] ×10<sup>5</sup></th><th>h<sub>2</sub> [–] ×10<sup>3</sup></th></tr></thead><tbody><tr><td></td><td>Doppler</td><td>5.1</td><td>0.23</td><td>41</td><td>0.58</td><td>–</td></tr><tr><td rowspan="2">No topography</td><td>Crossover</td><td>0.0017</td><td>0.0012</td><td>0.046</td><td>0.011</td><td>0.0018</td></tr><tr><td>X+D</td><td>0.0045</td><td>0.0012</td><td>0.036</td><td>0.0044</td><td>0.0049</td></tr><tr><td rowspan="2">DEM</td><td>Crossover</td><td>0.024</td><td>0.041</td><td>0.55</td><td>0.075</td><td>0.17</td></tr><tr><td>X+D</td><td>0.016</td><td>0.033</td><td>0.13</td><td>0.060</td><td>0.12</td></tr><tr><td rowspan="2">Small-scale (30 Hz)</td><td>Crossover</td><td>1.4</td><td>3.0</td><td>14</td><td>5.4</td><td>9.5</td></tr><tr><td>X+D</td><td>0.074</td><td>0.16</td><td>0.62</td><td>0.43</td><td>9.3</td></tr><tr><td rowspan="2">Small-scale (10 Hz)</td><td>Crossover</td><td>2.2</td><td>4.4</td><td>27</td><td>7.4</td><td>14</td></tr><tr><td>X+D</td><td>0.089</td><td>0.16</td><td>0.76</td><td>0.46</td><td>12</td></tr><tr><td rowspan="2">White noise σ = 12 m</td><td>Crossover</td><td>2.1</td><td>4.1</td><td>18</td><td>6.9</td><td>14</td></tr><tr><td>X+D</td><td>0.090</td><td>0.16</td><td>0.77</td><td>0.47</td><td>11</td></tr></tbody></table>

Regarding the true errors in Table 4, the benefit of a combined solution was also generally visible. When no topography was considered, the true errors in the estimated parameters did not necessarily improve. This may be due to the weight ratio between observation types, derived empirically without estimating rotation parameters in Section 5.2. The true error of $\delta_0$ increased for a few topography cases when combining crossovers with Doppler observations. Moreover, when considering small-scale topography and 10 Hz altimetry data, the true error on $h_2$ was larger for the combined solution than when considering crossovers only. However, in order to unambiguously analyze the true errors of this small set of parameters, one would need to repeat the experiment with different initial conditions, for instance by performing a Markov Chain Monte-Carlo experiment. Considering a 30 Hz altimeter sampling rate improved the formal errors by ~30% and up to 50% for $w_1$ with respect to a 10 Hz altimeter in the case of a crossovers only solution, and between 0% and 30% in the case of combinations of Doppler and crossovers. Differences between small-scale topography and pure white noise in the case of 10 Hz were more subtle. Indeed, as seen in Fig. 9 (center), the distribution of the crossovers was comparable, but not identical. The altimetry-only solution was worse in the case of small-scale topography, while the differences were similar when combining altimetry crossovers with Doppler observations.

#### 5.3.2. Inclusion of the gravity field coefficients and $k_2$

The same experiment was then repeated, but freely estimating the gravity field parameters and the $k_2$ Love number. Altimetry crossovers
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** True (filled bars) and formal (empty bars and solid line) errors on the estimate of rotation and orientation parameters, and of the Love number $h_2$, when considering either Doppler observations only, or when combining Doppler and crossover observations with different topography models.</figcaption></figure>

<caption>**Table 4**</caption>

True errors of the estimated rotation and orientation parameters and $h_2$ Love number, considering different observation types and different topography models for $\beta_{Earth} = 1^{\circ}$.

<table><caption>**Table 5**</caption><thead><tr><th>Topography</th><th>Obs. type</th><th>α<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>δ<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>w<sub>1</sub> [°/day]<br>×10<sup>9</sup></th><th>W<sub>f</sub> [°]<br>×10<sup>5</sup></th><th>h<sub>2</sub> [-]<br>×10<sup>3</sup></th></tr></thead><tbody><tr><td></td><td>Doppler</td><td>0.41</td><td>0.21</td><td>3.6</td><td>0.39</td><td>-</td></tr><tr><td rowspan="2">No topography</td><td>Crossover</td><td>0.15</td><td>0.08</td><td>0.014</td><td>0.17</td><td>0.0067</td></tr><tr><td>X+D</td><td>0.46</td><td>0.0087</td><td>3.7</td><td>0.014</td><td>0.0055</td></tr><tr><td rowspan="2">DEM</td><td>Crossover</td><td>0.41</td><td>0.59</td><td>2.3</td><td>0.37</td><td>0.7</td></tr><tr><td>X+D</td><td>0.14</td><td>0.34</td><td>0.86</td><td>0.069</td><td>0.14</td></tr><tr><td rowspan="2">Small-scale<br>(30 Hz)</td><td>Crossover</td><td>1.4</td><td>1.2</td><td>17</td><td>0.023</td><td>9.3</td></tr><tr><td>X+D</td><td>0.18</td><td>0.17</td><td>1.7</td><td>0.39</td><td>8.9</td></tr><tr><td rowspan="2">Small-scale<br>(10 Hz)</td><td>Crossover</td><td>2.5</td><td>7.7</td><td>8.1</td><td>0.4</td><td>0.5</td></tr><tr><td>X+D</td><td>0.24</td><td>0.23</td><td>1.9</td><td>0.14</td><td>5.6</td></tr><tr><td rowspan="2">White noise<br>σ = 12 m</td><td>Crossover</td><td>2.5</td><td>1.9</td><td>13</td><td>6.8</td><td>1.8</td></tr><tr><td>X+D</td><td>0.036</td><td>0.22</td><td>0.28</td><td>0.24</td><td>1.6</td></tr></tbody></table>

Formal errors of the estimated rotation and orientation parameters, and of the $h_2$ and $k_2$ Love numbers, estimated along with gravity field parameters when considering only Doppler observations and when combining observation types with different topography models.

<table><thead><tr><th>Topography<br>Obs. type</th><th>α<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>δ<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>w<sub>1</sub> [°/day]<br>×10<sup>6</sup></th><th>W<sub>f</sub> [°]<br>×10<sup>5</sup></th><th>h<sub>2</sub> [-]<br>×10<sup>3</sup></th><th>k<sub>2</sub> [-]<br>×10<sup>3</sup></th></tr></thead><tbody><tr><td>Doppler only</td><td>14</td><td>0.81</td><td>24</td><td>21</td><td>-</td><td>4.2</td></tr><tr><td>No topography</td><td>0.015</td><td>0.01</td><td>0.0013</td><td>0.043</td><td>0.015</td><td>0.0046</td></tr><tr><td>DEM</td><td>0.16</td><td>0.12</td><td>0.017</td><td>0.88</td><td>0.37</td><td>0.012</td></tr><tr><td>Small-scale (30 Hz)</td><td>1.2</td><td>0.53</td><td>0.77</td><td>7.5</td><td>9.8</td><td>0.19</td></tr><tr><td>Small-scale (10 Hz)</td><td>1.3</td><td>0.54</td><td>1</td><td>8.6</td><td>13</td><td>0.52</td></tr><tr><td>White noise</td><td>1.3</td><td>0.54</td><td>0.98</td><td>8.4</td><td>12</td><td>0.64</td></tr></tbody></table>

discrepancies alone did not allow for an estimation of the gravity field up to d/o 100. We therefore compared the Doppler only case (except for the $h_2$ Love number) and the different combinations between altimetry crossovers and Doppler observations, in terms of formal errors in Table 5 and of true errors in Table 6, which are also represented in Fig. 12.

As discussed in Section 4, the ground coverage available from Doppler observations was limited in the case of a near edge-on orbit, which can lead to an aliasing of the gravity field solution, unless properly constrained. Additionally, the rotation rate $w_1$ was highly

<caption>**Table 6**</caption>

True errors of the estimated rotation, and orientation parameters and $h_2$ and $k_2$ Love numbers, considering only Doppler observations and when combining observation types with different topography models, estimated with gravity field parameters.

<table><thead><tr><th>Topography<br>Obs. type</th><th>α<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>δ<sub>0</sub> [°]<br>×10<sup>4</sup></th><th>w<sub>1</sub> [°/day]<br>×10<sup>6</sup></th><th>W<sub>f</sub> [°]<br>×10<sup>5</sup></th><th>h<sub>2</sub> [-]<br>×10<sup>3</sup></th><th>k<sub>2</sub> [-]<br>×10<sup>3</sup></th></tr></thead><tbody><tr><td>Doppler only</td><td>22</td><td>1.7</td><td>23</td><td>11</td><td>-</td><td>6.5</td></tr><tr><td>No topography</td><td>0.12</td><td>0.083</td><td>0.0027</td><td>0.22</td><td>0.057</td><td>0.053</td></tr><tr><td>DEM</td><td>0.90</td><td>0.19</td><td>0.038</td><td>3.7</td><td>0.6</td><td>0.22</td></tr><tr><td>Small-scale 30 Hz</td><td>1.1</td><td>0.81</td><td>0.68</td><td>5.3</td><td>6.7</td><td>0.20</td></tr><tr><td>Small-scale 10 Hz</td><td>1.4</td><td>0.91</td><td>0.69</td><td>5.6</td><td>10</td><td>0.13</td></tr><tr><td>White noise</td><td>1.1</td><td>0.94</td><td>1.6</td><td>16</td><td>6.2</td><td>1.5</td></tr></tbody></table>

correlated with the non-zonal coefficients of the gravity field and its unconstrained estimation proved to be challenging (see Section 4). Nevertheless, we chose to analyze only unconstrained solutions in this section, to highlight the contribution of altimetry crossovers to make up for the depleted ground coverage, and to reduce correlations between parameters in this challenging scenario.

The estimation of all the rotation and orientation parameters improved in all the combinations considered with respect to a Doppler-only estimation, in terms of formal errors, but also in terms of true errors. The only exception was the true error of the main longitudinal libration amplitude $W_f$ when considering white noise only. However, as previously discussed, the true errors of a limited set of parameters of a single test are indicative but not conclusive. The formal errors of the constant term of the right ascension of the north pole $α_0$ were reduced by at least one order of magnitude, thanks to correlations between $α_0$ and $δ_0$ being reduced by at least a factor 4.8 (Fig. 13). The improvement of $δ_0$ was more moderate as the formal errors were reduced by 30% only in the small-scale topography and white noise cases alike.

Overall, the ideal cases without topography or with only DEM drastically reduced the errors in all the parameters considered, and also correlations between the rotation rate $w_1$ and the gravity field coefficients (Fig. 13 left). The correlation between $α_0$ and $δ_0$ was still significantly reduced when considering small-scale topography, from -0.63 to -0.08 (see Figs. 7 and 13, right). $w_1$ remained highly correlated with gravity field coefficients, showing that the benefits of altimetry are weakened when considering small-scale topography. However, lower correlations with respect to a Doppler-only solution were still noticeable. For instance the correlations between $w_1$ and $C_{21}$, $S_{21}$, or $C_{30}$ were reduced by a factor of at least 4, leading to lower correlations
<figure><img src="image_14.png" /><figcaption>Fig. 12. True (filled bars) and formal (empty bars and solid line) errors of the estimation of rotation, and orientation parameters and $h_2$ and $k_2$ Love numbers, considering different observation types and different topography models, estimated with gravity field parameters.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Correlations between the low-degree gravity field coefficients and rotation parameters estimated considering Doppler and crossovers observables with a DEM and without (left) or with (right) small-scale topography.</figcaption></figure>

between these coefficients and the other gravity field coefficients (see Fig. 13, right).

Formal errors on the $k_2$ Love number estimates were also reduced, by at least a factor 8 in the case of small-scale topography and 10 Hz altimetry. The gravity field solution significantly improved in the ideal cases (see Fig. 14). In the case of small-scale topography with a 10 Hz sampling rate, the gravity field solution was still improved with respect to a Doppler-only solution. Unlike the orbit differences discussed in Section 5.2, there were no large differences in the estimates of the geodetic parameters when considering small-scale topography or an equivalent white noise.

In these cases, gravity field and rotation parameters benefited most from a combination of altimetry crossovers and Doppler observations. The depleted ground coverage was at least partially mitigated, and the correlations arising from the estimation of rotation and orientation parameters were reduced. Orbit improvements discussed in Section 5.2 clearly contributed to the improved estimates of these parameters.

In this section, all parameters were estimated freely. However, constraining the rotation parameters as discussed in Section 4 would help further reducing correlations between the rotation rate $w_1$ and gravity field parameters.

# 6. Evaluation of the contribution of altimetry for different orbit configurations

This section discusses the added value of altimetry data on three Doppler-only fitted orbits representing the different geometries presented in Section 4. Only a 10 Hz altimeter sampling rate was considered with small-scale topography added on top of a DEM, and Callisto's rotation rate was not estimated, because of the high correlations with the gravity field parameters.

## 6.1. Improvement on the orbital elements

The relative weights between altimetry crossover discrepancies and Doppler observations derived following VCE as described in Section 3.2 are given in Table 7, along with the post-fit standard deviation of the crossover discrepancies, and the RMS of the orbit differences. The post-fit Doppler residuals and crossover discrepancies were overall very close to the values obtained after the Doppler-only fit. However, the orbits visibly improved with respect to the Doppler-only solution, with the global RMS being reduced by ~26% and up to 36% in the near face-on case (see Table 7). Our result that altimetry crossovers benefit the solution most for $\beta_{Earth}$ close to 90° agrees with the findings of [27].
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Difference (solid) and error (dashed) degree amplitude of gravity field solutions considering different topography models, using only Doppler observations (Dop) or a combination of crossovers and Doppler (X+D) observations.</figcaption></figure>

<caption>**Table 7**</caption>

Crossover discrepancy residuals and orbit differences for different orbits. The RMS of Doppler residuals is ~0.78 mHz in the three orbit cases with Doppler only and with or without adding altimetry crossovers.

<table><thead><tr><th rowspan="2">β<sub>Earth</sub><br>avg. [°]</th><th rowspan="2">Weight<br>ratio</th><th colspan="2">Stdev v [m]</th><th colspan="2">RMS orb. diff. [m]</th></tr><tr><th>Pre-fit</th><th>Post-fit</th><th>Pre-fit</th><th>Post-fit</th></tr></thead><tbody><tr><td>1</td><td>6.9 × 10<sup>-4</sup></td><td>8.33</td><td>8.32</td><td>1.95</td><td>1.39</td></tr><tr><td>72</td><td>6.5 × 10<sup>-4</sup></td><td>8.53</td><td>8.52</td><td>4.64</td><td>3.26</td></tr><tr><td>84</td><td>6.5 × 10<sup>-4</sup></td><td>8.57</td><td>8.56</td><td>6.99</td><td>4.58</td></tr></tbody></table>

As discussed in Section 4, and as shown in Fig. 5 for a Doppler-only solution, the osculating orbital elements are highly correlated with each other for both low and high $β_{Earth}$. However, correlations were notably reduced when the altimetry crossovers were combined with Doppler observations (see Fig. 15). For the near face-on orbit ($β_{Earth} = 84°$ in Fig. 15), the correlations between the osculating orbital elements were visibly reduced, but the correlations between $i, Ω$, and $u_0$ remained large. In case $β_{Earth} = 72°$, the improvement in correlations between orbital parameters was more moderate. However, the correlations between the semi-major axis $a$ and $i$, and between $Ω$ and $u_0$, slightly increased.

## 6.2. Estimation of the complete set of geodetic parameters

Doppler and altimetry crossovers observations were combined using VCE to estimate the gravity field coefficients along with the rotation parameters (except for the rotation rate) for the three orbit configurations. Since the gravity field solution still suffered from aliasing in case $β_{Earth} = 1°$, we used VCE to derive a Kaula constraint to regularize the gravity field solution to a level comparable with the more face-on orbit configurations, similarly to the Doppler-only solution (Section 4.2). The Kaula constant $K$ was estimated after combining the observations NEQs, resulting in $K = 0.00057$ (see Fig. 16).

The freely estimated rotation parameters improved in all cases with respect to a Doppler only solution (Table 8). The right ascension of the north pole $α_0$ significantly improved when adding crossovers to Doppler observations as the formal errors were reduced by one order of magnitude in case $β_{Earth} = 1°$, and were divided by at least a factor 7 in case of more face-on configurations. The improvement in declination $δ_0$ was more moderate, but still appreciable, around 30%. More significant improvements in the main libration amplitude $W_1$ and

<caption>**Table 8**</caption>

True errors $Δ$ and formal errors $f_{err}$ of the estimation of rotation, and orientation parameters and $h_2$ and $k_2$ Love numbers, with gravity field parameters.

<table><thead><tr><th>β<sub>Earth</sub></th><th>Obs</th><th>α<sub>0</sub> [°]<br>×10<sup>5</sup></th><th>δ<sub>0</sub> [°]<br>×10<sup>5</sup></th><th>W<sub>1</sub> [°]<br>×10<sup>5</sup></th><th>h<sub>2</sub> [-]<br>×10<sup>3</sup></th><th>k<sub>2</sub> [-]<br>×10<sup>3</sup></th></tr></thead><tbody><tr><td rowspan="3">1°</td><td>D</td><td>Δ<br>f<sub>err</sub></td><td>220<br>144</td><td>17.2<br>8.13</td><td>11<br>21.4</td><td>-<br>-</td><td>6.34<br>4.2</td></tr><tr><td>D + X</td><td>Δ<br>f<sub>err</sub></td><td>13.8<br>13.4</td><td>9.17<br>5.42</td><td>6.31<br>8.54</td><td>10.1<br>12.9</td><td>0.127<br>0.518</td></tr><tr><td>Kaula</td><td>Δ<br>f<sub>err</sub></td><td>151<br>122</td><td>14.2<br>7.29</td><td>2.12<br>12.6</td><td>-<br>-</td><td>0.436<br>3.42</td></tr><tr><td rowspan="3">72°</td><td>D</td><td>Δ<br>f<sub>err</sub></td><td>11.2<br>12</td><td>8.13<br>5.25</td><td>9.5<br>6.91</td><td>4.33<br>12.6</td><td>0.176<br>0.48</td></tr><tr><td>D + X</td><td>Δ<br>f<sub>err</sub></td><td>38.8<br>46.4</td><td>0.32<br>2.92</td><td>0.864<br>2.32</td><td>-<br>-</td><td>2.56<br>1.45</td></tr><tr><td>84°</td><td>Δ<br>f<sub>err</sub></td><td>2.78<br>5.89</td><td>0.67<br>2.22</td><td>0.0366<br>1.91</td><td>8.14<br>12.7</td><td>0.272<br>0.568</td></tr><tr><td rowspan="3">84°</td><td>D</td><td>Δ<br>f<sub>err</sub></td><td>3.48<br>39.7</td><td>0.738<br>2.54</td><td>4.04<br>1.82</td><td>-<br>-</td><td>0.564<br>0.333</td></tr><tr><td>D + X</td><td>Δ<br>f<sub>err</sub></td><td>1.13<br>5.7</td><td>1.5<br>1.77</td><td>3.03<br>1.38</td><td>29.6<br>13</td><td>0.473<br>0.189</td></tr></tbody></table>

Love number $k_2$ in the edge-on case were visible, with an improvement of ~50% and ~85%, respectively, in comparison with ~30% and ~40% in case $β_{Earth} = 72°$ and $β_{Earth} = 84°$. The formal errors on $h_2$ were overall 0.013, regardless of the $β_{Earth}$ angle. It can be noted that the other parameters also benefit from the Kaula constraint by means of a reduced aliasing of gravity field coefficients (see Table 8).

Moreover, correlations between rotation and gravity field parameters decreased with the addition of altimetry data (see Supplementary material). In the case of $β_{Earth} = 1°$, correlations between the gravity field parameters are still present (see Supplementary material), but mainly because of the poor coverage discussed in Section 4. The Kaula constraint helped reducing these residual correlations, thus bringing the formal errors of the low degree gravity field coefficients in the range of the combined solution for a face-on orbit.

## 7. Conclusion

This article presented the benefits of combining Doppler-only solutions with altimetry crossovers using the BSW software and the pyXover library. We proposed a novel implementation to combine these heterogeneous data sets at the normal equation level. We investigated the estimation of Callisto's orientation parameters, Love numbers, and gravity field parameters, which are key ingredients in providing additional constraints on the internal structure of the Jupiter moon. We simulated a 200 km orbiter around Callisto, tracked in X-band only by the Jiamusi ground station, part of the Chinese DSN. Assuming only a 3 months mission, the Doppler based solution was affected by large correlations between the rotation rate of Callisto and the gravity field parameters, all freely estimated.

Our results indicated that the addition of altimetry crossovers leads to a general improvement of the orbit, as well as of the uncertainties on the recovered rotation and gravity field parameters. All orbit geometries in this study benefited from the presence of an on-board altimeter, even though the benefit was more evident in terms of orbit accuracy for face-on orbits and in terms of rotation parameters for edge-on orbits. However, this study showed how such improvements are highly dependent on the ground truth found at destination, e.g., the moon's topography roughness. Moreover, the addition of altimetry data would allow for a consistent estimation of the Love number $h_2$ along with a gravity model, and also for a significant reduction of correlations among estimated parameters. We expect the benefits related to the addition of altimetry crossovers to increase for a longer mission, as the number of crossovers would drastically increase.
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Absolute correlations between initial osculating orbital elements averaged over 80 arcs for different $\beta_{Earth}$ fitted using Doppler observables and altimetry crossover observations with small-scale topography.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Difference degree amplitude of gravity field solutions considering different $\beta_{Earth}$, using only Doppler observations (dotted line) or a combination with altimetry crossovers (solid line).</figcaption></figure>

It is important to note that the relative benefit of altimetry may be reduced if a larger number of Doppler observations were available, e.g., with more ground stations actively tracking the probe, or if the noise on the Doppler observations were lower. On the other hand, the chosen Solar System configuration led to Doppler tracking signals propagating in a quiet solar plasma environment (cf. Fig. 1). The Doppler noise can be significantly larger for time periods during which the probe is closer to the Sun as seen from Earth, which would enhance the benefit of altimetry crossovers. Additionally, the contribution of altimetry is very sensitive to the interpolation error, which depends on the topography and on the geographical distribution of data, and which would vary with the celestial body and with a different orbit altitude. Altimetry contribution to the solution would also be enhanced by a more advanced weighting scheme, e.g., a latitude-dependent weight or setting a higher weight to crossovers corresponding to holes in the Doppler coverage. These would be interesting follow-up studies to this work.

In conclusion, these investigations provided guidelines for combining Doppler and altimetry crossovers datasets in the case of a planetary orbiter. It showed the benefits of such a combination when evaluating the performance of multiple mission scenarios in terms of orbit determination and the recovery of geodetic parameters. The provided methods and schemes can be readily applied to the analyses of future

missions towards the Jovian moons, as well as adapted to other celestial bodies. For instance, the methodology and tools developed in this work could be applied to the analysis of real data from forthcoming planetary orbiters such as BepiColombo and JUICE, both carrying a laser altimeter on board.

### CRediT authorship contribution statement

**William Desprats:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Conceptualization. **Stefano Bertone:** Writing – review & editing, Supervision, Software, Methodology, Conceptualization. **Daniel Arnold:** Writing – review & editing, Supervision, Investigation, Funding acquisition. **Martin Lasser:** Writing – review & editing, Software, Investigation. **Adrian Jäggi:** Writing – review & editing. **Michel Blanc:** Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

Calculations were performed on UBELIX, the HPC cluster at the University of Bern (http://www.id.unibe.ch/hpc). This study has been supported by the Swiss National Science Foundation (SNSF), Switzerland grant #200021_185056 "Callisto geodesy: A simulation study to support further space missions to the Jovian system". SB acknowledges support by NASA, USA under award number 80GSFC24M0006.

### Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.actaastro.2024.10.045.

### References

[1] J. Anderson, R. Jacobson, T. McElrath, W. Moore, G. Schubert, P. Thomas, Shape, mean radius, gravity field, and interior structure of Callisto, Icarus 153 (1) (2001) 157–161, http://dx.doi.org/10.1006/icar.2001.6664.

[2] M. Blanc, P.W. Bousquet, V. Dehant, B. Foing, M. Grande, L. Guo, A. Hutzler, J. Lasue, J. Lewis, M.-A. Perino, H. Rauer, Planetary Exploration Horizon 2061: A Long-Term Perspective for Planetary Exploration, Elsevier, 2023, http://dx.doi.org/10.1016/C2020-0-01626-9.
[3] O. Grasset, M.K. Dougherty, A. Coustenis, E.J. Bunce, C. Erd, D. Titov, M. Blanc, A. Coates, P. Drossart, L.N. Fletcher, H. Hussmann, R. Jaumann, N. Krupp, J.-P. Lebreton, O. Prieto-Ballesteros, P. Tortora, F. Tosi, T.V. Hoolt, JUpler ICy moons explorer (JUICE): An ESA mission to orbit Ganymede and to characterise the Jupiter system, Planet. Space Sci. 78 (2013) 1-21, http://dx.doi.org/10.1016/j.pss.2012.12.002.

[4] K. Khurana, M. Kivelson, D. Stevenson, G. Schubert, C. Russell, R. Walker, C. Polanskey, Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto, Nature 395 (6704) (1998) 777-780, http://dx.doi.org/10.1038/27394.

[5] D.E. Smith, F. Nimmo, K. Khurana, C.L. Johnson, M. Wieczorek, M.T. Zuber, C. Paty, A. Genova, E. Mazarico, L. Prockter, G.A. Neumann, J.E. Connerney, E.B. Bierhaus, S.J. Goossens, M.K. Barker, P.B. James, J. Head, Brown, Callisto: A guide to the origin of the Jupiter system [White paper], 2020, http://dx.doi.org/10.3847/25c2cfeb.8ef223f3.

[6] D.E. Smith, M.T. Zuber, R.M. Canup, T. Hurford Jr., B. Bierhaus, M.T. Bland, V.J. Bray, J.E. Connerney, A. Genova, C. Johnson, et al., MAGIC, A discovery proposal to the icy moon Callisto, in: AGU Fall Meeting Abstracts, Vol. 2019, 2019, pp. P34C-03.

[7] Y. Wei, Z. Yao, W. Wan, China's roadmap for planetary exploration, Nat. Astron. 2 (5) (2018) 346-348, http://dx.doi.org/10.1038/s41550-018-0456-6.

[8] M. Blanc, C. Wang, L. Li, M. Li, L. Wang, Y. Wang, Y. Wang, Y. Wang, Q. Zong, N. Andre, O. Mousis, D. Hestroffer, P. Vernazza, Gan De: Science objectives and mission scenarios for China's mission to the Jupiter System, in: EGU General Assembly Conference Abstracts, 2020, p. 20179, http://dx.doi.org/10.5194/egusphere-egu2020-20179.

[9] Q. Wang, Prospect of international cooperation in China lunar & planetary exploration programs, in: 73rd International Astronautical Congress (IAC 2022), 2022.

[10] M. Di Benedetto, P. Cappuccio, S. Molli, L. Federici, A. Zavoli, Analysis of 3GM Callisto gravity experiment of the JUICE mission, Planet. Space Sci. (2021) http://dx.doi.org/10.48550/ARXIV.2101.03401.

[11] P. Cappuccio, M. Di Benedetto, D. Durante, L. Iess, Callisto and Europa gravity measurements from JUICE 3GM experiment simulation, Planet. Sci. J. 3 (8) (2022) 199, http://dx.doi.org/10.3847/psj/ac83c4.

[12] H. Hussmann, K. Lingenauber, R. Kallenbach, K. Enya, N. Thomas, L.M. Lara, C. Althaus, H. Araki, T. Behnke, J.M. Castro-Marin, H. Eisenmenger, T. Gerber, M. Herranz de la Revilla, C. Hüttig, K. Ishibashi, J. Jiménez-Ortega, J. Kimura, M. Kobayashi, H.-G. Lötzke, A. Lichopoj, F. Lüdicke, I. Martínez-Navajas, H. Michaelis, N. Namiki, H. Noda, J. Oberst, S. Oshigami, J.P. Rodríguez García, J. Rodrigo, K. Rösner, A. Stark, G. Steinbrügge, P. Thabah, S. del Togno, K. Touhara, S. Villamil, B. Wendler, K. Wickhusen, K. Willner, The Ganymede laser altimeter (GALA): key objectives, instrument design, and performance, CEAS Space J. 11 (4) (2019) 381-390, http://dx.doi.org/10.1007/s12567-019-00282-8.

[13] G. Steinbrügge, A. Stark, H. Hussmann, P. Sohl, J. Oberst, Measuring tidal deformations by laser altimetry: A performance model for the JUICE Laser Altimeter, Planet. Space Sci. 117 (2015) 184-191, http://dx.doi.org/10.1016/j.pss.2015.06.013.

[14] P. Cappuccio, A. Hickey, D. Durante, M. Di Benedetto, L. Iess, F.D. Marchi, C. Plainaki, A. Milillo, A. Mura, Ganymede's gravity, tides and rotational state from JUICE'S 3GM experiment simulation, Planet. Space Sci. 187 (2020) 104902, http://dx.doi.org/10.1016/j.pss.2020.104902.

[15] G. Steinbrügge, T. Steinke, R. Thor, A. Stark, H. Hussmann, Measuring Ganymede's librations with laser altimetry, Geosciences 9 (7) (2019) 320, http://dx.doi.org/10.3390/geosciences9070302.

[16] C. Shum, B. Schutz, B. Tapley, B. Zhang, Altimeter crossover methods for precision orbit determination and the mapping of geophysical parameters, J. Astronaut. Sci. 38 (1990) 355-368.

[17] D.D. Rowlands, D.E. Pavlis, F.G. Lemoine, G.A. Neumann, S.B. Luthcke, The use of laser altimetry in the orbit and attitude determination of Mars Global Surveyor, Geophys. Res. Lett. 26 (9) (1999) 1191-1194, http://dx.doi.org/10.1029/1999gl090223.

[18] G.A. Neumann, D.D. Rowlands, F.G. Lemoine, D.E. Smith, M.T. Zuber, Crossover analysis of Mars Orbiter Laser Altimeter data, J. Geophys. Res.: Planets 106 (E10) (2001) 23753-23768, http://dx.doi.org/10.1029/2000je001381.

[19] F.G. Lemoine, D.E. Smith, D.D. Rowlands, M. Zuber, G. Neumann, D. Chinn, D. Pavlis, An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor, J. Geophys. Res.: Planets 106 (E10) (2001) 23359-23376, http://dx.doi.org/10.1029/2000je001426.

[20] D.D. Rowlands, F.G. Lemoine, D.S. Chinn, S.B. Luthcke, A simulation study of multi-beam altimetry for Lunar Reconnaissance Orbiter and other planetary missions, J. Geod. 83 (8) (2009) 709-721, http://dx.doi.org/10.1007/s00190-008-0285-y.

[21] E. Mazarico, D.D. Rowlands, G.A. Neumann, D.E. Smith, M.H. Torrence, F.G. Lemoine, M.T. Zuber, Orbit determination of the Lunar Reconnaissance Orbiter, J. Geod. 86 (3) (2012) 193-207, http://dx.doi.org/10.1007/s00190-011-0509-4.

[22] E. Mazarico, M.K. Barker, G.A. Neumann, M.T. Zuber, D.E. Smith, Detection of the lunar body tide by the Lunar Orbiter Laser Altimeter, Geophys. Res. Lett. 41 (7) (2014) 2282-2288, http://dx.doi.org/10.1002/2013gl059085.

[23] S. Goossens, K. Matsumoto, D.D. Rowlands, F.G. Lemoine, H. Noda, H. Araki, Orbit determination of the SELENE satellites using multi-satellite data types and evaluation of SELENE gravity field models, J. Geod. 85 (8) (2011) 487-504, http://dx.doi.org/10.1007/s00190-011-0446-2.

[24] E. Mazarico, G.A. Neumann, M.K. Barker, S. Goossens, D.E. Smith, M.T. Zuber, Orbit determination of the Lunar Reconnaissance Orbiter: Status after seven years, Planet. Space Sci. 162 (2018) 2-19, http://dx.doi.org/10.1016/j.pss.2017.10.004.

[25] M. Marabucci, The Determination of the Hermean Gravity and Topography from Radio Science and Laser Altimeter Data of the Mission BepiColombo (Ph.D. thesis), Sapienza Università di Roma, Facoltà di Ingegneria, 2013.

[26] A. HosseiniArani, BepiColombo Laser Altimeter (BELA) Performance Evaluation: from Laboratory Tests to Simulations of Flight Observations (Ph.D. thesis), Universität Bern, 2020.

[27] S. Villamil, D. Dirkx, A. Stark, H. Hussmann, Improvement of orbit determination using laser altimeter crossovers: JUICE mission case study, Acta Astronaut. 182 (2021) 587-598, http://dx.doi.org/10.1016/j.actaastro.2021.02.005.

[28] A. Genova, D.E. Smith, R. Canup, T. Hurford, S. Goossens, E. Mazarico, G.A. Neumann, M.T. Zuber, F. Nimmo, M. Wieczorek, E. Bierhaus, Geodetic investigations of the mission concept MAGIC to reveal Callisto's internal structure, Acta Astronaut. 195 (2022) 68-76, http://dx.doi.org/10.1016/j.actaastro.2022.02.013.

[29] W. Desprats, S. Bertone, D. Arnold, A. Jäggi, M. Blanc, Influence of low orbit design and strategies for gravity field recovery of Europa, in: D.O. Witasse (Ed.), Planet. Space Sci. 226 (2023) 105631, http://dx.doi.org/10.1016/j.pss.2022.105631.

[30] R. Dach, S. Lutz, P. Walser, P. Fridez, Bernese GNSS Software Version 5.2, University of Bern, Bern Open Publishing, 2015, http://dx.doi.org/10.7892/boris.72297.

[31] A. Jäggi, G. Beutler, U. Meyer, H. Bock, L. Mervart, The role of position information for the analysis of K-band data: Experiences from GRACE and GOCE for GRAIL gravity field recovery, in: N. Sneeuw, P. Novák, M. Crespi, F. Sansò (Eds.), VIII Hotine-Marussi Symposium on Mathematical Geodesy, Springer International Publishing, Cham, 2015, pp. 157-163, http://dx.doi.org/10.1007/1345_2015_63.

[32] D. Arnold, S. Bertone, A. Jäggi, G. Beutler, L. Mervart, GRAIL gravity field determination using the Celestial Mechanics Approach, Icarus 261 (2015) 182-192, http://dx.doi.org/10.1016/j.icarus.2015.08.015.

[33] S. Bertone, D. Arnold, V. Girardin, M. Lasser, U. Meyer, A. Jäggi, Assessing reduced-dynamic parametrizations for GRAIL orbit determination and the recovery of independent lunar gravity field solutions, Earth Space Sci. 8 (6) (2021) http://dx.doi.org/10.1029/2020ea001454.

[34] S. Bertone, W. Desprats, M.K. Barker, E. Mazarico, PyXover - A Python suite of altimetry analysis tools for planetary geodesy v0.2.0, 2024, http://dx.doi.org/10.5281/ZENODO.13963746.

[35] R.A. Jacobson, JUP310, 2013, URL https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/a_old_versions/jup310.bsp.

[36] W.M. Folkner, J.G. Williams, D.H. Boggs, R.S. Park, P. Kuchynka, The planetary and lunar ephemerides DE430 and DE431, Interplanet. Netw. Prog. Rep. 196 (1) (2014) 42-196.

[37] B.A. Archinal, C.H. Acton, M.F. A'Hearn, A. Conrad, G.J. Consolmagno, T. Duxbury, D. Hestroffer, J.L. Hilton, R.L. Kirk, S.A. Klioner, D. McCarthy, K. Meech, J. Oberst, J. Ping, P.K. Seidelmann, D.J. Tholen, P.C. Thomas, I.P. Williams, Report of the IAU working group on cartographic coordinates and rotational elements: 2015, Celest. Mech. Dyn. Astron. 130 (3) (2018) http://dx.doi.org/10.1007/s10569-017-9805-5.

[38] W.M. Kaula, Theory of Satellite Geodesy, Applications of Satellites to Geodesy, Blaisdell Publishing Company, Waltham, MA., USA, 1966, http://dx.doi.org/10.1063/1.3033941.

[39] F.G. Lemoine, S. Goossens, T.J. Sabaka, J.B. Nicholas, E. Mazarico, D.D. Rowlands, B.D. Loomis, D.S. Chinn, G.A. Neumann, D.E. Smith, M.T. Zuber, GRGM900C: A degree 900 lunar gravity model from GRAIL primary and extended mission data, Geophys. Res. Lett. 41 (10) (2014) 3382-3389, http://dx.doi.org/10.1002/2014GL060027.

[40] G. Petit, B. Luzum, IERS Conventions (2010), Technical Report, International Earth Rotation and Reference Systems Service, 2010.

[41] V. Lainey, L. Duriez, A. Vienne, Synthetic representation of the Galilean satellites' orbital motions from L1 ephemerides, Astron. Astrophys. 456 (2) (2006) 783-788, http://dx.doi.org/10.1051/0004-6361:20064941.

[42] E. Mazarico, D. Buccino, J. Castillo-Rogez, A.J. Dombard, A. Genova, H. Hussmann, W.S. Kiefer, J.I. Lunine, W.B. McKinnon, F. Nimmo, et al., The Europa Clipper gravity and radio science investigation, Space Sci. Rev. 219 (4) (2023) 30, http://dx.doi.org/10.1007/s11214-023-00972-0.

[43] G. He, M. Liu, X. Gao, X. Du, H. Zhou, H. Zhu, Chinese deep space stations: a brief review [antenna applications corner], IEEE Antennas Propag. Mag. 64 (1) (2022) 102-111.

[44] T.D. Moyer, Formulation for observed and computed values of Deep Space Network data types for navigation, 2003, http://dx.doi.org/10.1002/0471728470.
[45] I.I. Shapiro, Fourth test of general relativity, Phys. Rev. Lett. 13 (26) (1964) 789, http://dx.doi.org/10.1103/PhysRevLett.13.789.

[46] P.H. Phipps, Distribution of Plasma in the Io Plasma Torus from Radio Occultations During the Juno Epoch (Ph.D. thesis), Boston University, 2019.

[47] L.G. Casajus, A.I. Ermakov, M. Zannoni, J.T. Keane, D. Stevenson, D.R. Buccino, D. Durante, M. Parisi, R.S. Park, P. Tortora, S.J. Bolton, Gravity field of Ganymede after the Juno extended mission, Geophys. Res. Lett. 49 (24) (2022) http://dx.doi.org/10.1029/2022gl099475.

[48] S.W. Asmar, J.W. Armstrong, L. Iess, P. Tortora, Spacecraft Doppler tracking: Noise budget and accuracy achievable in precision radio science observations, Radio Sci. 40 (2) (2005) n/a-n/a, http://dx.doi.org/10.1029/2004rs003101.

[49] L. Iess, M. Di Benedetto, N. James, M. Mercolino, L. Simone, P. Tortora, Astra: Interdisciplinary study on enhancement of the end-to-end accuracy for spacecraft tracking techniques, Acta Astronaut. 94 (2) (2014) 699-707, http://dx.doi.org/10.1016/j.actaastro.2013.06.011.

[50] J.A. Barnes, A.R. Chi, L.S. Cutler, D.J. Healey, D.B. Leeson, T.E. McGunigal, J.A. Mullen, W.L. Smith, R.L. Sydnor, R.F.C. Vessot, G.M.R. Winkler, Characterization of frequency stability, IEEE Trans. Instrum. Meas. IM-20 (1971) 105-120, http://dx.doi.org/10.1109/tim.1971.5570702.

[51] L. Iess, M. Di Benedetto, M. Marabucci, P. Racioppa, et al., Improved Doppler tracking systems for deep space navigation, in: International Symposium on Space Flight Dynamics, ISSFD, 2012.

[52] S. Bertone, E. Mazarico, M.K. Barker, S. Goossens, T.J. Sabaka, G.A. Neumann, D.E. Smith, Deriving Mercury geodetic parameters with altimetric crossovers from the Mercury Laser Altimeter (MLA), J. Geophys. Res.: Planets 126 (4) (2021) http://dx.doi.org/10.1029/2020je006683.

[53] T. Van Hoolst, C. Jacobs, Mercury's tides and interior structure, J. Geophys. Res.: Planets 108 (E11) (2003) http://dx.doi.org/10.1029/2003JE002126.

[54] G. Steinbrügge, A. Stark, H. Hussmann, K. Wickhusen, J. Oberst, The performance of the BepiColombo Laser Altimeter (BELA) prior launch and prospects for Mercury orbit operations, Planet. Space Sci. 159 (2018) 84-92, http://dx.doi.org/10.1016/j.pss.2018.04.017.

[55] S. Goossens, E. Mazarico, Y. Ishihara, B. Archinal, L. Gaddis, Improving the geometry of Kaguya extended mission data through refined orbit determination using laser altimetry, Icarus 336 (2020) 113454, http://dx.doi.org/10.1016/j.icarus.2019.113454.

[56] K.J. Becker, M.S. Robinson, T.L. Becker, L.A. Weller, K.L. Edmundson, G.A. Neumann, M.E. Perry, S.C. Solomon, First global digital elevation model of Mercury, in: 47th Annual Lunar and Planetary Science Conference, Vol. 1903, 2016, p. 2959.

[57] N. Thomas, H. Hussmann, T. Spohn, L. Lara, U. Christensen, M. Affolter, T. Bandy, T. Beck, S. Chakraborty, U. Geissbühler, et al., The BepiColombo laser altimeter, Space Sci. Rev. 217 (2021) 1-62, http://dx.doi.org/10.1007/s11214-021-00794-y.

[58] G. Beutler, A. Jäggi, L. Mervart, U. Meyer, The Celestial Mechanics Approach: theoretical foundations, J. Geod. 84 (10) (2010) 605-624, http://dx.doi.org/10.1007/s00190-010-0401-7.

[59] G. Beutler, Methods of celestial mechanics. Vol. I: Physical, Mathematical, and Numerical Principles, Springer Berlin, Heidelberg, 2005, http://dx.doi.org/10.1007/b138225.

[60] J. Kusche, Noise variance estimation and optimal weight determination for GOCE gravity recovery, Adv. Geosci. 1 (2003) 81-85, http://dx.doi.org/10.5194/adgeo-1-81-2003, URL https://adgeo.copernicus.org/articles/1/81/2003/.

[61] L.J. Wood, Orbit determination singularities in the Doppler tracking of a planetary orbiter, J. Guid. Control Dyn. 9 (4) (1986) 485-494, http://dx.doi.org/10.2514/3.20136.

[62] R. Russell, S. Thurman, An analytic development of orbit determination for a distant, planetary orbiter, Telecommun. Data Acquis. Rep. (1989).

[63] C. Bonanno, A. Milani, Symmetries and rank deficiency in the orbit determination around another planet, Celest. Mech. Dyn. Astron. 83 (1/4) (2002) 17-33, http://dx.doi.org/10.1023/a:1020110226569.

[64] M. Lasser, Noise Modelling for GRACE Follow-On Observables in the Celestial Mechanics Approach (Ph.D. thesis), Universität Bern, 2023, http://dx.doi.org/10.48549/4127.

[65] G. Steinbrügge, D. Schroeder, M. Haynes, H. Hussmann, C. Grima, D. Blanken-ship, Assessing the potential for measuring Europa's tidal Love number h2 using radar sounder and topographic imager data, Earth Planet. Sci. Lett. 482 (2018) 334-341, http://dx.doi.org/10.1016/j.epsl.2017.11.028.