<figure><img src="image_1.png" /></figure>

Research Paper

# Coupled tidal tomography and thermal constraints for probing Mars viscosity profile

Alex Guinard *, Agnès Fienga, Anthony Mémin, Clément Ganino

Géoazur, CNRS, Observatoire de la Côte d'Azur, Université Côte d'Azur, Valbonne, France

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

**Keywords:**

Geophysics

Mars interior

Tides

solid body

## ABSTRACT

Computing the tidal deformations of Mars, we explored various Mars spherically symmetric internal structures with different types of interface between the mantle and the liquid core. By assessing their compatibility with a diverse set of geophysical observations we show that despite the very short periods of excitation, tidal deformation is very efficient to constrain Mars interior. We calculated densities and thicknesses for Martian lithosphere, mantle, core–mantle boundary layers and core and found them consistent with preexisting results from other methods. We also estimated new viscosities for these layers. We demonstrated that the geodetic records associated with thermal constraints are very sensitive to the presence of a 2-layered interface on the top of the liquid core in deep Martian mantle. This interface is composed by 2 layers of similar densities but very different viscosity and rheology: the layer on the top of the core is liquid (Newtonian, NBL) and the one at the base of the mantle, overlaying the liquid one, is an Andrade layer (ABL) with a viscosity in average 10 orders of magnitude greater than the Newtonian layer. Our results also indicate that the existence of this 2-layered interface significantly impacts the viscosity profiles of the mantle and the lithosphere. More precisely, models including the 2-layered interface do not display significant viscosity contrast between the mantle and the lithosphere, preventing mechanical decoupling between a lithosphere and the mantle immediately below. Such models are in favor of a stagnant lid regime that can be supported by the current absence of an Earth-like plate tectonics on Mars. Finally, in our results, the presence of liquid Newtonian layer at the top of the liquid core is incompatible with the existence of a solid inner core.

## 1. Introduction

The marsquake S1000a, detected by NASA InSight seismometer (Banerdt et al., 2020), represents a significant breakthrough, marking the first observation of a P-wave diffracted (hereafter labeled Pdiff) along the layer in between the lower mantle and the core also called the core–mantle boundary layer (CMB) (Horleston et al., 2022). On one hand, the majority of existing profiles, such as those presented in studies by Drilleau et al. (2022) or Duran et al. (2022), assume a mantle homogeneous in composition. However, these profiles show fits not matching or poorly matching with the differential measurements of arrival times between Pdiff and PP (P-wave reflected once at the planet’s surface). On the other hand, alternative profiles, like the one proposed by Durán et al. (2022), suggest the existence of a heterogeneous mantle. According to Samuel et al. (2022), these profiles require a significant reduction in velocity in the deep mantle to consistently explain the observations. Indeed, recent publications by Samuel et al. (2023) and Khan et al. (2023) highlight a thin layer of molten silicate above the core of Mars. This layer is believed to

result from the solidification of a primitive magma ocean, forming an iron-enriched basal layer containing heat-producing elements (e.g. Elkins-Tanton et al., 2003; Zeff and Williams, 2019).

Tidal forcing exerted on Mars induces a global long wavelength deformation of the planet. Due to the periodic nature of the orbit of planetary bodies, the induced tidal deformation on Mars is also periodic but with a slight delay that witnesses the anelastic tidal dissipation of the planet. As the tidal deformation and dissipation depend on the inner structure of Mars, it is also possible to probe Mars interior from tidal observations. Analysis of orbiter radio tracking have revealed the tidal change of the Martian gravity field (Konopliv et al., 2020). More specifically the non-dimensional tidal Love number $k_2$ relating linearly the degree-2 perturbation of the gravity potential to the external forcing potential, has been estimated. Furthermore, by considering the orbital interactions between Mars and its satellite Phobos, it was possible to detect the tidal dissipation induced by Phobos on Mars and to evaluate the tidal quality factor at the Phobos excitation frequency (Pou et al., 2022). In combining these two informations, it is then possible to

<footnote>* Corresponding author.
E-mail addresses: alex.guinard@geoazur.unice.fr (A. Guinard), agnes.fienga@oca.eu (A. Fienga).</footnote>
quantify the amplitude and the delay in the secular change of the gravity potential.

In this paper, we explore various spherically symmetric internal structures of Mars, by examining profiles that include or exclude a solid inner core. Simultaneously, we aim to detect the presence of a 2-layered interface in between the liquid core and the base of the mantle by assessing their compatibility with a diverse set of geophysical observations and thermal constraints. Section 2 outlines our semi-analytical approach for estimating tidal complex Love numbers for each Mars internal structure. Profiles consistent with observations are identified according to the procedure outlined in Section 2. The results of this statistical selection are presented in Section 3, where we also assess the compatibility of our profiles with existing thermal constraints to propose possible inner structures, the viscosity profiles being the most stringent aspect of this selection. In Section 4, we discuss the implications of these results and how tidal deformation, mass, moment of inertia and thermal constraints can help to decipher the Martian interior.

## 2. Tidal deformation and internal structure

### 2.1. Computing tidal deformation parameters

For a viscoelastic planetary body, the Love number $k_2$ depends on the period of the tidal forcing and the interior structure of the body. We use the software ALMA (plAnetary Love nuMbers cAlcula-tor) initially developed to compute Love numbers for a body subject to Heaviside time-history loading with applications to glacial iso-static adjustment studies (Spada and Boschi, 2006; Spada, 2008). With ALMA³, it has recently been extended to compute periodic tidal Love numbers (Melini et al., 2022). For all ALMA versions, the body is supposed to be spherically symmetric, homogeneous by layers, non-rotating and self-gravitating. This software numerically integrates the gravito-viscoelasticity equations in the frequency domain and subsequently employs a numerical Laplace inversion to retrieve Love numbers in the time domain using the Post-Widder formula (Post, 1930; Widder, 1934, 1941). As described in details in the paper of Briaud et al. (2023a), ALMA computes the Laplace-transformed solution of the equilibrium equations as follows

$$ \vec{x}(R, s) = f(s)[P_x W(s)J][P_b W(s)J]^{-1} \vec{b} \quad (1) $$

with $\vec{x}(R, s) = (u, v, \phi)^t$ where $u, v$ and $\phi$ are respectively the vertical and horizontal components of the displacement and the incremental potential. In Eq. (1), $R$ is the radius of the planet, $s$ is the Laplace variable, $f(s)$ is the Laplace transform of the forcing term, $W(s)$ is the matrix $(6 \times 6)$ which propagates the solution from the core radius to the surface, $P_x$ and $P_b$ are projection operations $(3 \times 6)$, $J$ is a matrix $(6 \times 3)$ which accounts for the boundary conditions between the core and the layer at the interface between the core and the mantle (see Section 2.4), and $\vec{b}$ expressing the tidal boundary conditions at the surface. The propagator $W$ is written as follows:

$$ W(s) = \prod_{j=L+1}^{1} Y_j(r_{j+1}, s) Y_j^{-1}(r_j, s) \quad (2) $$

where the index $i$ decreases from $j = L+1$ to $j = 1$, $r_j$ is the radius of each interface, $L$ is the number of layers, $r_1$ is the radius for the deepest interface and $r_{L+2} = R$. $Y(r, s)$ is the fundamental matrix $(6 \times 6)$ of the system of differential equations describing the radial part of the equilibrium and the Laplace equation (Spada and Boschi, 2006).

If the external forcing has a periodic dependence in time, the solution can be obtained by taking $f(s) = 1$ and by setting $s = i\omega$ in Eq. (1), where $\omega$ is the tidal forcing frequency and $i$ is the imaginary unit. In this case, the solution vector can be written as follows:

$$ \vec{x}(R, i\omega) = [P_x W(i\omega)J][P_b W(i\omega)J]^{-1} \vec{b} \quad (3) $$

The tidal Love numbers can then be obtained from the solution vector with the relation

$$ \begin{pmatrix} u(\omega) \\ v(\omega) \\ \phi(\omega) \end{pmatrix} = \phi_{\text{ext}} \begin{pmatrix} \frac{h(\omega)}{\gamma} \\ \frac{l(\omega)}{\gamma} \\ -(1 + k(\omega)) \end{pmatrix} \quad (4) $$

where $\phi_{\text{ext}}$ is the tidal-raising body potential and $\gamma$ is the surface gravitational acceleration. The Eq. (4) can be rewritten as follows

$$ \begin{pmatrix} h(\omega) \\ l(\omega) \\ k(\omega) \end{pmatrix} = \begin{pmatrix} \xi u \\ \xi v \\ -1 - \frac{\xi}{\gamma} \phi \end{pmatrix} \quad (5) $$

where $\xi = M/R$ is the ratio between mass $M$ and radius $R$.

This non-conventional strategy allows rapid and accurate computation of Love numbers for a spherically symmetric, non-rotating, incompressible and self-gravitating viscoelastic body. As a consequence, ALMA requires as inputs the parameters that describe the multi-layered 1-D rheological profile, namely radius, density, rigidity and viscosity for each layer. It also requires the frequency $\omega = \frac{2\pi}{T_f}$ of the periodic excitation $e^{i\omega t}$, where $t$ is time and $T_f$ is the period of the tidal forcing. In the frequency-domain, $k_2$ is complex. Due to energy dissipation, the response of the planet is delayed. This implies that the time-dependent $k_2$ follows the relation (Melini et al., 2022).

$$ k_2(t) = \|k_2(\omega)\| e^{i(\omega t - \epsilon)} \quad (6) $$

where $\epsilon$ is the phase lag that can be estimated with

$$ \tan(\epsilon) = - \frac{\operatorname{Im}(k_2(\omega))}{\operatorname{Re}(k_2(\omega))} \quad (7) $$

The ratio between the energy loss by dissipation and the total energy of the body is often given in term of quality factor $Q$ which also depends on the excitation frequency $\omega$ such as:

$$ Q(\omega) = - \frac{\|k_2(\omega)\|}{\operatorname{Im}(k_2(\omega))} \quad (8) $$

To compute Love numbers, ALMA also requires the rheological law for each layer that makes up the Martian profile. With the exception of Maxwell's rheology, the results are not very sensitive to the exact nature of the attenuation mechanisms that explain energy dissipation inside the planets, and in particular on Mars (Bagheri et al., 2019). In line with the study of Saliby et al. (2023) on the internal structure of Venus, we use the Andrade rheology whose complex compliance $J(\omega)$ defining the complex rigidity is given by:

$$ J(\omega) = \frac{1}{\mu} + \beta \frac{\Gamma(1 + \alpha)}{(i\omega)^\alpha} - \frac{i}{\eta\omega} \quad (9) $$

with $\beta$ is taken such as $\beta \approx \mu^{\alpha-1}/\eta^\alpha$ (Castillo-Rogez et al., 2011), $\Gamma$ is the Gamma function, $\mu$ is the rigidity, $\eta$ the viscosity and $\alpha$ determines the transient response duration in the primary creep. Accordingly to the work of Bagheri et al. (2019), $\alpha$ values that match with Martian geophysical constraints are in the interval $[0.2 : 0.4]$. Therefore we explored this parameter by considering for each tested profiles a value of $\alpha$ randomly selected from a uniform distribution between 0.2 and 0.4. The distribution of the $\alpha$ values for the selected models is given in Fig. A.1. Finally, in order to test the sensitivity of the models to the solid inner core characteristics, we use the Maxwell rheology as an alternative to a purely elastic body for the inner core (see Appendix B).

### 2.2. General inversion strategy

With the objective to use the tidal deformation as a tool to explore the inner structure of Mars we employed an inversion strategy as described by Briaud et al. (2023a,b) for the Moon and Saliby et al. (2023) for Venus. It consists in varying, using the random walk exploration method, the parameters required to estimate the tidal deformation. First, we use the observational constraints given in Table 1 to create
<table><caption>**Table 1**
Data used to constrain the internal structure of Mars.</caption><thead><tr><th>Data</th><th>Symbol</th><th>Value</th><th>3 &minus; &sigma;</th><th>Reference</th></tr></thead><tbody><tr><td>Mean radius (km)</td><td>R</td><td>3389.5</td><td>&plusmn;0.6</td><td>Seidelmann et al. (2002)</td></tr><tr><td>Total mass (kg)</td><td>M</td><td>6.417 &times; 10<sup>23</sup></td><td>&plusmn;8.943 &times; 10<sup>19</sup></td><td>Konopliv et al. (2016)</td></tr><tr><td>Normalized moment of inertia</td><td>C/MR<sup>2</sup></td><td>0.36379</td><td>&plusmn;0.0003</td><td>Konopliv et al. (2016)</td></tr><tr><td>Tidal Love number</td><td>k<sub>2</sub></td><td>0.174</td><td>&plusmn;0.024</td><td>Konopliv et al. (2020)</td></tr><tr><td>Tidal quality factor</td><td>Q</td><td>93</td><td>&plusmn;25.2</td><td>Pou et al. (2022)</td></tr><tr><td>Period of Phobos tides (h)</td><td>T<sub>s</sub></td><td>5.55</td><td></td><td></td></tr></tbody></table>

a database containing profiles that satisfy the total mass and radius of Mars as well as the moment of inertia (see Table 2). We then use ALMA to estimate the Love number $k_2(\omega)$ and the dissipation $Q(\omega)$ of Mars, at the period of Phobos tidal forcing. We subsequently compare these two estimated parameters to those derived from the observations (Table 1) computing the $\chi^2$ for each model from the database. The results of this selection are presented in Table 3 where we use the $3-\sigma$ of the $\chi^2$ distribution as a criterion for selecting profiles compatible with the observations.

## 2.3. Tidal observations

Geophysical Mars observations used as constraints for our inversion strategy are presented in Table 1 and correspond to: the mean radius R (Seidelmann et al., 2002), the total mass M (Konopliv et al., 2016), the normalized moment of inertia $C/MR^2$ (Konopliv et al., 2016), the tidal Love number $k_2$ (Konopliv et al., 2020) and the quality factor Q (Pou et al., 2022). Besides, Pou et al. (2022) showed that the difference between $k_2$ predicted at the excitation period of Phobos (5.55h) and that estimated from observations at the excitation period of the Sun (12.32h) is significantly lower than the error on the latter. It is also reasonable to assume that the observational uncertainties on these two values are of similar magnitude. Thus, as suggested by Pou et al. (2022), we consider $k_2$ and Q to be that induced by the Phobos excitation period. More about the Mars tidal parameters can be found in the review by Bagheri et al. (2022).

## 2.4. Profile configuration

To calculate the viscoelastic tidal Love numbers $k_2$ using ALMA, we need to set up 1D profiles that describe the hypothetical interior structure of Mars. The input parameters for ALMA are the radius, the density, the rigidity, the viscosity and the rheology for each layer. We consider three different configurations which are all composed of a crust (Cr), a lithosphere (Li), a mantle (Ma) and an uniform liquid core (LC). The first profile also includes a solid layer at the bottom of the mantle called Andrade basal layer (ABL) as we use the Andrade rheology for describing the deformation of this layer. The second profile has a liquid layer above the liquid core called Newton basal layer (NBL) similar to the model proposed by Khan et al. (2023). Finally, the last profile differs from the first two as it has a NBL overlain by a ABL between the mantle and the liquid core, as proposed by Samuel et al. (2023). These two new layers are named after the rheology that defines them (Newton for NBL and Andrade for ABL). We also explore profiles with a solid inner core (SIC) for these three different configurations and the results are presented in Appendix B. Hereafter we describe each layer and Table 2 shows the structure of the 1D profiles that we have considered and the exploration ranges of the parameters assumed for each layer.

**Crust.** The seismic studies of Kim et al. and Li et al. (2022) have allowed to give strong constraints on the Martian crust. Indeed, by combining seismology and gravimetry, Kim et al. managed to estimate the average thickness at 49 ± 7 km of the Martian crust thanks to the largest marsquake (Mw = 4.6) detected by the InSight mission. Li et al. (2022) used the data from the farthest marsquake from the lander

and concluded that the crust would be globally composed of two sub-layers with different propagation speeds rather than three as under the landing site. The choice of density is based on the study of Wieczorek et al. (2022) who evaluated three density profiles among which we will use the values of the crust with a uniform density equal to 2.9 g/cm$^3$. We use the elastic rheology for this layer.

**Lithosphere and mantle.** The input parameters concerning the lithosphere and the mantle are based on the work of Khan et al. (2021) and Drilleau et al. (2022). These studies show that Mars has a much thicker lithosphere than that of the Earth, about 500 km, which implies a reduction in the size of the mantle. However, recent research by Samuel et al. (2023) suggests that the thickness of the lithosphere is sensitive not only to thermal properties, but also to the composition of the deep mantle. The velocity profiles of Drilleau et al. (2022) revealed a slowing down of the S waves in this layer sometimes called LVZ (low velocity zone) in the literature. We consider the density of the lithosphere by exploring a range of values that are similar to the density profile of Plesa et al. (2021). The density of the mantle is deducted from the other layers to conserve the total mass of the planet. To follow the studies carried out by Bagheri et al. (2019), we use the Andrade rheology for the lithosphere and the mantle.

**CMB and liquid core.** Recent geophysical and geodynamic studies have proposed different structures at the interface between the core and the mantle of Mars. First, Samuel et al. (2021) proposed a model of the stratified Martian mantle deduced from geodynamic constraints. They show that the presence of a molten, iron-enriched silicate layer above the core is consistent with geophysical observations (seismic, geodetic, topographic, and gravity data). With this model, Samuel et al. (2023) proposed a more detailed configuration with a silicate layer subdivided into two parts: a partially molten upper layer and a completely molten lower layer. This suggests that these layers reduce the core size to 1650 ± 20 km, with a density of 6.5 g cm$^{-3}$. Subsequently, Khan et al. (2023), using InSight seismic data and other geodynamic and geochemical constraints, proposed an alternative with a fully molten silicate layer above the core, with a thickness of 150 ± 15 km and a density of 4.05 g cm$^{-3}$. This profile would induce a decrease of the core size to 1675 ± 30 km still compatible with Samuel et al. (2023). In comparison with the 1830 ± 30 km estimated in other studies such as Stähler et al. (2021), Samuel et al. (2023) suggest an increase to the density of the liquid core. We decide to explore these two types of profiles by choosing a wide range of density and viscosity values for the ABL, the NBL and the liquid core in order to explore as many structures as possible.

**Outer and inner core.** This paragraph only concerns the case of a model with a solid internal core, i.e. a model in which the core is subdivided in two parts: a liquid outer core and a solid inner core. The work of Hemingway and Driscoll (2021) aims at showing that the absence of a magnetic field on Mars does not necessarily imply the absence of a solid inner core. Indeed, thanks to simulations of thermal evolution of the planet, the study was able to demonstrate that the presence of a solid inner core is possible and this implies that the Martian dynamo will reactivate in the future as it has previously been able to reactivate for a short period according to the work of Lillis et al. (2005). Their simulation shows that the solid inner core could be about 400 km in radius and reach 7.7 g/cm$^{-3}$ in density. Moreover, recent seismic
<caption>**Table 2**</caption>

General parameters of 1D profiles for the interior of Mars. Values in square brackets [ ] indicate the range of parameters that vary randomly and uniformly between profiles. The stars * represent the values that are inferred from the other layers to preserve the size and total mass of Mars. For the α parameter of the Andrade rheology, we use a random value in the uniform distribution between 0.2 and 0.4 in agreement with the results of Bagheri et al. (2019). The distribution of the α parameter for the selected profiles is given in Appendix A (Fig. A.1).

<table><thead><tr><th>Layer</th><th></th><th>Crust</th><th>Lithosphere</th><th>Mantle</th><th>Andrade<br>Basal Layer</th><th>Newton<br>Basal Layer</th><th>Liquid<br>Core</th><th>Solid Inner<br>Core</th></tr><tr><th></th><th></th><th>Cr</th><th>Li</th><th>Ma</th><th>ABL</th><th>NBL</th><th>LC</th><th>SIC</th></tr></thead><tbody><tr><td>Radius</td><td>km</td><td>3,389.5</td><td>[3,319: 3,362]</td><td>*</td><td>a</td><td>b</td><td>[1,500: 2,000]</td><td>[0: 800]</td></tr><tr><td>Density</td><td>g/cm<sup>3</sup></td><td>[2.7: 3.1]</td><td>[3: 3.8]</td><td>*</td><td>a</td><td>b</td><td>[4: 10]</td><td>[4: 10]</td></tr><tr><td>Rigidity</td><td>GPa</td><td>[20: 25]</td><td>[20: 70]</td><td>[20: 130]</td><td>[20:180]</td><td>–</td><td>–</td><td>–</td></tr><tr><td>Viscosity</td><td>log<sub>10</sub>(Pa s)</td><td>–</td><td>[14: 30]</td><td>[18: 30]</td><td>[1: 30]</td><td>[1: 30]</td><td>[1: 30]</td><td>[20: 30]<sup>c</sup></td></tr><tr><td>Rheology</td><td></td><td>Elastic</td><td>Andrade</td><td>Andrade</td><td>Andrade</td><td>Newton</td><td>Newton</td><td>Elastic/Maxwell</td></tr></tbody></table>

a The radius of the ABL is a random value included in the last 100 kilometers of the mantle. Its density is a value randomly sampled between the density of the mantle and the density of the lower layer.

b The radius of the NBL is the radius of the LC, and the radius of the LC is randomly chosen in the first 200 kilometers. The density of the NBL is a value randomly sampled between the density of the ABL and the LC.

c The range of viscosity is used for the Maxwell case only.

studies by Irving et al. (2023) show that if a solid inner core exists then its radius does not exceed 750 km. We choose to explore a region around the values given by Khan et al. (2023) in the previous paragraph for the radius and the density of the outer core. For the inner core, we impose a maximum radius of 750 ± 50 km to follow the condition given by Irving et al. (2023). We use exploration intervals wide enough not to miss profiles compatible with the observations. We assume that the rheology of the inner core is either elastic or Maxwell and that of the outer core is Newtonian.

## 2.5. Viscosity profile sensitivity

The role of the viscosity of each layer is of primary importance when assessing the tidal dissipation and the tidal response delay of a planet. Hence we investigate the contribution of each layer in dissipating energy through tidal deformation and in bringing constraints on the viscosity structure of Mars. In order to assess the sensitivity of the tidal deformation (considering the imaginary and real parts of the $k_2$ Love number estimated with ALMA) to the viscoelastic behavior of the planet even for very rapid (few hours) excitation periods, we considered different possible viscosities for the lithosphere, the mantle, the ABL, the NBL and the fluid core, starting from possible 1D profiles (one per configuration) selected randomly among those meeting the geophysical constraints given in Table 1. We then plot the ratio $||k_2(\omega)||/Q$ and the value Re($k_2(\omega)$) computed for periods ranging from $10^{-4}$ to $10^8$ hours using the reference model (in black in Figs. 1–3). Fig. 1 present results for 1D profile with ABL only, Fig. 2 with NBL only and Fig. 3 with both ABL and NBL. We obtain that at the tidal excitation period of Phobos, Re($k_2(\omega)$) cannot discriminate the viscosity of the lithosphere and the mantle regardless the configuration of the profile (Figs. 1-B:D, 2-B:D and 3-B:D) unlike that of the ABL which covers an interval between $10^{15}$ and $10^{27}$ Pa s (Figs. 1-F and 3-F). Additionally, $||k_2(\omega)||/Q$ varies significantly depending on the viscosity of the mantle and of the ABL while the lithosphere seems to introduce more sensitivity to the viscosity for the profile with NBL only (Fig. 2-A) in comparison with the other two configurations (Figs. 1-A and 3-A). The presence of the NBL in the profiles drastically changes the sensitivity of the viscosity on the liquid core, in particular for models with both ABL and NBL. Figs. 1-G:H, 2-E:H and 3-G:J show indeed that the NBL behaves like the liquid core (in the case without NBL). The superposition in depth of two Newtonian layers makes the deepest layer less sensitive (or not at all) to viscosity variations, thus leading to a loss of information absorbed by the upper layer. Finally it is interesting to note that at 1 Hz ($10^{-4}$ hours), the $||k_2(\omega)||/Q$ values proposed in Figs. 1–3 are consistent with those estimated or expected with seismic measurements (i.e. Samuel

et al., 2023). Let us note that we also have made sensitivity tests for models with an elastic or a Maxwell solid inner core and models with 3 sublayers in the lithosphere. The results are given respectively in Appendices B and C.

# 3. Results

## 3.1. Geophysical constraints

Table 3 and Table B.1 give respectively the main characteristics of the profiles without and with SIC selected according to the mass, moment of inertia and tidal deformations (see Section 2). We also plotted their distributions in viscosity, density and thickness on Fig. 4 for profiles without SIC and Fig. B.1 for profiles with SIC.

We first present here the main differences between the results for the profiles with ABL-only or with NBL-only and profiles with both ABL and NBL. Overall, the radii and density parameters remain stable, regardless of the choice between the three configurations. As one can see on Fig. 5, the radii and densities obtained for the crust and the mantle are in agreement with literature (e.g. Khan et al., 2023; Samuel et al., 2023).

More generally, comparing intervals of exploration given in Table 2 with those of the selected models in Tables 3 and B.1, one can see that for both profiles without and with SIC, the rigidity is not well constrained as they present uniform distributions which sweep the entire exploration interval. On the contrary, on Fig. 4, we can see how the geophysical constraints are very useful for reducing the intervals of possible thicknesses and densities. This is particularly true for the LC density and thickness and for the density of the NBL. For the profiles with SIC (see Fig. B.1) the constraints are less stringent. Because we are considering viscoelastic rheologies of Mars and we use the quality factor Q obtained by Pou et al. (2022) as part of the profile selection, the viscosities of the different layers are also an interesting outcome of the selection as it can be seen in Table 3 and Figs. 4 and B.1. It is clearly visible for the viscosity of the profiles with NBL that the interval of possible viscosity is reduced by a factor 2 after the selection, ending up with a viscosity of the NBL systematically smaller than the one of the mantle and the one of the lithosphere. This is not the case for profiles with both ABL and NBL since 1 model out of 2 has a viscosity of the ABL lower than that of the NBL.

Finally, we can note that the presence of the NBL in the models implies a loss of information on the viscosity of the LC which is very little constrained if we compare to the case without NBL (ABL-only). Note that this loss of information is not visible when we consider only the median and extreme values of LC viscosity (Table 3), but becomes
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Behavior of the $||k_2(\omega)||/Q$ ratio and of the real part of $k_2$ as a function of the period for a profile with ABL only. Each curve corresponds to a different lithosphere (A-B), mantle (C-D), ABL (E-F), or core (H-I) viscosity. In black, the reference model considered and the red cross, the value and the error bars of $k_2/Q$ and of $\Re(k_2)$ corresponding to the period of excitation induced by Phobos at 5.55 h. The error bars in x-axis correspond to the values of Pou et al. (2022). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

obvious when we look at the entire distribution represented at Fig. 4 (3rd row). As it was predicted by the sensitivity analysis of Section 2.5, the NBL adopts the viscous signature of the LC, in agreement with the low sensitivity of the $||k_2(\omega)||/Q$ ratio to LC viscosity as shown in Figs. 1, 2 and 3. Most of the previous comments are also true for the profiles with SIC (see Fig. B.1).

## 3.2. Thermal constraints

### 3.2.1. Temperature estimation

We estimated the temperatures at the top of the mantle, at the top of the ABL and at the top of the NBL by considering the viscosity contrasts between the lithosphere and the mantle, between the mantle and the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Same as Fig. 1 for a profile with NBL only. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Same as Fig. 1 for a profile with both ABL and NBL. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4. Distribution in thickness, density and viscosity for profiles with ABL (1st column), with NBL (2nd column) and with ABL and NBL (3rd column) without solid inner core presented in Table 3. The gray bars in the background correspond to the exploration intervals of each parameter given in Table 2 before the $\chi^2$ filtering detailed in Section 2.</figcaption></figure>

<caption>Table 3</caption>

Radii at the top of the layer, densities and viscosities of profiles selected with the method described in Section 2. The values presented in this table represent the 0.5, 0.995 and 0.005 quantiles. Crust (Cr), Lithosphere (Li), Mantle (Ma), Andrade basal layer (ABL), Newton basal layer (NBL), Liquid core (LC)

<table><caption>Table 4</caption><thead><tr><th>Number of profile</th><th></th><th>Cr</th><th>Li</th><th>Ma</th><th>ABL</th><th>NBL</th><th>LC</th></tr></thead><tbody><tr><td rowspan="3">3,514</td><td colspan="8" style="text-align: center; font-weight: bold;">With ABL</td></tr><tr><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3,362</sup><sub>3,319</sub></td><td>2,820<sup>3,334</sup><sub>2,503</sub></td><td>1,698<sup>1,884</sup><sub>1,516</sub></td><td>-</td><td>1,645<sup>1,798</sup><sub>1,501</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3,110</sup><sub>2,70</sub></td><td>3.28<sup>3,54</sup><sub>3,00</sub></td><td>3.79<sup>4,98</sup><sub>3,40</sub></td><td>5.17<sup>7,68</sup><sub>3,64</sub></td><td>-</td><td>6.66<sup>8,26</sup><sub>2,01</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With NBL</td></tr><tr><td rowspan="3">1,792</td><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3,362</sup><sub>3,319</sub></td><td>2,890<sup>3,342</sup><sub>2,506</sub></td><td>-</td><td>1,692<sup>1,997</sup><sub>1,500</sub></td><td>1,579<sup>1,952</sup><sub>1,329</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3,110</sup><sub>2,70</sub></td><td>3.30<sup>3,55</sup><sub>3,01</sub></td><td>3.71<sup>4,66</sup><sub>3,47</sub></td><td>-</td><td>5.13<sup>7,67</sup><sub>3,62</sub></td><td>6.86<sup>9,30</sup><sub>3,33</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>-</td><td>20.75<sup>29,88</sup><sub>14,20</sub></td><td>21.63<sup>29,83</sup><sub>18,61</sub></td><td>-</td><td>6.96<sup>14,74</sup><sub>1,05</sub></td><td>12.78<sup>29,76</sup><sub>1,13</sub></td></tr><tr><td rowspan="3">3,471</td><td colspan="8" style="text-align: center; font-weight: bold;">With ABL and NBL</td></tr><tr><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3,362</sup><sub>3,319</sub></td><td>2,856<sup>3,340</sup><sub>2,505</sub></td><td>1,767<sup>2,065</sup><sub>1,525</sub></td><td>1,717<sup>1,996</sup><sub>1,502</sub></td><td>1,613<sup>1,960</sup><sub>1,329</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3,110</sup><sub>2,70</sub></td><td>3.29<sup>3,54</sup><sub>3,01</sub></td><td>3.76<sup>4,83</sup><sub>3,42</sub></td><td>4.49<sup>6,90</sup><sub>3,56</sub></td><td>5.21<sup>7,60</sup><sub>4,02</sub></td><td>6.54<sup>8,90</sup><sub>5,18</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Viscosity [log<sub>10</sub>(Pa s)]</td></tr><tr><td></td><td></td><td>-</td><td>20.87<sup>29,92</sup><sub>14,23</sub></td><td>22.08<sup>29,91</sup><sub>18,54</sub></td><td>9.52<sup>29,61</sup><sub>1,07</sub></td><td>9.23<sup>29,41</sup><sub>1,07</sub></td><td>13.70<sup>29,86</sup><sub>1,13</sub></td></tr></tbody></table>

Initialization values extracted from Samuel et al. (2023). The temperature of the top of the lithosphere is the average temperature taken at the bottom of the crust, for a thickness of the crust varying from 28 and 70 km.

<table><thead><tr><th>Profile</th><th>T lithosphere (K)</th><th>E* (kJ mol<sup>-1</sup>)</th><th>V* (cm<sup>3</sup> mol<sup>-1</sup>)</th><th>Distribution</th></tr></thead><tbody><tr><td>With ABL</td><td>750</td><td>320 ± 70</td><td>4.6 ± 2</td><td>HNormal</td></tr><tr><td>With NBL</td><td>750</td><td>[90 : 390]</td><td>[2.6 : 9.3]</td><td>Uniform</td></tr><tr><td>With ABL and NBL</td><td>750</td><td>120 ± 30</td><td>7.3 ± 2</td><td>HNormal</td></tr></tbody></table>

ABL or between the mantle and the NBL, and between the ABL and NBL (for the models with both ABL and NBL) with the hypothesis that all the layers are homogeneous. We use the equation for a convective mantle given by Nakada et al. (2012) such as:

$$ \eta(r) = \eta_L \exp \left[ - \frac{H^*}{R_g} \left( \frac{1}{T_L} - \frac{1}{T(r)} \right) \right] \qquad (10) $$

where $T_L$, $P_L$ and $\eta_L$ are respectively the upper layer temperature, pressure and viscosity values, $T(r)$ and $\eta(r)$ are respectively the temperature and viscosity of the layer considered. $R_g$ is the gas constant and $H^*$ is the activation enthalpy. Because the activation enthalpy ($H^*$) depends on the pressure ($P$) and varies with depth, we can rewrite Eq. (10) as a function of the activation energy of the upper layer ($H_L^*$) and the layer
<figure><img src="image_7.png" /><figcaption>Fig. 5. Density profiles from the profiles of the Table 3 following the geophysical constraints. Each point is the median of the density distribution of the layer considered and the error bars are given by the quantiles 0.995 and 0.005. The panels A and C represent densities of the profile with ABL-only and with ABL and NBL compared to Samuel et al. (2023). The panel B represents densities of the profile with NBL-only compared to Khan et al. (2023).</figcaption></figure>

considered ($H^*(r)$)

$$ \eta(r) = \eta_L \exp \left[ - \frac{1}{R_g} \left( \frac{H_L^*}{T_L} - \frac{H^*(r)}{T(r)} \right) \right] \quad (11) $$

or, equivalently

$$ \frac{T(r)}{T_L} = \frac{H^*(r)}{H_L^* + \ln\left(\frac{\eta(r)}{\eta_L}\right) R_g T_L} \quad (12) $$

Here, we set $H^* = E^* + PV^*$ with $E^*$ the activation energy and $V^*$ the activation volume given in the Table 4. Then Eq. (12) takes the following form

$$ \frac{T(r)}{T_L} = \frac{E^* + P(r)V^*}{E^* + P_L V^* + \ln\left(\frac{\eta(r)}{\eta_L}\right) R_g T_L} \quad (13) $$

The temperature at the top of the lithosphere (750 K) is deduced from Samuel et al. (2023) as the average temperature at the bottom of the crust with a thickness between 28 and 70 km. With Eq. (13), the viscosities presented in Table 3 and deduced from the method described in Section 2 and pressure estimated with Rivoldini et al. (2011), we computed the temperatures at the top of the considered layers (top mantle, top ABL, top NBL).

Following Samuel et al. (2023), we have used three sets of activation energy and volume: one set for models including only ABL, one set for models with both ABL and NBL and one set for NBL-only models (see Table 4).

Finally, in order to estimate the impact of the enthalpy uncertainty, we vary the activation energy and volume using normal distributions for both the ABL-only profiles and the ABL and NBL profiles and a uniform distribution for the ABL profiles based on uncertainties given by Samuel et al. (2023) and reported in Table 4. The intervals of exploration chosen for the ABL-only models follow those proposed for profiles without Basal Mantle Layer (BML) in Samuel et al. (2023). The $E^*$ and $V^*$ intervals for models with both ABL and NBL correspond to those from Samuel et al. (2023) with BML. Finally, for NBL-only models, we use large uniform intervals encompassing the former two. By design, these profiles match with the geophysical constraints described in Section 2.

### 3.2.2. Solidus and liquidus

At the CMB, several approaches are possible that are defined according to their thermal behavior. We used the solidus and liquidus profiles extracted from Samuel et al. (2023). These profiles was estimated

by Ruedas and Breuer (2017) for the solidus and by Duncan et al. (2018) for the liquidus for a homogeneous iron-rich mantle and have been modified in Samuel et al. (2021) for a heterogeneous one. We selected the models using the limits imposed by these profiles in the following manner. For ABL-only models, the profiles are selected in such a way that the temperature of the ABL is always smaller than the solidus. For models with only NBL, the selected profiles have the temperature of the NBL always greater than the liquidus. Finally, for models with both ABL and NBL, the profiles are such that the temperature of the NBL is greater than the one of the liquidus but also that the temperature of the ABL is greater than the solidus and smaller than the liquidus. In doing so, we selected models with gradients of temperatures at the interface between the mantle and the liquid core consistent with a partially melt layer.

### 3.2.3. Temperature inversion at the mantle-lithosphere boundary

An interesting feature of the temperature profiles presented in Fig. 6 is the temperatures at the top of the mantle. These temperatures have been obtained considering a lithosphere at 750 K as described in Section 3.2.1. However, some of the produced profiles have a mantle temperature below the lithosphere temperature. Such an inversion of temperature is difficult to interpret because the evolution of the internal temperature with depth in Mars interior is the outcome of the competition between, on one hand, residual heat and the heat produced by the radioactive decay of certain elements, primarily within the crust and mantle, and on the other hand, thermal dissipation through convection and conduction processes. Heat loss occurs through the surface, and as a result, the temperature within Mars necessarily increases with depth. In consequence, we decide to filter the profiles to those with a temperature of the mantle greater than the one of the lithosphere. The quantiles of the remaining profiles are represented with the red boxes of Fig. 6 and are given in Table 5.

Hereafter, we will employ the term “thermal filter”, which considers the solidus-liquidus constraints by Duncan et al. (2018) and detailed in Section 3.2.2, as well as the absence of temperature inversion between the mantle and lithosphere, as explained in Section 3.2.3. The distributions in thickness, density and viscosity for the profile selection considering the thermal filter are presented in Fig. 7 for models without SIC and in Appendix B for models with SIC (Fig. B.1).
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Temperature in Kelvin versus radius in km profiles for the models with ABL-only (A), with NBL-only (B) and both ABL and NBL (C) without inner core. Comparison with the temperature profiles of Samuel et al. (2023). The red dots give the median, 0.995 and 0.005 quantiles in accounting the thermal constraints. The red shaded area show the corresponding uncertainty zones. The black dashed lines indicate the limit of crystallization for the liquid core as given by Helffrich (2017). The orange dashed lines materialized the solidus (light orange) and liquidus (dark orange) both extracted from Samuel et al. (2023). The plain lines give the temperature versus depth profiles extracted from Samuel et al. (2023) for profiles with a solid basal layer (blue plain line in panel A) and for profiles with a liquid layer in contact with the core (black plain line in panels B and C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

**4. Discussion**

The use of thermal filters described in Section 3.2 induces in general a more detailed description of the distributions for all the parameters of the profiles, but with no systematic reduction of the possible intervals as one can see on Fig. 7 and Table 5. Some improvements and trends are however noticeable.

First, thanks to the limits imposed by the solidus and liquidus, it appears clearly that the existence of a solid inner core is only possible with ABL-only configuration.

Second, the use of the thermal filters reduces the interval of possible viscosities for the lithosphere for the 3 configurations. Furthermore, for the models with both ABL and NBL, the mean values of the viscosity for the lithosphere and the mantle are quite close, favoring the hypothesis of a stagnant lid regime, with no tectonics.

Third, in the models with both ABL and NBL, the contrast of viscosity between the ABL and the NBL is highly enhanced with a strong shift of ABL viscosity towards higher values when one accounts for the thermal filters (from $10^9$ to $10^{18}$ Pa s) and a more reduced interval of possible values (from 30 orders of magnitude to 14).

Finally, the addition of the thermal filter maintains the profiles in good agreement with those from Khan et al. (2023) and Samuel et al. (2023) (see Fig. A.3).

**4.1. The presence of a Newton basal layer prevents the existence of a solid core**

As seen in Section 3.2.2, on Fig. 6, are plotted the temperatures versus depth profiles estimated following the procedures described in Section 3.2. The red cross indicates the 0.5, 0.05 and 99.5 quantiles in temperatures and radii as given in Table 5 considering both the solidus-liquidus limit at the CMB and the constraint of the non inversion of the temperature gradient between the mantle and the lithosphere.

Helffrich (2017), proposed that the temperature of the CMB ($T_{CMB}$) is an indication of the existence of a liquid or solid core. For $T_{CMB}$ between 1300 K and 1450 K, it is indeed possible to have models with both a liquid core and a solid inner core whereas for $T_{CMB} > 1450$ K the core is entirely liquid. Consequently, only profiles with ABL-only obtained after filtering on the thermal constraints given in Section 3.2 follow the temperature conditions favoring the presence of a solid inner core (Fig. 6). On the other hand, profiles with a liquid layer in contact with LC (NBL-only and with both ABL and NBL) that match the solidus-liquidus thermal constraints, have temperatures of the CMB (by filtering greater than 1900 K) too high ($T_{CMB} > 1450$ K) for the existence of a solid inner core. We can then reject profiles with NBL only and profiles with both ABL and NBL with SIC and we detail the
<figure><img src="image_9.png" /><figcaption>Fig. 7. Comparisons of thickness, density and viscosity distributions deduced from geophysical and thermal constraints.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Comparisons of thickness, density and viscosity distributions deduced from geophysical and thermal constraints, including also the constraint imposed by an earth-like liquid core.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Cumulative histograms of viscosities for profiles with ABL (A,D), profiles with NBL (B,E) and profiles with ABL and NBL (C,F). (Top row A-C) Profiles without SIC meeting both geophysical and thermal constraints. (Bottom row D-F) Profiles without SIC meeting both geophysical and thermal constraints and with earth-like viscosity of the liquid core (smaller than $10^{10}$ Pa s). The lithosphere histogram is given in green, the mantle in blue, the ABL in yellow, the NBL in red and the LC in purple. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

results for profiles with ABL-only and with SIC obtained after thermal constraints in Table B.1.

## 4.2. Lithosphere and mantle viscosity contrast driven by both ABL and NBL models and possible link with stagnant lid regime

The use of the thermal filters allows the reduction of possible intervals for the viscosity of the lithosphere from about 16 orders of magnitude to 12 (see Table 5 and Fig. 7). But more importantly, it is noteworthy that in models with ABL only and with NBL only, the viscosity of the lithosphere is, for most of the models, above that of the mantle, as it is the case for the Earth where the viscosity of the lithosphere is about 2 orders of magnitude larger than the viscosity of the asthenosphere (Höink et al., 2011). But it is not the case for models with ABL and NBL where the viscosities of the lithosphere and mantle are almost equal. Indeed, as one can see on Fig. 9 presenting cumulative histograms of viscosities for the different layers and configurations, for a given quantile of the models (for example, 75%), the viscosity thresholds are significantly different for the mantle and the lithosphere for both ABL-only and NBL-only models. In our example, 75% of the ABL-only or NBL-only models have mantle viscosities smaller than $10^{20}$ Pa s and lithosphere viscosities smaller than $10^{25}$ Pa s. But for models with both ABL and NBL, the viscosity threshold at 75%, corresponding to $10^{20}$ Pa s, is the same for the lithosphere and the mantle.

This can also be seen on Fig. 10 showing histograms of viscosity ratios. In particular, Fig. 10-C shows the histograms of the log of the viscosity ratios between mantle and lithosphere. On these histograms, one

can see that in the case of both ABL and NBL models, the distribution of the $\eta_{Ma}/\eta_{Li}$ ratio is centered around 0 (equal viscosities). However, with ABL-only and NBL-only configuration, there is a clear offset of the distributions towards negative values (lithosphere viscosity being larger than the one of the mantle). So, the models with both ABL and NBL do not display significant viscosity contrast, preventing mechanical decoupling between a lithosphere and the mantle immediately below. Such models are in favor of a stagnant lid regime that can be supported by the current absence of a Earth-like plate tectonics on Mars. In such cases, the absence of viscosity contrast can be explained by the dryness of the Martian mantle. The deformation of the major mineral of terrestrial planet mantles, olivine, largely depends, indeed, on the temperature, and water content (Karato et al., 2008). If one supposes that the geotherm is quite similar between Earth and Mars, the significant difference in the viscosity of their mantle may indicate the extreme dryness of Martian mantle, as already proposed based on the chemical composition of Martian meteorites (Wanke and Dreibus, 1994). In the opposite, the models with higher viscosity contrasts (ABL-only and NBL-only) could be compatible with an Earth-like tectonics. But as Mars is clearly in stagnant-lid regime (Stern et al., 2018), this type of tectonic regime would not be related to the viscosity profile but to other reasons as for instance the limited tectonic forces at plate boundaries (Zhou et al., 2022). Note that in ABL and NBL configuration, the histograms present some gaps which are due to the thermal constraints imposed on the ABL and NBL layers which are stronger for this configuration than for those in ABL-only and NBL-only.
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Histograms of the log10 of the viscosity ratios between layers in contact for the three configurations of profiles: with ABL-only, with NBL-only and with both ABL and NBL. Figs. A:D were obtained after geophysical and thermal filtering and Figs. E:H correspond to models accounting also the LC constraint (LCC) imposed on the LC viscosity (see Section 4.4). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

In any case, an important conclusion is that the existence of a liquid basal layer below an Andrade layer at the base of the mantle can provide decisive information on the viscosity profiles for upper layers, and specifically for the lithosphere. The result remains true even when one considers multiple sublayers within the lithosphere (see Appendix C).

### 4.3. Viscosity profiles at the base of the mantle

Different patterns in the viscosity profiles can be identified for the ABL and NBL layers on Figs. 9 and 10.

First, as one can see on Fig. 10-D, for models with both ABL and NBL layers, the NBL (in contact with the liquid core) has a viscosity systematically smaller than the ABL one, by 10 orders of magnitude in average. This contrast is high but is compatible with the solidus and liquidus thresholds discussed in Section 3.2.2.

It can be explained by the wide range of possible activation values taken for computing the temperature contrast (Eq. (13) and Table 4). In comparison, on Figs. 6-B:C, one can see that for an ABL with a

temperature of about 2500 K, the contrast of temperature with the NBL deduced from Samuel et al. (2023) is of about 500 K. It is in between 50 K and 2000 K considering the minimum and the maximum of the temperature errorbars for ABL and NBL as estimated in this work.

Second, for models with a liquid layer at the CMB (NBL-only or both ABL and NBL), the viscosity of the layer in contact with the mantle is systematically lower than the one of the mantle (see Fig. 10-A). This is not the case for models with only an Andrade layer at the base of the mantle. As it is visible on Fig. 10-A, two trends are then possible: one with $\eta_{ABL} > \eta_{mantle}$ (with an average contrast of about 3 orders of magnitude) and with $\eta_{mantle} > \eta_{ABL}$ (with an average contrast of about 7 orders of magnitude). In comparison with the (Samuel et al., 2023) temperature profile on Fig. 6-A, it seems that the seismological profiles favor reduced temperature and consequently viscosity contrasts with hotter mantle temperatures.

Third, in comparison with the LC viscosity (Fig. 10-B), the ABL-only models present, in most of cases, viscosities greater than the one of the LC, which is what is naturally expected. However, in the case where a liquid layer stands on the top of the LC (NBL-only or both
<caption>**Table 5**</caption>

Radii (at the top of the layer), densities and viscosities of the profiles presented in Table 3 for which the temperatures for the top of the mantle, ABL and NBL were calculated for 300 different values of E* and V* determined randomly in the ranges given in Table 4. The values presented in this table represent the quantiles 0.5, 0.995 and 0.005 and follow both the geophysical (Section 2) and the thermal (Section 3.2) constraints (top) and geophysical, thermal and LC viscosity (Section 4.4) constraints (bottom).

<table><thead><tr><th>Number of profiles</th><th></th><th>Cr</th><th>Li</th><th>Ma</th><th>ABL</th><th>NBL</th><th>LC</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">With geophysical and thermal constraints</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With ABL</td></tr><tr><td rowspan="3">95,811</td><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3.362</sup><sub>3.319</sub></td><td>2,812<sup>3.331</sup><sub>2.503</sub></td><td>1,699<sup>1.878</sup><sub>1.518</sub></td><td>–</td><td>1,646<sup>1.798</sup><sub>1.501</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3.10</sup><sub>2.70</sub></td><td>3.28<sup>3.55</sup><sub>3.01</sub></td><td>3.78<sup>5.00</sup><sub>3.38</sub></td><td>5.16<sup>7.76</sup><sub>5.60</sub></td><td>–</td><td>6.66<sup>2.53</sup><sub>3.01</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>24.14<sup>29.92</sup><sub>17.86</sub></td><td>20.42<sup>29.23</sup><sub>18.46</sub></td><td>18.72<sup>29.72</sup><sub>2.68</sub></td><td>–</td><td>8.13<sup>28.48</sup><sub>1.07</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With NBL</td></tr><tr><td rowspan="3">13,574</td><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3.362</sup><sub>3.319</sub></td><td>2,890<sup>3.342</sup><sub>2.506</sub></td><td>–</td><td>1,692<sup>1.998</sup><sub>1.500</sub></td><td>1,579<sup>1.952</sup><sub>1.329</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3.10</sup><sub>2.70</sub></td><td>3.30<sup>3.55</sup><sub>3.01</sub></td><td>3.71<sup>4.66</sup><sub>3.47</sub></td><td>–</td><td>5.13<sup>7.67</sup><sub>3.62</sub></td><td>6.81<sup>9.31</sup><sub>5.33</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>22.94<sup>29.75</sup><sub>17.49</sub></td><td>20.43<sup>29.49</sup><sub>18.5</sub></td><td>–</td><td>8.65<sup>14.72</sup><sub>1.16</sub></td><td>12.05<sup>29.74</sup><sub>1.12</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With ABL and NBL</td></tr><tr><td rowspan="3">20,000</td><td>Radius [km]</td><td>3,389.5</td><td>3,345<sup>3.361</sup><sub>3.319</sub></td><td>2,824<sup>3.320</sup><sub>2.557</sub></td><td>1,745<sup>2.029</sup><sub>1.528</sub></td><td>1,706<sup>1.988</sup><sub>1.504</sub></td><td>1,580<sup>1.881</sup><sub>1.327</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.85<sup>3.09</sup><sub>2.70</sub></td><td>3.28<sup>3.52</sup><sub>3.02</sub></td><td>3.94<sup>4.44</sup><sub>3.47</sub></td><td>4.65<sup>7.02</sup><sub>3.69</sub></td><td>5.39<sup>7.60</sup><sub>4.21</sub></td><td>6.56<sup>8.89</sup><sub>5.52</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>21.61<sup>28.58</sup><sub>17.70</sub></td><td>21.64<sup>28.44</sup><sub>19.15</sub></td><td>18.25<sup>25.07</sup><sub>12.78</sub></td><td>8.43<sup>13.78</sup><sub>1.21</sub></td><td>12.28<sup>29.77</sup><sub>1.19</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With geophysical, thermal and LC viscosity constraints</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With ABL</td></tr><tr><td rowspan="3">63,024</td><td>Radius [km]</td><td>3,389.5</td><td>3,340<sup>3.362</sup><sub>3.319</sub></td><td>2,799<sup>3.338</sup><sub>2.503</sub></td><td>1,700<sup>1.877</sup><sub>1.519</sub></td><td>–</td><td>1,647<sup>1.797</sup><sub>1.500</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3.10</sup><sub>2.70</sub></td><td>3.27<sup>3.55</sup><sub>3.01</sub></td><td>3.79<sup>4.99</sup><sub>3.38</sub></td><td>5.15<sup>7.80</sup><sub>3.64</sub></td><td>–</td><td>6.66<sup>8.28</sup><sub>5.05</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>24.12<sup>29.94</sup><sub>17.97</sub></td><td>20.29<sup>28.69</sup><sub>18.40</sub></td><td>18.77<sup>29.72</sup><sub>4.38</sub></td><td>–</td><td>5.49<sup>9.98</sup><sub>1.04</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With NBL</td></tr><tr><td rowspan="3">5,424</td><td>Radius [km]</td><td>3,389.5</td><td>3,341<sup>3.362</sup><sub>3.319</sub></td><td>2,874<sup>3.346</sup><sub>2.504</sub></td><td>–</td><td>1,703<sup>1.998</sup><sub>1.501</sub></td><td>1,602<sup>1.952</sup><sub>1.332</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.92<sup>3.10</sup><sub>2.70</sub></td><td>3.30<sup>3.56</sup><sub>3.01</sub></td><td>3.73<sup>4.54</sup><sub>3.45</sub></td><td>–</td><td>5.05<sup>7.88</sup><sub>3.61</sub></td><td>6.75<sup>9.31</sup><sub>5.43</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>23.04<sup>29.83</sup><sub>17.52</sub></td><td>20.20<sup>28.97</sup><sub>18.55</sub></td><td>–</td><td>9.97<sup>14.41</sup><sub>1.13</sub></td><td>5.35<sup>9.92</sup><sub>1.06</sub></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">With ABL and NBL</td></tr><tr><td rowspan="3">7,710</td><td>Radius [km]</td><td>3,389.5</td><td>3,346<sup>3.361</sup><sub>3.321</sub></td><td>2,911<sup>3.285</sup><sub>2.537</sub></td><td>1,745<sup>1.983</sup><sub>1.545</sub></td><td>1,702<sup>1.959</sup><sub>1.520</sub></td><td>1,564<sup>1.833</sup><sub>1.327</sub></td></tr><tr><td>Density [g cm<sup>-3</sup>]</td><td>2.94<sup>3.09</sup><sub>2.71</sub></td><td>3.26<sup>3.52</sup><sub>3.02</sub></td><td>3.86<sup>4.44</sup><sub>3.47</sub></td><td>4.59<sup>6.41</sup><sub>3.69</sub></td><td>5.47<sup>7.49</sup><sub>4.21</sub></td><td>6.82<sup>8.89</sup><sub>5.70</sub></td></tr><tr><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>–</td><td>21.55<sup>28.58</sup><sub>18.13</sub></td><td>21.02<sup>27.72</sup><sub>19.65</sub></td><td>18.88<sup>25.07</sup><sub>13.90</sub></td><td>11.00<sup>13.78</sup><sub>1.22</sub></td><td>5.45<sup>9.87</sup><sub>1.19</sub></td></tr></tbody></table>

ABL and NBL), we obtain viscosity ratios showing that LC viscosity can be larger than the one of the overlain layer. As the thicknesses of the NBL is usually below 250 km (see Fig. 7), the differences in viscosities cannot be explained by the differences in pressure. However as proposed by Khan et al. (2023), the presence of heat producing elements could explain such a contrast of viscosities by an increase of the temperature, complemented by a small contribution of the tidal friction.

Finally, for models with both ABL and NBL layers, following Bierson and Nimmo (2016), the ABL viscosity can lead to a fraction of melt between 0 and 60%, consistent with the temperature constraints imposed by the solidus and liquidus and favoring partially melted layer. Further investigations are necessary on this aspect.

It is worth mentioning that, the viscosity distributions of the profiles with ABL-only without SIC (Fig. 7) are very similar to the one of the profiles with ABL-only with SIC (Fig. B.1). This gives credits to the hypothesis that the association of ABL and NBL is introducing more differences in the viscosity distributions than the introduction of the solid core.

### 4.4. Threshold on the liquid core viscosity

For the case of the Earth, the LC viscosity is expected to be smaller than $10^{10}$ Pa s. In our results, this condition is met by about 75% of

the models with ABL only, less than 50% for NBL only and about 33% for models with both NBL and ABL. As one can see on Figs. 9-D:F, 10-E:H and Table 5 (bottom part), if we consider only the models meeting the Earth-like LC viscosity constraint ($\eta_{LC} < 10^{10}$ Pa s), we do not see significant differences in the distributions of thicknesses and densities. The effect on the viscosity profiles is more visible by the improvement of the distribution of the viscosities but with no significant reduction of the possible intervals of values (see Fig. 8). The main conclusions driven previously remain the same, except that by imposing $\eta_{LC} < 10^{10}$ Pa s, the viscosity inversions between the LC and the upper layers (NBL, ABL and mantle) are marginals (see Fig. 10-F).

### 5. Conclusion

Despite the very short periods of excitation, tidal deformation is very efficient to constrain Mars internal structure. New estimations of densities and thicknesses were obtained in this work and give good agreement with previous studies. Combined with thermal constraints, this approach gives also new tools for estimating layer viscosities and thus, monitoring the Mars interior, and more precisely to discuss the existence of a BML.

A major result is the fact that the presence of a liquid layer on the top of the liquid core prevents the existence of a solid inner core. This result favors the scenario of Mars without solid inner core as it has been already proposed by Le Maistre et al. (2023).
A second interesting result is that the viscosity contrast between the lithosphere and the mantle can be used as a marker for the detection of a 2-layered interface between the liquid core (LC) and the mantle. This 2-layered interface would be composed by 2 layers of similar densities but very different viscosity and rheology: the layer on the top of the core is liquid (Newtonian, NBL) and the one at the base of the mantle, overlaying the liquid one, is an Andrade layer (ABL) with a viscosity in average 10 orders of magnitude higher than the Newtonian layer. What we have seen is that the models including the 2-layered interface do not display significant viscosity contrast between the mantle and the lithosphere, preventing mechanical decoupling between a lithosphere and the mantle immediately below. Such models are in favor of a stagnant lid regime that can be supported by the current absence of an Earth-like plate tectonics on Mars.

Finally, the presence of a liquid layer on the top of the LC imposes the mantle to have a viscosity greater than the viscosity of the layer immediately below (ABL or NBL). In the opposite, the presence of an Andrade layer (ABL) at the top of LC imposes a more viscous top layer. Inversion of viscosities between Andrade layers (ABL and mantle) or between Newton layers (NBL or LC) are possible in our models but more investigations are required to understand if these inversions can be explained by differences in chemical compositions or degrees of melt of the layers.

### CRediT authorship contribution statement

**Alex Guinard:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis.
**Agnès Fienga:** Writing – review & editing, Writing – original draft, Visualization, Software, Project administration, Methodology, Formal analysis.
**Anthony Mémin:** Writing – review & editing, Writing – original draft, Formal analysis.
**Clément Ganino:** Writing – review & editing, Writing – original draft, Formal analysis.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

No data was used for the research described in the article.

### Acknowledgments

The authors thank the reviewers for their constructive suggestions and A. Khan, H. Samuel and F. Nimmo for fruitful discussions. A.G. thanks Université Côte d'Azur EUR Spectrum for excellence grant and A.F and C.G. thank axe transverse *Terrestrial Planets* from Observatoire de la Côte d'Azur for support. This work was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence Université Côte d'Azur under reference number ANR-15-IDEX-01. The authors are grateful to the Université Côte d'Azur's Center for High-Performance Computing (OPAL infrastructure) for providing resources and support.

<figure><img src="image_13.png" /><figcaption>Fig. A.1. Distribution of the alpha values of the models in output from ALMA given in the Table 3 respecting the $\chi^2$ criterion described in Section 2.</figcaption></figure>

## Appendix A. Rheology and geophysical constraints

Profiles discussed in Section 4 were obtained in randomly sampling values for the $\alpha$ parameter of the Andrade rheology following a uniform distribution between 0.2 and 0.4 as proposed by Bagheri et al. (2019). Distribution of the $\alpha$ for the selected profiles is given in Fig. A.1. Profiles have been selected according to geophysical constraints as explained in Section 3. The distributions of normalized moment of inertia and the tidal Love number $k_2$ for the filtered profiles are given in Fig. A.2. The total mass is conserved by deriving the mantle density from the other randomly sampled density layers. We also provide with Fig. A.3, the density versus depth profiles obtained for models after thermal filtering.

## Appendix B. Results with a solid inner core and an andrade basal layer

Are given in this appendix the distribution in viscosity, density and thickness of the 6-layer profiles considering both the geophysical constraints and the thermal constraints on Fig. B.1. As explained in Section 4.1, only the case of 6-layer profiles with ABL-only is presented. By comparisons with Figs. 4 and 7, it is interesting to note that the introduction of the solid inner core induces a degradation of the thickness constraints on the liquid core. This is consistent with a weaker sensitivity of the geophysical constraints to the deeper layers. The fact that the interval of possible densities for the liquid core is similar with or without SIC is driven by the non-inversion of density between layers imposed to all the profiles and the small volume of the solid inner core (see Section 3 for discussion).

One can also consider the inner core as a Maxwell body instead of an elastic body. We implemented Maxwell SIC in our profiles and we have tested the sensitivity of the deformations (real and imaginary parts of the Love number $k_2$) in varying the viscosities of the inner core from $10^{20}$ to $10^{30}$ Pa s. Fig. B.2 presents the results showing no sensitivity on the SIC viscosity. The results obtained with an elastic inner core are also valid for a Maxwell SIC.
<figure><img src="image_14.png" /><figcaption>Fig. A.2. Distribution of the normalized moment of inertia and the tidal Love number $k_2$ values of the models with ABL and NBL in output from ALMA given in the Table 3 respecting the $\chi^2$ criterion described in Section 2.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. A.3. Same as Fig. 5 from the profiles of the Table 5 following both geophysical and thermal constraints.</figcaption></figure>

Table B.1

Radii, densities and viscosities rigidities of the profiles with ABL and with solid inner core presented for which the temperatures for the top of the mantle and the top of the CMB were calculated for 15 different values of $E^*$ and $V^*$ determined randomly from Gaussian manner in the ranges given in Table 4. The values presented in this table represent the quantiles 0.5, 0.995 and 0.005 and follow both the geophysical (Sec 2) and the thermal (Sec 3.2) constraints.

<table><thead><tr><th>Number of profiles</th><th>Cr</th><th>Li</th><th>Ma</th><th>ABL</th><th>LC</th><th>SIC</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center; font-weight: bold;">With ABL</td></tr><tr><td>5,721</td><td>Radius [km]</td><td>3,389.5</td><td>3,339<sup>3.362</sup><sub>3.319</sub></td><td>2,757<sup>3.316</sup><sub>2.502</sub></td><td>1,640<sup>1.876</sup><sub>1.372</sub></td><td>1,593<sup>1.797</sup><sub>1.351</sub></td><td>446<sup>7.88</sup><sub>8</sub></td></tr><tr><td></td><td>Density [g cm<sup>-3</sup>]</td><td>2.90<sup>3.10</sup><sub>2.70</sub></td><td>3.25<sup>1.53</sup><sub>3.00</sub></td><td>3.92<sup>4.99</sup><sub>3.44</sub></td><td>5.14<sup>8.27</sup><sub>3.72</sub></td><td>6.57<sup>9.34</sup><sub>4.97</sub></td><td>8.38<sup>9.99</sup><sub>5.36</sub></td></tr><tr><td></td><td>Viscosity [log<sub>10</sub>(Pa s)]</td><td>—</td><td>24.01<sup>29.90</sup><sub>17.97</sub></td><td>20.41<sup>28.76</sup><sub>18.46</sub></td><td>18.07<sup>29.84</sup><sub>2.54</sub></td><td>8.75<sup>28.98</sup><sub>1.07</sub></td><td>—</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Fig. B.1.** Density profile (top) and Distributions in thickness, density and viscosity (top) for the profiles with ABL and solid inner core, meeting the thermal constraints as given in Temperature profile (bottom). The intervals are given at 90% C.I. and can be found in Table B.1. The gray zones picture the exploration priors taken for each of the parameters after filtering with mass and moment of inertia.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. B.2. Same as Fig. 2 but in considering different range of viscosities for a Maxwell SIC.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. C.1. Distributions in thickness, density and viscosity (top) for the profiles with a lithosphere composed by 3 sublayers.</figcaption></figure>

## Appendix C. 3-Layer lithosphere

The impact of considering the lithosphere not as single layer but as 3 layers has been evaluated. The lithosphere has been subdivided in 3 layers of different densities, thicknesses and viscosities randomly sampled in such a way that the geophysical constraints are respected as described in Section 2. The obtained distributions for profiles meeting the geophysical requirements are given in Fig. C.1. As it is visible in this figure, these 3 layers do not show significant departure one from another indicating that the tidal deformations are not sensitive to the shallow layers (see Fig. C.2).

## Appendix D. Impact of viscosity gradient

Mars is hypothesized to be a single-plate planet, currently in a stagnant lid convection regime. Under these conditions, the viscosity of the cold lithosphere is expected to be significantly higher than that of the warmer underlying mantle. Additionally, in the deep mantle, the viscosity of a partially molten silicate layer, referred to as the ABL, is anticipated to be lower than that of the solid mantle above. Lastly, similar to Earth, the viscosity of the liquid core is expected to be much lower than that of the mantle. In this section, we present a revised and more constrained model, constructed in a manner consistent with Section 2, but with the added imposition of the following viscosity hierarchy: $\eta_{Li} > \eta_{Ma} > \eta_{ABL} > \eta_{NBL} > \eta_{LC}$. Fig. D.1, analogous to

Fig. 10, depicts the distributions of viscosity ratios between different contacting layers for three distinct configurations. Notable qualitative and quantitative differences emerge when compared to Fig. 10. Firstly, in Fig. D.1-C, we observe that, in configuration with both ABL and NBL, over 92% of the models exhibit an Earth-like viscosity ratio between the lithosphere and mantle (i.e., a ratio less than or equal to 2). In contrast, in either the ABL-only or NBL-only configurations, this percentage drops to 46%. Secondly, still in the configuration featuring both ABL and NBL, the viscosity ratio between these two layers remains within 5 orders of magnitude (see Fig. D.1-D), whereas in configurations without imposing a viscosity hierarchy, this ratio can range from 5 to 15 orders of magnitude (see Fig. 10-D). Thirdly, in the ABL-only configuration, the distribution of the viscosity ratio between the liquid core and the ABL shows a markedly different profile comparing Fig. 10-B to Fig. D.1-B. Specifically, while the average of this ratio is approximately -11 in the original manuscript, it rises to about -6 under the revised viscosity hierarchy constraints. Lastly, in the NBL-only configuration, the distribution of the viscosity ratio between the NBL and the mantle exhibits a much flatter curve, indicating a less constrained viscosity ratio. Thus, these imposed viscosity hierarchy constraints ($\eta_{Li} > \eta_{Ma} > \eta_{ABL} > \eta_{NBL} > \eta_{LC}$), if realistic, have significant implications for the expected viscosity ranges across Mars deep layers.
<figure><img src="image_19.png" /><figcaption>Fig. C.2. Same as Fig. 2 for a profile with 3 additional sublayers in the lithosphere. The deeper layers behave in the same way as the case where the lithosphere is not subdivided.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. D.1.** Histograms of the log10 of the viscosity ratios between layers in contact for the three configurations of profiles: with ABL-only, with NBL-only, with both ABL and NBL. Each panel was obtained after geophysical and thermal filtering. These distributions respect the criterion of decreasing viscosity with depth.</figcaption></figure>

## References

Bagheri, A., Efroimsky, M., Castillo-Rogez, J., et al., 2022. Adv. Geophys. 63, 231–320. doi:10.1016/bs.agph.2022.07.004.

Bagheri, A., Khan, A., Al-Attar, D., Crawford, O., Giardini, D., 2019. J. Geophys. Res: Planets 124 (11), 2703–2727.

Banerdt, W.B., Smrekar, S.E., Banfield, D., et al., 2020. Nat. Geosci. 13 (3), 183–189.

Bierson, C.J., Nimmo, F., 2016. J. Geophys. Res. (Planets) 121 (11), 2211–2224. http://dx.doi.org/10.1002/2016JE005005.

Briaud, A., Fienga, A., Melini, D., et al., 2023a. Icarus 115426.

Briaud, A., Ganino, C., Fienga, A., Mémin, A., Rambaux, N., 2023b. Nature 1–4.

Castillo-Rogez, J.C., Efroimsky, M., Lainey, V., 2011. J Geophys. Res. (Planets) 116 (E9), E09008. http://dx.doi.org/10.1029/2010JE003664.

Drilleau, M., Samuel, H., Garcia, R.F., et al., 2022. J Geophys. Res. Planets 127 (9), e2021JE007067.

Duncan, M.S., Schmerr, N.C., Bertka, C.M., Fei, Y., 2018. Geophys. Res. Lett 45 (19), 10,211–10,220. http://dx.doi.org/10.1029/2018GL078182.

Durán, C., Khan, A., Ceylan, S., et al., 2022. Geophys. Res. Lett. 49 (21), e2022GL100887.

Duran, C., Khan, A., Ceylan, S., et al., 2022. Phys. Earth Planet. Inter. 325, 106851.

Elkins-Tanton, L.T., Parmentier, E., Hess, P., 2003. Meteorit. Planet. Sci. 38 (12), 1753–1771.

Helffrich, G., 2017. Prog. Earth Planet. Sci. 4 (1), 24. http://dx.doi.org/10.1186/s40645-017-0139-4.

Hemingway, D.J., Driscoll, P.E., 2021. J Geophys. Res. Planets 126 (4), e2020JE006663.

Höink, T., Jellinek, A.M., Lenardic, A., 2011. Geochem. Geophys. Geosyst. 12 (10), http://dx.doi.org/10.1029/2011GC003698.

Horleston, A.C., Clinton, J.F., Ceylan, S., et al., 2022. Seism. Rec. 2 (2), 88–99.

Irving, J.C., Lekić, V., Durán, C., et al., 2023. Proc. Natl. Acad. Sci. 120 (18), e2217090120.

Karato, S.-I., Jung, H., Katayama, I., Skemer, P., 2008. Annu. Rev. Earth Planet. Sci. 36, 59–95. http://dx.doi.org/10.1146/annurev.earth.36.031207.124120.

Khan, A., Ceylan, S., van Driel, M., et al., 2021. Science 373 (6553), 434–438.

Khan, A., Huang, D., Durán, C., et al., 2023. Nature 622 (7984), 718–723.

Kim, D., Duran, C., Giardini, D., et al., 2023.

Konopliv, A.S., Park, R.S., Folkner, W.M., 2016. Icarus 274, 253–260.

Konopliv, A.S., Park, R.S., Rivoldini, A., et al., 2020. Geophys. Res. Lett. 47 (21), e2020GL090568.

Le Maistre, S., Rivoldini, A., Caldiero, A., et al., 2023. Nature 619 (7971), 733–737. http://dx.doi.org/10.1038/s41586-023-06150-0.

Li, J., Beghein, C., McLennan, S.M., et al., 2022. Nature Commun. 13 (1), 7950.

Lillis, R.J., Manga, M., Mitchell, D.L., Lin, R.P., Acuña, M.H., 2005. Lunar Planet. Sci. XXXVI, Part 12.

Melini, D., Saliby, C., Spada, G., 2022. Geophys. J. Int. 231 (3), 1502–1517. http://dx.doi.org/10.1093/gji/ggac263.

Nakada, M., Iriguchi, C., Karato, S.-i., 2012. Phys. Earth Planet. Inter. 208, 11–24.

Plesa, A.-C., Bozdağ, E., Rivoldini, A., et al., 2021. J Geophys. Res. Planets 126 (6), e2020JE006755.

Post, E.L., 1930. Trans. Amer. Math. Soc. 32 (4), 723–781.

Pou, L., Nimmo, F., Rivoldini, A., et al., 2022. J Geophys. Res. Planets 127 (11), e2022JE007291.

Rivoldini, A., Van Hoolst, T., Verhoeven, O., Mocquet, A., Dehant, V., 2011. Icarus 213 (2), 451–472.

Ruedas, T., Breuer, D., 2017. J Geophys. Res. Planets 122 (7), 1554–1579.

Saliby, C., Fienga, A., Briaud, A., Mémin, A., Herrera, C., 2023. Planetary and Space Sciences 231, 105677. http://dx.doi.org/10.1016/j.pss.2023.105677.
Samuel, H., Ballmer, M.D., Padovan, S., et al., 2021. J Geophys. Res. (Planets) 126 (4), e06613. http://dx.doi.org/10.1029/2020JE006613.

Samuel, H., Drilleau, M., Garcia, R., et al., 2022. Copernicus Meetings.

Samuel, H., Drilleau, M., Rivoldini, A., et al., 2023. Nature 622 (7984), 712–717. http://dx.doi.org/10.1038/s41586-023-06601-8.

Seidelmann, P.K., Abalakin, V.K., Bursa, M., et al., 2002. Celest. Mech. Dyn. Astron. 82 (1), 83–111. http://dx.doi.org/10.1023/A:1013939327465.

Spada, G., 2008. Comput. Geosci. 34 (6), 667–687.

Spada, G., Boschi, L., 2006. Geophys. J. Int. 166 (1), 309–321.

Stähler, S.C., Khan, A., Ceylan, S., et al., 2021. EGU General Assembly Conference Abstracts. pp. EGU21–13310.

Stern, R.J., Gerya, T., Tackley, P.J., 2018. Geosci. Front. 9 (1), 103–119. http://dx.doi.org/10.1016/j.gsf.2017.06.004, Lid Tectonics.

Wanke, H., Dreibus, G., 1994. Philosophical Transactions of the Royal Society of London Series A 349 (1690), 285–293. http://dx.doi.org/10.1098/rsta.1994.0132.

Widder, D.V., 1934. Trans. Amer. Math. Soc. 36 (1), 107–200.

Widder, D., 1941. Princeton Math. Ser. ..

Wieczorek, M.A., Broquet, A., Mclennan, S.M., et al., 2022. J Geophys. Res. Planets 127 (5), e2022JE007298.

Zeff, G., Williams, Q., 2019. Geophys. Res. Lett. 46 (20), 10997–11007.

Zhou, W.-Y., Olson, P.L., Shearer, C.K., et al., 2022. Earth Planet. Sci. Lett. 594, 117751. http://dx.doi.org/10.1016/j.epsl.2022.117751.