## 2 Functional and Stochastic Models for the Ranging Observations  
For the purpose of this study, let us consider the following variables. The three-element vector \(\mathbf{r}_{i}(t)=(x_{i}(t),y_{i}(t),z_{i}(t))^{\mathrm{T}}\) denotes the 3D coordinates identifying the position of object \(i\) (satellite's center of mass, or ground station reference point) at time \(t\) in the Earth-Centered Earth-Fixed (ECEF) reference frame. The ranging observable, denoted with \(\rho_{ij}(t)\), is the measured range between objects \(i\) and \(j\) at time \(t\),  
\[\rho_{ij}(t)=\|\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t-\Delta_{ij})\|+\varepsilon_{ ij}(t),\quad\varepsilon_{ij}(t)\sim\mathcal{N}(b_{ij}(t),\sigma_{ij}^{2}), \tag{1}\]  
where \(\Delta_{ij}\) is the signal flight time and \(\varepsilon_{ij}(t)\) is the measurement error (\(\sigma_{ij}^{2}\) denotes the measurement variance, and \(b_{ij}(t)\) is the measurement bias, generally time dependent). The latter is associated to a number of unmodeled effects: range biases due to residual clock offsets (we are here assuming one-way ranging based on the assumption of synchronized clocks at the two objects \(i\) and \(j\)), uncorrected relativistic effects, hardware delay calibration errors, satellite attitude modeling errors, as well as ground station coordinates uncertainties in GSR measurements.Functional Model (1) links the ranging observables with the unknown satellite coordinates. As a purely kinematic solution relying only on ISR and sparse GSR measurements is not feasible, satellite orbit is parametrized by a set of 15 parameters of a Galileo standard F/NAV navigation message [6]. We then denote the vector of unknown orbital parameters \(\mathbf{p}\) at time \(t_{oe}\) (time of ephemeris) for unknown \(m\) satellites as  
\[\mathbf{p}(t_{oe})\equiv\left[\mathbf{p}_{1}^{\mathrm{T}}(t_{oe}),\mathbf{p} _{2}^{\mathrm{T}}(t_{oe}),\ldots,\mathbf{p}_{m}^{\mathrm{T}}(t_{oe})\right]^{ \mathrm{T}}, \tag{2}\]  
where \(\mathbf{p}_{s}(t_{oe})\) is a 15-element vector containing the perturbed Keplerian parameters given in Table 1 (a), referred to satellite \(s\) at time of ephemeris \(t_{oe}\).  
Satellite Cartesian ECEF coordinates \(\mathbf{r}(t)\) are computed from the 15 F/NAV orbital parameters following the user algorithm (non-linear functions) summarized in Table 1 (b) [6]. The advantage of using an orbit model based on the 15 perturbed Keplerian parameterslies in its simplicity, which removes the need for a complex modeling of all forces acting on the satellite and their integration over time. Moreover, the data relay capabilities of the ISLs enable distributing ISR and GSR measurements across the whole constellation, potentially enabling autonomous orbit determination onboard each satellite. In this approach, the performance of orbit determination is determined by the accuracy of (a) the ranging observables, (b) the quality of the parametric model used to represent the actual satellite orbit and (c) on the length of the fitting interval (the longer the interval, the larger the errors, as it is shown in Section 5.1). It should be noted that such parametrized orbit has a very limited validity when extrapolating outside the fitting interval.