GoKawiilModelling the time-varying gravity field
The gravitational potential for Mars, U(r, λ, ϕ), can be expressed as a sum over 4π-normalized spherical harmonic coefficients C ℓm and S ℓm :
$$U(r,\lambda ,\phi )=\frac{{GM}}{r}\left[1+\mathop{\sum }\limits_{{\ell }=2}^{{\rm{\infty }}}{\left(\frac{{R}_{m}}{r}\right)}^{{\ell }}\mathop{\sum }\limits_{m=0}^{{\ell }}({C}_{{\ell }m}\cos m\lambda +{S}_{{\ell }m}\sin m\lambda ){P}_{{\ell }m}(\sin \phi )\right],$$ (1)
in which G is the gravitational constant, M is the mass of Mars, R m is a reference radius, P ℓm are fully normalized associated Legendre functions, r is the radial distance from the centre of the planet, λ is the planetocentric longitude and ϕ is the planetocentric latitude. The gravity field is expressed in the body-fixed frame of Mars (as defined in ref. 26), centred at the planet’s COM, so that the degree-1 coefficients vanish by definition. To describe the temporal variability of the Martian gravity field, periodic terms are introduced into the normalized coefficients. This formulation, implemented within the GEODYN II orbit determination framework55, allows the harmonic coefficients C ℓm and S ℓm to vary in time to account for periodic mass redistribution. This time-dependent expansion is written as:
$$\begin{array}{c}{C}_{{\ell }m}(t)={\bar{C}}_{{\ell }m}+\mathop{\sum }\limits_{k=1}^{3}\left[\Delta {C}_{{\ell }m}^{A,(k)}\cos ({\omega }^{(k)}t)+\Delta {C}_{{\ell }m}^{B,(k)}\sin ({\omega }^{(k)}t)\right],\\ {S}_{{\ell }m}(t)={\bar{S}}_{{\ell }m}+\mathop{\sum }\limits_{k=1}^{3}\left[\Delta {S}_{{\ell }m}^{A,(k)}\cos ({\omega }^{(k)}t)+\Delta {S}_{{\ell }m}^{B,(k)}\sin ({\omega }^{(k)}t)\right],\end{array}$$ (2)
in which \({\bar{C}}_{{\ell }m}\,{\rm{and}}\,{\bar{S}}_{{\ell }m}\) are static (that is, not changing with time) components of the Martian gravity field, ω(k) represents the angular frequency of the kth periodic component in which k = 1, 2, 3 respectively correspond to the Martian annual (or seasonal), semiannual and triannual periods, and t is time following the reference epoch (J2000). The amplitudes \(\Delta {C}_{{\ell }m}^{A,(k)}\), \(\Delta {S}_{{\ell }m}^{A,(k)}\), \(\Delta {C}_{{\ell }m}^{B,(k)}\) and \(\Delta {S}_{{\ell }m}^{B,(k)}\) represent the cosine and sine coefficients of the kth harmonic component, which together describe the temporal modulation of each spherical harmonic term at the modelled frequencies. For this work, we specifically extract coefficients at the Martian annual period for use in tidal tomography inversions (that is, \(\Delta {C}_{{\ell }m}^{A,B}\equiv \Delta {C}_{{\ell }m}^{A,B,(k=1)}\) and \(\Delta {S}_{{\ell }m}^{A,B}\equiv \Delta {S}_{{\ell }m}^{A,B,(k=1)}\), ω ≡ ω 1 = 1.058 × 10−7 rad s−1 and \(t=\frac{\omega }{2{\rm{\pi }}}n+{t}_{0}\) (n is an integer) are yearly intervals after the J2000 epoch t 0 ).
Gravity field inversion
We co-estimate parameters \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\) up to ℓ = 3 along with static coefficients \({\bar{C}}_{{\ell }m}\) and \({\bar{S}}_{{\ell }m}\) up to ℓ = 120 in equation (2) through our gravity inversion procedure. To do so, we reprocess Earth-based radiometric (X-band) Doppler tracking data acquired by the DSN from the MGS, ODY and MRO missions. The dataset spans approximately 16 years, corresponding to about one and a half solar cycles28. Following an initial step of correcting for non-gravitational effects on spacecraft acceleration, we minimize residuals between predicted spacecraft range rate values (that is, for an iteratively updated gravity field) and range rate observations. These steps are described in detail below.
Estimation of the gravitational coefficients was achieved through a batch least-squares analysis using orbital arcs of 2.5–8.0 days, with the arc length varying according to the mission phase. For each arc, partial derivatives of the Doppler observables were computed with respect to the estimated parameters and a global simultaneous inversion of all arcs and all parameters at once was performed to minimize the residuals between the observed and modelled tracking data. Our analysis represents an extension of previous analyses focused mainly on evaluating the temporal evolution of the zonal terms28,56. Results for the degree-3 and degree-2 zonal terms do not change beyond 2σ uncertainty when adding the non-zonal terms. Moreover, estimating only degree-3 time-varying terms or explicitly including non-zonal degree-2 terms does not substantially change degree-3 terms inferred from our analysis. As a further check, we performed a sensitivity analysis of post-fit residuals and empirical accelerations to the inclusion of the degree-3 time-variable gravity field using 2017 tracking data not used in the gravity field determination. The r.m.s. of fit for these independent test arcs improves only marginally (<1%) on inclusion of the degree-3 time-varying terms (Extended Data Fig. 7a). This is expected: long-period tidal parameters are constrained by stacking many years of data, so no individual arc is expected to show a marked improvement in fit. Similarly, the a posteriori amplitudes of the empirical accelerations change only marginally and without a systematic trend when the degree-3 time-varying terms are included (Extended Data Fig. 7b,c).
We correct for the contribution of the atmosphere on the Martian gravity field. To do this, we incorporate a time series of spherical harmonic coefficients representing the atmospheric load within the orbit determination procedure (see equation (4) in ref. 28). At each epoch, we expand the surface pressure P s (λ, ϕ, t) provided by the Mars Climate Database (MCD57) and include the solid-body response to this surface load through degree-dependent loading Love numbers \({k}_{{\ell }}^{{\prime} }\):
... continue reading