Modelling 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, Rm 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 t0).
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 Ps(λ, ϕ, 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} }\):
$$\begin{array}{c}\Delta {C}_{{\ell }m}^{{\rm{atm}}}(t)=\frac{3(1+{k}_{{\ell }}^{{\prime} })}{4{\rm{\pi }}\bar{\rho }g{R}_{m}(2{\ell }+1)}\oint {P}_{{\rm{s}}}(\lambda ,\phi ,t){P}_{{\ell }m}(\sin \phi )\cos (m\lambda ){\rm{d}}\Omega ,\\ \Delta {S}_{{\ell }m}^{{\rm{atm}}}(t)=\frac{3(1+{k}_{{\ell }}^{{\prime} })}{4{\rm{\pi }}\bar{\rho }g{R}_{m}(2{\ell }+1)}\oint {P}_{{\rm{s}}}(\lambda ,\phi ,t){P}_{{\ell }m}(\sin \phi )\sin (m\lambda ){\rm{d}}\Omega ,\end{array}$$
(3)
in which \(\bar{\rho }\) is Mars’s mean density, g is Mars’s mean surface gravity, Pℓm(sinϕ) are the same fully normalized associated Legendre functions used in equation (1), Ps(λ, ϕ, t) is the instantaneous surface pressure field and dΩ = cosϕdλdϕ is the element of solid angle on the unit sphere. The integrals project Ps(λ, ϕ, t) onto the cosine and sine basis of equation (2) and the factor \((1+{k}_{{\ell }}^{{\prime} })\) accounts for the deformation of the solid planet under this surface load. The loading Love numbers \({k}_{{\ell }}^{{\prime} }\) are computed for a spherically symmetric, self-gravitating 1D elastic Mars using the semi-analytic method described in ref. 58. The radially stratified reference interior used for this calculation is described in Extended Data Table 2 and yields \({k}_{2}^{{\prime} }=-0.14126\) and \({k}_{3}^{{\prime} }=-0.09109\). We note that atmospheric surface loads at ℓ = 2 and ℓ = 3 are comparable in magnitude, so mode coupling between harmonics (that is, ℓ = 2 to ℓ = 3) driven by lateral heterogeneity negligibly affect \({k}_{2}^{{\prime} }\) and \({k}_{3}^{{\prime} }\). For consistency, we also tested changes to the within-harmonic (diagonal) response of a laterally heterogeneous Martian interior with degree-1, order-(−1, 0, 1) shear modulus variations of 36%, 60%, 42% in the mantle (Fig. 2) using a finite-element code3,59,60,61,62,63,64,65. This model returns \({k}_{2}^{{\prime} }\) and \({k}_{3}^{{\prime} }\) values that agree with parameters from the 1D reference model to within <1%. We also test the impact of varying \({k}_{{\ell }}^{{\prime} }\) by 200% and 0% (relative to baseline values; ‘0%’ indicates a case with no deformational response to surface loads) on inversion results (Fig. 2) and find that these changes to assumed \({k}_{{\ell }}^{{\prime} }\) influence the median of recovered degree-1 shear modulus coefficient values by <10% (Extended Data Fig. 2).
The spherical harmonic expansions for the atmosphere of Mars in equation (3) are evaluated at 2-h intervals to capture diurnal and semidiurnal components and are included in the dynamical modelling of each spacecraft. The MCD and the Drag Temperature Model-Mars (DTM-Mars66) used for modelling non-gravitational forces (see below) apply to distinct altitude regimes (the lower atmosphere and the thermosphere, respectively) but share a common calibration against Mars atmospheric observations (and physical assumptions) and are therefore internally consistent for this analysis. The final retrieved static and time-varying gravity fields corrected for the atmosphere, given in equations (1) and (2) and Table 1 and shown as blue bars in Fig. 1, are therefore representative of the solid-body response only. Variations in the retrieved \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\) across these permutations remain stable on multi-year timescales (Extended Data Fig. 1a) and are generally smaller than the observational uncertainties reported in Fig. 1a, Table 1 and Extended Data Table 1 (Extended Data Fig. 1b); we correspondingly do not propagate these variations into the uncertainties for the atmosphere-corrected observations used as constraints on tidal tomography inversions. We do not correct for the seasonal condensation of atmospheric CO2 onto the surface (or associated deformational response) for non-zonal harmonics because related gravity perturbations, from both mass redistribution and surface deflection, are expected to be several orders of magnitude smaller than the observations in Table 1 or Extended Data Table 1 (refs. 58,67).
Past studies have tested the impact of several more modelling choices for the Martian atmosphere, including the incident solar flux (parameterized by means of the proxy quantity ‘extreme ultraviolet radiation’ or EUV) and the level of atmospheric dust opacity28,56. For example, Genova et al.28 find that, whereas the incident solar flux has a monotonic impact on gravity coefficients, dust opacity tends to amplify the impact of changes in incident solar flux on the time-variable gravity coefficients of Mars. Therefore, as an extra check, we reprocessed the atmospheric contributions from equation (3) to compare our baseline model (an average solar flux scenario, ‘Avg EUV’, F10.7 ≈ 130 sfu (ref. 57)) against endmember scenarios drawn from the MCD: (1) a cold-atmosphere case with minimum solar flux (F10.7 ≈ 70 sfu or solar flux units, for which 1 sfu = 10−22 W m−2 Hz−1 (ref. 57)) combined with dust-storm conditions (column-integrated visible optical depth τ ≳ 5 (ref. 68)) or ‘Cold min EUV dust’ and (2) a warm-atmosphere case with maximum solar flux (F10.7 ≈ 200 sfu (ref. 57)) combined with the same dust opacity or ‘Warm max EUV dust’. These differences between atmospheric models, when applied as corrections to the atmosphere-corrected observations shown in Fig. 1, propagate into ≤20% shifts in the inferred posterior distributions of degree-1 mantle shear modulus relative to our nominal inversion (Extended Data Fig. 2). The most substantial shifts (relative to the nominal case) occur for the ‘Cold min EUV dust’ scenario, although recovered posteriors for this case indicate a north–south pattern that still broadly aligns with the structure of the crustal dichotomy (degree-1 order-0 and degree-1 order-1 variations remain significant at 3σ and 2σ levels, respectively).
We correct for non-gravitational forces on each spacecraft with a detailed forward model of the non-conservative perturbations. To this end, we represent each spacecraft as a multi-panel body and assign thermo-optical properties to movable surfaces such as the high-gain antenna and solar arrays. We predict the thermospheric density along each orbit with DTM-Mars66 updated for the revised evolution of the main atmospheric constituents69. To absorb any residual mismodelling, we also co-estimate small periodic accelerations in the along-track and cross-track directions at frequencies tied to the orbital period, following ref. 28. We note that, although tidal forcing is theoretically stronger at diurnal frequencies than at the seasonal period1, the higher-frequency components of the time-variable gravity of Mars are also the most susceptible to residual non-gravitational errors. Estimated empirical accelerations therefore largely absorb diurnal tidal signals together with the perturbations they target28,56. For this reason, we do not attempt to recover diurnal tidal signals and restrict our analysis to the annual gravity field variations.
Modelling the tidal deformation of Mars
We directly model the impact of tidal deformation associated with Sun–Mars tides on coefficients ΔCℓm and ΔSℓm according to (superscript ‘p’ denotes ‘predicted’):
$$\Delta {C}_{{\ell }m}^{{\rm{p}}}-{\rm{i}}\Delta {S}_{{\ell }m}^{{\rm{p}}}={f}^{{\rm{p}}}(t)$$
(4)
in which:
$$\begin{array}{c}{f}^{{\rm{p}}}(t)=\frac{1}{{N}_{{\ell }m}}\mathop{\sum }\limits_{{{\ell }}^{{\prime} }=2}^{3}\mathop{\sum }\limits_{{m}^{{\prime} }=0}^{{{\ell }}^{{\prime} }}\frac{1}{2{{\ell }}^{{\prime} }+1}\frac{G{M}_{{\rm{S}}}}{{GM}}\frac{{R}^{{{\ell }}^{{\prime} }+1}}{{r}_{{\rm{S}}}^{{{\ell }}^{{\prime} }+1}}{P}_{{{\ell }}^{{\prime} }{m}^{{\prime} }}(\sin {\phi }_{S})\\ \,\times \left[({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },m}\cos ({m}^{{\prime} }{\lambda }_{{\rm{S}}})+{K}_{{{\ell }}^{{\prime} },-{m}^{{\prime} }}^{{\ell },m}\sin ({m}^{{\prime} }{\lambda }_{{\rm{S}}}))\right.\\ \,\left.-{\rm{i}}({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },-m}\cos ({m}^{{\prime} }{\lambda }_{{\rm{S}}})+{K}_{{{\ell }}^{{\prime} },-{m}^{{\prime} }}^{{\ell },-m}\sin ({m}^{{\prime} }{\lambda }_{{\rm{S}}}))\right].\end{array}$$
(5)
In equation (5), MS is the mass of the Sun and (λS, ϕS, rS) represent solar longitude, latitude and distance, respectively, in the body-fixed frame of Mars. We evaluate (λS, ϕS, rS) based on solar and Martian ephemerides extracted from NASA’s Planetary Data System evaluated over the mission durations of the MGS, ODY and MRO spacecraft. Note that equation(5 ) accounts for tides that arise from both Mars’s 0.0934 eccentricity and 25.2° spin-axis obliquity relative to the Sun through the impact of these orbital parameters on variations in λS, ϕS, rS. Nℓm is the normalization factor:
$${N}_{{\ell }m}=\sqrt{\frac{({\ell }-m)!\,(2-{\delta }_{0m})(2{\ell }+1)}{({\ell }+m)!}}$$
(6)
\({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },m}\) in equation (5) denotes ‘extended Love numbers’ (distinct from traditional Love numbers kℓm) that represents coupling between forcing at one harmonic (ℓ′, m′) and the gravitational response to this forcing at another harmonic (ℓ, m) for a given interior structure. We compute \({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },m}\) using the semi-analytic spectral method LOV3D (ref. 10) (see also refs. 70,71), which solves mass conservation, momentum and Poisson’s equations in the Fourier domain for a laterally heterogeneous body subject to tidal loading3,59,60,61,62. Given \({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },m}\), we compute \(\Delta {C}_{{\ell }m}^{A,B,p}\) and \(\Delta {S}_{{\ell }m}^{A,B,p}\) by evaluating fp(t) in equation (5) and then projecting this quantity into cosine and sine functions evaluated over the Martian annual period:
$$\begin{array}{c}\Delta {C}_{{\ell }m}^{A,{\rm{p}}}={\rm{\Re }}\,\left(\frac{\omega }{2{\rm{\pi }}}{\int }_{{t}_{{\rm{start}}}}^{{t}_{{\rm{end}}}}{f}^{{\rm{p}}}(t)\cos (\omega t){\rm{d}}t\right)\\ \Delta {C}_{{\ell }m}^{B,{\rm{p}}}={\rm{\Re }}\,\left(\frac{\omega }{2{\rm{\pi }}}{\int }_{{t}_{{\rm{start}}}}^{{t}_{{\rm{end}}}}{f}^{{\rm{p}}}(t)\sin (\omega t){\rm{d}}t\right)\\ \Delta {S}_{{\ell }m}^{A,{\rm{p}}}={\rm{\Im }}\,\left(-\frac{\omega }{2{\rm{\pi }}}{\int }_{{t}_{{\rm{start}}}}^{{t}_{{\rm{end}}}}{f}^{{\rm{p}}}(t)\cos (\omega t){\rm{d}}t\right)\\ \Delta {S}_{{\ell }m}^{B,{\rm{p}}}={\rm{\Im }}\,\left(-\frac{\omega }{2{\rm{\pi }}}{\int }_{{t}_{{\rm{start}}}}^{{t}_{{\rm{end}}}}{f}^{{\rm{p}}}(t)\sin (\omega t){\rm{d}}t\right).\end{array}$$
(7)
Note that equation (7) captures only the solid-body tidal response. However, zonal coefficients are also sensitive to seasonal CO2 exchange between the polar caps, which we model separately:
$$\Delta {C}_{{\ell }0}^{A,B,{\rm{p}}}=\Delta {C}_{{\ell }0}^{A,B,{\rm{p}},{\rm{solid}}}+\Delta {C}_{{\ell }0}^{A,B,{\rm{p}},{\rm{poles}}},$$
(8)
in which \(\Delta {C}_{{\ell }0}^{A,B,{\rm{p}},{\rm{solid}}}\) is computed from equation (7) and the polar-cap contribution \(\Delta {C}_{{\ell }0}^{A,B,{\rm{p}},{\rm{poles}}}\) is computed by representing the seasonal mass exchange as time-varying point masses ΔmN(t) and ΔmS(t) located at the north and south rotation poles, respectively. For point masses on the rotation axis, the corresponding (4π-normalized) coefficient is:
$$\Delta {C}_{{\ell }0}^{{\rm{poles}}}(t)=\frac{1}{M\sqrt{2{\ell }+1}}[\Delta {m}_{{\rm{N}}}(t)+{(-1)}^{{\ell }}\Delta {m}_{{\rm{S}}}(t)],$$
(9)
which, for the degree-2 and degree-3 zonal terms, reduces to
$$\Delta {C}_{20}^{{\rm{poles}}}(t)=\frac{\Delta {m}_{{\rm{N}}}(t)+\Delta {m}_{{\rm{S}}}(t)}{M\sqrt{5}},\,\Delta {C}_{30}^{{\rm{poles}}}(t)=\frac{\Delta {m}_{{\rm{N}}}(t)-\Delta {m}_{{\rm{S}}}(t)}{M\sqrt{7}}.$$
(10)
We use seasonal mass amplitudes of |ΔmN| = 6.2 × 1015 kg and |ΔmS| = 8.4 × 1015 kg (ref. 28), with ΔmN(t) and ΔmS(t) varying anti-phase over the Martian year (mass accumulates at one pole as it sublimates from the other). The harmonic components \(\Delta {C}_{{\ell }0}^{A,B,{\rm{p}},{\rm{poles}}}\) are then obtained by projecting \(\Delta {C}_{{\ell }0}^{{\rm{poles}}}(t)\) onto cos(ωt) and sin(ωt) following equation (7). The values of |ΔmN| and |ΔmS| from ref. 28 are empirically calibrated and implicitly absorb the solid-body loading response of the interior to seasonal polar-cap mass exchange. As an extra check of the sensitivity of the inclusion of zonal terms, we perform a tidal tomography inversion that entirely excludes zonal harmonics as constraints. The resulting posterior distributions for this scenario deviate by less than 15% from those of our nominal inversion (Extended Data Fig. 2).
Tidal tomography procedure
To constrain the structure of the Martian interior, we carry out a Bayesian inversion with Markov Chain Monte Carlo and a Metropolis sampling algorithm using PyMC (ref. 72). For our inversion, we vary elastic parameters relative to a reference Martian interior to fit observed \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\) (Table 1). Our detailed procedure is described below.
We consider the impact of 3D perturbations to the shear modulus of a 1D reference model on \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\). Our reference model is fixed (that is, not sampled) and is based on the analysis of ref. 29, which constrains Martian interior structure using the planet’s mean density, moment of inertia, seismic wave arrival times and the predicted composition of the Martian mantle (see Extended Data Table 2 for assumed mean shear modulus, bulk modulus and density values for each internal layer). We do not consider lateral changes in density or bulk modulus (the latter of which has a negligible effect on \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\); see ref. 9) for tidal tomography inversions (COM–COF offset is used only for a posteriori calculations described in Fig. 3). Because shear modulus (μ) is related to shear wave speed (Vs) and density (ρ) through \(\mu =\rho {V}_{{\rm{s}}}^{2}\), our approach effectively perturbs Vs while maintaining fixed ρ within each model layer. Note that our choice of reference model could substantially affect the non-coupled response of the interior of Mars (that is, \({K}_{{{\ell }}^{{\prime} },{m}^{{\prime} }}^{{\ell },m}\), in which ℓ, m = ℓ′, m′). However, from equation (5), forcing at ℓ = 3 harmonics is expected to be negligible (approximately three orders of magnitude smaller) compared with mode-coupled responses for these harmonics reported in Table 1 and so are ignored for inversions.
We consider a total of 30 parameters that describe 3D variations in shear modulus (that is, ℓ = 1–3 for the crust and mantle). These parameters are sampled as coefficients for spherical harmonic basis functions that comprise the ratio of the spatially variable shear modulus μ to that of the fixed reference model μref for each internal layer:
$$\begin{array}{l}\frac{\mu }{{\mu }_{{\rm{ref}}}}=\mathop{\sum }\limits_{{\ell }=1}^{3}\left[{\psi }_{{\ell }0}\,{\mu }_{{\ell }0}{Y}_{{\ell }0}(\theta ,\lambda )\right.\\ \,+\,\mathop{\sum }\limits_{m=1}^{{\ell }}\left.({\psi }_{{\ell }m}^{{\rm{c}}}\,{\mu }_{{\ell }m}^{{\rm{c}}}{Y}_{{\ell }m}^{{\rm{c}}}(\theta ,\lambda )+{\psi }_{{\ell }m}^{{\rm{s}}}\,{\mu }_{{\ell }m}^{{\rm{s}}}{Y}_{{\ell }m}^{{\rm{s}}}(\theta ,\lambda ))\right],\end{array}$$
(11)
in which the sampled coefficients are μℓ0 for m = 0 and \(({\mu }_{{\ell }m}^{{\rm{c}}},{\mu }_{{\ell }m}^{{\rm{s}}})\) for m ≥ 1. By normalizing 3D shear modulus variations by μref in equation (11), we implicitly account for the impact of variations in 1D shear modulus on mode-coupled responses such that our choice of reference model does not substantively affect inversion results. As an extra check, we have carried out tidal tomography inversions for which we sample the 1D shear and bulk moduli of the crust and mantle from a Gaussian distribution centred about our nominal reference model values (Extended Data Table 3) and with a standard deviation equal to 5% of these values (Extended Data Fig. 2). The resulting posterior distributions for recovered degree-1 variations in mantle shear modulus deviate minimally from those of our nominal inversion.
The real-form spherical harmonic basis functions (with θ the planetocentric co-latitude in the COM frame) are:
$$\begin{array}{c}{Y}_{{\ell }0}(\theta ,\lambda )=\sqrt{\frac{2{\ell }+1}{4{\rm{\pi }}}}{P}_{{\ell }0}(\cos \theta ),\\ {Y}_{{\ell }m}^{{\rm{c}}}(\theta ,\lambda )=\left\{\begin{array}{cc}\sqrt{\frac{(2{\ell }+1)}{2{\rm{\pi }}}\frac{({\ell }-m)!}{({\ell }+m)!}}{P}_{{\ell }m}(\cos \theta )\cos (m\lambda ),\, & m > 0,\\ 0,\, & m=0,\,\\ \,\end{array}\right.\\ {Y}_{{\ell }m}^{{\rm{s}}}(\theta ,\lambda )=\left\{\begin{array}{cc}\sqrt{\frac{(2{\ell }+1)}{2{\rm{\pi }}}\frac{({\ell }-m)!}{({\ell }+m)!}}{P}_{{\ell }m}(\cos \theta )\sin (m\lambda ),\, & m > 0,\\ 0,\, & m=0.\end{array}\right.\end{array}$$
(12)
The quantities ψℓ0, \({\psi }_{{\ell }m}^{{\rm{c}}}\) and \({\psi }_{{\ell }m}^{{\rm{s}}}\) indicate coefficients that normalize the right side of equation (11) so that an input \({\mu }_{{\ell }m},{\mu }_{{\ell }m}^{{\rm{c}}},\) or \({\mu }_{{\ell }m}^{{\rm{s}}}\) of 1 (assuming that all other sampled coefficients are zero) corresponds to a peak-to-peak variation in μ/μref of −1 to 1 (see, for example, equation (68) in ref. 9). Note that values presented in Fig. 2a–c indicate \({\mu }_{{\ell }m},{\mu }_{{\ell }m}^{{\rm{c}}},{\mu }_{{\ell }m}^{{\rm{s}}}\times 100{\rm{ \% }}\). To generate prior distributions, we parameterize \({\mu }_{{\ell }m},{\mu }_{{\ell }m}^{{\rm{c}}},{\mu }_{{\ell }m}^{{\rm{s}}}\) in equation (11) in terms of their base-10 logarithm (that is, \({\mu }_{{\ell }m}^{{\prime} },{\mu }_{{\ell }m}^{{{\rm{c}}}^{{\prime} }},{\mu }_{{\ell }m}^{{{\rm{s}}}^{{\prime} }}\)):
$$\begin{array}{l}{\mu }_{{\ell }m}^{{\prime} }\,=\,{\log }_{10}(1+{\mu }_{{\ell }m}),\\ \,{\mu }_{{\ell }m}^{{{\rm{c}}}^{{\prime} }}\,=\,{\log }_{10}(1+{\mu }_{{\ell }m}^{{\rm{c}}}),\\ \,{\mu }_{{\ell }m}^{{{\rm{s}}}^{{\prime} }}\,=\,{\log }_{10}(1+{\mu }_{{\ell }m}^{{\rm{s}}}).\end{array}$$
(13)
To generate ensembles of internal structure models for Markov chains, we sample uniform (that is, flat) prior probability distributions for \({\mu }_{{\ell }m}^{{\prime} },{\mu }_{{\ell }m}^{{{\rm{c}}}^{{\prime} }},{\mu }_{{\ell }m}^{{{\rm{s}}}^{{\prime} }}\) (that correspond to a range of \({\mu }_{{\ell }m},{\mu }_{{\ell }m}^{{\rm{c}}},{\mu }_{{\ell }m}^{{\rm{s}}}\times 100{\rm{ \% }}\) from −200% to 200%) and forward compute \(\Delta {C}_{{\ell }m}^{A,B,{\rm{p}}}\) and \(\Delta {S}_{{\ell }m}^{A,B,{\rm{p}}}\) using these values. Each ensemble consists of about 10,000 individual accepted model realizations (that is, approximately 500,000 samples total from 50 walkers). To speed up convergence, we use an adaptive sampling approach (the ‘tune’ functionality in PyMC) that dynamically adjusts step sizes based on the sensitivity of model outputs to input parameters72. We visually inspect Markov chains to discard initial burn-in steps (that is, typically the first roughly 10–20% of samples) and terminate inversions when parameter autocorrelation values are >0.99. Walker positions are updated on the basis of the likelihood function L:
$$\log L\propto -\frac{1}{2}{({\bf{X}}-{\bf{Y}})}^{{\rm{T}}}{{\boldsymbol{\Sigma }}}^{-1}({\bf{X}}-{\bf{Y}}),$$
(14)
in which X is the vector of observed \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\) in equation (2) (that is, parameters in Table 1 and Extended Data Table 1) and Y is the vector of model-predicted \(\Delta {C}_{{\ell }m}^{A,B,{\rm{p}}}\) and \(\Delta {S}_{{\ell }m}^{A,B,{\rm{p}}}\) in equation (7). For a comparison of model-predicted posteriors of \(\Delta {C}_{{\ell }m}^{A,B,{\rm{p}}}\) and \(\Delta {S}_{{\ell }m}^{A,B,{\rm{p}}}\) to observations, see Extended Data Fig. 3. Σ is a diagonal matrix (that is, each entry is independent) that incorporates observational covariances (that is, 15× formal uncertainties presented in Table 1). Other potential system constraints (for example, mean density, moment of inertia or quality factors) are not incorporated into vectors X or Y in equation (14).
Modelling temperature and iron content
In this section, we outline our methodology for the calculations shown in Fig. 3. We parameterize hemispheric variations in effective shear modulus in terms of the combined effect of temperature and composition on rheology:
$$\Delta \mu =\Delta T\frac{{\rm{d}}\mu }{{\rm{d}}T}+\Delta {X}_{{\rm{Fe}}}\frac{{\rm{d}}\mu }{{\rm{d}}{X}_{{\rm{Fe}}}},$$
(15)
in which ΔT is the hemispheric temperature contrast and ΔXFe is the hemispheric contrast in fayalite–forsterite mole fraction (Fa–Fo). To compute \(\frac{{\rm{d}}\mu }{{\rm{d}}T}\), we fit shear modulus versus forcing timescale data from Fig. 1a of ref. 34 (the ‘background-only’ Burgers model; note that all other calculations for tidal tomography inversions assume a purely elastic rheology) with power-law functions and extend these curves to the Martian annual forcing period (5.94 × 107 s; Extended Data Fig. 6) yielding:
$$\frac{{\rm{d}}\mu }{{\rm{d}}T}\approx -0.204\,{\rm{GPa}}\,{{\rm{K}}}^{-1}.$$
(16)
We note that several anelastic laws have been published (for example, refs. 73,74) but we have chosen that of ref. 34, as it has been used in several analogous applications for the Moon75, Mars76 and the Earth77,78. Although each anelastic law differs in details, the overall trend of shear modulus softening with lower frequency across all experimental laws is consistent. For compositional variations, we use \(\frac{{\rm{d}}\mu }{{\rm{d}}{X}_{{\rm{Fe}}}}=-0.3\,{\rm{GPa}}\,{\rm{per}}\,{\rm{ \% }}\) Fa–Fo based on experimental data presented in ref. 79.
We evaluate the impact of degree-1 variations in density on the COM–COF offset of Mars (that is, δzCOM-COF) by first assuming an offset vector of the form:
$$\delta {z}_{\text{COM-COF}}=\frac{1}{M}{\int }_{V}z\rho ({\bf{r}}){\rm{d}}V,$$
(17)
in which r is the position vector within Mars, ρ(r) is the local density, z is the component of r along the asymmetry axis, V is the volume of the planet and M is its mass. Converting equation (17) into spherical harmonics and integrating over the radial extent of the mantle, we obtain:
$$\delta {z}_{\text{COM-COF}}=\frac{1}{M}\sqrt{\frac{{\rm{\pi }}}{12}}\Delta \rho ({r}_{{\rm{Moho}}}^{4}-{r}_{{\rm{CMB}}}^{4}),$$
(18)
in which M = 6.4171 × 1023 kg is the mass of Mars, rMoho and rCMB are the radii of the Martian Moho and core–mantle boundary, respectively, and Δρ is the peak-to-peak amplitude of degree-1 variations in mantle density. We compute Δρ by assuming linear contributions from thermal expansion and iron enrichment:
$$\Delta \rho =-{\rho }_{{\rm{ref}}}\beta \Delta T+\Delta {X}_{{\rm{Fe}}}\frac{{\rm{d}}\rho }{{\rm{d}}{X}_{{\rm{Fe}}}},$$
(19)
in which ρref = 3,500 kg m−3 is a representative mantle density, β = 3 × 10−5 K−1 is the volume thermal expansion coefficient and \(\frac{{\rm{d}}\rho }{{\rm{d}}{X}_{{\rm{Fe}}}}=9.7\,{\rm{kg}}\,{{\rm{m}}}^{-3}\,{\rm{per}}\,{\rm{ \% }}\) Fa–Fo is the assumed density sensitivity to iron enrichment, following ref. 4. Equation (19) is formulated in terms of lateral density perturbations relative to a reference mantle, so the absolute mantle density and baseline iron content do not directly affect the calculation except insofar as they modify the adopted transfer function \(\frac{{\rm{d}}\rho }{{\rm{d}}{X}_{{\rm{Fe}}}}\).
For comparison with these calculations, we use COM–COF offsets from Table 4 of ref. 35, which gives Δx = 233 m, Δy = 1,428 m and Δz = 2,986 m, corresponding to a total offset magnitude of 3.3 km. Although Smith et al.35 do not report a formal uncertainty for this quantity, uncertainties are expected to be small (<0.02 km) relative to the 0–1.31-km range used here. We project the vector of the COM–COF offset onto the axis defined by the maximum shear modulus asymmetry in Fig. 2d (45° N, 138° W), which yields an upper-bound contribution of density variations in the mantle to COM–COF offset of 1.31 km. The lower bound of 0 km used in Fig. 3 corresponds to the opposite endmember case in which the observed COM–COF offset is explained entirely by crustal structure (for example, a case in which surface topography produces a COF shift of 1.31 km).
We note that the 0–1.31-km range is an approximate (that is, order-of-magnitude) estimate of the impact of mantle density on COM–COF offset and therefore does not encompass the full range of possible degree-1 density variations in the interior. For example, Moho relief can change the COM–COF offset by >50% if this effect is not negated by the impact of surface topography on the COM (for example, by means of Airy isostasy). Similarly, the depth dependence of \(\frac{{\rm{d}}\rho }{{\rm{d}}{X}_{{\rm{Fe}}}}\), including the impact of a possible density crossover near 600–700 km depth, could change the COM–COF offset by up to several hundred metres. Even so, these effects (and the 0–1.31-km range) are much smaller than the several tens of kilometres in COM–COF offset that would be required to explain the asymmetries shown in Fig. 3 purely in terms of compositional variations (for example, iron enrichment) in the mantle.
Martian melt ascent
Here we estimate the approximate timescale of ascent for melt for the mantle of Mars. McKenzie43 gives the characteristic timescale Λ for a small melt fraction to separate from a partially molten layer of thickness h:
$$\varLambda =\frac{h\eta \psi }{\kappa \Delta \rho g}$$
(20)
in which η is the melt viscosity, ψ is the melt fraction, Δρ is the density contrast between melt and matrix, g is the acceleration owing to gravity and κ is the permeability. For small melt fractions, the permeability is taken to be:
$$\kappa =\frac{{d}^{2}\,{\psi }^{n}}{C}$$
(21)
Here d is the grain size, n = 2 and C = 3,000. For the Martian mantle, we will take h = 300 km, Δρ = 100 kg m−3, g = 3.7 m s−2 and a grain size of 1 mm. If the areotherm crosses the solidus at 300 km depth, the temperature of the melts will be in the range 1,700–1,800 K. The resulting viscosity of Martian basaltic melts at these temperatures is roughly 1 Pa s (ref. 80). For a melt fraction of 0.02, the permeability is then about 10−13 m2 and the compaction timescale is approximately 3 Myr.
Ethics
This study did not involve human participants, human data or tissue, or animals, and therefore no ethical approval or informed consent was required. The work is based entirely on the analysis of publicly available spacecraft tracking and gravity data, and complies with all relevant institutional, national, and international guidelines and regulations.