RESEARCH ARTICLE
10.1029/2023JE007845
Key Points:• W
 e numerically simulated magmatismin the convecting mantle of the Moon
to understand its evolution• Magma is generated in the deep
mantle and ascends to the surface aspartially molten plumes driven bymelt buoyancy• Extension of partially molten regions
by the plumes causes the observed
expansion and active volcanism of theearly Moon
Supporting Information:

Supporting Information may be found in
the online version of this article.
Correspondence to:

K. U,
u-kenyo0822@g.ecc.u-tokyo.ac.jp
Citation:U, K., Kameyama, M., & Ogawa,M. (2023). The volcanic and radial
expansion/contraction history of theMoon simulated by numerical models
of magmatism in the convective mantle.
Journal of Geophysical Research:
Planets, 128, e2023JE007845. https://doi.
org/10.1029/2023JE007845
Received 29 MAR 2023
Accepted 22 AUG 2023

The Volcanic and Radial Expansion/Contraction History of theMoon Simulated by Numerical Models of Magmatism in theConvective Mantle
Ken'yo U1, Masanori Kameyama2, and Masaki Ogawa1

Department of Earth Sciences and Astronomy, The University of Tokyo, Meguro, Japan, 2Geodynamics Research Center,
Ehime University, Matsuyama, Japan
1

Abstract To understand the evolution of the Moon, we numerically modeled mantle convection andmagmatism in a two-dimensional polar rectangular mantle. Magmatism occurs as an upward permeable
flow of magma generated by decompression melting through the convecting matrix. The mantle is assumed
to be initially enriched in heat-producing elements (HPEs) and compositionally dense ilmenite-bearing
cumulates (IBC) at its base. Here, we newly show that magma generation and migration play a crucial role inthe calculated volcanic and radial expansion/contraction history. Magma is generated in the deep mantle byinternal heating for the first several hundred million years. A large volume of the generated magma ascends
to the surface as partially molten plumes driven by melt buoyancy; the magma generation and ascent causea volcanic activity and radial expansion of the Moon with the peak at 3.5–4 Gyr ago. Eventually, the Moon
begins to radially contract when the mantle solidifies by cooling from the surface boundary. As the mantle iscooled, the activity of partially molten plumes declines but continues for billions of years after the peak becausesome basal materials enriched in the dense IBC components hold HPEs. The calculated volcanic and radial
expansion/contraction history is consistent with the observed history of the Moon. Our simulations suggestthat a substantial fraction of the mantle was solid, and there was a basal layer enriched in HPEs and the IBC
components at the beginning of the history of the Moon.
Plain Language Summary We developed a numerical model of magmatism in the convecting
mantle to understand the volcanic and radial expansion/contraction history of the Moon. In the early period of
the calculated history, magma is generated in the deep mantle and ascends to the surface as partially molten
plumes driven by melt buoyancy. The plumes cause volcanic activity, and the extension of partially molten
regions by magma ascent causes radial expansion of the Moon. In its latter period, however, the Moon contracts
with time because partially molten regions solidify as they are cooled from the surface boundary. The activityof partially molten plumes declines but continues for billions of years because some materials that hostheat-producing elements (HPEs) are enriched in compositionally dense components and remain in the deep
mantle. The calculated history of radius change and volcanism is consistent with the observed lunar history.
Our simulations suggest that a substantial fraction of the mantle was solid, and a dense layer enriched in HPEs
developed at the base of the mantle at the beginning of the history of the Moon.
1. Introduction© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

U ET AL.

Understanding the mantle evolution of the Moon that is behind its observed history of volcanic activity andradius change has been a long-standing issue in studies of the interiors of terrestrial planets (e.g., Breuer &Moore, 2015; Kirk & Stevenson, 1989; Shearer et al., 2006; Solomon & Chaiken, 1976). The Moon expanded
globally by 0.5–5 km in its earlier history until around 3.8 Gyr ago as revealed by the gravity gradiometry data
(Andrews-Hanna et al., 2013, 2014; Liang & Andrews-Hanna, 2022; Sawada et al., 2016), and it then globallycontracted until today, as suggested from observations of tectonic features on the Moon (Frueh et al., 2023; Yueet al., 2017); some observations of fault scarps (thrust faults) suggest that the contraction for the past 100 Myr isaround 1 km or less (e.g., Clark et al., 2017; Klimczak, 2015; Matsuyama et al., 2021; van der Bogert et al., 2018;
Watters et al., 2010, 2015). The period when its radius reached the maximum coincides with that when the marevolcanism was active: mare volcanism became more active with time for the first several hundred million yearsof the lunar history, peaked at 3.5–3.8 Gyr ago, and then declined but continued until around 1.5 Gyr ago (e.g.,
Hiesinger et al., 2000, 2003; Morota et al., 2011; Whitten & Head, 2015). To clarify the mantle evolution that has

1 of 29

10.1029/2023JE007845

caused the observed features of the lunar history, we developed a two-dimensional (2-D) polar rectangular modelof the lunar mantle evolution where mantle convection and magmatism that transport heat, mass, and incompatible heat-producing elements (HPEs) are considered.
Various numerical models of mantle evolution have been advanced to account for the observed radius change,
in particular, the early expansion of the Moon (e.g., Shearer et al., 2006; Solomon, 1986). Classical sphericallysymmetric one-dimensional (1-D) models, where the radius change occurs only thermally, show that the early
expansion is reproduced in the models when the temperature in the deep mantle is initially 1100 K or less;these models suggest that subsequent temperature rise of more than 700 K should be caused by internal heatingin the deep mantle to account for the observed early expansion (Solomon & Chaiken, 1976). In the model of
Kirk and Stevenson (1989) that also takes volume change from compositional differentiation of the mantle into
account, the early expansion occurs even when the initial temperature of the deep mantle is as high as 1200 K.
Giant impact hypotheses for the origin of the Moon, however, suggest a much higher initial temperature for
the Moon (e.g., Canup, 2004; Cuk & Stewart, 2012; Lock et al., 2018; Pritchard & Stevenson, 2000; Rufu
et al., 2017; Stevenson, 1987). The mantle is expected to have been mostly molten to form the magma ocean after
the impact (e.g., Hosono et al., 2019; Newsom & Taylor, 1989). Even at the end of the mantle overturn, whichis expected to have occurred after solidification of the magma ocean, the temperature of the deep mantle is stillsuggested to have been as high as about 1800–1900 K (e.g., Alley & Parmentier, 1998; Boukaré et al., 2018;
Hess & Parmentier, 1995, 2001; H. Li et al., 2019). When such a high initial temperature is assumed, the early
expansion of the Moon is difficult to reproduce in the classical thermal history models (e.g., Solomon, 1986;
Solomon & Chaiken, 1976). The early expansion consistent with the observed history does not occur inthree-dimensional (3-D) spherical models where mantle convection is also considered. In the models of N. Zhang
et al. (2013a, 2013b), the early expansion does occur owing to internal heating, but the amplitude of the expansion
is much smaller than the observed value (e.g., Andrews-Hanna et al., 2013), and the period when the expansion
occurs is longer than 1 Gyr. In some models where the uppermost mantle is locally more enriched in HPEs on
the nearside as observed for the Procellarum KREEP terrane, or PKT (Laneuville et al., 2013; Wieczorek &Phillips, 2000), the expansion occurs only on the nearside and is too large to account for the observed expansion(e.g., Liang & Andrews-Hanna, 2022). To understand the observed early expansion of the Moon, U et al. (2022)
constructed a 1-D spherically symmetric model where volume change of the mantle by melting is considered inaddition to thermal expansion. The mantle expands by a few kilometers for the first several hundred million yearsof the calculated history when partially molten regions extend in the mantle by internal heating, suggesting thatmelting of the mantle played an important role in the early expansion of the Moon. Yet, mantle convection is notconsidered in this model.
Various numerical models have been advanced to account for the long-lasting volcanism of the Moon, too (e.g.,
Breuer & Moore, 2015; Shearer & Papike, 1999). Some earlier models show that partially molten regions persist
for billions of years in the upper mantle when the surface is covered with a crust enriched in HPEs or a blanketing
regolith layer (Konrad & Spohn, 1997; Spohn et al., 2001; Ziethe et al., 2009). In models where the uppermost
mantle in the nearside is more enriched in HPEs than that in the farside, the mantle has been partially molten
for more than 3 billion years (Laneuville et al., 2013, 2014, 2018). In these models, however, the distribution of
HPEs is spatially fixed or transported only by mantle convection; extraction of HPEs from the mantle to the crustby magmatism, which is known to reduce the activity of magmatism (e.g., Cassen & Reynolds, 1973; Cassen
et al., 1979; Ogawa, 2014), is not considered. In a model of Ogawa (2018a) where extraction of HPEs from the
mantle by migrating magma is considered, partially molten regions have been completely cooled, and magmatism stops within the first 2 Gyr of the calculated history; the end of this magmatism is too early to account for
the lunar mare volcanism (e.g., Hiesinger et al., 2003). It is also important to consider not only HPE-extraction
by migrating magma, but also the effects of structural evolution of the mantle in studies of the lunar volcanism(e.g., Hess & Parmentier, 1995; Stegman et al., 2003; N. Zhang et al., 2013a; W. B. Zhang et al., 2022; Zhong
et al., 2000). In the early period of the Moon, a compositionally dense layer, which is enriched in HPEs andilmenite-bearing cumulates (IBC) components, is suggested to have developed at the base of the mantle by crystal fractionation in the magma ocean and subsequent mantle overturn (e.g., Elkins-Tanton et al., 2011; Hess &Parmentier, 1995; Snyder et al., 1992). Earlier mantle convection models suggest that the basal layer becomes
thermally buoyant owing to internal heating and eventually migrates upward to the surface as hot plumes to cause
mare volcanism after around 4 Gyr ago (e.g., Hess & Parmentier, 1995; Zhong et al., 2000). Whether or not thebasal layer rises, however, depends on the compositional density contrast between the layer and the overlying
U ET AL.

2 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

mantle. These models assume conditions where the density contrast is low enough to allow the basal layer to riseby thermal convection (see Figure 1 in Le Bars and Davaille (2004)). The density contrast after mantle overturnis unclear and is influenced by the initial condition of overturn models (e.g., H. Li et al., 2019; Yu et al., 2019; W.
B. Zhang et al., 2022); further studies with various density content are necessary to fully understand to effects ofcompositional mantle structure on volcanic history.
To understand the mantle evolution in the Moon that is constrained by its volcanic and radial expansion/contraction history, we extend the 1-D spherically symmetric model we developed earlier (U et al., 2022) to a 2-D polar
rectangular model. In our previous model, we considered magma generation by internal heating and magma
migration that transports heat, HPE, and mass as well as volume change of the mantle by melting. In addition
to these effects, here in this study, we also include the effects of mantle convection. Heat and mass transport bymantle convection plays an important role in mantle evolution (e.g., Spohn et al., 2001; N. Zhang et al., 2017).
Upwelling flows of mantle convection generate magma by decompression melting, which also affects the mantle
evolution of the Moon (e.g., Ogawa, 2020). We also take the effects of compositional layering at the base of the
mantle formed by the magma ocean and mantle overturn on the subsequent mantle evolution into account. Thismodel is an extension of the one presented in Ogawa (2014, 2018a) in that we considered the volume changeof mantle by melting, more systematically studied the effects of the initial condition, and calculated in a 2-D polar
rectangle rather than a rectangle.

2. Model Description
A finite difference numerical code calculates the energy, mass, and momentum equations for mantle magmatism
and mantle convection in a two-dimensional polar rectangle
𝐴𝐴
𝐴𝐴 = [(𝑟𝑟𝑟 𝑟𝑟)| 385 km ≤ 𝑟𝑟 ≤ 1735 km, 0 ≤ 𝜃𝜃 ≤ 𝜋𝜋] ona mesh with 128 (radial direction) times 256 (lateral direction) mesh points under the Boussinesq approximation, where the inner and the outer radii correspond to the core and planetary radii, respectively, of the Moon(e.g., Viswanathan et al., 2019; Weber et al., 2011; Yan et al., 2015). The mantle contains incompatible HPEs
that decay with time. The mantle is considered a Newtonian fluid with a strongly T-dependent viscosity, andconvection is driven by thermal, composition, and melt buoyancy. The convecting materials are a binary eutectic
system between olivine-rich materials and IBC. Magmatism occurs as generation of basaltic magma enriched in
IBC materials and HPEs by decompression melting and internal heating together with upward permeable flow of
the generated magma through the matrix (McKenzie, 1984; Vander Kaaden et al., 2015); the permeable flow isdriven by the buoyancy of the melt. We assumed that the solid phase is always denser than the melt phase. Someearlier studies, however, suggest that IBC-rich magma can be denser than the coexisting matrix in the deep mantleof the Moon (e.g., Sakamaki et al., 2010; van Kan Parker et al., 2012; Xu et al., 2022). We also examined how
this melt-solid density inversion can influence our numerical results (Section 3.4). As mentioned there, we foundthat the density inversion does not affect the overall features of mantle evolution.
The crust of 35 km in thickness is placed on the top of the mantle (Wieczorek et al., 2013). The thermal diffusivityof the crust is about half that of the mantle, and hence the crust serves as a blanketing layer (Ziethe et al., 2009).
The temperature is fixed at Tsur = 270 K on the surface boundary, while the core is modeled as a heat bath of
uniform temperature; the vertical sidewalls are insulating. All of the boundaries are impermeable for both magma
and matrix and are shear stress-free. In the initial condition, we assume that the deep mantle is more enrichedin HPEs and the IBC components than the shallower mantle as earlier models of mantle differentiation by themagma ocean and subsequent mantle overturn suggest (e.g., Hess & Parmentier, 1995; Moriarty III et al., 2021;
Ringwood & Kesson, 1976). Several previous studies, however, suggest that a portion of the IBC materials
enriched in HPEs persist just beneath the crust even after mantle overturn (e.g., Schwinger & Breuer, 2022; Yu
et al., 2019; Zhao et al., 2019). Instead of explicitly simulating these remains of HPE-enriched materials, we
assume that the crust is uniformly enriched in HPEs.
2.1. The Properties of Materials
The convecting material is a binary eutectic system. The composition is written as AξB1−ξ; the mantle materials
are modeled as a mixture of olivine-rich materials (A) with a density of 3,300 kg m −3 (Elkins-Tanton et al., 2011)
and the IBC (B) with a density of 3,745 kg m −3 (Elkins-Tanton et al., 2011; Rapp & Draper, 2018; Shearer
et al., 2006; Snyder et al., 1992). The eutectic composition is A0.1B0.9 (ξe = 0.1) which corresponds to the basaltic
U ET AL.

3 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planetscomposition enriched in the IBC components (see Appendix A for the detailof the thermodynamic formulation).

Table 1The Meanings of the Symbols and Their Values
Symbol𝐴𝐴

Meaning

10.1029/2023JE007845

Value

The density ρ is written as𝜌𝜌 = (1 − 𝜙𝜙)𝜌𝜌s + 𝜙𝜙𝜙𝜙l .
(1)

Tsur

Surface temperature

270 K

T0

Solidus at the surface

1360 K

rp

Radius of the surface

1,735 kmrc

Radius of the core-mantle boundary (CMB)

385 kmrcrst

Radius of the Moho

1,700 kmρ0Reference density

3,300 kg m −3𝜌𝜌s = 𝜌𝜌0 [1 − 𝛼𝛼(𝑇𝑇 − 𝑇𝑇sur ) + 𝛽𝛽(1 − 𝜉𝜉s )],(2)κThermal diffusivity of the mantle

6 × 10 −7 m 2 s −1

κedd

Eddy diffusivity in largely molten region100κ at ϕ > 0.4

andwhere ρs is the density of solid phase, ρl the density of melt phase, and ϕ themelt-content. The densities depend on the temperature T and the content of
the end-member A in the solid phase ξs and that in the melt phase ξl as

}
{
Δ𝑣𝑣l
𝜌𝜌l = 𝜌𝜌0 1 − 𝛼𝛼(𝑇𝑇 − 𝑇𝑇sur ) + 𝛽𝛽(1 − 𝜉𝜉l ) −[1 + 𝛽𝛽(1 − 𝜉𝜉l )] ,
(3)
𝑣𝑣0

gsur

Gravitational acceleration at the surface

1.62 m s

gc

Gravitational acceleration at the CMB

0.55 m s −2ΔhLatent heat of melting

657 kJ kg −1

CpSpecific heat

1240 J K −1 kg −1

η0Reference viscosity

10 20 − 10 22 Pa s

ηmelt

Melt viscosity

ET

Sensitivity of viscosity to temperature

ϕ0

Reference melt-content

0.05ℎ = 𝐶𝐶𝑝𝑝 𝑇𝑇 + 𝜙𝜙Δℎ(1 + 𝐺𝐺),
(4)

𝐴𝐴𝝓𝝓0

Reference permeability

8.6 × 10 −15 m 2

αm

Thermal expansivity in the mantle

3 × 10 −5 K −1

αc

Thermal expansivity in the core

9 × 10 −5 K −1

where Cp is the specific heat; Δh the latent heat of melting; G the function which depends on Δvl/v0 (see Equation 10 below). Δvl/v0 expresses theamount of density reduction by melting as[]Δ𝜌𝜌zero − Δ𝜌𝜌∞Δ𝑣𝑣l
1
Δ𝜌𝜌∞ +
=
.
(5)
𝑣𝑣0
𝜌𝜌0
(𝑃𝑃 ∕𝜆𝜆 + 1)2−2

1 − 20 Pa s
11.3 × 10 −3 K −1

where ρ0 is the reference density; α the thermal expansivity; β = 0.135 aconstant estimated from the density of olivine-rich end-member A and thatof the IBC end-member B. The temperature is calculated from the energy
equation (Equation A3 in Appendix A) that describes the evolution of the“reduced” enthalpy h:Here, Δρ∞/ρ0 = 0.005 is the dimensionless density difference between solid and liquid phases when pressure𝑟𝑟
is infinity, and Δρzero/ρ0 = 0.22 that at zero-pressure.
𝐴𝐴
𝐴𝐴 = − ∫𝑟𝑟 𝜌𝜌0 𝑔𝑔𝑔𝑔𝑔𝑔 is the pressure; the gravitational accel𝑟𝑟 −𝑟𝑟𝑝𝑝

eration g depends on the depth
𝐴𝐴 as 𝐴𝐴 = 𝑔𝑔sur − (𝑔𝑔sur − 𝑔𝑔c ) 𝑟𝑟 𝑝𝑝−𝑟𝑟 (Garcia et al., 2011, 2012). We assumed the value𝑝𝑝

𝑐𝑐of λ = 16.42 GPa so that the solidus temperature is calculated from the Clausius-Clapeyron relationship (seeEquations 9 and 10 below) becomes close to that in the lunar mantle (Garcia et al., 2011, 2012; Katz et al., 2003).
Magma is generated by decompression melting and internal heating. The generated magma migrates upward as apermeable flow through the coexisting matrix driven by its buoyancy. Migrating magma transports heat, basaltic
materials, and HPEs. The difference between the velocity of magma u and that of matrix U is proportional to thedensity difference between them as𝑘𝑘𝝓𝝓
𝐮𝐮 − 𝐔𝐔 =𝑔𝑔(𝜌𝜌𝑠𝑠 − 𝜌𝜌𝑙𝑙 )𝐞𝐞𝐫𝐫 ,(6)𝜙𝜙𝜙𝜙melt

where kϕ is the permeability that depends on the melt-content ϕ as
( )3
𝜙𝜙
(7)
𝑘𝑘𝝓𝝓 = 𝑘𝑘𝝓𝝓0
𝜙𝜙0(McKenzie, 1984) and er is the unit vector in the radial direction. We assumed that ϕ0 = 0.05 and truncated ϕ at 0.4
for numerical reasons. The assumed range
𝐴𝐴 of 𝐴𝐴𝝓𝝓0 given in Table 1 is based on the earlier works of McKenzie (1984)
and Miller et al. (2014). In the top-most 150 km of the polar rectangular, we inserted
𝐴𝐴
𝐴𝐴′ = max(𝜙𝜙𝜙 𝜙𝜙𝑒𝑒 ) into the ϕin Equation 7 to mimic magma migration that occurs by crack propagation in the crust and the uppermost mantleof the Moon (Head & Wilson, 2017; Wilson & Head, 2017)[max(𝜙𝜙𝜙 𝜙𝜙𝑒𝑒 )]3
𝑘𝑘𝝓𝝓 = 𝑘𝑘𝝓𝝓0
,
(8)
𝜙𝜙30

U ET AL.

4 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 1. An illustration of the initial distributions of (a) bulk composition ξb, (b) internal heating rate q, (c) viscosity η at η0 = 10 20 Pa s and Tp = 1550 K, (d)temperature T, and (e) equivalent temperature h/Cp in the depth range shown by the rectangle in (d). In (d) and (e), the black dashed line is the solidus which is∗
calculated by Equation 9. Note that h/Cp = T holds for regions below the solidus (Equation 4). The meaning of variable parameters, 𝐴𝐴𝑇𝑇∗ , and Tp is described in
𝐴𝐴 𝐴𝐴crst
𝐴𝐴 l*,
Table 2.

where ϕe = 0.035 is an “effective porosity” for magma-migration by crack propagation (Kameyama et al., 1996).
The solidus temperature Tsolidus depends on the pressure as
𝑇𝑇solidus = 𝑇𝑇0 (1 + 𝐺𝐺),
(9)

where
𝑃𝑃

Δ𝑣𝑣l
1
𝐺𝐺 =
𝑑𝑑𝑑𝑑 𝑑(10)𝜌𝜌0 Δℎ ∫0 𝑣𝑣0Here, T0 is the solidus temperature at the surface (Katz et al., 2003).
The viscosity of the mantle η depends on temperature as𝜂𝜂 = 𝜂𝜂0 exp[𝐸𝐸𝑇𝑇 (𝑇𝑇ref − 𝑇𝑇 )],(11)where η0 is the reference viscosity; Tref = 1575 K is the reference temperature. ET is the sensitivity of viscosity to temperature; the adopted default value of ET is 11.3 × 10 −3𝐴𝐴K −1 (𝐴𝐴𝑇𝑇∗ = 6, see Section 2.3 below) which
implies that the viscosity decreases by a factor of 3 as temperature increases by 100 K (Figure 1c). This rangeof viscosity variation is appropriate for Newtonian rheology to mimic thermal convection in mantle materials(Dumoulin et al., 1999). In this study, we assumed that the viscosity of the mantle depends only on the temperature for simplicity. Some earlier studies, however, suggest that the viscosity depends on the contents of melt
Table 2The Variable Parameters of the Numerical Models
Variable parameter𝐴𝐴𝐴𝐴

U ET AL.

Meaning

Range of value

M*

Dimensionless reference permeability

5–100

Ra

Rayleigh number

𝐴𝐴𝑇𝑇∗

Dimensionless sensitivity of viscosity to temperature

Tp

The initial potential temperature at depths

1400–1700 Kl*

The thickness of the overturned layer after mantle overturn1/6.5–1/1.5∗
𝐴𝐴crstConcentration ratio of HPEs in the crust to the mantle

2.15 × 10 4–10 6
3–9

4–32

5 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

and IBC components (e.g., Dygert et al., 2016; Hirth & Kohlstedt, 2003; H. Li et al., 2019; Mei et al., 2002;
Scott & Kohlstedt, 2006; N. Zhang et al., 2017). We examined how the reduction of viscosity by melting and
IBC components can affect our numerical results in Section 3.3. As mentioned there, we confirmed that boththe ϕ-dependence and IBC-dependence of viscosity do not significantly affect the overall features of mantle
evolution although found that the ϕ-dependence suppresses “melt-fingers,” both of which we will describe in
Section 3 below.
In the estimate of the lunar radial expansion, we consider the effect of melting as well as thermal expansion/
contraction; the radius change ΔR is given by
[
)
]
𝜋𝜋 (
𝑟𝑟p
Δ𝑣𝑣l
1 𝛼𝛼c Δ𝑇𝑇c 2
𝛼𝛼
Δ𝑇𝑇
+
Δ𝑅𝑅
=
+
Δ𝜙𝜙
𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟
.
𝑟𝑟
m
(12)
c
∫𝑟𝑟c ∫0
𝑟𝑟p
2
𝑣𝑣0
In this equation, the first and second terms in the bracket of the right-hand side represent the volume changes of
the core and the mantle, respectively. Here, αc is the thermal expansivity of the core, ΔTc ≡ Tc(t) − Tc(t = 0) is thedeviation of the temperature Tc in the core at the elapsed time t from its initial value Tc(t = 0), and ΔT ≡ T(r, θ,t) − T(r, θ, t = 0) and Δϕ ≡ ϕ(r, θ, t) − ϕ(r, θ, t = 0) are the deviations in the temperature T and melt-content ϕin the mantle from their initial values, respectively. We assume that the initial condition in this model is the state
immediately after mantle overturn (see below), based on the fact that the tectonic features of expansion/contraction in the crust are not recorded during the solidification of the magma ocean (Elkins-Tanton & Bercovici, 2014).
2.2. The Initial Condition
The initial thermo-chemical state of the mantle is specified by the initial distributions of heating rate q, bulk
composition ξb, and “reduced” enthalpy h. Here, the bulk composition is calculated from the composition of thesolid phase ξs and that of the liquid phase ξl𝜉𝜉b = (1 − 𝜙𝜙)𝜉𝜉s + 𝜙𝜙𝜙𝜙l .
(13)Figure 1 shows an example of the initial distributions of these quantities (see Appendix B for the equations that
describe the distributions). This initial condition is motivated by earlier models of the mantle overturn that is
expected to have occurred after the solidification of the magma ocean (e.g., Alley & Parmentier, 1998; Boukaré
et al., 2018; Rapp & Draper, 2018; Snyder et al., 1992). In the last stage of crystal fractionation of the magma
ocean, a dense layer of the IBC enriched in HPEs develops at the top of the mantle. This layer is expected to sinkdown to the base of the mantle due to gravitational instability called mantle overturn (e.g., de Vries et al., 2010;
Hess & Parmentier, 1995; Ringwood & Kesson, 1976). To simulate the layering at the base of the mantle formedby the overturn, we assumed that the internal heating rate q and the content of the IBC components exponen[]tially increase with depth
𝐴𝐴 as exp −(𝑟𝑟 − 𝑟𝑟c )∕𝑙𝑙 where l represents the thickness of the overturned layer (Figures 1aand 1b). In Figure 1b, we adjusted the initial value of q so that its average in the entire calculated polar rectangle
becomes equal to the average value for the bulk silicate Moon. The internal heating rate in the topmost 35 km∗
is also higher than the average value for the mantle by a factor. We assumed this enriched layer to mimic
𝐴𝐴 of 𝐴𝐴crst
the crust and the KREEP (K, rare earth elements, and P-rich material) layer of the Moon (e.g., Ringwood &Kesson, 1976), which persist even after the mantle overturn (e.g., Zhao et al., 2019). As for the thermal initialcondition, we assumed h/Cp = Tp in the entire mantle except around the top and bottom boundaries (Figures 1dand 1e). Note that Tp is the actual temperature in the mid-mantle where materials are solid (Equation 4). For moredetails of the initial q, ξb, and h-distributions, see Appendix B.
2.3. The Parameter Values
𝐴𝐴

We carried out numerical experiments at various values of the non-dimensional reference permeability
𝑘𝑘𝝓𝝓 𝜌𝜌0 𝑔𝑔sur 𝐿𝐿
𝜌𝜌 𝛼𝛼 (Δℎ∕𝐶𝐶 )𝑔𝑔 𝐿𝐿3
𝐴𝐴 ∗ ≡ 0𝜅𝜅𝜅𝜅(see Equation A15) and Rayleigh number
𝐴𝐴
𝐴𝐴𝐴𝐴 ≡ 0 m 𝜅𝜅𝜅𝜅 𝑝𝑝 sur (see Equation A13), as well as
melt

0

the sensitivity of viscosity to temperature
𝐴𝐴
𝐴𝐴𝑇𝑇∗ = 𝐸𝐸𝑇𝑇 Δℎ∕𝐶𝐶𝑝𝑝 , the potential temperature Tp, the thickness of the∗overturned layer l* = l/L, and the ratio of the concentration of HPEs in the crust to that in the mantle, where𝐴𝐴
𝐴𝐴crst
L = rp − rc, is the depth of the mantle. The non-dimensional reference permeability M* is proportional
𝐴𝐴 to 𝐴𝐴𝝓𝝓0 ∕𝜂𝜂melt
and expresses how easily the upward migration of magma takes place. We varied the values of M* in the rangeof 5 ≤ M* ≤ 100, which corresponds to the change in the viscosity ηmelt of magma by a factor of 20 (1–20 Pa s).
U ET AL.

6 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 2. Snapshots of the distributions of (a) temperature T, (b) melt-content ϕ, (c) internal heating rate q, and (d)composition ξb calculated for Case Ref where M* = 100; Ra = 2.15 × 10 6, corresponding to η0 = 10 20
𝐴𝐴 Pa s; 𝐴𝐴𝑇𝑇∗ = 3;∗
Tp = 1550 K; l*𝐴𝐴= 4.5; 𝐴𝐴crst = 8. The elapsed times are indicated in the figure. In (d), the blue color stands for theilmenite-bearing cumulates components, while the red color for the olivine-rich end-member. The numbers [1]–[5]correspond with those of Figure 4.

On the other hand, the Rayleigh number Ra is inversely proportional to the reference viscosity and expresses how
readily mantle convection occurs. We estimated that this value is in the range of Ra = 2.15 × 10 4–10 6 corresponding to η0 = 10 22–10 20 Pa s), based on an earlier estimate of mantle rheology (Karato & Wu, 1993).
Other parameters listed in Table 2 are also varied to see how the numerical results depend on the thermal andcompositional structure of the mantle assumed in the initial conditions. The range of Tp is decided according
to that of initial temperature in Laneuville et al. (2013) and that of post-overturn temperature in some models(Boukaré et al., 2018; H. Li et al., 2019; Yu et al., 2019). l* is estimated from the density of the basal layer
which is enriched in the overturned materials obtained by the results of previous overturn models (e.g., H. Li
et al., 2019; Zhao et al., 2019). The thickness of the overturned layer is thinner and more concentrated in HPEs∗
and the IBC components with a lower l*. We also take the value
from the estimates in Konrad and
𝐴𝐴 of 𝐴𝐴crst
Spohn (1997) and Spohn et al. (2001).

3. Results
3.1. The Reference Case
Figures 2 and 3 show the reference case (Case Ref) calculated at the reference permeability M* = 100; theRayleigh number Ra = 2.15 × 10 6, corresponding to η0 = 10 20 Pa s; the sensitivity of viscosity on temperature
U ET AL.

7 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 3. (a) The horizontal averages of the temperature-distributions at various elapsed times for Case Ref. Also shown are
the horizontally averaged temperature at 4.4 Gyr for Case Ra2.15e4 (Ra = 2.15 × 10 4) and𝐴𝐴ET3 (𝐴𝐴𝑇𝑇∗ = 3). In (a), the gray andlight-gray areas are the temperature distributions in today's lunar mantle inferred by Khan et al. (2006, 2014), respectively;
the yellow areas are those estimated from the heat flux at the surface (Siegler & Smrekar, 2014; Siegler et al., 2022); thegray lines are those estimated under the assumption that the mantle consists of dry olivine (the solid line) and wet olivinewith 0.01 wt % H2O (the dashed line) by Karato (2013). Also shown are (b) the average temperature in the mantle (bulk)
and the temperature of the core; (c) the horizontal average of heat flux on the surface and the core-mantle boundary; (d) the
root-mean-square average of matrix-velocity in the mantle all plotted against time.

𝐴𝐴

𝐴𝐴𝑇𝑇∗ = 6; the potential temperature Tp = 1550 K; the initial thickness of the overturned layer l* = 1/4.5; the initial∗
crustal fraction of the 𝐴𝐴
HPEs 𝐴𝐴crst= 8 (Table 2).

3.1.1. Thermal and Structural Evolution of the Mantle
Figures 2 and 3 as well as the animation in (Movies S1–S5) illustrate how the mantle evolves dynamically by magmatism and mantle convection. The partially molten region in the uppermost mantle assumed in the initial condition
shrinks with time owing to conductive cooling from the surface boundary (Figure 2b for 0.36–0.68 Gyr). In contrast,
the temperature rises in the deep mantle that is enriched in HPEs in the initial condition, and magma is generated
there within the first 150 Myr (Figures 2a and 3a; see also Movies S1–S3). The distributions of melt-content, HPEs,
and bulk composition in the deep mantle are laterally uniform at the beginning of the calculated evolutionary historyof the mantle (Figures 2b–2d for 0.36 Gyr). However, finger-like structures, or “melt-fingers,” eventually develop
along the top of the partially molten region at a depth after around 0.4 Gyr. The fingertips are enriched in HPEs andthe IBC components because of their transport by upward migrating magma (Figure 2 for 0.40 Gyr). Following thedevelopment of melt-fingers, partially molten plumes develop and ascend often along the fingers (see the arrows in
Figure 2b) and induce peaks in the plot of rms-velocity in Figure 3d. Most of the partially molten plumes ascend to
the uppermost mantle and reach the depth levels as shallow as around 25 km by 0.7 Gyr (Figure 2b and Table 3; seealso magma flux in Figure S1 of the Supporting Information S1). These plumes in the uppermost mantle graduallysolidify as they are cooled from the surface boundary. Some partially molten plumes, however, solidify in the mid
and deep mantle (see the dashed circle in Figure 2c for 0.68 Gyr). The solidified materials enriched in HPEs andthe IBC components then sink to the deep mantle because of their compositionally induced negative buoyancy (seeFigure 2d; see also Movie S4 for around 0.60 and 0.84 Gyr). Internal heating by HPEs in these materials induces
further partial melting in the deep mantle. As a consequence, partially molten plumes develop until around 4 Gyr,although they become fainter with time as the HPEs decline after around 2 Gyr (2b and 3b–3d).
Some materials in the uppermost mantle are enriched in HPEs and the IBC components by partially molten
plumes (Figures 2c and 2d). These materials founder into the deep mantle owing to their compositionally induced
U ET AL.

8 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Table 3The Values of the Parameters Listed in Table 2 and Results

U ET AL.

Case #

M*

Ra𝐴𝐴

Ref

100

M50
M20
M5

∗
𝐴𝐴crst

ΔRpeak(km)tΔR (Gyr)

χGyr (km Gyr −))

Lϕ
(km)

4.5

8

3.00

0.66

−1.01

25

1550

4.5

8

4.44

0.70

−1.33

185

1550

4.5

8

6.62

0.89

−1.47

225

𝐴𝐴𝑇𝑇∗

Tp (K) 𝐴𝐴1/l*

2.15 × 10 6

6

1550

50

2.15 × 10

6

6

20

2.15 × 10 6

6

6

5

2.15 × 10

6

1550

4.5

8

8.95

1.73

−1.78

200

Ra2.15e4100

2.15 × 10 4

6

1550

4.5

8

3.84

1.04

−1.16

75

Ra2.15e5100

2.15 × 10

5

6

1550

4.5

8

3.80

1.05

−1.17

60

Ra7.15e5100

7.15 × 10 5

6

1550

4.5

8

3.53

0.83

−0.96

40

ET3

100

2.15 × 10

6

3

1550

4.5

8

3.06

0.78

−0.76

15

ET9

100

2.15 × 10 6

9

1550

4.5

8

3.30

0.66

−1.28

165

Tm1400100

2.15 × 10 6

6

1400

4.5

8

3.62

1.36

−0.57

30

Tm1475100

2.15 × 10 6

6

1475

4.5

8

3.34

1.13

−0.67

25

Tm1625100

2.15 × 10

6

6

1625

4.5

8

3.80

0.50

−1.00

30

Tm1700100

2.15 × 10 6

6

1700

4.5

8

4.23

0.42

−1.41

115

l1.5

100

2.15 × 10 6

6

1550

1.5

8

3.01

1.04

−1.17

190

l2.5

100

2.15 × 10 6

6

1550

2.5

8

3.66

0.86

−1.43

80

l3.5

100

2.15 × 10

6

6

1550

3.5

8

3.51

0.69

−0.97

170

l5.5

100

2.15 × 10 6

6

1550

5.5

8

2.92

0.85

−1.20

20

l6.5

100

2.15 × 10

6

6

1550

6.5

8

2.50

0.62

−0.93

15

l0

100

2.15 × 10 6

6

1550

0

8

0.22

1.45

−0.62

215

Fcrst4100

2.15 × 10 6

6

1550

4.5

4

3.85

0.69

−1.23

25

Fcrst16100

2.15 × 10 6

6

1550

4.5

16

1.69

1.07

−1.03

190

Fcrst32100

2.15 × 10 6

6

1550

4.5

32

−0.50

1.46

−0.95

395

M50-Ra2.15e450

2.15 × 10 4

6

1550

4.5

8

5.25

0.79

−0.83

205

M20-Ra2.15e4

20

2.15 × 10

4

6

1550

4.5

8

8.91

1.35

−0.82

305

M5-Ra2.15e4

5

2.15 × 10 4

6

1550

4.5

8

12.21

3.36

−0.37

285

M0-Ra2.15e4

02.15 × 10

4

6

1550

4.5

8

15.32

4.40

1.21

375

M50-l5.550

2.15 × 10 6

6

1550

5.5

8

4.00

0.73

−1.22

160

M20-l5.5

20

2.15 × 10

6

6

1550

5.5

8

7.09

0.83

−1.12

210

M5-l5.55

2.15 × 10 6

6

1550

5.5

8

9.00

1.80

−1.25

195

Tm1400-l5.5

100

2.15 × 10

6

6

1400

5.5

8

3.86

1.12

−0.10

155

Tm1700-l5.5

100

2.15 × 10 6

6

1700

5.5

8

4.09

0.35

−1.16

40

Tm1400-l1.5

100

2.15 × 10

6

6

1400

1.5

8

4.91

1.68

−0.84

210

Tm1700-l1.5

100

2.15 × 10 6

6

1700

1.5

8

−0.21

0.53

−1.08

105

Tm1400-l5.5

52.15 × 10

6

6

1400

5.5

8

9.83

3.40

−0.13

270

Ra2.15e4-l5.5

100

2.15 × 10 4

6

1550

5.5

8

3.44

1.09

−1.17

25

Ra2.15e5-l5.5

100

2.15 × 10

5

6

1550

5.5

8

3.15

1.10

−1.17

25

M50-Ra4-l5.550

2.15 × 10 4

6

1550

5.5

8

4.98

1.15

−1.17

175

M20-Ra4-l5.5

20

2.15 × 10 4

6

1550

5.5

8

9.43

1.62

−0.91

310

M5-Ra4-l5.5

52.15 × 10 4

6

1550

5.5

8

12.57

3.19

−0.37

290

M5-Ra5-l5.5

55

2.15 × 10

6

1550

5.5

8

12.46

3.02

−0.19

285

ET9-l5.5

100

2.15 × 10 6

9

1550

5.5

8

3.28

0.60

−0.87

85

l5.5-Fcrst16

100

2.15 × 10

6

1550

5.5

16

1.74

0.890

−0.85

185

6

9 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Table 3
Continued
Case #

M*

Ra𝐴𝐴

𝐴𝐴𝑇𝑇∗

Tp (K) 𝐴𝐴1/l*∗
𝐴𝐴crst

ΔRpeak(km)tΔR (Gyr)

χGyr (km Gyr −))

Lϕ
(km)

l5.5-Fcrst32

100

2.15 × 10 6

6

1550

5.5

32

−0.21

1.30

−0.91

315

ET3-l5.5

100

2.15 × 10 6

3

1550

5.5

8

2.99

0.80

−1.03

20

ET3-l1.5

6100

2.15 × 10

3

1550

1.5

8

4.00

1.11

−0.82

160

M5-ET3-l5.5

52.15 × 10 6

3

1550

5.5

8

10.40

2.34

−2.00

180

Ra4-ET3-l5.5

100

2.15 × 10

4

3

1550

5.5

8

3.88

1.12

−1.67

25

Ra4-ET9-l5.5

100

2.15 × 10 4

9

1550

5.5

8

3.00

0.66

−1.17

150

ET9-l1.5

100

2.15 × 10 6

9

1550

1.5

8

2.19

0.85

−1.39

160

5

2.15 × 10 6

9

1550

5.5

8

8.75

1.62

−1.91

185

Tm14-l5.5-Fcrst16

100

2.15 × 10

6

6

1400

5.5

16

2.46

1.06

0.07

225

Tm17-l5.5-Fcrst16

100

2.15 × 10 6

6

1700

5.5

16

2.93

0.45

−1.10

50

ET3-l5.5-Fcrst16

100

2.15 × 10

6

3

1550

5.5

16

1.96

0.83

−0.94

50

ET9-l5.5-Fcrst16

100

2.15 × 10 6

9

1550

5.5

16

1.80

0.80

−1.27

195

Ra4-l5.5-Fcrst16

100

2.15 × 10 4

6

1550

5.5

16

2.14

1.46

−1.10

100

Ra5-l5.5-Fcrst16

100

2.15 × 10 5

6

1550

5.5

16

2.06

1.58

−0.86

115

l1.5-Fcrst16

100

2.15 × 10

6

6

1550

1.5

16

0.83

1.45

−0.97

90

l3.5-Fcrst16

100

2.15 × 10 6

6

1550

3.5

16

1.88

0.98

−0.69

185

M5-l5.5-Fcrst16

5

2.15 × 10 6

6

1550

5.5

16

4.37

1.79

−0.71

240

M20-l5.5-Fcrst16

20

2.15 × 10 6

6

1550

5.5

16

3.63

1.07

−0.60

260

M50-l5.5-Fcrst16

50

2.15 × 10

6

6

1550

5.5

16

2.65

0.80

−1.05

205

No-conv-HPEs-Mass-tr

100

2.15 × 10 2

6

1550

4.5

8

−0.65

1.11

−0.71

195

No-HPEs-Mass-tr

100

2.15 × 10

6

6

1550

4.5

8

−4.12

0.49

−0.80

20

No-Mass-tr

100

2.15 × 10 6

6

1550

4.5

8

−5.01

0.55

−1.30

15

vl-change100

2.15 × 10

6

6

1550

4.5

8

3.67

0.70

−1.08

65

melt-dep100

2.15 × 10 6

6

1550

4.5

8

2.57

0.57

−0.71

35

melt-dep-Ra4

1002.15 × 10

4

6

1550

4.5

8

6.84

0.91

−0.28

305

IBC-dep100

2.15 × 10 6

6

1550

4.5

8

3.20

0.75

−0.91

20

rl-100 km

100

2.15 × 10

6

6

1550

4.5

8

2.45

0.68

−1.23

85

beta0202 (β = 0.202)

100

2.15 × 10 6

6

1550

4.5

8

5.21

0.70

−0.85

100

beta0067 (β = 0.067)

100

2.15 × 10 6

6

1550

4.5

8

−0.60

2.38

−0.85

20

M5-ET9-l5.5

Note. ΔRpeak and tΔR stand for the magnitude and timing of the peak of radial expansion, respectively; χGyr the contraction
rate for the past 1 Gyr obtained by least squares method; Lϕ the shallowest depth levels of partially molten regions other thanthat in the initial state.

negative buoyancy (see the arrow in Figure 2d). The return flow of this foundering produces magma that is notenriched in HPEs and the IBC components by decompression melting (Figures 2b–2d for 1.44 Gyr; see also
Movies S2–S4 for around 1.28 and 1.44 Gyr).
We will describe more about melt-fingers and partially molten plumes in Section 3.1.3, below.
3.1.2. Radius Change
Figure 4 shows how the radius changes with time owing to the mantle evolution shown in Figures 2 and 3. The
radius of the Moon changes with time for two reasons, thermal expansion/contraction, and melting of the mantle(see Equation 12), as indicated by the blue and red lines in the figure, respectively. The blue line in Figure 4 showsthat the Moon thermally expands by 0.4 km in the early history as the mantle is heated up and then contracts until
the end of the calculation as the mantle is cooled (Figure 3b). On the other hand, the red line indicates the Moon
U ET AL.

10 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

expands by 3 km for the first 0.7 Gyr owing to widespread partial meltingof the mantle caused by melt-fingers and partially molten plumes (the peak[3] in Figure 4). After that, the Moon gradually contracts with time as the
mantle solidifies. Note that the contraction is not monotonous: slight expansion occurs several times when partially molten plumes develop (see [4] and[5] in Figures 2 and 4). As a whole, the black line in Figure 4 indicates that
the Moon radially expands by about 3 km during the first 0.7 Gyr and thencontracts at the rate of around −1.0 km Gyr −1 until the end of the calculated
history (Table 3).
3.1.3. Melt-Fingers and Partially Molten Plumes

Figure 4. The history of radius change of planet calculated in the referencecase (Case Ref). The radius change ΔR is defined by Equation 12; the blue
and red lines indicate the contribution of thermal expansion/contraction andmelting, respectively, to the total radius change (the black line). The numbers[1]–[5] correspond with those of Figure 2.

To understand why melt-fingers grow from the partially molten region at
depths around 0.4 Gyr in Figures 2a and 2b, we delineate the velocity field of
the matrix around a nascent finger in the frame for 0.38 Gyr of Figure 5 (seethe bump indicated by the yellow arrow). This frame shows that the matrix
expands around the fingertip to let it grow upward. To show why this expansion occurs, we decompose the flow field U into the components driven by
volume change of the matrix U vc and that driven by the buoyancy force U by as
(14)𝐔𝐔 = 𝐔𝐔𝐯𝐯𝐯𝐯 + 𝐔𝐔𝐛𝐛𝐛𝐛 ∶we calculated U vc from Equations 6, A1 and A2 with Ra = 0 and the originalvalue of M*, while U by from these equations with M* = 0 and the originalvalue of Ra. We confirmed that the flow field U in Figure 5b for 0.38 Gyris close to that of U vc, implying that the matrix expansion around the fingertip is caused by the injection of upward migrating magma; upward migration of magma with respect to the coexisting matrix is the cause of growthof melt-fingers. (We recalculated melt-fingers on a mesh with twice higher
resolution and confirmed that fingers still grow; melt-fingers are not an artifact of numerical instability.)
After 100 Myr from the development of melt-fingers (Figure 5 for 0.48 Gyr),some stems of the melt-fingers become thicker with time as the fingertips
migrate upward from the partially molten layer at depth (see the yellow dashedcircle in the figure). In the thickened area, both the matrix and the melt migrate
upward (Figure 5). We decomposed the flow field U (see Equation 14) andfound that U around the stem is dominated by Uby, implying that the thickened
area is a partially molten plume driven by the buoyancy of the melt.
3.2. The Occurrence of Melt-Fingers and Partially Molten Plumes
To understand under what condition melt-fingers and partially molten
plumes observed in the reference case grow, we calculated the model at various values of the reference permeability M*, and the Rayleigh number Ra.

Figure 5. Snapshots of the (a) melt-content ϕ and (b) matrix-velocity field in
the reference case (Case Ref). In (b), the red and blue colors show the regions
where the convective flow points upward and downward, respectively; thearrows express the direction of convective flow but not its magnitude. The
contour lines show the ϕ-distribution with the contour interval of 0.1 startingfrom 0.

U ET AL.

The reference permeability M* influences the growth of melt-fingers; thesnapshots of Case M5 (M* = 5 and other parameters fixed at their default
values) illustrate that the melt-fingers do not occur at the lower reference
permeability of this case (Figure 6a); we confirmed that U vc is negligibly
small in this case. Instead of melt-fingers, partially molten plumes that arebroader than those observed in the reference case rise to the depth levels of
200 km by around 1.7 Gyr (Figure 6a for 1.72 Gyr and Table 3). The plumes
laterally extend and form a continuous layer of partially molten materials(Figure 6a). The partially molten layer then solidifies upon cooling fromthe surface boundary; the layer, however, still remains in the mid-mantle at
4.4 Gyr (Figure 6a). We observed that melt-fingers grow only in the cases
calculated at M* ≥ 50 (see Figure S2 in Supporting Information S1).
11 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 6. Snapshots of the distributions of melt-content ϕ calculated in cases with (a) a reduced reference permeability(Case M5) and (b) a reduced Rayleigh number (Case Ra2.15e4). The numbers [1]–[5] correspond with those of Figures 14a
and 14e.

At M* = 100, we also calculated several models at various values of Ra in the range from 2.15 × 10 2 to 2.15 × 10 6
(Table 3); we found that melt-fingers develop regardless of the values of Ra, implying that buoyancy does not play
any role in the growth of melt-fingers (Figure 6b). This result reinforces the above conclusion that melt-fingers
develop because of matrix expansion around fingertips (Section 3.1.3).
We further calculated Case No-conv-HPEs-Mass-tr where the convection is sluggish (Ra = 2.15 × 10 2) and transport of HPEs and the IBC components is negligible: the partition coefficient of HPEs D is 1 (see Equation A11 inAppendix A), and the model starts from a compositionally uniform mantle with the bulk composition of ξb = ξe.
We found that melt-fingers develop even in this case, implying that the development of melt-fingers does notdepend on these parameters.
In contrast to melt-fingers, the growth of partially molten plumes depends on the Rayleigh number. At
Ra = 2.15 × 10 4 that is lower than Ra of the reference case by a factor of 100 (Case Ra2.15e4), only melt-fingers
develop, and partially molten plumes do not grow (Figure 6b); we confirmed that U by is negligibly small inthis case. Melt-fingers then grow upward to the depth level of around 75 km by around 1.0 Gyr (Figure 6b andTable 3) and then expand laterally, to make the most part of the mantle partially molten (Figure 6b for 1.04 Gyr).
After that, the partially molten region shrinks monotonously with time owing to conductive cooling from thesurface boundary (Figure 6b for 2.80 Gyr). We found that partially molten plumes grow only in the cases calculated at Ra ≥ 7.15 × 10 5 (see Figure S3 in Supporting Information S1). This result reinforces that partially molten
plumes are driven by their buoyancy.
The melt-fingers observed here are different from the pencil-shaped upwellings reported by Whitehead andLuther (1975) despite the similarity of shape. These upwellings develop when there is a viscosity contrast between
the upwellings and the surrounding fluid, whereas our melt-fingers grow even when there is no viscosity contrast.
3.3. The Effects of More Sophisticated Rheology
In the cases described above, we assumed that the viscosity of the mantle depends only on the temperature (seeEquation 11) for simplicity. However, some earlier studies show that the viscosity depends on the contents ofmelt and the IBC components. In this section, we considered the reduction of viscosity by melting (e.g., Hirth &Kohlstedt, 2003; Mei et al., 2002; Scott & Kohlstedt, 2006) and the IBC components (e.g., Dygert et al., 2016; H.
Li et al., 2019; N. Zhang et al., 2017). The viscosity depends on the melt-content ϕ and the ratio of the olivinerich materials to the IBC components Fol = (ξb − ξe)/(1 − ξe) as well as the temperature as][𝜂𝜂 = [𝐹𝐹ol 𝜂𝜂0 + (1 − 𝐹𝐹ol )𝜂𝜂IBC ]exp 𝐸𝐸𝑇𝑇 (𝑇𝑇ref − 𝑇𝑇 ) − 𝛼𝛼𝝓𝝓 𝜙𝜙 ,
(15)U ET AL.

12 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 7. The same as Figure 2 but for Case melt-dep, where viscosity reduction by melting is considered.

where, αϕ = 26 is the melt-weakening factor (Mei et al., 2002), and the viscosity reduction by melting is truncatedat its value at ϕ = 0.4; ηIBC = 7.5 × 10 18 Pa s is the viscosity of the IBC components (Dygert et al., 2016). Notethat Fol depends on the composition ξb, which changes with time by magmatism. We calculated two cases wherethe ϕ-dependent viscosity (Fol = 1 is assumed in Equation 15; Case melt-dep) and the IBC-dependent viscosity(αϕ = 0 is assumed in Equation 15; Case IBC-dep) are considered.
Figures 7 and 8 where the results of Case melt-dep are presented show that the effects of ϕ-dependent viscosity
do not significantly affect the overall features of mantle evolution observed for Case Ref. The partially molten
region in the deep mantle vertically extends at the beginning of the calculated history of the mantle (Figure 7 for
0.32 Gyr). Partially molten plumes then develop at the top of the partially molten region (Figure 7 for 0.36 Gyr)and rise to the uppermost mantle. The plume activity continues over billions of years (Figure 7 for 0.64–2.40 Gyr).
In this case, the Moon expands by 2.6 km for the first several million years and then contracts at the rate of around−0.7 km Gyr −1 until the end of the calculated history (Figure 8). These overall features are the same as thoseobserved for Case Ref.
However, melt-fingers observed for Case Ref do not grow in this case. To clarify the reason for the absence ofmelt-fingers, we decomposed the flow field U around the partially molten bump in Figure 9b into the components
driven by volume change of the matrix U vc and that driven by buoyancy force U by as we did in Section 3.1.3 ofthe main article and found that U by ≫ U vc holds. U by represents the rotation of matrix within the bump observedin the figure. This rotation is enabled by the reduction in viscosity along the head of the bump that accommodates the strong shear caused by the rotational flow. The strong buoyancy-driven flow U by overshadows U vc thatinduces a melt-finger.
U ET AL.

13 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

In contrast, we found no conspicuous effects of the IBC-dependent viscosity
on numerical results in Case IBC-dep: the overall features of the mantle evolution observed in the reference case were not affected by the IBC-dependent
viscosity, and both melt-fingers and partially molten plumes developed as in
the reference case. The effects of the IBC-dependent viscosity are negligible
in our study because the viscosity variation by this dependence is only anorder of magnitude at most for the variation in the content of the IBC components realized in the numerical models.
3.4. The Effects of Density Inversion Between Melts and Matrix
The density of the solid phase with ξs = 1 is always denser than that of themelt phase with ξl = 0.1 in the cases described so far (see Equations 2 and 3,
and the red line in Figure 10a). This occurred because we estimated Δvl/
v0 in Equation 3 from the solidus curve of lunar mantle materials (Garcia
et al., 2011, 2012; Katz et al., 2003). (see Equations 9 and 10 and Figure 10b
as well as the assumed values of Δρ∞/ρ0 = 0.005, Δρzero/ρ0 = 0.22, and
λ = 16.42 GPa in Equation 5.) However, some studies suggest that IBC-rich
Figure 8. The same as Figure 4 but for Case melt-dep.
magma can be denser than the coexisting matrix in the deep mantle ofthe Moon (e.g., Sakamaki et al., 2010; van Kan Parker et al., 2012; Xuet al., 2022). To see how this possible density inversion affects our numericalresults, we calculated the two additional cases (vl-change and beta0202). In Case vl-change, which is calculated
at ρ∞/ρ0 = 0.002, Δρzero/ρ0 = 0.35, and λ = 4.56 GPa (see Equation 5), melt phase is denser than that of thecoexisting at radius r less than around 800 km (see the blue line and P1 in Figure 10a). Note that this change
in parameter values slightly raises the solidus temperature as shown in Figure 10b. To see if this higher solidus
temperature affects our numerical results in Case vl-change, we also calculated Case beta0202 (see the black line
in Figure 10a), where β is 0.202 (the density of IBC (end-member B) is 3,965 kg m −3) but the default values are
assumed for Δρ∞/ρ0, Δρzero/ρ0, and λ.
Figures 11 and 12 where the results of Case vl-change are presented show that the effects of density inversion
between melts and coexisting matrix do not significantly affect the overall features of mantle evolution observedfor Case Ref (see also Movies S6–S10). For the first few hundred million years, magma generated in the deep
mantle sinks to the core-mantle boundary (CMB), making the mantle more enriched in HPEs and the IBC
components at its base (Figure 11 for 0.52 Gyr). Partially molten plumes, however, eventually grow along the
top of the partially molten region and ascend to carry the enriched materials at the base of the mantle to theuppermost mantle over billions of years (Figure 11 for 0.64–2.16 Gyr). The Moon radially expands by 3.7 km
during this time (see Figure 12). The plume activity then declines with timeowing to cooling of the mantle (Figure 11 for 2.16–4.00 Gyr) and the Moon
accordingly contracts at a rate of 1 km Gyr −1. These overall features are
the same as those observed in Case Ref. (Note, however, that melt-fingers
observed in Case Ref do not grow in this case because of the density inversion. The absence of melt-fingers does not affect the overall evolutionary
history.)
We found that the overall features of mantle evolution calculated in Case
beta0202 are also the same as those calculated in Case vl-change. The slightdifference in the solidus curve shown in Figure 10b does not influence thenumerical results.

Figure 9. The same as Figure 5 but for Case melt-dep. In (b), the red andblue colors show the regions where the convective flow points upward anddownward, respectively; the arrows express the direction of convective flow
but not its magnitude. The contour lines show the ϕ-distribution with thecontour interval of 0.2 starting from 0.

U ET AL.

It may look counterintuitive that the partially molten regions in the deep
mantle eventually migrate upward as partially molten plumes despite thedensity inversion between melts and the coexisting matrix there. This density
inversion does not influence the overall upwelling flow, because melting of
the mantle always causes a decrease in the bulk density (see Movie S10);
the density inversion just lets magma migrate downward within upwelling
partially molten plumes in the deep mantle.
14 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 10. An illustration of (a) the density of melt phase with ξl = 0.1, and (b) the solidus assumed in Case vl-change where the density inversion between solid andmelt phases is considered. In (a), the density difference between the solid phase and the melt phase is generally consistent with that suggested by some experimental
studies (van Kan Parker et al., 2012; Xu et al., 2022). In (b), the black dashed line is the solidus curve for the lunar mantle materials (Garcia et al., 2011, 2012; Katz
et al., 2003).

Figure 11. The same as Figure 2 but for Case vl-change, where the magma generated in the deep mantle is denser than thecoexisting matrix.

U ET AL.

15 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

3.5. The Effect of the Initial Distribution of Mantle Composition
To see how plume magmatism continues for billions of years in the referencecase (Figure 2), we further calculated Case No-Mass-tr where we assumed
a compositionally uniform mantle with the bulk composition ξb = ξe in thewhole mantle (Table 3). In this case, most of partially molten plumes ascend
to the uppermost mantle by around 0.8 Gyr, and plume activity declines afterthat (Figure 13). The early decline of plume activity results from the early
extraction of HPEs from the deep mantle where magma is mostly generated.
In contrast, HPE-extraction by early plume magmatism is not so efficient in
the reference case where the deep mantle is assumed to be enriched in thedense IBC components in the initial condition (see Section 3.1.1). The duration time of plume activity depends on the initial distribution of the dense
IBC components.
3.6. Dependence of Radial Expansion on Model Parameters

Figure 12. The same as Figure 4 but for Case vl-change.

We further carried out numerical experiments to show how the radial
expansion-history depends on the model parameters: the reference permeability M*; the Rayleigh number Ra; the sensitivity of viscosity to temper𝐴𝐴ature 𝐴𝐴𝑇𝑇∗ ; the potential temperature Tp; the thickness of the overturned layer∗
l*; the initial concentration ratio of the HPEs in the crust to the mantle
𝐴𝐴𝐴𝐴crst(Figure 14 and Table 3).

Among the parameters, the reference permeability M* influences the expansion-history most strongly(Figure 14a); a lower M* leads to a later peak of expansion with a larger amplitude. The larger expansion is due
to a more slowly migrating magma that retains HPEs in the deep mantle for a longer period, and it allows moreheat to build up in the mantle (see Figure 6a and Table 3).
∗
The initial ratio of HPE-concentration in the crust to that in the mantle
also substantially influences the𝐴𝐴𝐴𝐴crst∗magnitude of radial expansion (Figure 14b); a higher𝐴𝐴𝐴𝐴crst reduces radial expansion because the mantle is more
depleted in HPEs and less magma is generated in the deep mantle.

The potential temperature Tp affects the timing and, to some extent, the amplitude of radial expansion as shown in
Figure 14c. A higher Tp leads to an earlier peak of expansion because it implies an earlier generation of partiallymolten regions in the deep mantle and an earlier extension of the regions into the uppermost mantle. A higher Tp
also leads to a slightly larger amplitude of the early expansion at Tp above 1550 K (Figure 14c). At this high Tp,
the magma is generated immediately in the deep mantle and spreads into the uppermost mantle before most of an
initial partially molten layer in the uppermost mantle is not completely cooled. The contraction due to solidification of the initial melt is smaller, and the earlier expansion becomes larger.
The thickness of the overturned layer l* has an effect on the beginning of radial expansion (Figure 14d). A higher
1/l* where the deep mantle is more enriched in HPEs induces an earlier beginning of the expansion. Note thatin Case l0 where HPEs and the composition ξb are uniform in the whole mantle (1/l* = 0), the amplitude of

Figure 13. Snapshots of the distributions of internal heating rate q calculated in a case where the composition ξb is kept
uniform with the eutectic composition ξe (Case No-Mass-tr). The contour lines show the distribution of melt-content with
contour level of 0 and 0.1.

U ET AL.

16 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 14. Plots of radius change against time calculated at various values of (a) the reference permeability M*, (b) the initial crustal fraction of the heat-producing
elements
𝐴𝐴
𝐴𝐴crst ∗, (c) the potential temperature Tp, (d) the thickness of the overturned layer l*, (e) the Rayleigh number Ra, and (f) the sensitivity of viscosity totemperature
𝐴𝐴
𝐴𝐴𝑇𝑇∗ . The dotted lines show the reference case (Case Ref) presented in Figures 2–4. In (a) and (e), the numbers [1]–[5] correspond with those of Figures 6aand 6b, respectively.

radial expansion is smaller than in cases where the overturned layer is considered (Figure 14d). This is because
enough magma generation does not occur to cause the early expansion in Case l0. (Case l0 is different from Case
No-HPEs-Mass-tr where the initial HPE-distribution is preserved throughout the calculated history (the sameas Case No-conv-HPEs-Mass-tr, but Ra = 2.15 × 10 6). HPE-distribution is assumed to be uniform in the initialcondition but becomes heterogeneous by magma-transport in Case l0.)
U ET AL.

17 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Figure 15. An illustration of the thermal and structural history of the lunar mantle inferred from our numerical study (CaseRef). The yellow color stands for melting.

Conversely, although the horizontal averages of the temperature-distributions depend on the values of Ra
𝐴𝐴 and 𝐴𝐴𝑇𝑇∗(Figure 3a), the dependences of expansion-history on the Rayleigh number Ra and the sensitivity of viscosity totemperature
𝐴𝐴
𝐴𝐴𝑇𝑇∗ are negligible as shown in Figures 14e and 14f.

4. Discussions
Figure 15 illustrates how the mantle evolves in the reference case shown in Figures 2–5. The partially molten region
in the uppermost mantle shrinks with time, whereas that in the deep mantle expands on the earliest stage of the calculated history (Figure 15a). Melt-fingers then develop along the top of the partially molten region at depth, and thefingers extend upward (Figure 15b). After the growth of those, partially molten plumes driven by melt buoyancy riseto the uppermost mantle to cause plume magmatism (Figure 15c). The Moon expands in the stage of Figures 15a–15c
because of melting of the mantle. These partially molten plumes transport HPEs and the IBC components from thedeep mantle to the uppermost mantle. The plume magmatism then declines with time after around 2 Gyr of the calculated history because of the depletion of HPEs in the deep mantle (Figure 15d). In the latter period of the calculated
history, the Moon gradually contracts by cooling and solidification of the partially molten mantle. The reference casefits in best with the observed history of the Moon among the models calculated here as we will discuss in Section 4.2.
4.1. Comparisons With Earlier Models
4.1.1. The Mantle Evolution Caused by Magmatism and Convection
A comparison with earlier classical 1-D models of lunar thermal history in the literature shows the crucial role that heattransport by migrating magma and mantle convection plays in the reference case. In earlier 1-D models where onlyinternal heating and thermal conduction are considered, the temperature in the deep mantle monotonously increases
to the solidus temperature due to internal heating, while the lithosphere monotonously thickens with time owing tocooling from the surface boundary throughout the calculated history (e.g., Solomon & Chaiken, 1976; Solomon &Toksöz, 1973; Toksőz & Solomon, 1973; Wood, 1972). These models show that the deep mantle becomes extensively
molten and the lithosphere becomes as thick as around 600 km at present. In our reference case where heat transport
by migrating magma and mantle convection is also considered, however, the temperature is below the solidus in mostpart of the deep mantle, and the thickness of the lithosphere is by less than around 350 km at 4.4 Gyr (Figure 3a).
A comparison with earlier 2- or 3-D models of the lunar thermal history where mantle convection also is considered, on the other hand, shows the crucial roles that magma plays in heat transport in the convecting mantle. The
vigor of thermally driven mantle convection is controlled by the distance from the threshold for the onset of thermal convection on the plane of the Rayleigh number Rad and the viscosity contrast ηtop/ηbot shown in Figure 16(Yanagisawa et al., 2016). As shown in the figure, our reference case (Case Ref) is calculated under the condition
that thermally driven mantle convection is more sluggish than that of the earlier models of lunar mantle convection
Z-E100, Z-E200 (N. Zhang et al., 2013a), and Zi-M650 (Ziethe et al., 2009). The lithosphere at 4.4 Gyr in ourreference case is (see the black line in Figure 3a), however, thinner than that in these models (see e.g., Figure 3c in
N. Zhang et al. (2013a)). Even in the models where thermally driven mantle convection does not occur at all (CaseRa2.15e4), the lithosphere (see the red line in Figure 3a) is substantially thinner than that in the earlier models KS-S
(Konrad & Spohn, 1997; Spohn et al., 2001) and Ln (Laneuville et al., 2013). These differences arise because, in
our models, the uppermost mantle is kept hot by melt-fingers and partially molten plumes rather than totally solid
plumes. Heat transport by melt-fingers and partially molten plumes is an essential part of our models. (Note that eventhe lithosphere of 1-D spherically symmetric mantle of an earlier model where heat transport by migrating magma
U ET AL.

18 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

is considered (U et al., 2022) is already thinner than that of the earlier models(Konrad & Spohn, 1997; Laneuville et al., 2013; Spohn et al., 2001). For a morequantitative estimate of the contribution of heat transport by melt-fingers andpartially molten plumes to thinning of the lithosphere, further calculations in a
3-D spherical mantle are necessary.)

Figure 16. Plots of the critical Rayleigh number for onset of thermal
convection in the 3-D spherical mantle heated from the core-mantle boundary(CMB), taken from Figure 8b in Yanagisawa et al. (2016). The viscosity is
assumed to depend on the temperature and the sensitivity of the dependence
is measured by the viscosity contrast between the surface boundary ηtop andthe CMB ηbot. Z-E100 and Z-E200 (the yellow circles) correspond to the case
H50E100v5e20 and H50E200v5e20 of N. Zhang et al. (2013a); the brown
arrow (Ln) by Laneuville et al. (2013); the red circle (Zi-M650) by the case
M650 of Ziethe et al. (2009); the blue rounded rectangle (KS-S) by Konrad
and Spohn (1997) and Spohn et al. (2001). The black circles (Ref, Ra2.15e4,
and ET3) also correspond to the Case Ref, Ra2.15e4, and ET3. Note that somemodels in this figure use the Arrhenius law for the temperature dependence,
rather than the exponential (Frank-Kamenetskii) dependence considered inour study. The difference is not important in the convecting part of the mantle
beneath the stagnant lid.

In order to understand the evolution of the compositionally stratified mantle
predicted from the hypotheses of the magma ocean and mantle overturn, masstransport by migrating magma and magma-driven convective flow of the mantle
is indispensable. The lunar mantle is expected to have been compositionallystratified with a layer enriched in the compositionally dense IBC components
at the base of the mantle after putative crystal fractionation in the magma
ocean and subsequent mantle overturn (e.g., Elkins-Tanton et al., 2011; Hess &Parmentier, 1995; Snyder et al., 1992). Some mantle convection models suggestthat the basal layer eventually rises as upwelling plumes owing to its thermal
buoyancy as the layer is heated by HPEs (e.g., Stegman et al., 2003; N. Zhang
et al., 2017; W. B. Zhang et al., 2022; Zhong et al., 2000). In these models,
however, the excess compositional density of the basal layer with respect to theoverlying olivine-rich mantle is less than that suggested by some recent models
of lunar mantle overturn (e.g., Yu et al., 2019; Zhao et al., 2019), greater than
160 kg m−3 (see Figure 6 in H. Li et al. (2019)). At this large compositional
density contrast, the basal layer is convectively stable and does not ascend as
solid plumes by thermal buoyancy alone, as inferred from earlier laboratory
experiments of mantle convection (see Figure 1 in Le Bars and Davaille (2004)).
Although the compositional density contrast in our reference case (l* = 4.5) isaround 180 kg m−3 (Figure 1), a large fraction of the IBC-rich materials in thelayer is extracted and transported to the uppermost mantle by melt-fingers andby partially molten plumes driven by melt buoyancy (Figure 2); mass transport by migrating magma and magma-driven convective flow of the mantle are
crucial for understanding the structural evolution of the mantle.

Our model also shows that a careful modeling of HPE transport is important
to correctly understand the thermal and volcanic history of lunar mantle. The
lunar volcanism has continued until 1–2 Gyr ago (e.g., Hiesinger et al., 2000;
Whitten & Head, 2015), and some earlier mantle convection models
conclude that the volcanism has continued for such a long period because
partially molten regions remained in the cooling mantle for billions of years(Konrad & Spohn, 1997; Spohn et al., 2001; Ziethe et al., 2009). In models
where the uppermost mantle is locally more enriched in HPEs, the partiallymolten regions observed in there persist for more than 3 Gyr (Laneuvilleet al., 2013, 2014, 2018). Volcanism, however, extracts HPEs from the mantle leading to the decline of subsequent volcanic activity (Cassen & Reynolds, 1973; Cassen et al., 1979; Ogawa, 2014, 2018a). In particular,
Ogawa (2018a) suggests that partially molten regions disappear within 2 Gyr since the beginning of the calculated
history, too short to be a model of the lunar mare volcanism, owing to the extraction of HPEs by magmatism. In ourmodels, in contrast, magmatism continues for a much longer time despite that HPE transport by magma is considered (Figures 2 and 3) because a compositionally dense IBC-enriched layer is assumed at the base of the mantle in
the initial condition. When HPE- and IBC-enriched materials in the basal layer are transported by melt-fingers andpartially molten plumes, the magma often solidifies on the way to the surface and sinks again to the deep mantle(Figure 2; see also Movies S1–S5). The initial mantle stratification is necessary for magmatism to continue long.
Indeed, in Case No-Mass-tr where the composition ξb is kept uniform (ξb = ξe), the mantle becomes completely
solid, and magmatism declines much earlier than that in the reference case (Figure 13).
Because of the explicit implementation of a model of magma generation and migration into that of mantle
convection, our models allow us to infer volcanic history directly from the calculated history of mantle melting. Earlier models have discussed volcanic history based on the calculated distribution of partially molten
regions, especially the depth of the top of the regions and the calculated rate of magma generation in the

U ET AL.

19 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

mantle (e.g., Laneuville et al., 2018; Solomon & Toksöz, 1973; Spohn et al., 2001; U et al., 2022; Wieczorek& Phillips, 2000; Wood, 1972). These models assume that the magma ascends from partially molten regions as
deep as 200–800 km to the crust. However, it is not clear if magma eruption is well correlated with the depth
and magma generation rate. Magma may not be able to make its way to the surface when the stress state in thelithosphere is horizontally compressive and hence when the Moon is contracting (Solomon, 1986; Solomon& Head, 1979). In our models, in contrast, magma rises directly to the base of the crust by melt-fingers andpartially molten plumes (Figures 2 and 3). Several issues, however, still remain on magma migration in the crustand the uppermost mantle. It is unclear how mare basaltic magma can ascend through the crust that is not denser
than the magma (e.g., Head & Wilson, 1992). The detail of formation of dikes in the uppermost mantle is also
important for understanding the magma ascent (e.g., Head & Wilson, 2017; Wilson & Head, 2003, 2017). Theseissues may be important for also constructing a more refined model of thermal history of the Moon as noted by
Lourenço et al. (2018): the thermal history can depend on the ratio of extrusive to intrusive volcanism, which
can substantially depend on the porosity, thickness, and density of the crust (e.g., Head & Wilson, 2020; Morota
et al., 2009; Solomon, 1975; Taguchi et al., 2017). In future studies, a more refined modeling of magma migration through the crust and the mantle is needed. (As we have already described in Section 2 above, the possible
density inversion between the solid phase and the melt phase in the mantle does not affect the overall features ofmantle evolution; see Section 3.4.)
For a more realistic simulation of the evolution of lunar mantle, it is essential to extend the model to a 3-D
spherical shell (e.g., Laneuville et al., 2013; N. Zhang et al., 2017). 2-D annular models of mantle convection
tend to predict a higher average temperature in the mantle than 3-D spherical models do, especially when thecore size is small, as is the case for the Moon (Guerrero et al., 2018). Although it is computationally challenging, modeling magma generation and migration in a 3-D spherical mantle is a promising avenue for future
research.
The crustal dichotomy between nearside and farside is also a long-standing issue in studies of mantle dynamics
in the Moon (e.g., Cho et al., 2012; Jolliff et al., 2000; Lawrence et al., 2000). Some studies suggest that the lunar
dichotomy is caused by an exogenous agent, such as the South Pole-Aitken impact (e.g., Arai et al., 2008; Joneset al., 2022; N. Zhang et al., 2022). In the future study of 3-D spherical mantle, it is important to assume theinitial condition in which the uppermost mantle in the nearside is more enriched in HPEs than that in the farside(Laneuville et al., 2013, 2014, 2018; Wieczorek & Phillips, 2000).
4.1.2. The Radial Expansion/Contraction
Our reference case shows that the volume change of the mantle by melting is a key for understanding the radiuschange of the mantle. To cause the observed early expansion of the Moon, classical earlier models suggest thatthe initial temperature in the deep mantle was less than 1200 Kirk & Stevenson 1989; Solomon, 1986; Solomon& Chaiken, 1976). This upper limit is, however, substantially lower than that expected from earlier models ofmantle overturn that start from giant-impact hypotheses in the literature (e.g., Canup, 2004; Lock et al., 2018;
Pritchard & Stevenson, 2000). These overturn models suggest that the initial temperature in the deep mantle isapproximately 1800–1900 K (e.g., Boukaré et al., 2018; Elkins-Tanton et al., 2011; Hess & Parmentier, 1995;
H. Li et al., 2019). When such a high initial temperature is assumed, the early expansion occurs only nearside(Laneuville et al., 2013) or is much smaller than the observed expansion (N. Zhang et al., 2013a). In a model
where the blanket effect in the crust is taken into account (Ziethe et al., 2009; N. Zhang et al., 2013b), thermal
expansion of more than 1 km does occur but continues for longer than 1 Gyr, too long to account for the radial
expansion of the Moon (e.g., Andrews-Hanna et al., 2013); the thermal expansion continues for such a long
period even in our reference case, as shown by the blue line in Figure 4. In our model, however, global expansion
occurs by a few kilometers within the first 0.7 Gyr of the calculated history because of melting of the mantle. The
radial expansion by melting occurs earlier and is larger than the thermal expansion (Figure 4), suggesting that
mantle melting dominates the radial expansion/contraction history of the Moon.
The radial expansion/contraction history depends on the spatial dimensionality. In a spherically symmetric model
where the distribution of melt-content depends only on the radial coordinate r, the amplitude of radial expansion
is around 1.2 km (see Figure 2 in U et al. (2022)). In contrast, the amplitude in our 2-D model is around 3 km
(Figure 4), substantially larger than that in the 1-D model. To predict more quantitatively the radial expansion/contraction history of the Moon, it is necessary to develop a model in a 3-D spherical shell (Laneuvilleet al., 2013; N. Zhang et al., 2013a, 2013b).
U ET AL.

20 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

4.2. Comparison With the Observed Features of the Moon
4.2.1. The Radial Expansion/Contraction
The Moon globally expands by a few kilometers for the first several hundred million years in our referencecase owing to the volume change caused by melting of the mantle (Figure 4). The timing and amplitudeof radial expansion are consistent with those of the Moon inferred from its gravity field in the literature(Andrews-Hanna et al., 2013, 2014; Liang & Andrews-Hanna, 2022; Sawada et al., 2016). After the radial
expansion, the Moon begins to radially contract at around 1 Gyr after the start of the calculation (Figure 4).
The timing of contraction is consistent with the beginning of compressive tectonics on the Moon (Frueh
et al., 2023; Yue et al., 2017). The Moon in our model contracts at a rate of approximately −1.0 km Gyr −1
in the past billion years (Table 3), which is also consistent with the estimates obtained from observationsof fault scarps on the Moon (e.g., Clark et al., 2017; Klimczak, 2015; van der Bogert et al., 2018; Watters
et al., 2010). As a whole, the calculated history of radius change of the reference case is consistent with thatof the Moon.
Among our models presented in Figure 14, the radius changes calculated at M* ≥ 50, Tp ≈ 1600 K,
𝐴𝐴 and 𝐴𝐴crst ∗ ≤ 16
are consistent with the observed one. When the reference permeability M* is lower than 50, the amplitude of radial
expansion is much larger than the estimate for the Moon because generated magma in the deep mantle stays there
for a long time (Figures 6 and 14a). In the cases where the initial concentration ratio of the HPEs in the crust to themantle𝐴𝐴
𝐴𝐴crst ∗ is larger than 16 (the average heating rate in the mantle is less than 7.6 pW kg −1 in the initial condition), sufficient partially molten regions cannot develop to cause the early expansion (Figure 14b). Our models also
show that the timing and amplitude of radial expansion strongly depend on the presence of overturned IBC components enriched layer l* (Figure 14d and Table 3). In Case l0 where the overturned layer is not assumed (1/l* = 0),
indeed, the timing and amplitude of calculated radial expansion are later and smaller than the observed expansion.
Our models suggest that a substantial fraction of the mantle should have been solid at the beginning of the lunar
history for the observed early expansion to take place (Figures 1a and 14c). This result urges us to rethink the
lunar formation process. Many of previous simulations of the lunar formation suggest that most materials were hot
enough to evaporate immediately after the Moon-forming giant impact (e.g., Canup, 2014; Lock et al., 2018). A
substantial fraction of the mantle may have been still molten at the beginning of the lunar evolutionary history after
the mantle overturn if the Moon was formed so hot. The more solid mantle after the mantle overturn we suggest here
would be more consistent with the “immediate origin” model of the Moon (Kegerreis et al., 2022) where the giant
impact is simulated at a resolution more than one order of magnitude higher than that employed in earlier studies;
this model suggests that the outer material of the Moon is heated to at least 4000 K by the impact, but the deeper
interior is only a few hundred Kelvin warmer than the pre-impact temperature, assumed to be around 2000 K. (Notethat our preferred initial potential temperature of around 1600 K (Figure 14c) is consistent with the temperature
after the mantle overturn in the earlier models of H. Li et al. (2019) and Yu et al. (2019) where the temperature of
the pre-overturn mantle is assumed to be around 2000 K close to the value suggested by Kegerreis et al. (2022)).
4.2.2. The Volcanic Activity
Our reference case is also consistent with the observed history of mare volcanism of the Moon. For the first 0.4 Gyrof the calculated history, the partially molten regions in the deep mantle grow only slowly, while the partiallymolten region in the uppermost mantle shrinks over time (Figure 2). This early stage is likely to correspond to
the period during which mare volcanism on the Moon was not so active (e.g., Hiesinger et al., 2003; Whitten &Head, 2015). The growth of melt-fingers and subsequent partially molten plumes observed in Figure 2b accountfor the lunar mare volcanism that became active after around 4 Gyr ago and peaked at 3.5–3.8 Gyr ago, andthen gradually declined over a period of billions of years (e.g., Hiesinger et al., 2000; Morota et al., 2011); thecalculated activity of partially molten plumes is indeed peaked at around 3.7 Gyr ago and then gradually declines
with time. These overall features of the calculated volcanic history are not affected by the detail of the rheology(Section 3.3) and possible density inversion between melts and matrix (Section 3.4). In contrast, in a model
where plumes do not appear (Case Ra2.15e4 calculated at lower Rayleigh number Ra = 2.15 × 10 4), the volcanic
activity develops only during the first several million years by melt-fingers (Figure 6b), suggesting that partiallymolten plumes play an important role in the volcanic history of the Moon.
Note that magma not enriched in HPEs is generated in the latter period of the calculated history (Figure 2 for
1.44 Gyr). This volcanism is caused by a return flow of a foundering material enriched in IBC components (seeU ET AL.

21 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

the arrow in Figure 2d). This may account for the volcanism of HPE-depleted young basalts (Che et al., 2021;
Q.-L. Li et al., 2021; Su et al., 2022).
4.2.3. The Temperature Profile
The depth-profile of the horizontally averaged temperature calculated at 4.4 Gyr in our reference case shows thatthe lithosphere develops as a thermal boundary layer of the convective mantle (Figure 3a), while the temperature
profile suggested for the present Moon (e.g., Karato, 2013; Khan et al., 2006, 2014; Sonett, 1982) is closer to athermal diffusion profile. The difference is not large in the shallow mantle (until the depth level of around 200 km),but the calculated temperature in the mid-mantle is considerably higher than that in the present Moon (Figure 3a).
As a consequence, the mantle is partially molten in the mid-mantle (from around r = 1,100–1,300 km), whereas
seismic evidence suggests that the partially molten region occurs only at the base of the mantle (e.g., Latham
et al., 1973; Nakamura et al., 1973; Tan & Harada, 2021; Weber et al., 2011). A temperature profile consistent
with the observed one was not obtained at other parameter values (Figure 3a). This difficulty may be a consequence of the assumed 2-D polar rectangular geometry of the convecting vessel and calls for further numericalcalculations in a 3-D spherical shell where the mid-mantle tends to be more strongly cooled and mantle convection is more sluggish (Guerrero et al., 2018).

5. Conclusions
To understand the volcanic and radial expansion/contraction history of the Moon, we developed a 2-D polar
rectangular numerical model of mantle evolution illustrated in Figure 15. The internally heated mantle of themodel evolves by the transport of heat, mass, and HPEs by mantle convection and migrating magma that is generated by decompression melting and internal heating.
Our simulations show that magma generation and migration play a crucial role in the volcanic and radial expansion/contraction history of the Moon. Magma is generated in the deep mantle by internal heating and eventuallyascends to the surface as partially molten plumes driven by melt buoyancy for the first several hundred million
years (Figures 2 and 15). This stage is likely to correspond to the period during which mare volcanism became
active after 4 Gyr ago with the peak at 3.5–3.8 Gyr ago (e.g., Whitten & Head, 2015). Subsequent magma
ascents by partially molten plumes decline with time but continue for billions of years after the peak becausesome materials that host HPEs are enriched in the IBC and remain in the deep mantle by their negative buoyancy (Figure 2). This activity accounts for the lunar mare volcanism that gradually declined after the peak (e.g.,
Hiesinger et al., 2003). The model which accounts for the observed mare volcanism is also consistent with theradial expansion/contraction history of the Moon, which globally expanded in its earlier history until around
3.8 Gyr ago and then contracted with time (e.g., Andrews-Hanna et al., 2013; Frueh et al., 2023). In our model,
the Moon expands by a few kilometers for the first several hundred million years and then contracts over time
(Figure 4). The lunar expansion is due to the extension of partially molten regions by partially molten plumes
that extract magma generated in the deep mantle (Figures 6a and 14); the subsequent contraction is caused bysolidification of the regions due to cooling from the surface boundary. The early expansion by mantle melting
suggested here implies that a substantial fraction of the mantle should have been solid, and there was a layer
enriched in HPEs and the IBC components at the base of the mantle in the Moon at the beginning of its history(Figures 13, 14c, and 14d and Table 3). In order to construct a more realistic thermal history model of the Moon,
it is necessary to extend our model to a 3-D spherical shell geometry and to introduce a lateral heterogeneity inthe initial thermal and compositional condition in future work.

Appendix A: The Basic Equations
In this section, we describe the basic equations that are not referred in Section 2.1. The continuity equation is∇ ⋅ 𝐔𝐔 = −∇ ⋅ [𝜙𝜙(𝐮𝐮 − 𝐔𝐔)].
(A1)The momentum equation for mantle convection is
[ (
)]
= 0.
−∇𝑃𝑃 + 𝜌𝜌𝜌𝜌 + ∇ ⋅ 𝜂𝜂 ∇𝐔𝐔+𝑡𝑡 ∇𝐔𝐔
(A2)Where the superscript t means transpose of a matrix.
U ET AL.

22 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Migration of magma is calculated in the energy equation (Katz, 2008),
𝜕𝜕(𝜌𝜌0 ℎ)
Δ𝑣𝑣l+ ∇ ⋅ (𝜌𝜌0 ℎ𝐔𝐔) = − ∇ ⋅ [𝜌𝜌0 ℎl 𝜙𝜙(𝐮𝐮 − 𝐔𝐔)] −𝜌𝜌0 𝑔𝑔𝑔𝑔𝑔𝑔𝑟𝑟 + ∇ ⋅ (𝑘𝑘∇𝑇𝑇 )
𝜕𝜕𝜕𝜕
𝑣𝑣0
(A3)
+ 𝜌𝜌0 𝑞𝑞 + ∇ ⋅ [𝜅𝜅edd ∇(𝜌𝜌0 ℎ)],where hl = h(ϕ = 1). k = ρ0Cpκ is the thermal conductivity, and κ is the thermal diffusivity. We assumed that thematrix disintegrates and that a strong turbulent diffusion occurs with the eddy diffusivity of κedd = 50κ in a largely
molten region with ϕ > 0.4 (Kameyama et al., 1996; Ogawa, 2020). κedd is assumed to gradually increase withincreasing ϕ as ϕ 3 (Ogawa, 2018b; U et al., 2022). We also assumed that the crustal thermal diffusivity is 0.48
times smaller than the mantle thermal diffusivity, taking into account the blanket effect of the crust and regolithlayers (Ziethe et al., 2009).
We calculate the phase diagram of the binary eutectic materials AξB1−ξ as well as the temperature, and melt-content
from the chemical potential defined for the materials. In the solid phase, the chemical potential μ solid is𝜇𝜇𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 = −𝐶𝐶𝑝𝑝 𝑇𝑇 ∕𝜎𝜎𝐴𝐴 − 𝑇𝑇 (𝑆𝑆0 + 𝐶𝐶𝑝𝑝 ln 𝑇𝑇 ∕𝜎𝜎𝐴𝐴 ) + 𝑃𝑃 ∕𝜌𝜌0
(A4)for both of the end-members A and B where S0 is an arbitrary constant. In the liquid state where the melt behaves
ideal solution in this model, the chemical potential μ liquid of the end-members i (i = A, B) is
𝜇𝜇𝑖𝑖𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 = 𝜇𝜇𝑖𝑖𝑙𝑙0 + 𝑅𝑅𝑅𝑅 ln 𝜉𝜉𝑖𝑖𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 ,
(A5)

and
(
)
𝜇𝜇𝑖𝑖𝑙𝑙0 = 𝜇𝜇𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 + Δℎ∕𝜎𝜎𝑖𝑖 1 + 𝐺𝐺 − 𝑇𝑇 ∕𝑇𝑇𝑖𝑖𝑙𝑙0 .
(A6)Here, R is the gas constant 8.3 J mol −1 K −1; σi the molar mass of each end-members (A = 140.69 g mol −1;
B = 151.71 g mol
𝐴𝐴 −1); 𝐴𝐴𝑖𝑖𝑙𝑙0 the dry solidus of the end member i at zero-pressure.
The core is regarded as a heat bath of a uniform temperature Tc that changes with time as
𝑑𝑑𝑑𝑑
𝑐𝑐𝑝𝑝c 𝜌𝜌core 𝑉𝑉core c = −𝑆𝑆𝑆𝑆 𝑆
(A7)
𝑑𝑑𝑑𝑑

where cpc = 675 J K −1 kg −1 (Ziethe et al., 2009) is the specific heat of the core, ρcore = 6,200 kg m −3 (Kronrod
et al., 2022) the core density, Vcore the volume of the core, and S the surface area of the CMB. The heat flux at the
CMB f is calculated from
𝜋𝜋 (
)
1
𝜕𝜕𝜕𝜕
𝑓𝑓 = −
𝑑𝑑𝑑𝑑𝑑
𝑘𝑘
(A8)
𝜋𝜋 ∫0
𝜕𝜕𝜕𝜕 𝑟𝑟=𝑟𝑟c

In this study, we neglect the internal heating in the core.
The heating rate q changes with time as( )𝑡𝑡
𝑞𝑞 = 𝑞𝑞tr∗ 𝑞𝑞0 exp − ,
(A9)
𝜏𝜏

where q0 = 14.7 pW kg −1 is the average initial heating rate at 4.4 Gyr ago estimated from the total amount of
HPEs in the current Moon (see Table 2 in U et al. (2022)); τ the decay time of HPEs. We approximate thisvalue as τ = 1.5 Gyr, an average of the decay times of 235U and 40K (Kameyama et al., 1996; Ogawa, 2020). The
non-dimensional𝐴𝐴value 𝐴𝐴tr∗ changes with migrating magma as
(
)
[
)
]
(
𝜕𝜕𝜕𝜕tr∗(A10)+ ∇ ⋅ 𝑞𝑞tr∗ 𝐔𝐔 = −∇ ⋅ 𝑞𝑞l∗ 𝜙𝜙(𝐮𝐮 − 𝐔𝐔) + ∇ ⋅ 𝜅𝜅edd ∇𝑞𝑞tr∗ .
𝜕𝜕𝜕𝜕
𝐴𝐴Here, 𝐴𝐴l∗ is the internal heating rate in the melt as
𝐷𝐷𝐷𝐷tr∗,𝑞𝑞l∗ =
(A11)(𝐷𝐷 − 1)𝜙𝜙 + 1

U ET AL.

23 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

and, D = 100 is the partition coefficient of HPEs between the solid phase and the melt phase.
The bulk composition ξb also changes with time owing to the mass transport by melt and matrix as𝜕𝜕𝜕𝜕b+ ∇ ⋅ (𝜉𝜉b 𝐔𝐔) = −∇ ⋅ [𝜉𝜉l 𝜙𝜙(𝐮𝐮 − 𝐔𝐔)] + ∇ ⋅ [𝜅𝜅edd ∇(𝜉𝜉b )].
(A12)
𝜕𝜕𝜕𝜕The basic equations are converted into their non-dimensional forms using the length scale L = rp − rc, the temperature scale Δh/Cp, and times scale L 2/κ. The momentum equation in its non-dimensional form is
[ (
)]
= 0.
−∇𝑃𝑃 ∗ + 𝑅𝑅𝑅𝑅𝑅𝑅∗ 𝐞𝐞𝐫𝐫 + ∇ ⋅ 𝜂𝜂 ∗ ∇𝐔𝐔∗ +𝑡𝑡 ∇𝐔𝐔∗
(A13)
where
[ ( ∗
)]
𝜂𝜂 ∗ = exp 𝐸𝐸𝑇𝑇∗ 𝑇𝑇𝑟𝑟𝑟𝑟𝑟𝑟− 𝑇𝑇 ∗ .
(A14)The non-dimensional relative velocity u* − U* is written as
{
}
2
Δ𝑣𝑣l
∗
∗
∗ ∗ 𝜙𝜙−𝛽𝛽(𝜉𝜉−𝜉𝜉𝐞𝐞𝐫𝐫 .
+
𝛽𝛽(1
−
𝜉𝜉
𝐮𝐮
−
𝐔𝐔
=
−𝑀𝑀
𝑔𝑔
[1
)]
)
l
s
l
(A15)
𝑣𝑣0
𝜙𝜙30
The non-dimensional energy equation is written as
[
]
Δ𝑣𝑣
𝜕𝜕𝜕∗
+ ∇ ⋅ (ℎ∗ 𝐔𝐔∗ ) = − ∇ ⋅ ℎ∗l 𝜙𝜙(𝐮𝐮∗ − 𝐔𝐔∗ ) − 𝑁𝑁 ∗ 𝑔𝑔 ∗ l 𝜙𝜙𝜙𝜙∗𝑟𝑟 + ∇ ⋅ (𝜅𝜅 ∗ ∇𝑇𝑇 ∗ )𝜕𝜕𝜕𝜕
𝑣𝑣0(A16)( ∗)+ 𝑞𝑞 ∗ + ∇ ⋅ 𝜅𝜅edd∇ℎ∗ ,where 𝐴∗ = 𝑇𝑇 ∗ + 𝜙𝜙(1 + 𝐺𝐺), and N* ≡ gsurL/Δh.
𝐴𝐴

Appendix B: Details of the Initial Condition
In this section, we describe the details of the initial condition. The initial distribution of the temperature is
obtained from that of the “reduced” enthalpy (see Equation 4), which isℎ = min(ℎsur , ℎmantle ),
(B1)

where
)]
[
(
𝑟𝑟
,
ℎsur = 𝐶𝐶𝑝𝑝 𝑇𝑇sur + 𝛿𝛿crst 1 −
(B2)
𝑟𝑟𝑝𝑝
⎧
if 𝑟𝑟 𝑟 𝑟𝑟lp⎪𝐶𝐶𝑝𝑝 𝑇𝑇
]
[
ℎ
=
.
(B3)
mantle()( 𝑟𝑟l −𝑟𝑟 )2⎨
if 𝑟𝑟 𝑟 𝑟𝑟l⎪𝐶𝐶𝑝𝑝 𝑇𝑇p + 𝑇𝑇c − 𝑇𝑇p 𝑟𝑟 −𝑟𝑟l c⎩Here, δcrst and rl are constants arbitrarily chosen. We assumed δcrst = 79.5 × 10 3 K to mimic the temperature-increase
with depth in the topmost 35 km (i.e., the crust) of the Moon, and rl = 550 km which is estimated from earlier
numerical models of the post-overturn stratification (e.g., Boukaré et al., 2018; Mitchell, 2021). Tp is the potential
temperature; Tc = 1875 K the initial temperature of the core. The value of Tc is based on the assumption that thetemperature of the core is higher than the mantle at 4.4 Gyr ago (Alley & Parmentier, 1998; Boukaré et al., 2018;
Maurice et al., 2020; Morbidelli et al., 2018). (Note that rl is not well-constrained from earlier modeling studies.
To see how rl affects our numerical results, we calculated Case rl-100 km where rl is 100 km and described thiscase in Text S3 of the Supporting Information S1. As discussed there, we confirmed that the reduction in rl does
not influence the overall features of mantle evolution observed in the reference case although found that the Moon
contracts by −2 km for the first 100 Myr. The early contraction is caused by a rapid cooling of the core (see Equation 12), which occurs because the reduction in rl enhances the heat flow at the CMB.)∫𝑟𝑟 p 𝑞𝑞𝑞𝑞𝑞𝑞𝑞𝑞 = ∫𝑟𝑟 p 𝑞𝑞0 𝑟𝑟𝑟𝑟𝑟𝑟, where
We assumed that the initial distribution of internal heat source q satisfies
𝐴𝐴
c
c
−1
q0 = 14.7 pW kg is the average value of internal heating rate at 4.4 Gyr ago (Laneuville et al., 2018; U
𝑟𝑟

U ET AL.

𝑟𝑟

24 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845∗
times higher than that in the mantle
et al., 2022). The total initial heating rate in the topmost 35 km 𝐴𝐴
qcrst is 𝐴𝐴crst(see Table 2):𝑟𝑟∫𝑟𝑟 p 𝑞𝑞crst 𝑟𝑟𝑟𝑟𝑟𝑟
crst∗
𝐹𝐹crst = 𝑟𝑟crst,
(B4)∫𝑟𝑟 𝑞𝑞m 𝑟𝑟𝑟𝑟𝑟𝑟cwhere,
𝑟𝑟crst𝑟𝑟p

𝑟𝑟p

𝑞𝑞crst 𝑟𝑟𝑟𝑟𝑟𝑟𝑟
𝑞𝑞m 𝑟𝑟𝑟𝑟𝑟𝑟 =𝑞𝑞0 𝑟𝑟𝑟𝑟𝑟𝑟 −(B5)∫𝑟𝑟crst∫𝑟𝑟c∫𝑟𝑟c

We assumed the initial distribution of internal heating rate in the mantle qm to be
𝑞𝑞m = 𝑞𝑞ol + 𝑞𝑞add ,(B6)where qol = 2.77 pW kg −1 is the internal heating rate of the olivine-rich materials with ξb = 1 (Yu et al., 2019),
while qadd is that of additional IBC components that increases with depth as[]𝑟𝑟 − 𝑟𝑟c
𝑞𝑞add = Δ𝑞𝑞 exp −
.
(B7)
𝑙𝑙
Here, l is the thickness of the basal layer which is enriched in the overturned materials (see Table 2); Δq is aconstant which is calculated from𝑟𝑟crst𝑟𝑟crst𝑟𝑟crst𝑞𝑞add 𝑟𝑟𝑟𝑟𝑟𝑟𝑟
𝑞𝑞ol 𝑟𝑟𝑟𝑟𝑟𝑟 +𝑞𝑞m 𝑟𝑟𝑟𝑟𝑟𝑟 =
(B8)
∫𝑟𝑟c
∫𝑟𝑟c
∫𝑟𝑟c

We calculated the initial distribution of the bulk composition ξb in the mantle from that of the internal heatingrate. We assumed that the IBC components with the composition ξe = 0.1 are 7.5 times more enriched in HPEs
than the bulk silicate Moon q0 (Hess & Parmentier, 1995) and that the initial content of the IBC components is
proportional to qadd as
1 − 𝜉𝜉b
𝑞𝑞add
=
.
(B9)
1 − 𝜉𝜉e
7.5𝑞𝑞0In the topmost 35 km, where we assumed ξb = 1 for simplicity.
Acknowledgments
The authors would like to extend theirsincere appreciation to T. Yanagisawa
and T. Miyagoshi at JAMSTEC for theirconstructive comments. The authors
also thank the editor L. Montési, and theanonymous reviewers for their insightfulcomments and suggestions, which helped
improve the quality of this paper. This
work was supported by JST SPRING,
Grant JP-MJSP2108 of Japan. This work
was also supported by: the Joint Usage/
Research Center PRIUS at Ehime University, the Earth Simulator of Japan Agency
for Marine-Earth Science and Technology(JAMSTEC), “Exploratory Challenge on
Post-K Computer” (Elucidation of theBirth of Exoplanets [Second Earth] andthe Environmental Variations of Planetsin the Solar System), and “Program for
Promoting Research on the Supercomputer Fugaku” (Toward a unified view of
the universe: from large scale structures
to planets). Animations and some figures
were drawn with the ParaView by Sandia
National Laboratory, Kitware Inc., and
Los Alamos National Laboratory.

U ET AL.

Data Availability Statement
The original data used to produce Table 3 and figures, the plots of Figures 1, 3, 4, and 14 as well as the numericalcode used to construct this work are found at U et al. (2023).

References
Alley, K. M., & Parmentier, E. M. (1998). Numerical experiments on thermal convection in a chemically stratified viscous fluid heated frombelow: Implications for a model of lunar evolution. Physics of the Earth and Planetary Interiors, 108(1), 15–32. https://doi.org/10.1016/
S0031-9201(98)00096-X
Andrews-Hanna, J. C., Asmar, S. W., Head, J. W., III, Kiefer, W. S., Konopliv, A. S., Lemoine, F. G., et al. (2013). Ancient igneous intrusions andearly expansion of the Moon revealed by GRAIL gravity gradiometry. Science, 339(6120), 675–678. https://doi.org/10.1126/science.1231753
Andrews-Hanna, J. C., Besserer, J., Head, J. W., III., Howett, C. J. A., Kiefer, W. S., Lucey, P. J., et al. (2014). Structure and evolution of the lunar
procellarum region as revealed by grail gravity data. Nature, 514(7520), 68–71. https://doi.org/10.1038/nature13697
Arai, T., Terada, H., Yamaguchi, A., & Ohtake, M. (2008). A new model of lunar crust: Asymmetry in crustal composition and evolution. Earth
Planets and Space, 60(4), 433–444. https://doi.org/10.1186/BF03352808
Boukaré, C.-E., Parmentier, E. M., & Parman, S. W. (2018). Timing of mantle overturn during magma ocean solidification. Earth and Planetary
Science Letters, 491(1), 216–225. https://doi.org/10.1016/j.epsl.2018.03.037
Breuer, D., & Moore, W. (2015). 10.08—Dynamics and thermal history of the terrestrial planets, the Moon, and Io. In G. Schubert (Ed.), Treatiseon geophysics (2nd ed., pp. 255–305). https://doi.org/10.1016/B978-0-444-53802-4.00173-1
Canup, R. M. (2004). Simulations of a late lunar-forming impact. Icarus, 168(2), 433–456. https://doi.org/10.1016/j.icarus.2003.09.028
Canup, R. M. (2014). Lunar-forming impacts: Processes and alternatives. Philosophical Transactions of the Royal Society A: Mathematical,
Physical and Engineering Sciences, 372(2024), 20130175. https://doi.org/10.1098/rsta.2013.0175
Cassen, R., & Reynolds, R. T. (1973). Role of convection in the Moon. Journal of Geophysical Research, 78(17), 3203–3215. https://doi.
org/10.1029/JB078i017p03203

25 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Cassen, R., Reynolds, R. T., Graziani, F., Summers, A., McNellis, J., & Blalock, L. (1979). Convection and lunar thermal history. Physics of the
Earth and Planetary Interiors, 19(2), 183–196. https://doi.org/10.1016/0031-9201(79)90082-7
Che, X., Nemchin, A., Liu, T., Long, D., Norman, M. D., Joy, K. H., et al. (2021). Age and composition of young basalts on the Moon, measured
from samples returned by Chang’e-5. Science, 374(6569), 887–890. https://doi.org/10.1126/science.abl7957
Cho, Y., Morota, T., Haruyama, J., Yasui, M., Hirata, N., & Sugita, S. (2012). Young mare volcanism in the Orientale region contemporary with
the Procellarum KREEP terrane (PKT) volcanism peak period around 2 billion years ago. Geophysical Research Letters, 39(11), L11203.
https://doi.org/10.1029/2012GL051838
Clark, J. M., Hurtado, J. D., Jr., Hiesinger, H., van der Bogert, C. H., & Bernhardt, H. (2017). Investigation of newly discovered lobate scarps:
Implications for the tectonic and thermal evolution of the Moon. Icarus, 298, 78–88. https://doi.org/10.1016/j.icarus.2017.08.017
Cuk, M., & Stewart, S. T. (2012). Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning. Science,
338(1047), 1047–1052. https://doi.org/10.1126/science.1225542
de Vries, J., van den Berg, A., & van Westrenen, W. (2010). Formation and evolution of a lunar core from ilmenite-rich magma ocean cumulates.
Earth and Planetary Science Letters, 292(1–2), 139–147. https://doi.org/10.1016/j.epsl.2010.01.029
Dumoulin, C., Doin, M.-P., & Fleitout, L. (1999). Heat transport in stagnant lid convection with temperature- and pressure-dependent Newtonian
or non-Newtonian rheology. Journal of Geophysical Research, 104(B6), 12759–12777. https://doi.org/10.1029/1999JB900110
Dygert, N., Hirth, G., & Liang, Y. (2016). A flow law for ilmenite in dislocation creep: Implications for lunar cumulate mantle overturn. Geophysical Research Letters, 43(2), 532–540. https://doi.org/10.1002/2015GL066546
Elkins-Tanton, L., & Bercovici, D. (2014). Contraction or expansion of the Moon’s crust during magma ocean freezing? Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372(2024), 20130240. https://doi.org/10.1098/rsta.2013.0240
Elkins-Tanton, L., Burgess, S., & Yin, Q.-Z. (2011). The lunar magma ocean: Reconciling the solidification process with lunar petrology andgeochronology. Earth and Planetary Science Letters, 304(3–4), 326–336. https://doi.org/10.1016/j.epsl.2011.02.004
Frueh, T., Hiesinger, H., Bogert, C. H., Clark, J. D., Watters, T. R., & Schmedemann, N. (2023). Timing and origin of compressional tectonism inmare tranquillitatis. Journal of Geophysical Research: Planets, 128, e2022JE007533. https://doi.org/10.1029/2022JE0075
Garcia, R. F., Gagnepain-Beyneix, J., Chevrot, S., & Lognonné, P. (2011). Very preliminary reference Moon model. Journal of GeophysicalResearch, 188(1–2), 96–113. https://doi.org/10.1016/j.pepi.2011.06.015
Garcia, R. F., Gagnepain-Beyneix, J., Chevrot, S., & Lognonné, P. (2012). Erratum to “very preliminary reference Moon model”. Physics of the
Earth and Planetary Interiors, 202–203, 89–91. https://doi.org/10.1016/j.pepi.2012.03.009
Guerrero, J. M., Lowman, J. P., Deschamps, F., & Tackley, P. J. (2018). The influence of curvature on convection in a temperature-dependent
viscosity fluid: Implications for the 2-D and 3-D modeling of Moon’s. Journal of Geophysical Research: Planets, 123(7), 1863–1880. https://doi.org/10.1029/2017JE005497
Head, J. W., & Wilson, L. (1992). Lunar mare volcanism: Stratigraphy, eruption conditions, and the evolution of secondary crusts. Geochimica et
Cosmochimica Acta, 56(6), 2155–2175. https://doi.org/10.1016/0016-7037(92)90183-J
Head, J. W., & Wilson, L. (2017). Generation, ascent and eruption of magma on the Moon: New insights into source depths, magma supply,
intrusions and effusive/explosive eruptions (part 2: Predicted emplacement processes and observations). Icarus, 283, 176–223. https://doi.
org/10.1016/j.icarus.2016.05.031
Head, J. W., & Wilson, L. (2020). Magmatic intrusion-related processes in the upper lunar crust: The role of country rock porosity/permeability
in magmatic percolation and thermal annealing, and implications for gravity signatures. Planetary and Space Science, 180, 104765. https://doi.org/10.1016/j.pss.2019.104765
Hess, P. C., & Parmentier, E. M. (1995). A model for the thermal and chemical evolution of the Moon’s interior: Implications for the onset of marevolcanism. Earth and Planetary Science Letters, 134(3–4), 501–514. https://doi.org/10.1016/0012-821X(95)00138-3
Hess, P. C., & Parmentier, E. M. (2001). Thermal evolution of a thicker KREEP liquid layer. Journal of Geophysical Research, 106(E11),28023–28032. https://doi.org/10.1029/2000JE001416
Hiesinger, H., Head, J. W., III., Wolf, U., Jaumann, R., & Neukum, G. (2003). Ages and stratigraphy of mare basalts in Oceanus Procellarum,
Mare Nubium, Mare Cognitum, and Mare Insularum. Journal of Geophysical Research, 108(E7), 5065. https://doi.org/10.1029/2002JE001985
Hiesinger, H., Jaumann, R., Neukum, G., & Head, J. W. (2000). Ages and stratigraphy of mare basalts in Oceanus Procellarum, Mare Nubium,
Mare Cognitum, and Mare Insularum. Journal of Geophysical Research, 105(E12), 29239–29275. https://doi.org/10.1029/2000JE001244
Hirth, G., & Kohlstedt, D. L. (2003). Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. Geophysical Monograph Series, 138, 83–105. https://doi.org/10.1029/138GM06
Hosono, N., Karato, S.-I., Makino, J., & Saitoh, T. R. (2019). Terrestrial magma ocean origin of the Moon. Nature Geoscience, 12(6), 418–423.
https://doi.org/10.1038/s41561-019-0354-2
Jolliff, B. L., Gillis, J. J., Haskin, L., Korotev, R. L., & Wieczorek, M. A. (2000). Major lunar crustal terranes: Surface expressions and crustalmantle origins. Journal of Geophysical Research, 105(E2), 4197–4216. https://doi.org/10.1029/1999JE001103
Jones, M. J., Evans, A. J., Johnson, B. C., Weller, M. B., Andrews-Hanna, J. C., Tikoo, S. M., & Keane, J. T. (2022). A South Pole–Aitken impact
origin of the lunar compositional asymmetry. Science Advances, 8(14). https://doi.org/10.1126/sciadv.abm8475
Kameyama, M., Fujimoto, H., & Ogawa, M. (1996). A thermo-chemical regime in the upper mantle in the early Earth inferred from a numerical model of magma-migration in a convecting upper mantle. Physics of the Earth and Planetary Interiors, 94(3–4), 187–215. https://doi.
org/10.1016/0031-9201(95)03102-2
Karato, S. (2013). Geophysical constraints on the water content of the lunar mantle and its implications for the origin of the Moon. Earth andPlanetary Science Letters, 384, 144–153. https://doi.org/10.1016/j.epsl.2013.10.001
Karato, S., & Wu, P. (1993). Rheology of the upper mantle: A synthesis. Science, 260(5109), 771–778. https://doi.org/10.1126/science.260.5109.771
Katz, R. F. (2008). Magma dynamics with the enthalpy method: Benchmark solutions and magmatic focusing at mid-ocean ridges. Geochemistry,
Geophysics, Geosystems, 49(12), 2099–2121. https://doi.org/10.1093/petrology/egn058
Katz, R. F., Spiegelman, M., & Langmuir, C. H. (2003). A new parameterization of hydrous mantle melting. Geochemistry, Geophysics, Geosystems, 4(9), 1073. https://doi.org/10.1029/2002GC000433
Kegerreis, J. A., Ruiz-Bonilla, S., Eke, V. R., Massey, R. J., Sandnes, T. D., & Teodoro, L. F. A. (2022). Immediate origin of the Moon as a
post-impact satellite. The Astrophysical Journal Letters, 937(2), L40. https://doi.org/10.3847/2041-8213/ac8d96
Khan, A., Connolly, J. A. D., Pommier, A., & Noir, J. (2014). Geophysical evidence for melt in the deep lunar interior and implications for lunar
evolution. Journal of Geophysical Research: Planets, 119(10), 2197–2221. https://doi.org/10.1002/2014JE004661
Khan, A., Maclennan, J., Taylor, S. R., & Connolly, J. A. D. (2006). Are the Earth and the Moon compositionally alike? Inferences on lunar
composition and implications for lunar origin and evolution from geophysical modeling. Journal of Geophysical Research, 111(E5), E05005.
https://doi.org/10.1029/2005JE002608

U ET AL.

26 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Kirk, R. L., & Stevenson, D. J. (1989). The competition between thermal contraction and differentiation in the stress history of the Moon. Journalof Geophysical Research, 94(B9), 12133–12144. https://doi.org/10.1029/JB094iB09p12133
Klimczak, C. (2015). Limits on the brittle strength of planetary lithospheres undergoing global contraction. Journal of Geophysical Research:
Planets, 120(12), 2135–2151. https://doi.org/10.1002/2015JE004851
Konrad, W., & Spohn, T. (1997). Thermal history of the Moon: Implications for an early core dynamo and post-accretional magmatism. Advancesin Space Research, 19(10), 1511–1521. https://doi.org/10.1016/S0273-1177(97)00364-5
Kronrod, E., Matsumoto, K., Kuskov, O. L., Kronrod, V., Yamada, R., & Kamata, S. (2022). Towards geochemical alternatives to geophysicalmodels of the internal structure of the lunar mantle and core. Advances in Space Research, 69(7), 2798–2824. https://doi.org/10.1016/j.
asr.2022.01.012
Laneuville, M., Taylor, J., & Wieczorek, M. A. (2018). Distribution of radioactive heat sources and thermal history of the Moon. Journal ofGeophysical Research: Planets, 123(12), 3144–3166. https://doi.org/10.1029/2018JE005742
Laneuville, M., Wieczorek, M. A., Breuer, D., Aubert, J., & Ruckriemen, T. (2014). A long-lived lunar dynamo powered by core crystallization.
Earth and Planetary Science Letters, 401, 251–260. https://doi.org/10.1016/j.epsl.2014.05.057
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. https://doi.org/10.1002/jgre.20103
Latham, G., Dorman, J., Duennebier, F., Ewing, M., Lammlein, D., & Nakamura, Y. (1973). Moonquakes, meteoroids, and the state of the lunar
interior. In Proceedings of the lunar science conference (Vol. 4, pp. 2515–2527).
Lawrence, D. J., Feldman, W. C., Barraclough, B. L., Binder, A. B., Elphic, R. C., Maurice, S., et al. (2000). Small-area thorium features on the
lunar surface. Journal of Geophysical Research, 105(E8), 20307–20331. https://doi.org/10.1029/1999JE001177
Le Bars, M., & Davaille, A. (2004). Whole layer convection in a heterogeneous planetary mantle. Journal of Geophysical Research, 109(B3),
B03403. https://doi.org/10.1029/2003JB002617
Li, H., Zhang, N., Ling, Y., Wu, B., Dygert, N. K., Huang, J., & Parmentier, E. M. (2019). Lunar cumulate mantle overturn: A model constrainedby ilmenite rheology. Journal of Geophysical Research: Planets, 124(5), 1357–1378. https://doi.org/10.1029/2018JE005905
Li, Q.-L., Zhou, Q., Liu, Y., Xiao, Z., Lin, Y., Li, J.-H., et al. (2021). Two-billion-year-old volcanism on the Moon from Chang’e-5 basalts. Nature,600(7887), 54–58. https://doi.org/10.1038/s41586-021-04100-2
Liang, W., & Andrews-Hanna, J. C. (2022). Probing the source of ancient linear gravity anomalies on the Moon. Icarus, 380, 114978. https://doi.
org/10.1016/j.icarus.2022.114978
Lock, S., Stewart, S., Petaev, M., Leinhardt, z., Mace, M., Jacobsen, S., & Cuk, M. (2018). The origin of the Moon within a terrestrial synestia.
Journal of Geophysical Research: Planets, 123(4), 910–951. https://doi.org/10.1002/2017JE005333
Lourenço, D. L., Rozel, A. B., Gerya, T., & Tackley, P. J. (2018). Efficient cooling of rocky planets by intrusive magmatism. Journal of Geophysical Research, 11(5), 322–327. https://doi.org/10.1038/s41561-018-0094-8
Matsuyama, L., Keane, J. T., Trinh, A., Beuthe, M., & Watters, T. R. (2021). Global tectonic patterns of the Moon. Geophysical Research Letters,
358(114202), 114202. https://doi.org/10.1016/j.icarus.2020.114202
Maurice, M., Tosi, N., Schwinger, S., Breuer, D., & Kleine, T. (2020). A long-lived magma ocean on a young Moon. Science Advances, 6(28).
https://doi.org/10.1126/sciadv.aba8949
McKenzie, D. (1984). The generation and compaction of partially molten rock. Journal of Petrology, 25(3), 713–765. https://doi.org/10.1093/
petrology/25.3.713
Mei, S., Bai, W., Hiraga, T., & Kohlstedt, D. L. (2002). Influence of melt on the creep behavior of olivine–basalt aggregates under hydrous conditions. Earth and Planetary Science Letters, 201(3–4), 491–507. https://doi.org/10.1016/S0012-821X(02)00745-8
Miller, K. J., Zhu, W.-L., Montési, L. G. J., & Gaetani, G. A. (2014). Experimental quantification of permeability of partially molten mantle rock.
Earth and Planetary Science Letters, 388(15), 273–282. https://doi.org/10.1016/j.epsl.2013.12.003
Mitchell, R. N. (2021). Chang’e-5 reveals the Moon’s secrets to a longer life. The Innovation, 2(4), 100177. https://doi.org/10.1016/j.xinn.2021.100177
Morbidelli, A., Nesvorny, D., Laurenz, V., Marchi, S., Rubie, D. C., Elkins-Tanton, L., et al. (2018). The timeline of the lunar bombardment:
Revisited. Icarus, 305, 262–276. https://doi.org/10.1016/j.icarus.2017.12.046
Moriarty, D. P., III., Dygert, N., Valencia, S. N., Watkins, R. N., & Petro, N. E. (2021). The search for lunar mantle rocks exposed on the surfaceof the Moon. Nature Communications, 12(4659), 4659. https://doi.org/10.1038/s41467-021-24626-3
Morota, T., Haruyama, J., Honda, C., Ohtake, M., Yokota, Y., Kimura, J., et al. (2009). Mare volcanism in the lunar farside Moscoviense
region: Implication for lateral variation in magma production of the Moon. Geophysical Research Letters, 36(21), L21202. https://doi.
org/10.1029/2009GL040472
Morota, T., Haruyama, J., Ohtake, M., Matsunaga, T., Kawamura, T., Yokota, Y., et al. (2011). Timing and duration of mare volcanism in thecentral region of the northern farside of the Moon. Earth Planets and Space, 63(1), 5–11. https://doi.org/10.5047/eps.2010.02.009
Nakamura, Y., Lammlein, D., Latham, G., Ewing, M., Dorman, J., Press, F., & Toksöz, N. (1973). New seismic data on the state of the deep lunar
interior. Science, 181(4094), 49–51. https://doi.org/10.1126/science.181.4094.49
Newsom, H. E., & Taylor, S. R. (1989). Geochemical implications of the formation of the Moon by a single giant impact. Nature, 338(6210),
29–34. https://doi.org/10.1038/338029a0
Ogawa, M. (2014). A positive feedback between magmatism and mantle upwelling in terrestrial planets: Implications for the Moon. Journal ofGeophysical Research: Planets, 119(11), 2317–2330. https://doi.org/10.1002/2014JE004717
Ogawa, M. (2018a). The effects of magmatic redistribution of heat producing elements on the lunar mantle evolution inferred from numericalmodels that start from various initial states. Planetary and Space Science, 151, 43–55. https://doi.org/10.1016/j.pss.2017.10.015
Ogawa, M. (2018b). Magmatic differentiation and convective stirring of the mantle in early planets: The effects of the magmatism-mantle
upwelling feedback. Geophysical Journal International, 215(3), 2144–2155. https://doi.org/10.1093/gji/ggy413
Ogawa, M. (2020). Magmatic differentiation and convective stirring of the mantle in early planets—2: Effects of the properties of mantle materials. Geophysical Journal International, 220(2), 1409–1420. https://doi.org/10.1093/gji/ggz499
Pritchard, M. E., & Stevenson, D. J. (2000). Thermal aspects of a lunar origin by giant impact. In R. M. Canup, & K. Righter (Eds.), Origin of the
Earth and Moon (pp. 179–196). University of Arizona Press. https://doi.org/10.2307/j.ctv1v7zdrp
Rapp, J. F., & Draper, D. S. (2018). Fractional crystallization of the lunar magma ocean: Updating the dominant paradigm. Meteoritics & Planetary Science, 53(7), 1432–1455. https://doi.org/10.1111/maps.13086
Ringwood, A. E., & Kesson, S. E. (1976). A dynamic model for mare basalt petrogenesis. In Proceedings of the lunar science conference (Vol.
7, pp. 1697–1722).
Rufu, R., Aharonson, O., & Perets, H. B. (2017). A multiple-impact origin for the Moon. Nature Geoscience, 10(2), 89–94. https://doi.
org/10.1038/NGEO2866

U ET AL.

27 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Sakamaki, T., Ohtani, E., Urakawa, S., Suzuki, S., Katayama, Y., & Zhao, D. (2010). Density of high-Ti basalt magma at high pressure and originof heterogeneities in the lunar mantle. Earth and Planetary Science Letters, 299(3–4), 285–289. https://doi.org/10.1016/j.epsl.2010.09.007
Sawada, N., Morota, T., Kato, S., Ishihara, Y., & Hiramatsu, Y. (2016). Constraints on timing and magnitude of early global expansion of
the Moon from topographic features in linear gravity anomaly areas. Geophysical Research Letters, 43(10), 4865–4870. https://doi.
org/10.1002/2016GL068966
Schwinger, S., & Breuer, D. (2022). Employing magma ocean crystallization models to constrain structure and composition of the lunar interior.
Physics of the Earth and Planetary Interiors, 322, 106831. https://doi.org/10.1016/j.pepi.2021.106831
Scott, T., & Kohlstedt, D. L. (2006). The effect of large melt fraction on the deformation behavior of peridotite. Geophysical Monograph Series,
246(3–4), 177–187. https://doi.org/10.1016/j.epsl.2006.04.027
Shearer, C. K., Hess, P. C., Wieczorek, M. A., Pritchard, M. E., Parmentier, E. M., Borg, L. E., & Wiechert, U. (2006). Thermal and magmatic
evolution of the Moon. Reviews in Mineralogy and Geochemistry, 60(1), 365–518. https://doi.org/10.2138/rmg.2006.60.4
Shearer, C. K., & Papike, J. (1999). Magmatic evolution of the Moon. American Mineralogist, 84(1), 1469–1494. https://doi.org/10.2138/
am-1999-1001
Siegler, M. A., & Smrekar, S. E. (2014). Lunar heat flow: Regional prospective of the apollo landing sites. Journal of Geophysical Research:
Planets, 119(1), 47–63. https://doi.org/10.1002/2013JE004453
Siegler, M. A., Warren, P., Franco, K. L., Paige, D., Feng, J., & White, M. (2022). Lunar heat flow: Global predictions and reduced heat flux.
Journal of Geophysical Research: Planets, 127(9), e2022JE007182. https://doi.org/10.1029/2022JE007182
Snyder, G. A., Taylor, L. A., & Neal, C. R. (1992). A chemical model for generating the sources of mare basalts: Combined equilibrium and fractional crystallization of the lunar magmasphere. Geochimica et Cosmochimica Acta, 56(10), 3809–3823. https://doi.
org/10.1016/0016-7037(92)90172-F
Solomon, S. C. (1975). Mare volcanism and lunar crustal structure. In Proceedings of the lunar science conference (Vol. 6, pp. 1021–1042).
Solomon, S. C. (1986). On the early thermal state of the Moon (pp. 435–452). Lunar and Planetary Institute.
Solomon, S. C., & Chaiken, J. (1976). Thermal expansion and thermal stress in the Moon and terrestrial planets: Clues to early thermal history.
In Proceedings of the lunar science conference (Vol. 7, pp. 3229–3243).
Solomon, S. C., & Head, J. W. (1979). Vertical movement in mare basins: Relation to mare emplacement, basin tectonics, and lunar thermal
history. Journal of Geophysical Research, 84(B4), 1667–1682. https://doi.org/10.1029/JB084iB04p01667
Solomon, S. C., & Toksöz, M. N. (1973). Internal constitution and evolution of the Moon. Physics of the Earth and Planetary Interiors, 7(1),15–38. https://doi.org/10.1016/0031-9201(73)90037-X
Sonett, C. P. (1982). Electromagnetic induction in the Moon. Reviews of Geophysics, 20(3), 411. https://doi.org/10.1029/RG020i003p00411
Spohn, T., Konrad, W., Breuer, D., & Ziethe, R. (2001). The longevity of lunar volcanism: Implications of thermal evolution calculations with 2D
and 3D mantle convection models. Icarus, 149(1), 54–65. https://doi.org/10.1006/icar.2000.6514
Stegman, D. R., Jellinek, A. M., Zatman, S. A., Baumgardner, J. R., & Richards, M. A. (2003). An early lunar core dynamo driven by thermochemical mantle convection. Nature, 421(6919), 143–146. https://doi.org/10.1038/nature01267
Stevenson, D. J. (1987). Origin of the Moon-the collision hypothesis. Annual Review of Earth and Planetary Sciences, 15(1), 271–315. https://doi.org/10.1146/annurev.ea.15.050187.001415
Su, B., Yuan, J., Chen, Y., Yang, W., Mitchell, R. N., Hui, H., et al. (2022). Fusible mantle cumulates trigger young mare volcanism on the cooling
Moon. Science Advances, 8(42). https://doi.org/10.1126/sciadv.abn2103
Taguchi, M., Morota, T., & Kato, S. (2017). Lateral heterogeneity of lunar volcanic activity according to volumes of mare basalts in the farsidebasins. Journal of Geophysical Research: Planets, 122(7), 1505–1521. https://doi.org/10.1002/2016JE005246
Tan, Y., & Harada, Y. (2021). Tidal constraints on the low-viscosity zone of the Moon. Icarus, 365, 114361. https://doi.org/10.1016/j.
icarus.2021.114361
Toksőz, M. N., & Solomon, S. C. (1973). Thermal history and evolution of the Moon. The Moon, 7(3–4), 251–278. https://doi.org/10.1007/bf00564634
U, K., Hasumi, H., & Ogawa, M. (2022). Effects of magma-generation and migration on the expansion and contraction history of the Moon. Earth
Planets and Space, 74(78), 78. https://doi.org/10.1186/s40623-022-01631-4
U, K., Kameyama, M., & Ogawa, M. (2023). Numerical models of magmatism in the convective mantle of the Moon [Dataset]. Figshare. https://doi.org/10.6084/m9.figshare.22297921.v1
van der Bogert, C. H., Clark, J. D., Hiesinger, H., Banks, M. E., Watters, T. R., & Robinson, M. S. (2018). How old are lunar lobate scarps? 1.
Seismic resetting of crater size-frequency distributions. Icarus, 306, 225–242. https://doi.org/10.1016/j.icarus.2018.01.019
Vander Kaaden, K. E., Agee, C. B., & McCubbin, F. M. (2015). Density and compressibility of the molten lunar picritic glasses: Implications forthe roles of ti and Fe in the structures of silicate melts. Geochimica et Cosmochimica Acta, 149, 1–20. https://doi.org/10.1016/j.gca.2014.10.029
van Kan Parker, M., Sanloup, C., Sator, N., Guillot, B., Tronche, E. J., Perrillat, P., et al. (2012). Neutral buoyancy of titanium-rich melts in the
deep lunar interior. Nature Geoscience, 5(3), 186–189. https://doi.org/10.1038/ngeo1402
Viswanathan, V., Rambaux, N., Fienga, A., Laskar, J., & Gastineau, M. (2019). Observational constraint on the radius and oblateness of the lunar
core-mantle boundary. Geophysical Research Letters, 46(13), 7295–7303. https://doi.org/10.1029/2019GL082677
Watters, T. R., Robinson, M. S., Beyer, R. A., Banks, M. E., Bell, J. F., Pritchard, M. E., et al. (2010). Evidence of recent thrust faulting on theMoon revealed by the lunar reconnaissance orbiter camera. Science, 329(5994), 936–940. https://doi.org/10.1126/science.1189590
Watters, T. R., Robinson, M. S., Collins, G. C., Banks, M. E., Daud, K., Williams, N. R., & Selvans, M. M. (2015). Global thrust faulting on theMoon and the influence of tidal stresses. Geology, 43(10), 851–854. https://doi.org/10.1130/G37120.1
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.
https://doi.org/10.1126/science.1199375
Whitehead, J. A. J., & Luther, D. S. (1975). The generation and compaction of partially molten rock. Journal of Geophysical Research, 80(5),705–717. https://doi.org/10.1029/JB080i005p00705
Whitten, J. L., & Head, J. W. (2015). Lunar cryptomaria: Physical characteristics, distribution, and implications for ancient volcanism. Icarus,
247, 150–171. https://doi.org/10.1016/j.icarus.2014.09.031
Wieczorek, M. A., Neumann, G. A., Nimmo, F., Kiefer, W. S., Taylor, G. J., Melosh, R. J., et al. (2013). The crust of the Moon as seen by GRAIL.
Science, 339(6029), 671–675. https://doi.org/10.1126/science.1199375
Wieczorek, M. A., & Phillips, R. J. (2000). The “procellarum KREEP terrane”: Implications for mare volcanism and lunar evolution. Journal ofGeophysical Research, 105(E8), 20417–20430. https://doi.org/10.1029/1999JE001092
Wilson, L., & Head, J. W. (2003). Deep generation of magmatic gas on the Moon and implications for pyroclastic eruptions. Icarus, 30(12).
https://doi.org/10.1029/2002GL016082
Wilson, L., & Head, J. W. (2017). Generation, ascent and eruption of magma on the Moon: New insights into source depths, magma supply,
intrusions and effusive/explosive eruptions (part 1: Theory). Icarus, 283, 146–175. https://doi.org/10.1016/j.icarus.2015.12.039

U ET AL.

28 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets

10.1029/2023JE007845

Wood, J. A. (1972). Thermal history and early magmatism in the Moon. Icarus, 16(2), 229–240. https://doi.org/10.1016/0019-1035(72)90070-X
Xu, M., Jing, Z., Van Orman, J. A., Yu, T., & Wang, Y. (2022). Experimental evidence supporting an overturned iron-titanium-rich melt layer inthe deep lunar interior. Geophysical Research Letters, 49(13), e2022GL099066. https://doi.org/10.1029/2022GL099066
Yan, J., Xu, L., Li, F., Matsumoto, K., Rodriguez, J. A. P., Miyamoto, H., & Dohm, J. M. (2015). Lunar core structure investigation: Implicationof GRAIL gravity field model. Advances in Space Research, 55(6), 1721–1727. https://doi.org/10.1016/j.asr.2014.12.038
Yanagisawa, T., Kameyama, M., & Ogawa, M. (2016). Numerical studies on convective stability and flow pattern in three-dimensional sphericalmantle of terrestrial planets. Geophysical Journal International, 206(3), 1526–1538. https://doi.org/10.1093/gji/ggw226
Yu, S., Tosi, N., Schwinger, S., Maurice, M., Breuer, D., & Xiao, L. (2019). Overturn of ilmenite-bearing cumulates in a rheologically weak lunar
mantle. Journal of Geophysical Research: Planets, 124(2), 418–436. https://doi.org/10.1029/2018je005739
Yue, Z., Michael, G. G., Di, K., & Liu, J. (2017). Global survey of lunar wrinkle ridge formation times. Earth and Planetary Science Letters, 477,
14–20. https://doi.org/10.1016/j.epsl.2017.07.048
Zhang, N., Ding, M., Zhu, M.-H., Li, H., Li, H., & Yue, Z. (2022a). Lunar compositional asymmetry explained by mantle overturn following theSouth Pole–Aitken impact. Nature Geoscience, 15(1), 37–41. https://doi.org/10.1038/s41561-021-00872-4
Zhang, N., Dygert, N., Liang, Y., & Parmentier, E. (2017). The effect of ilmenite viscosity on the dynamics and evolution of an overturned lunar
cumulate mantle. Geophysical Research Letters, 44(13), 6543–6552. https://doi.org/10.1002/2017GL073702
Zhang, N., Parmentier, E., & Liang, Y. (2013a). A 3-D numerical study of the thermal evolution of the Moon after cumulate mantle overturn:
The importance of rheology and core solidification. Journal of Geophysical Research: Planets, 118(9), 1789–1804. https://doi.org/10.1002/
jgre.20121
Zhang, N., Parmentier, E., & Liang, Y. (2013b). Effects of lunar cumulate mantle overturn and megaregolith on the expansion and contraction
history of the Moon. Geophysical Research Letters, 40(19), 5019–5023. https://doi.org/10.1002/grl.50988
Zhang, W. B., Zhang, N., & Li, H. Y. (2022b). Abundances of lunar heat-producing elements constrained by a 3-D numerical model oftitanium-rich basaltic eruption. Chinese Journal of Geophysics, 65(1), 119–136. https://doi.org/10.6038/cjg2022P0753
Zhao, Y., de Vries, J., van den Berg, A. P., Jacobs, M. H. G., & van Westrenen, W. (2019). The participation of ilmenite-bearing cumulates in lunar
mantle overturn. Earth and Planetary Science Letters, 511(1), 1–11. https://doi.org/10.1016/j.epsl.2019.01.022
Zhong, S., Parmentier, E. M., & Zuber, M. T. (2000). A dynamic origin for the global asymmetry of lunar mare basalts. Earth and Planetary
Science Letters, 177(3–4), 131–140. https://doi.org/10.1016/S0012-821X(00)00041-8
Ziethe, R., Seiferlin, K., & Hiesinger, H. (2009). Duration and extent of lunar volcanism: Comparison of 3D convection models to mare basalt
ages. Planetary and Space Science, 57(7), 784–796. https://doi.org/10.1016/j.pss.2009.02.002

U ET AL.

29 of 29

21699100, 2023, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JE007845 by Cochrane Japan, Wiley Online Library on [20/09/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Journal of Geophysical Research: Planets