Assessing Reduced‐Dynamic Parametrizations for GRAIL Orbit Determination and the Recovery of Independent Lunar Gravity Field Solutions

Orbit determination of probes visiting Solar System bodies is currently the main source of our knowledge about their internal structure, inferred from the estimate of their gravity field and rotational state. Nongravitational forces acting on the spacecraft need to be accurately included in the dynamical modeling (either explicitly or in the form of empirical parameters) not to degrade the solution and its geophysical interpretation. In this study, we present our recovery of NASA GRAIL orbits and our lunar gravity field solutions up to degree and order 350. We propose a systematic approach to select an optimal parametrization with empirical accelerations and pseudo‐stochastic pulses, by checking their impact against orbit overlaps or, in the case of GRAIL, the very precise inter‐satellite link. We discuss how parametrization choices may differ depending on whether the goal is limited to orbit reconstruction or if it also includes the solution of gravity field coefficients. We validate our setup for planetary geodesy by iterating extended lunar gravity field solutions from pre‐GRAIL gravity fields, and we discuss the impact of empirical parametrization on the interpretation of gravity solutions and of their error bars.

about the spacecraft properties or attitude modeling are available. The choice of these parametrizations often results from the past experience of the scientist and it constitutes one of the main differences among the solutions provided by different groups. However, their impact on the solution of both orbit and geophysical parameters is usually not analyzed in a systematic way.
The NASA GRAIL mission  is the first planetary mission equipped with a sub-μm/s precision inter-satellite link, operating in Ka-band to measure variations in the relative distance change of the two probes, GRAIL-A and GRAIL-B (Klipstein et al., 2013). Orbit information is then available over the whole orbit of the probes around the Moon, thus allowing for a highly accurate recovery of the lunar gravity field on both sides of the Moon (see, e.g., Goossens et al., 2020;Konopliv et al., 2014;Lemoine et al., 2014). Moreover, the Ka-Band Range Rate (KBRR) link also allows for an independent check of the quality of the recovered orbits for both probes.
In this study, we present the planetary extension within the research version of the Bernese GNSS Software (BSW, Dach et al., 2015). We use it to provide independent solutions for the orbits of GRAIL-A and GRAIL-B and for the gravity field of the Moon up to degree and order (d/o) 350 in its spherical harmonic expansion based on data from the GRAIL Primary Mission (PM, March 1, 2012to April 19, 2012. We present a thorough analysis of our choice of dynamical and empirical orbit modeling, and we study the impact of this so-called "reduced-dynamics" approach (Jäggi et al., 2011) on both orbit and gravity field recovery. The former can be used as guideline for the processing of Doppler data, while the latter focuses on the sensitivity of gravity field coefficients to the parametrization of the solution, not to be overlooked when inferring geophysical properties of the body of interest.
Also, when appropriate, we provide details about Doppler processing in the BSW and specifics on a limited number of critical steps in the processing of deep-space Doppler data which have proved difficult to find in the literature.

GRAIL Orbit Determination in the Bernese (GNSS) Software
Orbit determination (OD) for the two GRAIL satellites is mainly based on Doppler observations by the Deep Space Network (DSN) for its absolute positioning and on KBRR provided by the Lunar Gravity Ranging System (LGRS, Klipstein et al., 2014) for the high accuracy measurement of the inter-satellite relative velocity. Moreover, in the specific case of a Moon orbiter, the latter is assuring a continuous tracking even when the probes are not visible from Earth (i.e., when flying over the far side of the Moon).
The research version of the BSW can readily deal with KBRR data as an heritage of previous GRACE (e.g., Meyer et al., 2016, 03) and GRAIL (e.g., Arnold et al., 2015) activities at the Astronomical Institute of the University of Bern (AIUB). We extended the framework of BSW processing to include nongravitational force models, reference frames and an accurate modeling of Doppler observables (following the guidelines of Moyer, 2003) to comply with standards of OD software currently used in planetary OD and geodesy (e.g, Evans et al., 2018;GRGS, 2018;Pavlis et al., 2013;Serra et al., 2019).

Doppler observations O
 from Orbit Determination Files (ODF) are imported into our internal format and eventually converted to the desired integration time T C , that is, the interval over which the phase of the tracking signal is accumulated. Orbit integration from a priori initial elements, a background force model and an accurate modeling of light propagation are used to compute Doppler observables C  , as detailed in Section 2.1. These are in turn used to set up the observation equations required for the least squares recovery of arc-specific orbital and global (e.g., gravity field coefficients) parameters.
For the computation of GRAIL orbits, we follow a procedure similar to Jäggi et al. (2008) (application to GRACE) and Arnold et al. (2015) (GRAIL solution based on dynamic positions by the GRAIL science team), by performing the following steps: • Set up of two Normal Equation (NEQ) systems based on GRAIL-A and GRAIL-B Doppler observations as part of a reduced-dynamic OD procedure. Per orbit arc we estimate six initial osculating elements, empirical accelerations (in each of the three directions), as well as pseudo-stochastic pulses in the same direction as the accelerations. This limited number of empirical orbit parameters is meant to handle model deficiencies. Pulses falling in intervals with no Doppler coverage are tightly constrained to 0 in this phase to ensure robustness. This applies to both one-way and two-way Doppler observations. • Set up of one NEQ based on KBRR observations using the same orbit parametrization. In addition to the parameters specified in the first step, specific Ka-band parameters may be set up (e.g., a time bias). Due to the presence of the orbit parameters of both satellites, the normal equations set up in this step are singular when used alone. • Doppler and KBRR NEQs are combined by setting an appropriate weighting ratio, for example, based on the relative accuracy of the observations or determined by Variance Component Estimation (VCE, Kusche, 2003). • Solution of the least squares problem by numerically efficient algorithms, for example, LAPACK (Anderson et al., 1999) DOPTRF function or, for smaller problems, SYMINV2 (Rutishauser, 1963).
The orbits emerging from this procedure then serve as a priori orbits for the setup of the generalized OD and the solution of NEQs containing both orbit and gravity field parameters.

A Priori Orbit and Forces Modeling
We use the positions and velocities provided by the GRAIL science team as initial conditions for each daily arc and perform orbit integration for each GRAIL probe.
As detailed in Arnold et al. (2015), the Celestial Mechanics Approach (CMA, Beutler et al., 2010) uses the equations of motion in the inertial system expressed as where GM M denotes the gravity parameter of the Moon,  r is the selenocentric position of the probe and  f collects all perturbing accelerations. Dots indicate derivatives with respect to time. The second-order differential Equations 1 require six initial or boundary conditions for a particular solution. In the framework of the CMA the satellite's motion is described as an initial value problem. The initial conditions  0 ( ) r t and   0 ( ) r t at an initial epoch t 0 are defined by six Keplerian osculating elements. q 1 , …, q d , parametrizing the perturbing accelerations, are either arc-specific orbit parameters like, for example, scaling factors of nongravitational or empirical accelerations, or general parameters such as gravity field coefficients. Let us denote the 6 + d parameters (initial conditions and q i ) collectively as p i .
The perturbing forces appearing on the right-hand side of Equation 1 are given by where V denotes the lunar gravity potential in the body-fixed reference frame, i f T is a rotation matrix relating the body-fixed with the inertial system,  b a are the third-body accelerations,  t a denote accelerations due to the tidal deformation of the Moon,  r a are relativistic corrections,  n a nongravitational accelerations and  e a empirical accelerations. Details about integration settings and force models are provided in Table 1.
The initial orbital elements and, possibly, dynamical and stochastic parameters are then fitted to the Doppler data by using a classical least-square procedure for each arc. The full model of Doppler data  ( ( ), ) r t f  is linearized around values based on the a priori orbit  0 ( ) f are the unknowns of a least-squares adjustment for all orbit parameters and Doppler-specific parameters (e.g., clock parameters), respectively.

Table 1 Setup and Models Used for GRAIL Orbit Determination
We model the impact of these forces on the GRAIL probes represented by the 28-plates macromodel given in Fahnestock et al. (2012) with optical properties as detailed in Wirnsberger et al. (2019).
Remaining mismodelings are accounted for by empirical accelerations and pseudo-stochastic pulses (Jäggi et al., 2006) estimated in the orbital frame (radial, R, along-track, S, and cross-track, W) directions. For the scope of this study, it is sufficient to mention that we use empirical accelerations as constant or 1 cycle-per-revolution (1-cpr) sine and cosine terms, with the same amplitude estimated for the whole arc duration. Concerning pseudo-stochastic pulses (instantaneous velocity changes imposed on the probe at fixed time intervals) a different amplitude is estimated for each pulse. A difference between pulses and empirical accelerations is then that between one pulse and the next the orbit is fully dynamical, only determined by explicitly modeled forces. The estimated amplitude of pulses is constrained to zero by user-determined statistical bounds (which in the following we characterize as "loose" or "tight" with respect to the a priori measurement error); pulses falling in time slots with few or no observations are tightly constrained to zero to make the system more robust.

Tracking Data: Retrieval and Pre-Processing
Both one-way X-band Doppler (1W) and two-way S-band Doppler (2W) tracking measurements of GRAIL, as well as KBRR data, are available on NASA's Planetary Data System (PDS) Geosciences Node. We convert Doppler data from the binary Orbit Determination File (ODF) format TRK-2-18 (Shin, 2008) to ASCII and import them to daily tracking files for each observation type (and each satellite), and to a set of daily ramp table files for each station. We then accumulate Doppler data over 10 s integration time (from the original 1 s) to reduce noise and computational burden. In this study, we only use 2W Doppler data along with inter-satellite KBRR, although similar results can be achieved by using 1W Doppler data by including additional parameters related to the frequency of the on-board clock (Moyer, 2003). Table 2 shows statistics about the distribution of tracking data along the PM phase, eventually divided by tracking station and site.
Based on GRAIL orbits integrated with the force model presented in sections 2.1 and 2.2, we compute light-times between the tracking station and the probes, including all required time scale transformations and relativistic corrections, that is, Shapiro time delays. Tracking station coordinates at January 1, 2012 are extrapolated using coordinates and velocities from Folkner et al. (2014). Precise light times expressed in TDB (Barycentric Dynamical Time, Bretagnon & Brumberg, 2003) then lead to "computed" Doppler frequencies by introducing a reference frequency REF  , as detailed in Moyer (2003), Equations 13-54. For documentation purposes, Hz, with high-part (Hi) and low-part (Lo) stored in columns 21-22 of the ODF file. Computing the light-time between tracking-station and probe is a sensitive step regarding numerical precision: in particular, using 8-bit floats for the reception and emission epochs results in numerical noise exceeding 10 mHz for a Moon probe. Using quadruple precision floats, or splitting each epoch in seconds and fraction of seconds within the light-time computation algorithm, allows to keep numerical precision compatible with sub-mHz level computed Doppler. The difference of model and observations provides Doppler residuals, which we use as basis for screening the observations. A mix of threshold-based screening (for large outliers) and the manual detection of time spans to be excluded from the computations are used. All identified observations, and those below 15° elevation (where atmospheric noise has a stronger impact), are then flagged for exclusion. Table 2 shows the statistics of excluded observations for each tracking station, satellite, and observation type in a typical run. A similar procedure is used for KBRR data, plus data at specific locations, for example, those affected by systematic mismodeling of SRP, are flagged by ad-hoc tools. Exclusion flags are eventually updated at each iterations, based on newly estimated orbits and/or force model. In the end, our analysis is supported by BERTONE ET AL. Note. Around 5.5% of measurements are excluded in an average GRAIL run (the exact number can change depending on the setup and a priori modeling). Canberra DSN site (DSS-34 and DSS-45) performs slightly worse than average. ∼1.5 × 10 6 KBRR measurements collected almost continuously (with just a few "technical" gaps, e.g., due to maneuvers) over the PM with a sampling of 5 s.

Orbit Results and Quality Assessment
Our approach to OD is based on a mix of direct modeling and empirical parametrization using both accelerations and stochastic terms, as explained in Section 2. Although possible within the BSW, we do not apply any scaling factor to our models of nongravitational forces. Empirical parametrizations are regularly applied for Earth satellites, where dense GNSS data support their determination (Jäggi et al., 2006). For planetary probes, however, one mostly relies on Doppler data providing a much weaker constraint. We use both orbit overlaps and the very accurate measurements of the inter-satellite KBRR link to evaluate orbits based on an extensive set of empirical parametrizations. We show that the use of empirical parameters can partly replace a detailed modeling of nongravitational forces for planetary geodesy, but also that parametrizations should be carefully chosen, especially when only constrained by sparse Doppler observations. It goes without saying that improved Doppler residuals do not necessarily correspond to improved orbits. Throughout this section, when not differently specified, we use the high-quality GRAIL gravity field GRG-M900C, truncated at d/o 350 (for computation time reasons), as background field, implying that we are not confronted with large background model errors on the far-side of the Moon.

Doppler-Based Solutions
Two-way S-band Doppler observations are the main source of information for GRAIL absolute positioning. We first integrate GRAIL-A and GRAIL-B orbits with a force model including only gravitational forces acting on the satellites, i.e., the extended gravity potential of the Moon, J 2 of Earth and third-body forces of point masses by other Solar System planets and the Sun. We fit to Doppler data a set of initial Keplerian elements for each daily arc, and, optionally, a set of constant and 1-cpr empirical accelerations in R, S, and W directions. Based on these orbits, we compute range-rates between the two probes and compare them with KBRR observations. The Root Mean Square Error (RMSE) of the KBRR residuals over the whole PM constitutes our evaluation criteria for a set of parametrizations, including empirical accelerations and stochastic parameters. Figure 1 shows the result of such tests over the whole parameter space when only explicitly modeling gravitational forces, left, and when also modeling non-gravitational forces, right. Color bar values are centered around the RMSE of the purely dynamical case, so that green cells show an improvement due to empirical and/or stochastic parameters.
BERTONE ET AL.
10.1029/2020EA001454 6 of 21 One notices that in both cases a full parametrization with constant and 1-cpr accelerations in all directions would be detrimental when using Doppler data. An explicit modeling of nongravitational forces reduces the overall residuals by ∼60% even in a purely dynamical processing, as expected in the GRAIL case, when an accurate macro-model and attitude information are available. Still, estimating a few empirical accelerations on top of modeling nongravitational forces allows to further improve the results, but only marginally. Moreover, we notice that empirical accelerations in R and S directions have a similar impact on the orbit. However, when aiming at gravity recovery, we would generally avoid estimating constant radial accelerations since they are strongly correlated with accelerations due to the GM of the central body.
Doppler observations from Earth take place between two "extreme" geometries for a lunar probe on a polar orbit, each repeating every two weeks over the three months of the PM phase: • The "face-on" geometry, when the orbital plane is perpendicular to the line of sight. In this particular configuration, the full orbit of the probe around the Moon is visible from Earth. However, the line-ofsight velocity of the probe with respect to the observer (and hence to the Doppler frequency shift) is close to 0, which limits the quality of OD; • The "edge-on" geometry, when the orbital plane is parallel to the line of sight. A large portion of the orbit is then not visible from Earth (up to 80%, in our analysis). The along-track component of the orbit is well constrained by observations, while the cross-track component is difficult to determine, that is, its inclination i and the right ascension of its ascending node Ω.
In Figure 2, we compare the value of KBRR residuals in these scenarios. Orbits are best determined during edge-on days, showing that geometry of Doppler observations plays a stronger role than sparser data coverage along the orbit. We also find that most parametrizations with empirical accelerations do not improve the determination during face-on arcs (i.e., most top half cells are dark red, indicating a larger RMSE than the reference solution without empirical terms). This applies to both cases, with and without the modeling of nongravitational forces, showing that observation geometry also strongly limits the quality of estimated empirical parameters.
As mentioned in Section 2.2, we base our OD setup on both empirical accelerations, whose amplitude is estimated only once over each arc (i.e., resulting in a maximum of nine additional parameters per arc), and on a larger set of pseudo-stochastic pulses with a fixed spacing. Based on the results shown in Figure 1 relatively loose (lower than the a priori σ associated to observations) to tightly constrained (ratio of the constraint to 0 and the a priori σ associated to observations smaller than one). Because of computational time, we do not explicitly model nongravitational forces for this experiment.
A loose constraining would generally degrade the results. Figure 3 compares the chosen set of orbit solutions over the whole PM. It shows that, when only Doppler tracking is available, adding stochastic pulses on top of empirical accelerations only marginally improves the orbits. We also show (on the right) that face-on orbits are generally the worst reconstructed, independently from the parametrization. Nevertheless, we notice that the estimation of pulses in W is more stable during these days, as expected.
To summarize, we presented a systematic approach to the choice of the optimal empirical orbit parametrization for Doppler-based recovery. Our results highlighted that different choices might be best suited for different orbital configurations and a priori knowledge of, for example, nongravitational forces. Still, for practical reason, a globally best-performing parametrization can be chosen, as we do here. The choice of an adapted parametrization using empirical terms significantly improves the quality of Doppler-based orbits, especially when no modeling of nongravitational forces is available. Different from the GNSS-based orbit recovery of Earth-satellites, the estimate of local pseudo-stochastic parameters only marginally improves our Doppler-based solutions (and actually degrades it in most cases). When independent measurements of the orbit quality, e.g., range or Delta-Differential One-Way Ranging (Δ-DOR), are not available, orbit overlaps can be used for a systematic evaluation, as shown in Section 3.2.

Impact of KBRR Data
In the case of GRAIL, orbits are further improved based on accurate and continuous inter-satellite KBRR measurements. We extend the analysis presented in Section 3.1 to the combined processing of Doppler and KBRR data. For these preliminary analysis, the two datasets are weighted with a fixed ratio of 1:10 4 in favor BERTONE ET AL.
As KBRR data enter the OD process in a combined least squares fit with Doppler data, we cannot use KBRR residuals to independently evaluate the orbit quality (more estimated parameters will automatically result in lower residuals). We then recur to the analysis of orbit overlaps: GRAIL-A/B orbits are propagated for 25h from the daily set of estimated arc-parameters, resulting in overlapping orbit position batches of 1h length. When applied to the Doppler-only scenario, the orbit overlaps analysis led to results consistent with those shown in Figures 1-3. Figure 4 then shows the RMSE of overlaps differences over the full PM for the same sets of modeling and parametrizations presented in Section 3.1. As expected, an improved modeling of nongravitational forces results in improved orbit overlaps for the purely dynamical case. We confirm the results of Section 3.1 and identify an optimal parametrization using a constant acceleration in the along-track direction and 1-cpr accelerations in either radial or along-track direction. These parameters provide some freedom to absorb eventual mismodelings in the radial direction. We also explore the results of parametrizations based on stochastic pulses. Highly accurate and continuous KBRR data allow for a more robust estimate of stochastic pulses, so that several scenarios using stochastic pulses result in a significant reduction of orbit overlaps residuals. Figure 5 indicates that an almost free estimate of pulses degrades nevertheless the orbits. Weak constraints in all directions result in improved orbits with respect to orbits parametrized with empirical accelerations only (here represented by tight constraints applied in all directions). Also, estimating pulses every 30 min is a good compromise allowing to absorb orbit mismodelings and at the same time limit potential correlations with low-degree coefficients of the gravity field to be recovered.

Gravity Field Determination
The orbits determined in the combined Doppler and KBRR OD serve as a priori information for a common orbit and gravity field estimation based on the accumulation of daily normal equations over the whole PM. We first validate the orbit parametrization identified in Section 3 and presented in Table 3 for the recovery of the global lunar gravity field up to d/o 350 in spherical harmonic expansion. Then, we provide independent solutions for lunar gravity, which will serve as basis for discussing the impact of empirical parametrizations on gravity recovery in Section 5.
A classical least squares adjustment is used by extending Equation 3 to a generalized orbit improvement of the previously computed a priori orbits  0 ( ) j r t (j = 1, 2 denotes the satellites). Therefore, beside arc-and satellite-specific parameters with a priori values p 0ij , with i = 1, …, 6 + d, we also set up corrections to the spherical harmonic coefficients representing the central body gravity field as global "solve-for" parameters. Observation residuals and partial derivatives with respect to these parameters are computed to set-up daily Normal Equations (NEQs), which are progressively stacked over the full PM. The final NEQ, where all arc parameters have been pre-eliminated (see, e.g., Beutler, 2005), is then passed to the solver, providing both updated values and formal errors of global parameters.

Comparing Parametrizations for Gravity Field Solutions
In Section 3, we presented a systematic analysis of the best performing parametrizations for orbit recovery in terms of empirical accelerations and stochastic pulses. In Figure 6, we compare gravity field results for a subset of parametrizations identified in Figure 5 in terms of difference degree amplitudes Δ l , defined as where Δ lm C and Δ lm S are the differences between spherical harmonics coefficients of degree l and order m of a given solution and corresponding coefficients of the reference solution. In the following, the reference solution is GRGM900C .
While stochastic pulses improve the solution in most cases (when compared with a solution only including empirical accelerations-green curve), we get slightly lower differences when tightening pulses constraints in the radial direction. This is expected, because of correlations with gravity field coefficients. In the following, we hence parametrize orbits by a constant acceleration in S, 1-cpr in R and along-track and cross-track pulses every 30 min with a "medium" constraint to 0. Note. The choice of this set of arc parameters is detailed and motivated in Section 3, Figure 4, and by the gravity field solutions shown in Figure 6.

Optimal Weighting of Observables and Mission Phases
Data weighting is an important factor when dealing with a heterogeneous mix of observations, such as in the case of the GRAIL mission. In the first part of this study, we used a fixed weighting of 1:10 4 in favor of KBRR. In reality, as also shown in Lemoine et al. (2013), the consistency of KBRR with the orbits heavily depends on the truncation of the background gravity field solution. Also, the appropriate weighting would change with the altitude of the periapse of the GRAIL orbits over the different phases of the PM (see, e.g. . For our gravity processing, we thus consider three sub-phases within GRAIL PM in 2012, based on changes in periapse altitude of the two satellites: 1. 02 March-10 March (periapse ∼20 km); 2. 11 March-10 May (∼40-50 km); 3. 11 May-29 May (∼20 km).
We use Variance Component Estimation (VCE, see e.g., Kusche, 2003) to determine weights to assign to each observable and phase. Doppler and KBRR NEQs are hence accumulated separately, resulting in six NEQs containing all arc-related and gravity parameters. Because of computational limitations (KBRR is singular with respect to orbit parameters, so that all arc parameters have to be explicitly kept in the NEQ throughout the process), we cannot directly apply this technique to our final setup. We hence compute VCE weights for a reduced setup, with orbits parametrized by only six keplerian elements per arc and by empirical accelerations in along-track (constant) and radial (1-cpr) directions. Figure 7 compares several solution strategies, including VCE-derived weights, in terms of their difference degree amplitudes. For the "VCE, all NEQs" solution, we overweight KBRR NEQs for each phase relatively to Doppler NEQs by a factor 2.8 × 10 2 , 3.6 × 10 2 , and 1.5 × 10 2 , respectively. Differences from the fixed weights used in Section 3.2 are expected, as we base our VCE analysis on a preliminary lunar gravity solution of d/o 300 iterated from SGM150J. While more adapted to our work in Section 4.3, such fields are less consistent with KBRR data than GRGM900C (used as background field in Section 3), thus resulting in lower weights for KBRR (as already noted by, e.g., Lemoine et al., 2013), especially for phases including low altitude data. Estimated VCE weights were compared and validated against the GEODYN II software (Pavlis & Nicholas, 2017) used by Figure 6. Difference degree amplitudes with respect to GRGM900C of lunar gravity field solutions resulting from using several parametrizations for pulses in R, S, and W directions, either loose (l), medium (m), or tight (t)-defined by the ratio of weights given to the constraint with the a priori σ associated with the observations. Most solutions are superposed and differences are mostly visible in the formal errors (dashed lines). Lemoine et al. (2013Lemoine et al. ( , 2014 to compute several reference GRAIL solutions. Letting VCE determine relative weights between both observables and phases ("VCE, all NEQs") results in a clear improvement against fixed weights, especially at low degrees. Another possible approach would be to first combine Doppler and KBRR data in each NEQ with a fixed 1: 10 4 weight, and then estimate relative weights for the three different sub-phases. One would then have the option to "pre-eliminate" arc-parameters before accumulating NEQs over each of the sub-phases, significantly reducing the computational burden. Unfortunately, this approach does not significantly modify "FIX" weights and only allows for a marginal improvement over the baseline scenario, as shown in Figure 7, "VCE, phases only".

Lunar Gravity Field Solution AIUB-GRL350
We analyze two possible paths toward a high-resolution lunar gravity field solution based on our processing. One option is to use a priori information from a GRAIL-based solution, for example, GRGM900C  to set-up the extended OD processing. Using the full GRGM900C, we re-estimate all coefficients up to d/o 350 to assess the "best-case scenario" expected for our resulting field at this resolution. Figure 8 compares difference degree amplitudes of such solution, which we call AIUB-GRL350B, with GRAIL-based solutions up to a similar d/o by other groups, namely GrazGLM420 (Wirnsberger et al., 2019) and GL420 (Zuber, Smith, Lehman, et al., 2013). We get a very good consistency with the reference solution GRG-M900C, especially at low resolutions (up to d/o 75).
The second option is to start the iterative solution process from a pre-GRAIL gravity field, in order to obtain a truly independent solution, which we call AIUB-GRL350A. Figure 9 shows the progression and improvements during our iterative process from the SELENE-based field SGM150J (Goossens et al., 2011). Several iterations were necessary at each resolution, in order to let coefficient values converge before further enlarging the parameter space (e.g., as clearly visible between iterations 1 and 3, both up to d/o 200).
Including days from "phase 3" (11-May to 29-May) right from the start resulted in a degraded solution because of the low (∼20 km) GRAIL periapsis and limited parameter space, insufficient to accommodate the strong gravity signal. We first introduced phase 3 days when expanding the gravity field solution to d/o 300 after "iteration 3". The "iteration 6" solution in Figure 10 shows the significant improvement resulting from BERTONE ET AL.
10.1029/2020EA001454 12 of 21 Figure 7. Difference degree amplitudes with respect to GRGM900C of lunar gravity field solutions computed with: a fixed 1: 10 4 relative weighting between Doppler and KBRR data (FIX, green); with VCE determined weights for each of the three phases (but still fixed between observables, blue), and with six estimated weights (one per each NEQ, purple). Dashed lines represent the formal errors associated to each solution.
including these data when enlarging the parameter space to accommodate the additional signal. Even then, VCE estimated weights for this data were significantly lower than others phases (see Section 4.2). After 9 iterations, we converge on solution AIUB-GRL350A. As shown in Figure 9, AIUB-GRL350A shows a very good consistency with GRGM900C up to d/o 50, where the missing information in terms beyond degree 350, and subsequent omission error, cause increased differences. Further improvements would require an extended parameter space, which is out of our current computational capabilities and of the scope of this work.
BERTONE ET AL.

10.1029/2020EA001454
13 of 21  Figure 9. Difference degree amplitudes of orbit and gravity field improvement iterations from SGM150J. The reference field GRGM900C is also shown for comparison with coefficients differences.
On top of difference degree amplitudes (which assume GRGM900C as a reference solution), we also evaluate the improvement along the iterative process by checking KBRR residuals and by independent metrics such as the correlation of our gravity fields with topography-induced gravity. The daily RMSE of post-fit KBRR residuals for some notable iteration steps is shown in Figure 10. It shows that our gravity field solutions get increasingly consistent with the very accurate KBRR data. Furthermore, the solutions show a decreasing correlation with orbit geometry while the far-side gravity field improves (mainly visible up to iteration 3). Larger KBRR residuals result for the first and last phases of the PM, when GRAIL orbits are lower and the resolution of the estimated field is insufficient. While further iterations to larger fields would likely further flatten these differences, they are out of the scope of this work due to computational limitations.
We compute correlations of our gravity solutions with the Lunar Orbiter Laser Altimeter (LOLA, Smith et al., 2010) topography-induced gravity by following the procedure outlined by Wieczorek and Phillips (1998). Correlations shown in Figure 11 improve at high resolutions when expanding the parameters space over subsequent iterations, which is expected on the Moon since a decrease in correlation with topography is mainly due to the omission error in the gravity field solution. Correlations of AIUB-GRL350B do not drop below 0.93 up to d/o 350, thus almost coinciding with its a priori field GRGM900C (whose coefficients above d/o 350 are kept fixed to their a priori values).
Finally, we analyze differences between the free-air gravity anomalies (Heiskanen & Moritz, 1967) of our two largest solutions, AIUB-GRL350A/B. Figure 12 (top) shows gravity anomalies derived from AI-UB-GRL350B using a Moon reference radius of 1,738 km: most terrain features and details are clearly visible on both sides of the Moon. Differences between the two solutions are below 10% (see Figure 12, bottom) and are localized along a limited number of GRAIL tracks over areas on the near side where the satellites were lowest (see, e.g., Figure 2 of Lemoine et al., 2014), and in cratered areas of the far-side.

How Much Can We Rely on Empirical Parameters?
As detailed information about the shape, attitude and optical properties of planetary probes is not often available, a certain degree of empirical modeling of nongravitational forces is usually a necessity.
In the case of GRAIL, we can use KBRR residuals to evaluate the orbit degradation due to a purely empirical parametrization. Figure 13 shows several hours of KBRR residuals based on our AIUB-GRL350B lunar field and by several explicit or empirical modelings of nongravitational forces. While residuals are on the level of few μm/s for all configurations, the best performing "purely empirical" modeling identified in Section 3 shows large deviations at both the epochs of stochastic pulses and light-shadow transitions. The explicit modeling of radiation pressure and satellite eclipses results in reduced systematic signatures and improved residuals below the μm/s level. Still, Figure 14 shows that high quality lunar gravity field solutions are achievable with an appropriate stochastic parametrization, even when no explicit modeling of nongravitational forces is used. No a priori information on the satellite geometrical and optical properties is required to achieve this result. As a comparison, the best performing setup identified in Section 3, including detailed information of nongravitational forces acting on the TWO satellites, "only" improves the former solution (in terms of consistency with GRGM900C) and its formal errors by around one order of magnitude.

Impact of Empirical Parametrizations on Gravity Field Solutions
At the same time, other geodetic aspects will also be affected by the choice of empirical parametrization. We show in Figure 15 the variability of the low-degree coefficients of lunar gravity field solutions presented in section 4.1 (i.e., only differing in their empirical parametrization) in units of formal errors associated to GRGM900C solution.
Differences between these solutions (including GL420, for reference) are significantly larger than the formal error bars associated with GRAIL coefficients, especially for C21, S21, and degree 3 coefficients. However, we also notice that even when considering the impact of different parametrizations, degree 3 coefficients BERTONE ET AL.

10.1029/2020EA001454
15 of 21 Figure 11. Correlations between the gravity field induced by LOLA topography and several of our lunar gravity field solutions.
from LLR (see, e.g., Viswanathan et al., 2018, values included in Figure 15) are significantly different from GRAIL solutions, as highlighted by recent analysis.
As degree 2 coefficients are used for frame definition purposes, the chosen parametrization also impacts geophysical analysis depending, for example, on the pole position. As an example, we compute different positions of the dynamical lunar pole derived from degree 2 coefficients from the solutions in Figure 15. We use both the approach outlined in Wahr et al. (2015) (and references therein) and the rotation to the Principal Axes frame defined by the diagonalization of the tensor of inertia, getting consistent results. We show in Figure 16 that differences among GRAIL solutions result in different coordinates for the gravitational position of the lunar pole (as a comparison, LLR derived orientation is accurate at ∼10 mas level). The impact of the orbit parametrization should then be adequately taken into account when interpreting gravity field results for geophysical studies, for example, when comparing GRAIL estimates with those of independent techniques such as LLR (Viswanathan et al., 2019). This is best done by evaluating the spread of multiple solutions as different and independent as possible (i.e., based on different software, parametrization, or combination of datasets), rather than based on provided formal errors. Similar procedures are common practice in Earth geodesy, for instance in the analysis of the Earth gravity field and of its temporal variations (Jäggi et al., 2019

Conclusions
At AIUB, OD capabilities from Doppler deep-space tracking have been recently developed in the framework of the Bernese GNSS Software following the guidelines of Moyer (2003) and the most recent conventions for planetary reference frames and ephemerides. This allows to apply the tools and expertise of the BSW group to deep-space and planetary-borne OD and geodetic studies. We present our solutions for the lunar gravity field based on tracking data from the GRAIL mission. Orbiting the Moon in 2012, data from these twin probes significantly increased our knowledge of the Moon's gravity and of its interior (see, e.g., Chappaz et al., 2017;Zuber, 2014). We detail our modeling of two-way S-band Doppler and our approach to get an optimal combination of Doppler and KBRR data.
Orbit improvement often involves the estimate of some sort of empirical parameters, especially for planetary probes whose geometrical and optical properties are often badly known, and for which accelerometer measurements are often not available. In this study, we thus outlined a procedure to choose an appropriate parametrization for orbit and gravity recovery, by a systematic analysis of results on both levels. We used accurate GRAIL KBRR data to evaluate the influence of empirical parameters adapted to the common scenario of Doppler-based OD, using a combination of empirical accelerations and stochastic pulses. We conclude that, although a combination of explicit modeling and empirical parametrization gives the best results, a purely empirical modeling without additional information gives satisfactory results if data coverage is sufficient (especially concerning stochastic pulses). For practical reasons, in this work we chose to apply the same parametrization to all arcs of the GRAIL PM. However, based on Figure 2, we note that future analysis could adopt different empirical parametrizations to each arc, depending on tracking geometry and predominant forces acting on the spacecraft. One can see light-shadow transitions and the impact of pulses (their epochs are marked by the vertical black dotted lines) in the purely empirical solution (red curve). As discontinuities due to eclipses are modeled (green and red curves), the amplitude of stochastic pulses are significantly reduced to the sub-μm/s level. Figure 14. Difference degree amplitudes for two gravity solutions estimated up to d/o 200 (higher degree coefficients are fixed to GRGM900C), with (green) and without (red) the explicit modeling of nongravitational forces (solar and planetary radiation pressure), compared to GRGM900C. Difference degree amplitudes of GL420 are also shown as reference. Degree 2 has been omitted to improve the visualization. As all mission scenarios are different, the setup resulting from our analysis will not necessarily fit another probe, and the proposed algorithm should rather be applied as a pre-processing step when approaching a new dataset.
We then extended our procedure to the combined Doppler-and KBRR-based OD, showing how continuous KBRR help constraining empirical parameters and the overall orbit recovery. We identified a set of optimal parametrizations by evaluating orbit overlaps and tested their impact on gravity field recovery. Using a refined parametrization and VCE-derived weights for each observation over different phases of GRAIL PM, we proved the ability of the BSW to recover highly resolved gravity fields in the planetary environment. We presented d/o 350 solutions based on both GRAIL GRGM900C  field, or fully independent and iterated from the pre-GRAIL SELENE field SGM150J (Goossens et al., 2011). Evaluation in terms of difference-degree amplitudes, post-fit KBRR residuals, correlations with LOLA topography and free-air gravity anomalies prove our solutions to be comparable to those released by the GRAIL science team and by other groups at similar spatial resolutions.
Finally, based on our results and analyses, we discussed the impact of explicit modeling and empirical parameters on gravity results, showing that different parametrizations can result in differences significantly larger than the error bars associated with, in this case, GRAIL solutions. Especially for low-degree coefficients, which are used in the definition of dynamical reference frames and pole coordinates, under-estimating error bars can erroneously influence the geodetic interpretation of released solutions. We thus stress the importance of assessing the uncertainty associated with gravity field solutions from the dispersion of independently estimated coefficients, when available, rather than from the provided (rescaled) formal errors. BERTONE ET AL.

10.1029/2020EA001454
18 of 21 Figure 15. We highlight here the sensitivity of low degree coefficients to the different empirical parametrizations used by most, if not all, planetary geodesy groups. Differences from GRGM900C of a set of gravity field coefficients resulting from 3 different orbit parametrizations using empirical accelerations and pulses are shown in units of formal errors of GRGM900C coefficients. We include GRAIL and LLR solutions by other groups, for reference. Figure 16. Displacement of the lunar pole (arcseconds, degrees) computed from degree 2 gravity coefficients (see Wahr et al., 2015, and references therein) just differing in the empirical parametrization of orbits. By definition, the LLR pole occupies the origin of the plot. Point sizes indicate the dispersion expected from formal errors on gravity coefficients, as comparison.