Tidal Heating Anomaly Zone at the core-mantle boundary of the Moon

Abstract

Young volcanism revealed by the Chang’e samples suggests an additional heat source for the generation of young lunar basalts. While radiogenic heating, accretion energy, and impact heating provided insufficient heat, early lunar tidal heating may have played a significant role. The magnitude of tidal heating is mainly controlled by the internal viscosity and orbital parameters. The recently observed ultralow-viscosity zone (ULVZ) at the lunar core-mantle boundary (CMB) indicates spatial viscosity heterogeneity, highlighting the need to reconstruct the Moon’s internal distribution of tidal heating. Using the available orbital configuration at ca. 3.2 Ga, we present a new tidal heating model revealing two symmetric enhanced tidal heating anomaly zones (THAZ) at the CMB whose heat production is an order of magnitude higher than the contemporaneous radiogenic heating. The THAZ would generate lunar mantle plumes from the CMB which transferred deep-seated heat to shallow regions. This heat delivery process is capable of providing the requisite thermal source for the formation of high-Ti basalts (e.g., mare basalt samples from Apollo 11 and Apollo 17). The additional heat provided by the THAZ is also consistent with a long-lived lunar dynamo lasting ca. 500 million years.

Graphical abstract

Keywords

Tidal heating / ultralow-viscosity zone / lunar thermal evolution / core-mantle boundary / lunar volcanism

1 Introduction

The discovery of young lunar basalts (~2 Ga) by the Chang’e-5 mission requires additional heat sources to sustain a prolonged volcanism (Li et al., 2021). Conventional mechanisms like radiogenic (Laneuville et al., 2018) and impact heating (Nesvorný et al., 2023) are inadequate to explain this activity (Carlson, 2021). Specifically, heat production from radiogenic elements (e.g., K, Th, U) within KREEP materials decays exponentially; by ~2.0 Ga, it was too weak to independently drive widespread mantle melting. Furthermore, large-scale impact events capable of generating significant heat had largely ceased by this epoch (Yue et al. 2022). Given the rapid overall cooling of the relatively small Moon, an additional dynamic internal heat source is necessary. Tidal heating emerges as a potential alternative. Early in the formation of the Earth-Moon system, a high orbital eccentricity (Meyer et al., 2010; Tian et al., 2017) would amplify the tidal dissipation (Ćuk and Stewart, 2012), which generated substantially greater internal heating compared to that of today. Previous models, however, oversimplified lunar structure by treating the entire Moon as a rigid body and neglecting its internal viscosity variations (Heller et al., 2021; Sleep et al., 2014). Recent numerical simulations of the lunar tidal response have revealed an ultralow-viscosity zone (ULVZ) at the lunar core-mantle boundary indicative of strong spatial heterogeneity in mantle rheology (Harada et al., 2014). Since the intensity of tidal heating depends nonlinearly on the viscosity structure, this heterogeneity in viscosity indicates the potential effect of tidal heating in the early thermal evolution of the Moon.

2 Methods

2.1 Tidal heating calculation

Considering the variation of the lunar orbit eccentricity in the early history of the Moon, the eccentricity-driven tides of the Moon were calculated. For a spherically symmetric, non-rotating body, the instantaneous spherically symmetric elastic deformation caused by tidal forces can be represented as a vector y with six components (y₁₋₆) (Sabadini, 2004; Takeuchi and Saito, 1972; Tobie et al., 2005):

(1)\frac{{\text{d}y}_i}{\text{d}r}=\mathit{A}\cdot \mathit{y},

The matrix A represents the fundamental matrix for each layer. The vector y is composed of the six components: 1) Radial displacement (y₁): elastic deformation along the radial direction, directly proportional to the tidal potential amplitude; 2) Tangential displacement (y₂): shear deformation perpendicular to the radial direction, dependent on the spherical harmonic degree l of tidal forcing; 3) Radial stress (y₃): normal stress component balancing tidal pressure; 4) Tangential stress (y₄): shear stress generated by differential displacement across layers, critical for viscous dissipation; 5) Gravitational potential perturbation (y₅): tidal potential correction due to mass redistribution, governed by Poisson’s equation; 6) Tidal stress potential (y₆): non-hydrostatic component derived from the divergence of the displacement field. The specific expression of the radial functions (which represent instantaneous elastic spherical deformations induced by tidal forces for a spherically symmetric non-rotating satellite model):

(2)\frac{\mathrm{d}{y}_1}{\mathrm{d}r}=-\frac{2\lambda }{\left(\lambda +2\mu \right)r}{y}_1+\frac{1}{\left(\lambda +2\mu \right)}{y}_3+\frac{\lambda l\left(l+1\right)}{\left(\lambda +2\mu \right)r}{y}_2,

(3)\frac{\mathrm{d}{y}_2}{\mathrm{d}r}=-\frac{1}{r}{y}_1+\frac{1}{r}{y}_2+\frac{1}{\mu }{y}_4,

(4)\begin{array}{rl}\frac{\mathrm{d}{y}_3}{\mathrm{d}r}=& \left[-{\omega }^2\rho r^2-4\rho gr+\frac{4\mu \left(3\lambda +2\mu \right)}{\left(\lambda +2\mu \right)}\right]\frac{1}{r^2}{y}_1-\frac{4\mu }{\left(\lambda +2\mu \right)r}{y}_3+\\ & \left[l\left(l+1\right)\rho gr-\frac{2\mu \left(3\lambda +2\mu \right)l\left(l+1\right)}{\left(\lambda +2\mu \right)}\right]\frac{1}{r^2}{y}_2+\\ & \frac{l\left(l+1\right)}{r}{y}_4+\frac{\rho }{r}{y}_6,\end{array}

(5)\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{y}_4}{\mathrm{d}r}}=& \left[\rho gr-\frac{2\mu \left(3\lambda +2\mu \right)}{\left(\lambda +2\mu \right)}\right]\frac{1}{r^2}{y}_1-\frac{\lambda }{\left(\lambda +2\mu \right)r}{y}_3+\\ & \left[-\rho {\omega }^2{r}^2+\frac{2\mu }{\left(\lambda +2\mu \right)}\left(\lambda \left(2l^2+2l-1\right)+2\mu \left(2l^2+l-1\right)\right)\right]\\ & \frac{1}{r^2}{y}_2-\frac{3}{r}{y}_4+\frac{\rho }{r}{y}_5,\end{array}

(6)\frac{\mathrm{d}{y}_5}{\mathrm{d}r}=-4\pi G\rho y_1-\frac{l+1}{r}{y}_5+y_6,

(7)\frac{\mathrm{d}{y}_6}{\mathrm{d}r}=-\frac{4\pi G\rho l\left(l+1\right)}{r}{y}_1+\frac{4\pi G\rho l\left(l+1\right)}{r}{y}_2+\frac{l-1}{r}{y}_6.

In these equations, λ and μ represent Lamé constants. The density and gravitational acceleration are denoted by \rho and ɡ, respectively. The deformation is driven by an external tidal force with an angular velocity ω. The angular order of the deformation is given by the spherical harmonic degree l. Finally, G represents the universal gravitational constant.

The resolution of radial functions using the propagation matrix technique (Sabadini, 2004; Takeuchi and Saito, 1972), satisfies the fulfilment of boundary conditions for forced oscillations:

y_1\left(0\right)=0,y_2\left(0\right)=0,y_5\left(0\right)=0,

(8)y_3\left(R_s\right)=0,y_4\left(R_s\right)=0,y_6\left(R_s\right)=\frac{2l+1}{R_s},

For small inclinations and at first-order in eccentricity, the tidal potential function can be represented as follows (Kaula, 1964):

(9)\begin{array}{rl}\mathit{\Phi }=& R_{\text{s}}^2{\omega }^{2}e[-\frac{3}{2}{P}_2^0\left(\cos \theta \right)\cos \omega t+\frac{1}{4}{P}_2^2\left(\cos \theta \right)\times \\ & \left(\text{3cos}\omega t\text{cos2}\varphi \text{+4sin}\omega t\text{sin2}\varphi \right)],\end{array}

where R_s is the radius, e is the eccentricity, θ is the co-latitude, ϕ is the longitude, and P_l^m is the associated Legendre polynomial of degree l and order m. Finally, the stress tensor σ_ij and strain tensor ε_ij are computed for each point within the celestial body.

For a viscoelastic body, described by a Maxwell rheology, the deformation involves both elastic storage and viscous dissipation. The frequency-dependent complex shear modulus \tilde{\mu }(\omega ) is given by

(10)\tilde{\mu }\left(\omega \right)=\frac{{\mu }_{\text{E}}{\omega }^2{\tau }^2}{1+{\omega }^2{\tau }^2}+i\frac{\omega \tau }{1+{\omega }^2{\tau }^2},

where \tau =\eta /{\mu }_{\text{E}} is the Maxwell relaxation time. Although the lunar mantle behaves as a stiff elastic body at the tidal forcing frequency (where \omega \tau \gg 1), we utilize the Maxwell model to explicitly calculate the viscous phase lag that generates heat. Following Nakada and Karato (2012), we assume that steady-state creep dominates the dissipation mechanism at the high temperatures relevant to the deep mantle. This approach captures the temperature-dependent energy dissipation governing the long-term thermal evolution without introducing unconstrained parameters associated with transient creep. Based on the complex power dissipation (Beuthe, 2013; Henning and Hurford, 2014; Tobie et al., 2005), the tidal heating rate at each point is obtained by integrating the product of the complex stress and strain tensors over an orbital period:

(11){\dot{E}}_t=\frac{\omega }{2}\left[\text{Im(}{\tilde{\sigma }}_{ij}\text{)Re(}{\tilde{\epsilon }}_{ij}\text{)-Re(}{\tilde{\sigma }}_{ij}\text{)Im(}{\tilde{\epsilon }}_{ij})\right],

{\tilde{\sigma }}_{ij}and {\tilde{\epsilon }}_{ij} are the components of the complex stress and complex strain tensors, respectively.

2.2 Thermal numerical simulation

To simulate the thermo-chemical evolution of spherical bodies, our numerical models are constructed using the mantle convection code CitcomS. The dynamics are governed by the fundamental conservation laws of mass, momentum, and energy. We assume the material to be incompressible and apply the Boussinesq approximation, under which density variations are neglected except in the buoyancy term. Consequently, the conservation equations for mass, momentum, and energy are simplified as follows:

(12)u_{i,i}=0,

(13){\sigma }_{ij,j}+\mathrm{R}\mathrm{a}\left(T-BC\right){\hat{e}}_r=0,

(14)\frac{\mathrm{\partial }T}{\mathrm{\partial }t}+u_i{T}_{,i}={\left(\kappa T_{,i}\right)}_{,i}+{\gamma }_{\text{rad}}+{\gamma }_{\text{tidal}}.

To track the dynamic evolution of the dense compositional layer, the advection equation for the compositional field C is explicitly included as:

(15)\frac{\mathrm{\partial }C}{\mathrm{\partial }t}+u_i{C}_{,i}=0,

where, u_i, σ_ij, Τ, and к represent velocity, Cauchy stress tensor, temperature, thermal diffusivity, respectively. The Cauchy stress tensor σ_ij, including both pressure and viscous stresses is defined as

(16){\sigma }_{ij}=-P{\delta }_{ij}+2\eta {\epsilon }_{ij},

where P is pressure, δ_ij is the Kronecker delta, η is viscosity, and ε_ij is the strain-rate tensor, given by

(17){\epsilon }_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i}).

In the momentum equation, {\hat{e}}_r denotes the unit vector in the radial direction, reflecting the spherical geometry. The energy conservation equation incorporates two distinct heat sources: γ_rad represents the radiogenic heating rate and γ_tidal explicitly denotes the spatially varying tidal heating rate calculated from the viscoelastic model.

The buoyancy term in the momentum equation is derived from a linearized equation of state, which couples the fluid density ρ to both temperature and composition to account for thermal and compositional buoyancy:

(18)\rho ={\rho }_0(1-\alpha \mathrm{\Delta }T{T}^{{}^{\prime }}+\frac{\mathrm{\Delta }{\rho }_c}{{\rho }_0}C),

where T^{{}^{\prime }} and C are the nondimensional temperature and compositional fields, respectively, and ΔT is the characteristic temperature difference between the top and bottom of the lunar mantle. ρ₀ is the reference density of the lunar mantle, α is the thermal expansivity, and Δρ_c is the intrinsic density contrast between the distinct chemical components.

The dynamic evolution of the system, primarily characterized by the temperature field T and velocity field u_i, is fundamentally governed by two dimensionless parameters: the thermal Rayleigh number Ra and the buoyancy number B. The Rayleigh number dictates the vigor of thermal convection and is defined as

(19)\text{Ra=}\frac{{\rho }_0\text{g}\alpha \mathrm{\Delta }T{R}_m^3}{\kappa {\eta }_0},

where g is the gravitational acceleration, R_m is the lunar mantle thickness, and η₀ is the reference viscosity. The buoyancy number B represents the ratio of the intrinsic chemical density contrast to the maximum thermal density anomaly, defined as

(20)B=\frac{\mathrm{\Delta }{\rho }_c}{{\rho }_0\alpha \mathrm{\Delta }T}.

In our spherical model geometry, the total planetary radius R_s is set to 1690 km, which comprises a fluid core radius R_c of 390 km and a solid mantle thickness R_m of 1300 km. The viscosity in the model is temperature-dependent:

(21)\eta ={\eta }_0\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{E}{R{T}_{\dim }}-\frac{E}{R{T}_{\text{ref}}}\right),

where η is the viscosity evaluated at the dimensional absolute temperature T_dim. E is the activation energy and R is the gas constant. The reference temperature T_ref is set to 1573 K, which reasonably represents the characteristic interior temperature of the early lunar mantle following the solidification of the Lunar Magma Ocean, constrained to remain below the peridotite solidus (Zhang et al., 2013). The temperature field evolves and changes the viscosity of the lunar mantle. The initial temperature structure of the model is determined based on the solidus line of olivine (Vlaar et al., 1994). By considering different surface temperatures, we can delineate the temperature profile of the early magma ocean (Ziethe et al., 2009). Furthermore, the thickness of the lunar crust formed by surface cooling is determined according to the solidus line of olivine. This crustal layer is treated as an elastic plate, whereas the underlying lunar mantle behaves as a viscous medium. The model is initialized with a spherically symmetric geometry, where the solid mantle spans from the core-mantle boundary (R_c = 390 km) to the lunar surface (R_s = 1690 km). The initial thermal structure is constrained by the cooling profile of the magma ocean. Based on this, the initial viscosity structure of the 1300-km-thick mantle is parameterized into several depth-dependent layers. From the core-mantle boundary upwards, the reference viscosity is defined as 2.8×10²⁰ Pa·s (0–812.5 km, serving as a low-viscosity baseline that incorporates ULVZ characteristics), 1.8×10²² Pa·s (812.5–1079.0 km), 4.5×10²⁴ Pa·s (1079.0–1137.5 km), and 8.5×10²⁴ Pa·s (1137.5–1300.0 km). The rigid outer layer, encompassing the lunar crust and upper lithosphere, corresponds to the uppermost layer of our model (depths of 0–162.5 km from the surface), exhibiting the highest structural viscosity of 8.5×10²⁴ Pa·s. For the boundary conditions, mechanical boundaries are set to free-slip at both the surface and the CMB to allow unhindered lateral flow. Thermal boundary conditions are strictly isothermal at these interfaces, with the surface temperature fixed at 273 K and the CMB temperature set according to the initial thermal state of the core.

Our tidal heating models assume uniform material properties in each layer, but real systems have strong lateral viscosity variations due to convection-driven temperature changes. The models also do not fully capture radial viscosity shifts from evolving convection. To address this, we adjusted tidal heating at each timestep by applying a correction factor f to each element’s heating rate (Han and Showman, 2008; Roberts and Nimmo, 2008; Sotin et al., 2002):

(22)f=\left[\frac{\omega \eta /\mu }{1+{\left(\omega \eta /\mu \right)}^2}\right]/\left[\frac{\omega {\eta }_0/\mu }{1+{\left(\omega {\eta }_0/\mu \right)}^2}\right],

where ω is orbital angular frequency, μ is shear modulus and η₀ is reference viscosity. At each timestep of the convection calculation, the governing equations are solved, and the viscosity for each element is updated according to Eq. (21) based on the current temperature field. This parameterized approach allows for an efficient first-order coupling of tidal heating with mantle convection. However, this is a simplification. As shown by Běhounková et al. (2010), more rigorous fully coupled viscoelastic models predict locally amplified dissipation in regions of extreme thermal gradient (e.g., within plumes). Consequently, our calculated heating rates in these regions should be viewed as conservative estimates.

Using CitcomS to simulate the thermal evolution, the entire model is divided into 12 spherical caps. To balance the significant computational cost of simulating more than 1.5 billion years of global thermal-chemical evolution with high-frequency tidal heating updates, a resolution of 16×16×16 elements per cap is utilized. While this resolution successfully captures the first-order global plume dynamics, we acknowledge that future high-resolution simulations will be necessary to fully resolve the finer-scale thermochemical entrainment and localized rheological feedbacks. Within this grid framework, the spatially varying tidal heating rate γ (r, θ, φ) in Eq. (14) is mapped to each element by assigning the precomputed value at its geometric center. This procedure allows the spatially heterogeneous tidal heating field to be incorporated into the thermal evolution calculation at each timestep.

2.3 Radioactive heating rates for the lunar mantle and crust

Our study primarily focuses on long-lived radioactive elements (half-life {\text{t}}_{1/2}>{10}^9 years), while short-lived isotopes such as ²⁶Al (half-life years) decayed within a few million years after their formation (Larsen et al., 2011). Some lunar samples indicate that the concentrations of uranium and thorium on the Moon are similar to those on Earth, though the potassium concentration is relatively low (Fig. 1) (Taylor, 1982; Taylor and Wieczorek, 2014). According to the model provided by Iriyama (1973), the concentrations of key radioactive elements in the Moon are as follows: ²³⁸U=4.15×10⁻⁸ g/g; ²³⁵U=0.0303×10⁻⁸ g/g; ²³²Th=16.6×10⁻⁸ g/g; ⁴⁰K=1.48×10⁻⁸ g/g. The calculated radioactive heat production (Fig. 1) suggests that the early Moon’s internal heat from radioactive decay could have exceeded 2 TW. The temporal evolution of radioactive heating is dynamically coupled to the redistribution history of heat-producing elements in the lunar interior. As the magma ocean cooled, differentiation resulted in the formation of the lunar mantle and crust (Taylor, 1982). Taylor et al. (2006) suggest that the elements in the lunar highland crust are significantly more enriched than the overall composition of the Moon. Therefore, the concentrations of radioactive elements in the lunar mantle and crust must be calculated separately. The K, U, and Th abundances in the highland crust represent 60%–70% of the Moon’s total values. Using the heating generation rates provided by Rybach (1988), the total radioactive heat generation rate is given by

(23)A={10}^{-5}\rho \left(9.52C_{\text{U}}+2.56C_{\text{Th}}+3.48C_{\text{K}}\right)

where A is the total radioactive heat generation rate in μWm^‒3, ρ is the density in kg m^‒3. The concentrations are given in weight-ppm for uranium C_u and thorium C_Th, and in weight % for potassium C_K (Clauser, 2011).

2.4 Initialization of the ilmenite-bearing cumulate (IBC) layer

To investigate the potential role of late-stage lunar magma ocean crystallization products in our tidal heating model, we incorporate an ilmenite-bearing cumulate (IBC) layer. This layer represents dense, titanium-rich material that is hypothesized to have formed during the final stages of magma ocean solidification and subsequently undergone mantle overturn, potentially sinking toward the core-mantle boundary. In our simulations, this process is modeled by initializing a compositionally distinct layer with a density ρ_IBC of 3390 kg·m⁻³ and a thickness of 300 km adjacent to the CMB (Zhang et al., 2013). Its rheological properties are consistent with the surrounding mantle, but its higher density provides the initial buoyancy contrast that drives its evolution. To initiate mantle convection, a small standard thermal perturbation of 1% is applied at the CMB as an initial condition. However, the substantial thermal anomaly required to mobilize the dense IBC layer is not pre-prescribed, but rather dynamically generated by the continuous tidal heating. The resulting dynamical interaction between the dense layer and the tidal-driven thermal upwellings is a key focus of our analysis, and its overturn process is visualized in Supplementary Movie S1.

3 Results

3.1 Construction of a tidal heating model

Our model considers tidal dissipation solely within the mantle because geophysical evidence confirms a predominantly fluid lunar core (radius: 250–340 km) (Briaud et al., 2023; Weber et al., 2011). In such a fluid core, negligible shear rigidity results in insignificant tidal energy dissipation; thus, core contribution to thermal budgets is negligible. The solidus temperature of olivine-rich rocks provides a key constraint on the thermal structure of the lunar mantle following the cooling of the magma ocean. As the lunar crust continued to solidify, the mantle material beneath it gradually crystallized as it cooled. The peridotite solidus serves as an indicator of the melting state of the mantle. Based on the cooling temperature profile of the early lunar magma ocean, we can infer the existence of a thermal boundary layer several tens of kilometers thick at the lunar surface. Beneath this layer, extending to depths of several hundred to approximately a thousand kilometers, the mantle temperature increases with depth and is predominantly composed of olivine. The deeper mantle layer, extending to the lunar core, follows an adiabatic gradient (Konrad and Spohn, 1997; Hirschmann, 2000).

The viscoelastic structure of the lunar mantle affects the magnitude of the spatial tidal heating rate. Nimmo et al. (2012) showed that the viscoelastic structure of the lunar mantle (particularly grain size and thermal profile) governs the amplitude of tidal heating by reconciling observed tidal dissipation Q and Love numbers without invoking mantle melting. Figure 2 presents the relationship between the reference viscosity η₀ and the tidal heating rates. Within the Moon’s viscosity range (10²⁰ to 10²¹ Pa·s) (Cassen et al., 1979), the magnitude of tidal heating rate decreases with increasing viscosity. This is due to the increased viscosity limiting the Moon’s interior deformation under Earth’s gravitational tidal forces. Therefore, the regions within the Moon that experience higher heat production are those where the viscosity is lower.

Besides the viscosity, the early tidal heating rate is highly dependent on the factors of orbital period T and orbital eccentricity e. The construction of the 3D tidal heating model has two steps: 1) calculation of the spatial heating rate; 2) numerical simulation of the thermal evolution using the updated heating rates. With the specific viscosity structure, the tidal heating rate was then derived based on the orbital eccentricity and period (Fig. 3). Informed by recent orbital evolution models (Daher et al., 2021), we selected eccentricity values of e=0.05, representing the present-day low, and e=0.5, representing an estimated early high, for our comparative analysis. We selected two contrasting orbital eccen-tricities with orbital periods ranging from 1 to 27 solar days. When the orbital period exceeds this range, the resulting tidal heating rate becomes computationally insignificant due to the exponential decrease in tidal dissipation efficiency with period. Elevated orbital eccentricity induces amplified variations in lunar orbital radius during the evolution, creating intensified differential gravitational stresses that drive enhanced viscoelastic deformation and subsequent tidal dissipation heating. Conversely, the orbital period affects the rate of change of potential energy at each point on the Moon, so as the period increases with the Moon’s rotation around Earth, the tidal heating rate decreases. Changing the orbital eccentricity and period would only affect the magnitude of the tidal heating rate. In other words, the distribution of the tidal heating within the mantle is independent of the selections of these orbital parameters. With the calculated spatial heating rate, a simulation of the thermal evolution of lunar mantle was developed. The lunar mantle parameter values used in the model are given in Table 1.

3.2 Tidal Heating Anomaly Zone (THAZ) model

The presence of ULVZ at the CMB (Khan et al., 2014) is the fundamental driver for the formation of THAZ. It has an extremely low viscosity, estimated at ca. 3×10¹⁶ Pa·s (Harada et al., 2014; Matsumoto et al., 2015; Briaud et al., 2023), and its outer boundary can extend to depths of 560−580 km (Harada et al., 2016; Tan and Harada, 2021). The ULVZ likely originated from IBC downwelling during the crystallization of the lunar magma ocean (Xu et al., 2022), a process that occurred within a relatively short period of less than a few tens of Myr (Li et al., 2019; Zhao et al., 2019). This short process ensures that the ULVZ was established early in the lunar history, soon after the differentiation of the lunar mantle (van Kan Parker et al., 2012). During the early period of the moon, the rheological anomaly contributed from the ULVZ critically enhanced the tidal dissipation within the solid mantle, concentrating the tidal heating at the CMB, which gave rise to the formation of THAZ (Fig. 4).

The early orbital parameters in our models are derived from geologic observations of tidal laminae with known depositional ages (Williams, 1990; Kvale et al., 1999; Williams, 2000). At 3.2 Ga, the Earth-Moon distance was about 70% of the present value, and the orbital period was approximately 16 present solar days ({\mathit{T}}_{\text{lunar}}={\displaystyle \frac{\mathit{l}{\mathit{y}}_1\times {\mathit{t}}_{\mathit{s}}}{24}}≈16, where ly₁ is the number of sidereal days in the sidereal month; t_s is the duration of solar day) (Eulenfeld and Heubeck, 2023). The orbital eccentricity is constrained between 0.05 and 0.5, whose lower bound is taken from the current value and the maximum eccentricity of 0.5 was inferred from early Earth-Moon orbital resonance (Garrick-Bethell et al., 2006). The maximum heating rate at the CMB was estimated to range between 6.2\times {10}^{-9} and 6.2\times {10}^{-7}\text{W}/{\text{m}}^3. The tidal heating is not only strong but also highly localized. To quantify its dominance over other heat sources, we compare the heating rates within the high-heat anomaly zone (defined as where tidal heating exceeds 50% of its peak) at ~3.2 Ga (Table 2). As shown in Table 2, within this concentrated zone constituting only ~2.5% of the mantle, the integrated tidal heating power for a high eccentricity (150 GW) exceeds the local radiogenic heating power (34.9 GW) by a factor of ~4. The peak local tidal heat flux surpasses the average mantle radiogenic flux by two orders of magnitude, underscoring its role as the dominant local heat source.

The calculated spatial distribution of heating rates was then incorporated into numeric simulation for the thermal evolution. The simulation reveals that tidal heating within the lunar mantle exhibits a highly non-uniform distribution under the specified orbital parameters (e=0.5, T=16 days). The volumetric heat production rate shows a sharp peak at the CMB, with the strongest dissipation concentrated within the equatorial region of this interface. With strong tidal heating (e=0.5), the lunar cooling would be slowed down significantly, resulting in the persistence of THAZ for over 1000 Myr (Fig. 5).

4 Discussion

The heating from THAZ alters the thermal structure of the lunar mantle, resulting in a substantially elevated temperature at the CMB compared to the initial profile, which would cause partial melting in the mantle. Figure 6 illustrates the evolution of mantle temperature relative to the solidus of peridotite in the lunar interior. At a high eccentricity (e=0.5), the temperatures in the lower mantle, gradually raised by tidal heating, exceed the solidus of peridotite, causing partial melting. In comparison, even at the lower bound of tidal heating (e=0.05), such heating was sufficient to raise mantle temperature close to that of the peridotite solidus. Even if this process alone did not initiate a widespread partial melting, it significantly pre-heated the mantle. The whole evolution is given in Appendix Fig. 1. When coupled with other heat sources likely present in the early Moon (e.g., radiogenic heating), this additional thermal energy would progressively elevate the mantle temperature above the solidus, ultimately leading to partial melting.

The significant temperature increases by up to 400 K in our model is primarily attributed to tidal heating, with limited contribution from the radiogenic decay during the Moon’s early stage. Our calculation reveals that the radioactive heat production rate in the lunar mantle was mostly diminished by ca. 3.2 Ga, with only 6×10⁻⁹ W/m³ (Fig. 7). This value is up to two orders of magnitude lower than the contemporary tidal heating production rate at the CMB (10⁻⁸ to 10⁻⁷ W/m³, Fig. 3). Although the early differentiated lunar crust has a larger radioactive heat generation than that of the mantle, we focus on the difference of tidal heating production between the CMB and the average bulk mantle.

THAZ model provides a sustained heat source for the long-lived lifetime of the core dynamo. Paleomagnetic evidence indicates that the dynamo persisted for > 500 Myr (Garrick-Bethell et al., 2009; Shea et al., 2012), requiring prolonged thermal support to maintain core convection. Our THAZ model provides an interpretation for this phenomenon: tidal dissipation concentrated at the CMB efficiently maintained an elevated local temperature. By reducing effective core-to-mantle heating transport, this process preserved the thermal gradient driving compositional convection within the core. This process not only delayed lunar core cooling but also enhanced convection within the core. It exerts a significant influence on both sustaining the dynamo action and amplifying the local magnetic field strength. Subsequent tidal heating diminished rapidly due to lunar orbital recession, becoming insufficient to power the dynamo after ca. 3.2 Ga. This result is consistent with the timing of the observed rapid decline in the paleomagnetic field (Strauss et al., 2021).

Recent studies suggest that the lunar high-Ti basalts originated as partial melts of overturned IBC but were extensively modified by reactive flow in the lunar mantle (Klaver et al., 2024). This IBC overturn occurred during the late stages of lunar magma ocean (LMO) cooling, driven by density contrasts resulting from compositional differences (Li et al., 2019; Zhao et al., 2019). However, the process lacks a sustained heat source to promote partial melting, which would enhance buoyancy to overcome the density contrast between the overturned IBC and the surrounding mantle. The tidal heating proposed in our THAZ model could fill the gap by providing the required, additional heat source. The ca. 3.7 Ga mare basalt samples of Apollo 11 and Apollo 17 demonstrate basaltic volcanism at that time. Their geochemical features of having high Ti content suggest the potential involvement of IBC in their formation process. The deep-seated heat in the THAZ model could drive mantle upwellings that transport the overturned IBC upward to facilitate the formation of high-Ti basalts. We propose that this magmatism was driven by the intense tidal heating characteristic of the early Earth-Moon evolution. Specifically, the continuous, localized tidal heating produces a localized temperature excess of up to ~200 K at the CMB. This intense heat provides the critical thermal buoyancy needed to overcome the compositional density of the IBC layer. Furthermore, the vigorous mantle plumes generated by this heating are sufficiently strong to reach the shallow mantle. These large-scale upwellings provide a viable mechanism for the massive ancient mare basalt eruptions, which represent the lunar equivalents of terrestrial Large Igneous Provinces (LIPs). Constrained by independent geological proxies that indicate a close lunar orbit at ~3.2 Ga, our model demonstrates that tidal dissipation (under high eccentricity) remained sufficient to mobilize the IBC during this period. The subsequent cessation of High-Ti volcanism is consistent with the secular decay of tidal forces, which eventually became too weak to drive the overturn of the dense IBC layer. Supplementary Video S1 details the tidal heating-driven re-overturn and upward migration of the IBC (density of IBC-rich material, {\rho }_{\text{IBC}}=\text{3390 kg}\cdot {\text{m}}^{-3}) from the CMB.

In addition, our model is capable of explaining the incomplete overturn of IBC. The ca. 2.0 Ga Chang’e-5 basalt has low Mg# and negative anomalies in Ti and Ta (Wang et al., 2024), in sharp contrast to the features of peridotite-derived low-Ti basalts. Rather, the young basalts most likely formed through partial melting of a shallow (< 100 km) IBC pyroxenite source. The evidence implies that some IBC resided at the shallow part of the lunar mantle. The vigorous thermal upwelling from the deep mantle shown in our THAZ model is consistent with these observations. These upwellings impeded the subsidence of IBC, despite the presence of substantial compositional density contrasts favoring sinking. We emphasize that these plumes were transient features sustained by intense early tidal heating (> 3.2 Ga); their dissipation following orbital recession explains the absence of present-day geophysical signatures. Therefore, while the intense early THAZ-driven heating gradually decayed, the trapped fertile IBC materials permanently remained in the shallow mantle. This incomplete overturn preserved the essential source reservoirs needed for the younger volcanism at ~2.0 Ga, as documented by the Chang’e-5 samples. Future high-resolution seismic and geochemical surveys may yet reveal deep thermochemical remnants of this ancient epoch.

Our model predicts equatorial upwellings, contrasting with the polar-amplified patterns of Enceladus due to distinct mechanical boundary conditions. While a subsurface ocean decouples Enceladus’s shell to favor polar heating, the Moon’s solid-solid interface promotes equatorial dissipation (Weber et al., 2011), an effect amplified here by the deep ULVZ. However, such equatorially symmetric upwellings do not directly explain the observed hemispheric asymmetry of lunar mare basalts (concentrated on the nearside). The surface distribution of volcanism is primarily governed by a first-order heterogeneity in crustal thickness, likely established by early asymmetric impact bombardment (e.g., the South Pole–Aitken basin), which acted as a filter for magma eruption. Our model elucidates a potent deep-mantle heat source and its dynamical consequences, while the surficial volcanic record reflects the interplay between such deep processes and pre-existing lithospheric structure.

While our model provides a coherent framework, it is important to acknowledge the methodological simplifications. The coupling between tidal heating and mantle convection employs a parameterized scheme (Roberts and Nimmo, 2008; Nakada and Karato, 2012). Such parameterized schemes may underestimate the localized enhancement of tidal dissipation within regions of sharp rheological contrast, such as the hot upwelling plumes predicted by our model (Běhounková et al., 2010). Their fully coupled models demonstrate that the strong temperature dependence of viscosity can create positive feedback, leading to peak dissipation rates potentially higher than our estimates. Therefore, this implies that our results provide a conservative lower bound on the intensity of tidal heating.

Furthermore, we explicitly acknowledge that our current thermal evolution model does not include the latent heat of melting or magmatic heat advection (heat-pipe effect). We recognize that these processes would act as efficient cooling mechanisms that compete with tidal heating once widespread melting begins. Quantifying these cooling effects requires a fully coupled two-phase flow model, which is beyond the scope of this study. Nevertheless, our results identify tidal dissipation as a potent candidate energy source. Unlike radiogenic heating, it is physically concentrated at the CMB. This study, therefore, establishes the energetic potential of tidal heating to initiate the overturn of the IBC layer, providing a theoretical basis for future comprehensive simulations to rigorously evaluate the thermal balance between tidal heat production and magmatic heat loss.

5 Conclusions

In conclusion, young lunar volcanism requires an internal heat source beyond conventional radiogenic heating. Our numerical simulations show that early tidal dissipation was strongly localized by the ULVZ at the CMB, generating two symmetric heating zones (THAZ). This intense, localized heat fundamentally altered the mantle’s thermal structure. Crucially, it preserved the thermal gradient needed to sustain a long-lived lunar dynamo for over 500 million years. Furthermore, the THAZ generated vigorous mantle upwellings that mobilized the dense IBC layer. This process provides a unified mechanism for lunar volcanism. It explains both ancient high-Ti mare basalts and the shallow IBC reservoirs that sourced the young Chang’e-5 volcanism. As the Moon receded from Earth, this THAZ-driven heating decayed rapidly. This decay explains the subsequent cessation of widespread high-Ti volcanism and the decline of the lunar magnetic field at ~3.2 Ga. Although our model gives a conservative estimate, it establishes tidal heating as a key driver of early lunar dynamics. This work provides a new perspective for future thermo-magmatic models and seismic surveys to explore this ancient epoch.

6 Appendix Figure

References

[1] Běhounková , M. , Tobie , G , Choblet , G , Čadek , O (2010). Coupling mantle convection and tidal dissipation: Applications to Enceladus and Earth‐like planets.Journal of Geophysical Research: Planets, 115(E9):

[2] Beuthe M (2013). Spatial patterns of tidal heating.Icarus, 223(1): 308–329

[3] Briaud A, Fienga A, Melini D, Rambaux N, Mémin A, Spada G, Saliby C, Hussmann H, Stark A, Viswanathan V, Baguet D (2023). Constraints on the lunar core viscosity from tidal deformation.Icarus, 394: 115426

[4] Carlson R W (2021). Robotic sample return reveals lunar secrets.Nature Publishing Group UK London,

[5] Cassen P, Reynolds R T, Graziani F, Summers A, McNellis J, Blalock L. (1979). Convection and lunar thermal history.Phys Earth Planet Inter, 19(2): 183–196

[6] Clauser C (2011). Radiogenic Heat Production of Rocks. In Encyclopedia of Solid Earth Geophysics, pp. 1018–1024 doi: 10.1007/978-3-030-58631-7_74

[7] Ćuk M and Stewart S T (2012). Making the Moon from a fast-spinning Earth: a giant impact followed by resonant despinning. science 338:1047–1052 doi: 10.1126/science.1225542

[8] Daher H, Arbic BK, Williams JG, Ansong JK, Boggs DH, Müller M, Schindelegger M, Austermann J, Cornuelle BD, Crawford EB, Fringer OB, Lau HCP, Lock SJ, Maloof AC, Menemenlis D, Mitrovica JX, Green JAM, Huber M (2021). Long‐term Earth‐Moon evolution with high‐level orbit and ocean tide models.Journal of Geophysical Research: Planets, 126(12): e2021JE006875

[9] Eulenfeld T, Heubeck C (2023). (). Constraints on Moon’s orbit 3.2 billion years ago from tidal bundle data..Journal of Geophysical Research: Planets, 128(1): e2022JE007466

[10] Garrick-Bethell I, Weiss B P, Shuster D L, Buz J (2009). Early lunar magnetism.Science, 323(5912): 356–359

[11] Garrick-Bethell I, Wisdom J, Zuber M. T (2006). Evidence for a past high-eccentricity lunar orbit.Science, 313(5787): 652–655

[12] Han L, Showman A P (2008). (). Implications of shear heating and fracture zones for ridge formation on Europa.Geophys Res Lett, 35(3):

[13] Harada Y, Goossens S, Matsumoto K, Yan J, Ping J, Noda H, Haruyama J (2014). Strong tidal heating in an ultralow-viscosity zone at the core–mantle boundary of the Moon.Nat Geosci, 7(8): 569–572

[14] Harada Y, Goossens S, Matsumoto K, Yan J, Ping J, Noda H, Haruyama J (2016). The deep lunar interior with a low-viscosity zone: Revised constraints from recent geodetic parameters on the tidal response of the Moon.Icarus, 276: 96–101

[15] Heller R, Duda J P, Winkler M, Reitner J, Gizon L (2021). Habitability of the early Earth: liquid water under a faint young Sun facilitated by strong tidal heating due to a closer Moon.Palaontol Z, 95(4): 563–575

[16] Henning W G, Hurford T (2014). Tidal heating in multilayered terrestrial exoplanets.Astrophys J, 789(1): 30

[17] Hirschmann M M (2000). (). Mantle solidus: Experimental constraints and the effects of peridotite composition..Geochem Geophys Geosyst, 1(10):

[18] Hirt C, Featherstone W (2012). A 1.5 km-resolution gravity field model of the Moon.Earth Planet Sci Lett, 329: 22–30

[19] IRIYAMA J (1973). Thermal structure within the convecting moon.Journal of Physics of the Earth, 21(2): 77–96

[20] Karato S I, Wu P (1993). Rheology of the upper mantle: A synthesis.Science, 260(5109): 771–778

[21] Kaula W M (1964). Tidal dissipation by solid friction and the resulting orbital evolution.Rev Geophys, 2(4): 661–685

[22] Khan A, Connolly J, Pommier A, Noir J (2014). Geophysical evidence for melt in the deep lunar interior and implications for lunar evolution.J Geophys Res Planets, 119(10): 2197–2221

[23] Klaver M, Klemme S, Liu X N, Hin R C, Coath C D, Anand M, Lissenberg C J, Berndt J, Elliott T (2024). Titanium-rich basaltic melts on the Moon modulated by reactive flow processes.Nat Geosci, 17(2): 118–123

[24] Konrad W, Spohn T (1997). Thermal history of the Moon: Implications for an early core dynamo and post-accertional magmatism.Adv Space Res, 19(10): 1511–1521

[25] Kuskov O, Kronrod V (1998). Constitution of the Moon: 5. Constraints on composition, density, temperature, and radius of a core.Phys Earth Planet Inter, 107(4): 285–306

[26] Kvale E P, Johnson H W, Sonett C P, Archer A W, Zawistoski A (1999). Calculating lunar retreat rates using tidal rhythmites.J Sediment Res, 69(6): 1154–1168

[27] Laneuville M, Taylor J, Wieczorek M (2018). Distribution of radioactive heat sources and thermal history of the Moon.J Geophys Res Planets, 123(12): 3144–3166

[28] Laneuville M, Wieczorek M A, Breuer D, Tosi N (2013). Asymmetric thermal evolution of the Moon.Journal of Geophysical Research: Planets, 118(7): 1435–1452

[29] Larsen K K, Trinquier A, Paton C, Schiller M, Wielandt D, Ivanova M A, Connelly J N, Nordlund Å, Krot A N, Bizzarro M (2011). Evidence for magnesium isotope heterogeneity in the solar protoplanetary disk.Astrophys J Lett, 735(2): L37

[30] Li H, Zhang N, Liang Y, Wu B, Dygert N J, Huang J, Parmentier E M (2019). Lunar cumulate mantle overturn: a model constrained by ilmenite rheology.J Geophys Res Planets, 124(5): 1357–1378

[31] Li Q L, Zhou Q, Liu Y, Xiao Z, Lin Y, Li J H, Ma H X, Tang G Q, Guo S, Tang X, Yuan J Y, Li J, Wu F Y, Ouyang Z, Li C, Li X H (2021). Two-billion-year-old volcanism on the Moon from Chang’e-5 basalts.Nature, 600(7887): 54–58

[32] Matsumoto K, Yamada R, Kikuchi F, Kamata S, Ishihara Y, Iwata T, Hanada H, Sasaki S (2015). Internal structure of the Moon inferred from Apollo seismic data and selenodetic data from GRAIL and LLR.Geophys Res Lett, 42(18): 7351–7358

[33] Meyer J, Elkins-Tanton L, Wisdom J (2010). Coupled thermal–orbital evolution of the early Moon.Icarus, 208(1): 1–10

[34] Nakada M, Karato S I (2012). Low viscosity of the bottom of the Earth’s mantle inferred from the analysis of Chandler wobble and tidal deformation.Phys Earth Planet Inter, 192: 68–80

[35] Nesvorný D, Roig F V, Vokrouhlický D, Bottke W F, Marchi S, Morbidelli A, Deienno R (2023). Early bombardment of the moon: Connecting the lunar crater record to the terrestrial planet formation.Icarus, 399: 115545

[36] Nimmo , F. , Faul , U. H., Garnero , E. J. (2012). Dissipation at tidal and seismic frequencies in a melt‐free Moon.Journal of Geophysical Research: Planets, 117(E9): E09005

[37] Nimmo F, Kleine T, Morbidelli A (2024). Tidally driven remelting around 4.35 billion years ago indicates the Moon is old.Nature, 636(8043): 598–602

[38] Pitre L, Sparasci F, Risegari L, Guianvarc’h C, Martin C, Himbert M E, Plimmer M D, Allard A, Marty B, Giuliano Albo P A, Gao B, Moldover M R, Mehl J B (2017). New measurement of the Boltzmann constant k by acoustic thermometry of helium-4 gas.Metrologia, 54(6): 856–873

[39] Roberts J H, Nimmo F (2008). Tidal heating and the long-term stability of a subsurface ocean on Enceladus.Icarus, 194(2): 675–689

[40] Rybach L (1988). Determination of heat production rate.Handbook of terrestrial heat-flow density determination, 4: 125–142

[41] Sabadini, R (2004). Applications of Normal Mode Relaxation Theory to Solid Earth Geophysics. In V Hotine-Marussi Symposium on Mathematical Geodesy: Matera, Italy June 17–21, 2003 (pp. 262-263). Berlin, Heidelberg: Springer Berlin Heidelberg.

[42] Shea E K, Weiss B P, Cassata W S, Shuster D L, Tikoo S M, Gattacceca J, Grove T L, Fuller M D (2012). A long-lived lunar core dynamo.Science, 335(6067): 453–456

[43] Sleep N H, Zahnle K J, Lupu R E (2014). Terrestrial aftermath of the Moon-forming impact.Philos Trans A Math Phys Eng Sci, 372(2024): 20130172

[44] Sotin , C , Head III, J W, Tobie , G (2002). Europa: Tidal heating of upwelling thermal plumes and the origin of lenticulae and chaos melting.Geophysical Research Letters, 29(8): 74–1

[45] Strauss B. , Tikoo S. , Gross J. , et al. (2021). Constraining the decline of the lunar dynamo field at≈ 3.1 Ga through paleomagnetic analyses of Apollo 12 mare basalts.Journal of Geophysical Research: Planets, 126(3): e2020JE006715

[46] Takeuchi H, Saito M (1972). Seismic Surface Waves: Chap. 5 of Methods in Computational Physics Volume 11, edited by BA Bolt. Academic Press, New York.

[47] Tan Y, Harada Y (2021). Tidal constraints on the low-viscosity zone of the Moon.Icarus, 365: 114361

[48] Taylor G J, Wieczorek M A (2014). Lunar bulk chemical composition: a post-Gravity Recovery and Interior Laboratory reassessment.Philos Trans A Math Phys Eng Sci, 372(2024): 20130242

[49] Taylor S R (1982). Lunar and terrestrial crusts: a constrast in origin and evolution.Phys Earth Planet Inter, 29(3−4): 233–241

[50] Taylor S R, Taylor G J, Taylor L A (2006). The moon: a Taylor perspective.Geochim Cosmochim Acta, 70(24): 5904–5918

[51] Tian Z, Wisdom J, Elkins-Tanton L (2017). Coupled orbital-thermal evolution of the early Earth-Moon system with a fast-spinning Earth.Icarus, 281: 90–102

[52] Tobie G, Mocquet A, Sotin C (2005). Tidal dissipation within large icy satellites: Applications to Europa and Titan.Icarus, 177(2): 534–549

[53] van Kan Parker M, Sanloup C, Sator N, Guillot B, Tronche E J, Perrillat J P, Mezouar M, Rai N, van Westrenen W (2012). Neutral buoyancy of titanium-rich melts in the deep lunar interior.Nat Geosci, 5(3): 186–189

[54] Vlaar N, Van Keken P, Van den Berg A (1994). Cooling of the Earth in the Archaean: consequences of pressure-release melting in a hotter mantle.Earth Planet Sci Lett, 121(1−2): 1–18

[55] Wang C, Xu Y G, Zhang L, Chen Z, Xia X, Lin M, Guo F (2024). A shallow (< 100 km) ilmenite-bearing pyroxenitic source for young lunar volcanism.Earth Planet Sci Lett, 639: 118770

[56] Weber , R C, Lin , P Y, Garnero , E J, Williams , Q , Lognonné , P (2011). Seismic detection of the lunar core.science, 331(6015): 309–312

[57] Williams G E (1990). Tidal rhythmites: key to the history of the Earth’s rotation and the lunar orbit.J Phys Earth, 38(6): 475–491

[58] Williams G E (2000). Geological constraints on the Precambrian history of Earth’s rotation and the Moon’s orbit.Rev Geophys, 38(1): 37–59

[59] Xu , M , Jing , Z , Van Orman, J A, Yu , T , Wang , Y (2022). Experimental evidence supporting an overturned iron‐titanium‐rich melt layer in the deep lunar interior.Geophysical Research Letters, 49(13): e2022GL099066

[60] Yu S, Tosi N, Schwinger S, Maurice M, Breuer D, Xiao L (2019). Overturn of ilmenite‐bearing cumulates in a rheologically weak lunar mantle.J Geophys Res Planets, 124(2): 418–436

[61] Yue Z, Di K, Wan W, Liu Z, Gou S, Liu B, Peng M, Wang Y, Jia M, Liu J, Ouyang Z (2022). Updated lunar cratering chronology model with the radiometric age of Chang’e-5 samples.Nat Astron, 6(5): 541–545

[62] Zhang N, Parmentier E, Liang Y (2013). A 3‐D numerical study of the thermal evolution of the Moon after cumulate mantle overturn: The importance of rheology and core solidification.J Geophys Res Planets, 118(9): 1789–1804

[63] Zhao Y, De Vries J, van den Berg A, Jacobs M H G, van Westrenen W (2019). The participation of ilmenite-bearing cumulates in lunar mantle overturn.Earth Planet Sci Lett, 511: 1–11

[64] Ziethe R, Seiferlin K, Hiesinger H (2009). Duration and extent of lunar volcanism: Comparison of 3D convection models to mare basalt ages.Planet Space Sci, 57(7): 784–796

RIGHTS & PERMISSIONS

The Author(s). Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)