Tidal tomography reveals a thermal anomaly beneath Mars’s crustal dichotomy

Mars exhibits a prominent dichotomy in crustal structure. Whereas low-lying plains cover the northern hemisphere of Mars, the southern hemisphere is more topographically elevated, rugged and densely cratered¹⁶. These asymmetries potentially accompany an approximately 25 km mean variation in crustal thickness¹¹ (or a roughly 200 kg m⁻³ variation in crustal density^17,18), as well as spatial differences in crustal magnetization¹⁹ and seismic wave attenuation²⁰. The northern lowlands may also have held surface water before the loss of Mars’s putative oceans to space²¹. Understanding the origin of the crustal dichotomy of Mars therefore constrains the planet’s ancient hydrology and the timing of surface conditions capable of supporting life²².

The origin of the hemispheric crustal dichotomy of Mars is widely debated. Some studies suggest that a giant impact disrupted the primordial mantle of the planet and either thickened the crust across the southern hemisphere^23,24 or excavated the northern hemisphere²⁵. Other studies suggest that mantle convection over one hemisphere drives subsidence and resurfacing of the northern crust or thickening and uplift of the southern crust^12,13. These formation hypotheses for the crustal dichotomy predict long-lived structures within the Martian deep interior that may persist into the present day^12,13,23,24. For example, crustal thickening over the southern hemisphere could enhance thermal insulation and radiogenic heating of the deeper interior, resulting in a structurally weak underlying mantle^14,15. Here we investigate the character of these potential heterogeneities at depth by analysing the gravitational response of Mars to seasonal (that is, 687 Earth days) tidal interactions with the Sun.

The Martian gravity field can be expressed in terms of spherical harmonic coefficients of degree ℓ and order m (C_ℓm and S_ℓm), with the full wavelength at a given ℓ usually defined as roughly 2πR/ℓ (R = 3,396 km)²⁶. Temporal variations can be modelled as cyclic perturbations to these coefficients and separated into in-phase (or ‘cosine’, A) and (a quarter cycle) out-of-phase (or ‘sine’, B) components over a given period (that is, \(\Delta {C}_{{\ell }m}^{A,B}\) and \(\Delta {S}_{{\ell }m}^{A,B}\)) (Methods). For a spherically symmetric planet, forcing at a given degree and order excites deformation only at the same degree and order. In this case, tidal forcing acts almost entirely at degree-2 and will contribute negligibly to the gravity field of Mars at higher degrees (that is, \(\Delta {C}_{3m}^{A,B},\Delta {S}_{3m}^{A,B}\approx 0\)). However, lateral heterogeneity in the internal structure of Mars (for example, north–south variations in shear modulus) can couple with degree-2 forcing to produce substantial degree-3 signals such that \(\Delta {C}_{3m}^{A,B},\Delta {S}_{3m}^{A,B}\ne 0\) (ref. ⁸). Moreover, although static gravity fields carry some sensitivity to deep-seated density structure (especially at long wavelengths), their sensitivity is weighted towards structure in the shallower interior at all degrees (equation (4) of ref. ²⁷). By contrast, the low-order non-zonal (that is, m ≠ 0) time-variable gravity field of Mars is most sensitive to lateral variations in the shear modulus of the mantle, as these heterogeneities generate larger mass redistributions than comparable perturbations within shallower layers^4,9. Measuring time-variable gravity signals (for example, by analysing the trajectory of orbiting spacecraft) correspondingly provides a direct means to constrain the extent of Mars’s deep-seated heterogeneities⁴.

To recover the seasonal degree-3 tidal signals of Mars, we analyse Earth-based X-band Doppler tracking data acquired by the NASA Deep Space Network (DSN) from several Martian orbiters. We use radio science data from MGS, ODY and MRO spanning a total of 16 years. Our data processing procedure follows the approach used to determine the time-varying zonal (m = 0) gravity field of Mars in ref. ²⁸ but is extended in this study to include non-zonal ℓ = 2 and ℓ = 3 time-varying gravity field coefficients over the seasonal (or annual) 687-day period.

Recovered non-zonal degree-3 gravity field coefficient perturbations are presented in Table 1 (for zonal and degree-2 coefficients, see Extended Data Table 1). To account for modelling uncertainties, particularly those arising from unmodelled non-gravitational accelerations, we inflate the formal uncertainties of these coefficients by a factor of 15 (see, for example, ref. ⁴) in Table 1 and in all subsequent analysis. As in ref. ²⁸, we directly model the effects of the atmosphere by integrating the impact of the Mars general circulation model on gravity coefficients (Methods and Extended Data Figs. 1 and 2).

We recover non-zonal ℓ = 3 perturbations to the Martian gravity field that deviate substantially from values predicted for a spherically symmetric Mars with an atmosphere (Fig. 1a). For example, the cosine coefficient of the ℓ = 3, m = 1 gravity field perturbation (that is, \(\Delta {C}_{31}^{A}\)) deviates approximately 300% from expectations for a spherically symmetric Mars with atmospheric loading (with >99.99% confidence, the difference between the green and red bars in Fig. 1a), suggesting coupling between tidal deformation and laterally heterogeneous internal structure in the mantle of Mars (Fig. 1b). To constrain the nature of these asymmetries, we perform a Bayesian, Markov chain Monte Carlo inversion (Methods). Inversions incorporate both degree-2 and degree-3 (zonal and non-zonal) time-variable gravity fields of Mars as constraints (Table 1, Extended Data Table 1 and Methods).

Our model space consists of lateral shear modulus variations imposed onto a 1D reference interior constrained by Mars’s degree-2 tidal Love number, moment of inertia and seismic travel-time data²⁹ (Extended Data Table 2). We expand these variations in a spherical harmonic basis truncated at degree ℓ = 3 across two internal layers: the crust (0–50 km depth) and the mantle (50–1,560 km depth). Lateral crustal thickness and density variations on Mars negligibly (<0.3%) affect the time-variable gravity field and are therefore ignored for inversions. Heterogeneities in the core (including a potential solid inner layer³⁰) also minimally affect time-variable gravity fields (Fig. 1b). We consequently prescribe a uniform elastic structure below 1,560 km depth. We use the semi-analytical spectral code LOV3D (ref. ¹⁰) to forward compute coupling between degree-2 forcing and lateral heterogeneity for candidate interior structures (Methods). Note that the zonal gravity field variations of Mars are also sensitive to seasonal CO₂ condensation over the polar caps³¹. We accordingly model m = 0 terms as the sum of the solid-body tide and a seasonal mass exchange of 6.2 × 10¹⁵ kg and 8.4 × 10¹⁵ kg over the northern and southern polar caps, respectively²⁸ (Methods).

Our results indicate that degree-1 (ℓ = 1) shear modulus structure, which reflects a hemispheric pattern, enhances the amplitude of ℓ = 3 coefficients for Martian interiors subject to seasonal tidal forcing at ℓ = 2 (Fig. 1b). For example, the ℓ = 2, m = 1 harmonic associated with the obliquity tide of Mars interacts with north–south (that is, ℓ = 1, m = 0) shear modulus structure (lower/higher shear modulus values in the mid-latitudes of the southern/northern portions of the meridional hemisphere that spans −90° E to 90° E longitude) to regionally increase/decrease outward radial deformation. The resulting mass displacement yields an ℓ = 3, m = 1 gravity signature that enhances the amplitude of \(\Delta {C}_{31}^{A}\). Interaction between other hemispherical variations such as east–west (ℓ = 1, m = −1) and meridional (ℓ = 1, m = 1) asymmetries and the full degree-2 Martian tide (that is, components arising from both obliquity and eccentricity of Mars) is more complex¹⁰ but produces deviation of several degree-3 gravity coefficients from zero (Extended Data Fig. 3). By contrast, the degree-2 time-variable gravity field does not measurably deviate from predictions for a spherically symmetric interior (Extended Data Fig. 3), suggesting minimal internal heterogeneity that strongly couples with ℓ = 2 forcing¹⁰ (for example, ℓ = 2 patterns; Extended Data Fig. 4).

Our inversions predict statistically significant degree-1 variations in the shear modulus of the Martian mantle. Specifically, we resolve north–south, meridional and east–west patterns with peak-to-peak amplitudes (and 3σ uncertainties) of 60 ± 40%, 42 ± 38% and 36 ± 34%, respectively (Fig. 2a–c), as well as an overall variation of 81 ± 60% (or >20%). Inversions fit all time-variable gravity field constraints to within 3σ uncertainty (Extended Data Fig. 3) and do not recover any statistically significant lateral variations in shear modulus for the crust of Mars (Extended Data Table 3 and Extended Data Fig. 4a) or for degrees higher than ℓ = 1 in the mantle (Extended Data Table 3 and Extended Data Fig. 4b). Note that heterogeneity amplitudes represent averages over the entire mantle depth range. This averaging is necessary to produce statistically significant results, because the sensitivity of degree-3 time-variable gravity to 3D structure varies gradually with depth (Fig. 1b), resulting in substantial non-uniqueness between recovered heterogeneity amplitudes when the mantle is subdivided into several layers (Extended Data Fig. 5). Even so, our inversions do not necessarily preclude more vertically localized structure within the mantle of Mars (for example, at the base of the lithosphere¹⁴ or near the core–mantle boundary³²).

The hemispheric mantle shear modulus variation of Mars broadly matches the spatial pattern of the crustal dichotomy (Fig. 2d). For example, combining maximum likelihood solutions for recovered degree-1 order-1, order-0 and order-−1 variations (Fig. 2a–c) results in a region of increased shear modulus beneath the northern lowlands (centred at 45° N, 138° W over Vastitas Borealis; Fig. 2d) and a corresponding reduction in shear modulus beneath the southern highlands (centred at 45° S, 42° E near Hellas basin; Fig. 2d). The zero-value contour of the mantle shear modulus variation of Mars also tracks the boundary of the crustal dichotomy, including its southward (or northward) deflection over the anti-meridional (or meridional) hemisphere. The correspondence between surface geology and mantle stiffness variations weakens between 40° W and 139° W near the volcanic Tharsis rise region, potentially because of overprinting of the dichotomy boundary by this structure¹⁶. The lack of anomalies observed directly beneath Tharsis (Fig. 2d) also suggests that present-day internal structures associated with this feature are small in amplitude (that is, compared with asymmetries in Fig. 2), shallow or restricted to a relatively small (ℓ > 3) lateral region within the interior³³.

Our modelling suggests that the anomalous degree-3 gravity field of Mars arises owing to a substantial, >20% lateral variation in the mantle’s effective shear modulus at the seasonal timescale. What could generate such a large contrast? At a given pressure, the rigidity of mantle rock depends on both temperature and composition. At seismic timescales, crystalline temperature differences would need to exceed an unrealistically large >1,000 K to produce observed shear modulus variations⁴. However, when extended to the Martian annual period (Methods and Extended Data Fig. 6), the effective shear modulus of olivine becomes approximately 15–20 times more sensitive to temperature³⁴. As such, recovered shear modulus variations can mostly be explained by a hemispheric temperature anomaly of approximately 200–400 ºC in the present-day Martian mantle (Fig. 3 and Methods). Thermal modelling also suggests that the insulating effect of the thicker crust in the southern highlands can result in a hemispheric temperature contrast >200 ºC (Fig. 3f of ref. ¹⁵), consistent with our results. The centre-of-mass–centre-of-figure (COM–COF) offset of Mars also limits the potential compositional component (for example, iron, water content) of asymmetries because, in isolation, such differences (and their associated density structure) would need to produce an offset that is approximately 50 times larger than observations³⁵ to account for the inferred shear modulus variations. Even so, our models allow for up to 5% iron enrichment (that is, variation in forsterite–fayalite mole fraction or ΔFo–Fa) within the mantle of the southern highlands (Fig. 3).

Mantle heterogeneities may preserve signatures of, and thus provide independent constraints on, geodynamic processes that influence the evolution of the crustal dichotomy. We consider three such candidate processes: a giant impact, spontaneous degree-1 upwelling (that persists into the present day) and enhanced thermal insulation or radiogenic heating of the southern highlands mantle by a thick overlying crust. Spontaneous degree-1 upwelling or thermal insulation would produce a mainly thermal mantle anomaly^12,13,14,15 as well as potential secondary compositional signatures through decompression melting below the lithosphere^12,13,14,15. However, enhanced melt extraction in these scenarios would deplete the southern highlands mantle in iron^12,36, opposite to the observations (Fig. 3). By contrast, a giant impact is (in isolation) not expected to result in a long-term secular variation in temperature of several hundred kelvin over a given hemisphere³⁷. A hybrid scenario could resolve this mismatch. For example, following the formation of the northern lowland crust from the Borealis impact, the northern hemisphere mantle is first depleted in iron³⁶, leaving the southern mantle relatively enriched. The resultant thicker southern crust could also trigger preferential thermal insulation¹⁵ or degree-1 upwelling³⁸ below the highlands of Mars (the latter would be favoured absent variations in crustal thickness or radiogenic material abundance¹⁸), increasing the present-day temperature of this region (Figs. 3 and 4).

Temperature variations may also influence the dissipative properties of the mantle of Mars. For example, comparisons of low-frequency marsquakes detected by the InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) lander indicate substantially lower quality factors (Q ≈ 500) for events originating in the nearby Terra Cimmeria region of the southern highlands than for events from Cerberus Fossae in the northern lowlands (Q ≈ 800–2,000)²⁰. These differences in attenuation can be attributed to lateral variations in the temperature of the mantle of a few hundred kelvin (Fig. 4 of ref. ²⁰). A warm southern highlands mantle might also result in localized melting¹⁴, which could reduce the 200–400 ºC temperature anomaly required to explain gravity field observations (Fig. 3). However, the Q values inferred for the southern highlands remain well above those expected for partially molten olivine (Q ≈ 100 or lower)²⁰. Even so, future electromagnetic sounding measurements (for example, ref. ³⁹) could resolve whether the inferred thermal anomaly in the southern highlands mantle is associated with isolated magma pockets or with a more continuous molten layer at depth (for example, refs. ^14,32).

The preferential magnetization of the southern highlands crust¹⁹ may also reflect the effects of a laterally heterogeneous interior. For example, upwelling beneath the southern highlands may have heated the overlying crust above Curie temperatures while Mars’s ancient dynamo was active (4.5–4.1 billion years ago)⁴⁰. These temperatures may have subsequently fallen below Curie temperatures before the cessation of the dynamo, allowing affected portions of the southern highlands to acquire their present-day remanent magnetization^13,41. Alternatively, increased heat flow at the core–mantle boundary (and in the overlying mantle) may have regionally enhanced convection in the outer core of Mars, resulting in a stronger dynamo below the southern highlands⁴². Inferred asymmetries in mantle structure may therefore constitute a present-day remnant of geodynamic processes associated with the generation of Mars’s early dynamo (Fig. 4).

Although melt is not required by our results (Fig. 3), some thermal models predict that higher temperatures in the southern highlands mantle would increase magma production in this region¹⁴. Such pockets of deep-seated melt are expected to ascend and erupt onto the surface over timescales of several Myr (ref. ⁴³) (Methods). Yet recent (tens of Myr) volcanic activity is confined to Cerberus Fossae in the north⁴⁴. We propose two potential solutions to this apparent paradox. The first is that the thicker (and potentially lower-density⁴⁵) crust in the southern highlands may present a barrier to melt ascent. Under this scenario, melt would stall in the mid-crust and generate local intrusive features (Fig. 4). These intrusions would be difficult to detect in surface images but may be resolved from very-high-resolution (ℓ > 300) measurements of the static gravity field of Mars⁴⁶. The second solution is that melt ascent and eruption requires a background extensional stress field⁴⁷. Tectonic features in the southern highlands are dominantly compressional (for example, ref. ⁴⁸), as expected from the slow cooling of Mars. By contrast, Cerberus Fossae is a region undergoing localized extension⁴⁹. In any case, the petrology of volcanic rocks sampled from different locations over the crust of Mars could constrain the spatial differences in mantle potential temperature implied by our results⁵⁰.

Tidal tomography—or the inference of lateral variations in interior structure from a body’s time-variable response to tidal forcing—has been used to investigate the mantle of Earth² and, more recently, the Moon⁴. Here we use nearly two decades of precise spacecraft tracking data to extend this approach to Mars. Future missions with dedicated gravity experiments (similar to GRACE for Earth⁵¹ or GRAIL for the Moon⁵²) could resolve similar (or finer) variations in Martian mantle structure over much shorter mission durations⁵³. In the future, these techniques may also be applied to other planetary bodies that exhibit pronounced low-order asymmetries in structure, including Ganymede, Mercury, Io and Enceladus. Because tidal tomography largely relies on remote measurements, it represents a powerful tool for future missions seeking to characterize the interiors of planetary bodies without landed spacecraft.

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:

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. ²⁶), 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 framework⁵⁵, 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:

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.058 × 10⁻⁷ rad s⁻¹ and \(t=\frac{\omega }{2{\rm{\pi }}}n+{t}_{0}\) (n is an integer) are yearly intervals after the J2000 epoch t₀).

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 cycles²⁸. 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 terms^28,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. ²⁸). At each epoch, we expand the surface pressure P_s(λ, ϕ, t) provided by the Mars Climate Database (MCD⁵⁷) and include the solid-body response to this surface load through degree-dependent loading Love numbers \({k}_{{\ell }}^{{\prime} }\):

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), P_s(λ, ϕ, 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 P_s(λ, ϕ, 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. ⁵⁸. 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 code^{3,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-Mars⁶⁶) 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 CO₂ 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 opacity^28,56. For example, Genova et al.²⁸ 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’, F_10.7 ≈ 130 sfu (ref. ⁵⁷)) against endmember scenarios drawn from the MCD: (1) a cold-atmosphere case with minimum solar flux (F_10.7 ≈ 70 sfu or solar flux units, for which 1 sfu = 10⁻²² W m⁻² Hz⁻¹ (ref. ⁵⁷)) combined with dust-storm conditions (column-integrated visible optical depth τ ≳ 5 (ref. ⁶⁸)) or ‘Cold min EUV dust’ and (2) a warm-atmosphere case with maximum solar flux (F_10.7 ≈ 200 sfu (ref. ⁵⁷)) 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-Mars⁶⁶ updated for the revised evolution of the main atmospheric constituents⁶⁹. 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. ²⁸. We note that, although tidal forcing is theoretically stronger at diurnal frequencies than at the seasonal period¹, 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 target^28,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’):

in which:

In equation (5), M_S is the mass of the Sun and (λ_S, ϕ_S, r_S) represent solar longitude, latitude and distance, respectively, in the body-fixed frame of Mars. We evaluate (λ_S, ϕ_S, r_S) 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, r_S. N_ℓm is the normalization factor:

\({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. ¹⁰) (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 loading^{3,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 f^p(t) in equation (5) and then projecting this quantity into cosine and sine functions evaluated over the Martian annual period:

Note that equation (7) captures only the solid-body tidal response. However, zonal coefficients are also sensitive to seasonal CO₂ exchange between the polar caps, which we model separately:

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 Δm_N(t) and Δm_S(t) located at the north and south rotation poles, respectively. For point masses on the rotation axis, the corresponding (4π-normalized) coefficient is:

which, for the degree-2 and degree-3 zonal terms, reduces to

We use seasonal mass amplitudes of |Δm_N| = 6.2 × 10¹⁵ kg and |Δm_S| = 8.4 × 10¹⁵ kg (ref. ²⁸), with Δm_N(t) and Δm_S(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 |Δm_N| and |Δm_S| from ref. ²⁸ 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. ⁷²). 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. ²⁹, 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. ⁹) 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 (V_s) and density (ρ) through \(\mu =\rho {V}_{{\rm{s}}}^{2}\), our approach effectively perturbs V_s 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:

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:

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. ⁹). 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} }}\)):

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 parameters⁷². 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:

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:

in which ΔT is the hemispheric temperature contrast and ΔX_Fe 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. ³⁴ (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 × 10⁷ s; Extended Data Fig. 6) yielding:

We note that several anelastic laws have been published (for example, refs. ^73,74) but we have chosen that of ref. ³⁴, as it has been used in several analogous applications for the Moon⁷⁵, Mars⁷⁶ and the Earth^77,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. ⁷⁹.

We evaluate the impact of degree-1 variations in density on the COM–COF offset of Mars (that is, δz_COM-COF) by first assuming an offset vector of the form:

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:

in which M = 6.4171 × 10²³ kg is the mass of Mars, r_Moho and r_CMB 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:

in which ρ_ref = 3,500 kg m⁻³ is a representative mantle density, β = 3 × 10⁻⁵ K⁻¹ 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. ⁴. 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. ³⁵, 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.³⁵ 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. McKenzie⁴³ gives the characteristic timescale Λ for a small melt fraction to separate from a partially molten layer of thickness h:

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:

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⁻³, g = 3.7 m s⁻² 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. ⁸⁰). For a melt fraction of 0.02, the permeability is then about 10⁻¹³ m² 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.

1. Hoolst, T. V., Dehant, V., Roosbeek, F. & Lognonné, P. Tidally induced surface displacements, external potential variations, and gravity variations on Mars. Icarus 161, 281–296 (2003).

   Article  ADS  Google Scholar 

2. Lau, H. C. P. et al. Tidal tomography constrains Earth’s deep-mantle buoyancy. Nature 551, 321–326 (2017).

   Article  ADS  CAS  PubMed  Google Scholar 

3. Berne, A., Simons, M., Keane, J. T., Leonard, E. J. & Park, R. S. Jet activity on Enceladus linked to tidally driven strike-slip motion along tiger stripes. Nat. Geosci. 17, 385–391 (2024).

   Article  ADS  CAS  Google Scholar 

4. Park, R. S. et al. Thermal asymmetry in the Moon’s mantle inferred from monthly tidal response. Nature 641, 1188–1192 (2025).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

5. Smith, D. E. et al. The gravity field of Mars: results from Mars Global Surveyor. Science 286, 94–97 (1999).

   Article  ADS  CAS  PubMed  Google Scholar 

6. Saunders, R. S. et al. 2001 Mars Odyssey mission summary. Space Sci. Rev. 110, 1–36 (2004).

   Article  ADS  CAS  Google Scholar 

7. Zuber, M. T. et al. Mars Reconnaissance Orbiter Radio Science Gravity Investigation. J. Geophys. Res. Planets 112, E05S07 (2007).

   Article  Google Scholar 

8. Zhong, S., Qin, C., A, G. & Wahr, J. Can tidal tomography be used to unravel the long-wavelength structure of the lunar interior. Geophys. Res. Lett. 39, L24201 (2012).

   Article  Google Scholar 

9. Qin, C., Zhong, S. & Wahr, J. Elastic tidal response of a laterally heterogeneous planet: a complete perturbation formulation. Geophys. J. Int. 207, 89–110 (2016).

   Article  ADS  Google Scholar 

10. Rovira-Navarro, M., Matsuyama, I. & Berne, A. A spectral method to compute the tides of laterally heterogeneous bodies. Planet Sci. J. 5, 129 (2024).

    Article  Google Scholar 

11. Wieczorek, M. A. & Zuber, M. T. Thickness of the Martian crust: improved constraints from geoid-to-topography ratios. J. Geophys. Res. Planets 109, E01009 (2004).

    Article  ADS  Google Scholar 

12. Zhong, S. & Zuber, M. T. Degree-1 mantle convection and the crustal dichotomy on Mars. Earth Planet. Sci. Lett. 189, 75–84 (2001).

    Article  ADS  CAS  Google Scholar 

13. Roberts, J. H. & Zhong, S. Degree-1 convection in the Martian mantle and the origin of the hemispheric dichotomy. J. Geophys. Res. Planets 111, E06013 (2006).

    Article  ADS  Google Scholar 

14. Gibet, V. B., Michaut, C., Wieczorek, M. & Lognonné, P. A positive feedback between crustal thickness and melt extraction for the origin of the Martian dichotomy. J. Geophys. Res. Planets 127, e2022JE007472 (2022).

    Article  ADS  Google Scholar 

15. Plesa, A. C. et al. The thermal state and interior structure of Mars. Geophys. Res. Lett. 45, 12198–12209 (2018).

    Article  ADS  Google Scholar 

16. Andrews-Hanna, J. C., Zuber, M. T. & Banerdt, W. B. The Borealis basin and the origin of the martian crustal dichotomy. Nature 453, 1212–1215 (2008).

    Article  ADS  CAS  PubMed  Google Scholar 

17. Kim, D. et al. Global crustal thickness revealed by surface waves orbiting Mars. Geophys. Res. Lett. 50, e2023GL103482 (2023).

    Article  ADS  Google Scholar 

18. Wieczorek, M. A. et al. Insight constraints on the global character of the Martian crust. J. Geophys. Res. Planets 127, e2022JE007298 (2022).

    Article  ADS  Google Scholar 

19. Lillis, R. J., Frey, H. V. & Manga, M. Rapid decrease in Martian crustal magnetization in the Noachian era: implications for the dynamo and climate of early Mars. Geophys. Res. Lett. 35, L14203 (2008).

    Article  ADS  Google Scholar 

20. Sun, X. & Tkalčić, H. Constraints on the origin of the Martian dichotomy from southern highlands marsquakes. Geophys. Res. Lett. 51, e2024GL110921 (2024).

    Google Scholar 

21. Jakosky, B. M. et al. Loss of the Martian atmosphere to space: present-day loss rates determined from MAVEN observations and integrated loss through time. Icarus 315, 146–157 (2018).

    Article  ADS  CAS  Google Scholar 

22. Hurowitz, J. A. et al. Redox-driven mineral and organic associations in Jezero Crater, Mars. Nature 645, 332–336 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

23. Reese, C. C., Orth, C. P. & Solomatov, V. S. Impact megadomes and the origin of the martian crustal dichotomy. Icarus 213, 433–442 (2011).

    Article  ADS  Google Scholar 

24. Cheng, K. W. et al. Mars’s crustal and volcanic structure explained by southern giant impact and resulting mantle depletion. Geophys. Res. Lett. 51, e2023GL105910 (2024).

    Article  ADS  Google Scholar 

25. Marinova, M. M., Aharonson, O. & Asphaug, E. Mega-impact formation of the Mars hemispheric dichotomy. Nature 453, 1216–1219 (2008).

    Article  ADS  CAS  PubMed  Google Scholar 

26. Konopliv, A. S. et al. Mars high resolution gravity fields from MRO, Mars seasonal gravity, and other dynamical parameters. Icarus 211, 401–428 (2011).

    Article  ADS  Google Scholar 

27. Hemingway, D. J. & Mittal, T. Enceladus’s ice shell structure as a window on internal heat production. Icarus 332, 111–131 (2019).

    Article  ADS  Google Scholar 

28. Genova, A. et al. Seasonal and static gravity field of Mars from MGS, Mars Odyssey and MRO radio science. Icarus 272, 228–245 (2016).

    Article  ADS  Google Scholar 

29. Stähler, S. C. et al. Seismic detection of the martian core. Science 373, 443–448 (2021).

    Article  ADS  PubMed  Google Scholar 

30. Bi, H. et al. Seismic detection of a 600-km solid inner core in Mars. Nature 645, 67–72 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

31. Smith, D. E., Zuber, M. T., Haberle, R. M., Rowlands, D. D. & Murphy, J. R. The Mars seasonal CO₂ cycle and the time variation of the gravity field: a general circulation model simulation. J. Geophys. Res. 104, 1885–1896 (1999).

    Article  ADS  CAS  Google Scholar 

32. Samuel, H. et al. Geophysical evidence for an enriched molten silicate layer above Mars’s core. Nature 622, 712–717 (2023).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

33. Broquet, A. & Andrews-Hanna, J. C. Geophysical evidence for an active mantle plume underneath Elysium Planitia on Mars. Nat. Astron. 7, 160–169 (2023).

    Article  ADS  Google Scholar 

34. Jackson, I. & Faul, U. H. Grainsize-sensitive viscoelastic relaxation in olivine: towards a robust laboratory-based model for seismological application. Phys. Earth Planet. Inter. 183, 151–163 (2010).

    Article  ADS  Google Scholar 

35. Smith, D. E. et al. Mars Orbiter Laser Altimeter: experiment summary after the first year of global mapping of Mars. J. Geophys. Res. Planets 106, 23689–23722 (2001).

    Article  ADS  Google Scholar 

36. Shaw, D. M. Trace element fractionation during anatexis. Geochim. Cosmochim. Acta 34, 237–243 (1970).

    Article  ADS  CAS  Google Scholar 

37. Ruedas, T. & Breuer, D. On the relative importance of thermal and chemical buoyancy in regular and impact-induced melting in a Mars-like planet. J. Geophys. Res. Planets 122, 1554–1579 (2017).

    Article  ADS  CAS  Google Scholar 

38. Citron, R. I., Manga, M. & Tan, E. A hybrid origin of the Martian crustal dichotomy: degree-1 convection antipodal to a giant impact. Earth Planet. Sci. Lett. 491, 58–66 (2018).

    Article  ADS  CAS  Google Scholar 

39. Civet, F. & Tarits, P. Electrical conductivity of the mantle of Mars from MGS magnetic observations. Earth Planets Space 66, 85 (2014).

    Article  ADS  Google Scholar 

40. Ruiz, J. The very early thermal state of Terra Cimmeria: implications for magnetic carriers in the crust of Mars. Icarus 203, 173–179 (2009).

    Article  Google Scholar 

41. Nimmo, F. & Stevenson, D. J. Influence of early plate tectonics on the thermal evolution and magnetic field of Mars. J. Geophys. Res. Planets 105, 11969–11979 (2000).

    Article  ADS  Google Scholar 

42. Yan, C. et al. An ancient Martian dynamo driven by hemispheric heating: effect of thermal boundary conditions. Planet. Sci. J. 4, 11 (2023).

    Article  Google Scholar 

43. McKenzie, D. Some remarks on the movement of small melt fractions in the mantle. Earth Planet. Sci. Lett. 95, 53–72 (1989).

    Article  ADS  Google Scholar 

44. Berman, D. C. & Hartmann, W. K. Recent fluvial, volcanic, and tectonic activity on the Cerberus plains of Mars. Icarus 159, 1–17 (2002).

    Article  ADS  Google Scholar 

45. Phillips, M. S. et al. Widespread ancient anorthosites in the lower crust of Mars. Commun. Earth Environ. 6, 1026 (2025).

    Article  Google Scholar 

46. Losarcos, J. M. & Dombard, A. J. Lava tubes on Earth, the Moon, and Mars: detection, evolution, and exploration potential. Space Sci. Rev. 221, 127 (2025).

    Article  ADS  CAS  Google Scholar 

47. Solomon, S. C. On volcanism and thermal tectonics on one-plate planets. Geophys. Res. Lett. 5, 461–464 (1978).

    Article  ADS  Google Scholar 

48. Andrews-Hanna, J. C. & Broquet, A. The history of global strain and geodynamics on Mars. Icarus 395, 115476 (2023).

    Article  Google Scholar 

49. Stähler, S. C. et al. Tectonics of Cerberus Fossae unveiled by marsquakes. Nat. Astron. 6, 1376–1386 (2022).

    Article  ADS  Google Scholar 

50. Baratoux, D., Toplis, M. J., Monnereau, M. & Gasnault, O. Thermal history of Mars inferred from orbital geochemistry of volcanic provinces. Nature 472, 338–341 (2011).

    Article  ADS  CAS  PubMed  Google Scholar 

51. Wahr, J., Swenson, S., Zlotnicki, V. & Velicogna, I. Time-variable gravity from GRACE: first results. Geophys. Res. Lett. 31, L11501 (2004).

    Article  ADS  Google Scholar 

52. Park, R. S., Vaughan, A. T. et al. Gravity Imaging Radio Observer (GIRO) for planetary science and mission opportunities. Planet. Sci. J. 6, 127 (2025).

    Article  Google Scholar 

53. Keane, J. T. et al. Next-generation planetary geodesy: results from the 2021 Keck Institute for Space Studies Workshops. In Proc. 53rd Lunar and Planetary Science Conference, LPI Contribution No. 2678, id. 1622 https://ui.adsabs.harvard.edu/abs/2022LPICo2678.1622K (2022).

54. Smith, D. E. & MOLA Science Team. MGS-M-MOLA-5-IEGDR-L3-V1.0: IEG100_A ASCII topography/radius table. Initial Experiment Gridded Data Record (IEGDR), 1° × 1° grid, ASCII table format. NASA Planetary Data System (PDS), Geosciences Node https://pds-geosciences.wustl.edu/mgs/mgs-m-mola-5-iegdr-l3-v1/mgsl_2xxx/data/ieg100_a.tab (1999).

55. Pavlis, D. E. & Nicholas, J. B. GEODYN II system description (vols. 1–5), contractor report. SGT Inc. https://earth.gsfc.nasa.gov/geo/data/geodyn-documentation (2017).

56. Konopliv, A. S., Park, R. S. & Folkner, W. M. An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data. Icarus 274, 253–260 (2016).

    Article  ADS  Google Scholar 

57. Millour, E. et al. The Mars Climate Database (version 5.3). Technical report; http://www-mars.lmd.jussieu.fr.

58. Wagner, N. L., James, P. B., Ermakov, A. I. & Sori, M. M. Evaluating the use of seasonal surface displacements and time-variable gravity to constrain the interior of Mars. J. Geophys. Res. Planets 129, e2023JE008053 (2024).

    Article  ADS  Google Scholar 

59. Berne, A., Simons, M., Keane, J. & Park, R. S. Inferring the mean thickness of the outer ice shell of Enceladus from diurnal crustal deformation. J. Geophys. Res. Planets 128, e2022JE007712 (2023).

    Article  ADS  Google Scholar 

60. Berne, A., Simons, M., Keane, J. & Park, R. S. Using tidally-driven elastic strains to infer regional variations in crustal thickness at Enceladus. Geophys. Res. Lett. 50, 311–318 (2023).

    Article  Google Scholar 

61. Spitale, J. N. et al. Curtain-based maps of eruptive activity in Enceladus’s south-polar terrain at 15 Cassini epochs. Planet. Sci. J. 6, 67 (2025).

    Article  Google Scholar 

62. Bagheri, A. et al. Exploring the interior structure and mode of tidal heating in Enceladus. Planet. Sci. J. 6, 245 (2025).

    Article  Google Scholar 

63. Berne, A. C. Tidal Dynamics of Laterally Heterogeneous Planetary Bodies. PhD thesis, California Institute of Technology https://resolver.caltech.edu/CaltechTHESIS:08022025-042237605 (2026).

64. Park, R. S. et al. The global shape, gravity field, and libration of Enceladus. J. Geophys. Res. Planets 129, e2023JE008054 (2024).

    Article  ADS  Google Scholar 

65. Meyer, C. R. et al. Implications of shallow-shell models for topographic relaxation on icy satellites. J. Geophys. Res. Planets 130, e2025JE009247 (2025).

    Article  ADS  Google Scholar 

66. Bruinsma, S. & Lemoine, F. G. A preliminary semiempirical thermosphere model of Mars: DTM-Mars. J. Geophys. Res. Planets 107, 5085 (2002).

    Article  ADS  Google Scholar 

67. Sanchez, B. V., Rowlands, D. D. & Haberle, R. M. Variations of Mars gravitational field based on the NASA/Ames general circulation model. J. Geophys. Res. Planets 111, E06010 (2006).

    Article  ADS  Google Scholar 

68. Madeleine, J.-B., Forget, F., Millour, E., Montabone, L. & Wolff, M. J. Revisiting the radiative impact of dust on Mars using the LMD Global Climate Model. J. Geophys. Res. Planets 116, E11010 (2011).

    Article  ADS  Google Scholar 

69. Genova, A. et al. Long-term variability of CO₂ and O in the Mars upper atmosphere from MRO radio science data. J. Geophys. Res. Planets 120, 849–868 (2015).

    Article  ADS  CAS  Google Scholar 

70. Rovira-Navarro, M. et al. Prospects of using tidal tomography to constrain Ganymede’s interior. Geophys. Res. Lett. 52, e2025GL114708 (2025).

    Article  ADS  Google Scholar 

71. Rovira-Navarro, M. et al. Using tidal tomography to unveil the 3D structure of moons. In Proc. EPSC-DPS Joint Meeting 2025 id. EPSC-DPS2025-49 (Copernicus Meetings, 2025).

72. Patil, A., Huard, D. & Fonnesbeck, C. J. PyMC: Bayesian stochastic modelling in Python. J. Stat. Softw. 35, 1–81 (2010).

    Article  PubMed  PubMed Central  Google Scholar 

73. Yamauchi, H. & Takei, Y. Polycrystal anelasticity at near-solidus temperatures. J. Geophys. Res. Solid Earth 121, 7790–7820 (2016).

    Article  ADS  Google Scholar 

74. McCarthy, C., Takei, Y. & Hiraga, T. Experimental study of attenuation and dispersion over a broad frequency range: 2. The universal scaling of polycrystalline materials. J. Geophys. Res. Solid Earth 116, B09207 (2011).

    ADS  Google Scholar 

75. Nimmo, F., Faul, U. H. & Garnero, E. J. Dissipation at tidal and seismic frequencies in a melt-free Moon. J. Geophys. Res. Planets 117, E09005 (2012).

    Article  ADS  Google Scholar 

76. Nimmo, F. & Faul, U. H. Dissipation at tidal and seismic frequencies in a melt-free, anhydrous Mars. J. Geophys. Res. Planets 118, 2558–2569 (2013).

    Article  ADS  CAS  Google Scholar 

77. Dalton, C. A. & Faul, U. H. The oceanic and cratonic upper mantle: clues from joint interpretation of global velocity and attenuation models. Lithos 120, 160–172 (2010).

    Article  ADS  CAS  Google Scholar 

78. Lau, H. C. P. & Faul, U. H. Anelasticity from seismic to tidal timescales: theory and observations. Earth Planet. Sci. Lett. 508, 18–29 (2019).

    Article  ADS  CAS  Google Scholar 

79. Chung, D. H. Effects of iron/magnesium ratio on P- and S-wave velocities in olivine. J. Geophys. Res. 75, 7353–7361 (1970).

    Article  ADS  CAS  Google Scholar 

80. Chevrel, M. O., Baratoux, D., Hess, K. U. & Dingwell, D. B. Viscous flow behavior of tholeiitic and alkaline Fe-rich martian basalts. Geochim. Cosmochim. Acta 124, 348–365 (2014).

    Article  ADS  CAS  Google Scholar 

81. Berne, A., Matsuyama, I. & Rovira-Navarro, M. Example code for “Tidal tomography reveals a thermal anomaly beneath Mars’s crustal dichotomy”. Zenodo https://doi.org/10.5281/zenodo.20823472 (2026).

Download references

Authors and Affiliations

1. Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ, USA

   Alexander Berne, Isamu Matsuyama & Angela Marusiak

2. Department of Earth, Environmental and Planetary Sciences, Brown University, Providence, RI, USA

   Nicholas Wagner & Harriet C. P. Lau

3. Planetary Geology, Geophysics and Geochemistry Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA

   Sander Goossens

4. Department of Mechanical and Aerospace Engineering, Sapienza University of Rome, Rome, Italy

   Antonio Genova

5. Department of Geoscience and Engineering, Delft University of Technology, Delft, The Netherlands

   Marc Rovira-Navarro

6. Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA

   Amirhossein Bagheri

7. Department of Earth, Planetary, and Space Sciences, University of California, Los Angeles (UCLA), Los Angeles, CA, USA

   Chuan Qin

8. Institute for Geophysics, Jackson School of Geosciences, University of Texas at Austin, Austin, TX, USA

   Douglas Hemingway

9. Institute of Earth Sciences, Academia Sinica, Taipei, Taiwan

   Kar Wai Cheng

10. Department of Physics, University of Colorado Boulder, Boulder, CO, USA

    Shijie Zhong

11. Department of Earth and Planetary Sciences, University of California, Santa Cruz, Santa Cruz, CA, USA

    Francis Nimmo

Contributions

A. Berne conceived and designed this study, under the supervision of I.M. A. Berne drafted the article and constructed Figs. 1–4. N.W. provided expertise on modelling the atmosphere of Mars to produce coefficients presented in Table 1. A. Berne, S.G. and A.G. analysed the spacecraft data to determine a model for the static and time-varying gravity field of Mars. M.R.-N., A. Bagheri and C.Q. developed numerical models used for predicting tidal deformation for laterally heterogeneous bodies in the study. K.W.C. and S.Z. provided expertise on the extent that formation scenarios for the crustal dichotomy can be constrained from results presented in this study. A. Bagheri, F.N., A.M. and H.C.P.L. provided insight into the relationship between partial melt and the dissipative properties of the mantle of Mars. D.H. provided expertise on the relationship between partial melt and magnetism in this study. F.N. provided expertise on the relationship between tectonics, volcanism and the presence of mantle partial melt for Mars. All authors discussed the results of the study and commented on the manuscript at each stage of revision.

Corresponding author

Correspondence to Alexander Berne.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Gabriel Tobie and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions