Effects of magma-generation and migration on the expansion and contraction history of the Moon

Geological and geodetic observations of the Moon from spacecraft revealed that it expanded by a few km for the first several hundred million years and then contracted later. The period when the planet expanded most coincides with that when the mare volcanism of the Moon was active. Given the high initial temperature of the deep mantle inferred from the giant impact and mantle overturn hypotheses of the Moon, the observed early expansion is difficult to account for by thermal expansion only. To understand the observed radial change of the Moon, we numerically calculated the thermal evolution of a one-dimensional spherically symmetric mantle caused by transport of heat, mass, and incompatible heat-producing elements (HPEs) by migration of magma that is generated by internal heating. The mantle is assumed to be enriched in HPEs at its base in the initial condition. The calculated mantle expands for the first several hundred million years by melting of the deep mantle and upward migration of the generated magma to the uppermost mantle; the top of the partially molten region rises to the depth level of around 300 km, which is shallow enough to generate mare basalts of the Moon. The migrating magma, however, extracts HPEs from the deep interior, and the planet then contracts gradually by cooling and solidification of the partially molten mantle. We obtained a thermal history model that is consistent with the observed history of radial change of the Moon when the initial mid-mantle temperature TM≈1600K\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{{\text{M}}} \approx 1600 \,{\text{K}}$$\end{document} and the initial ratio of the concentration of HPEs in the crust to that of the mantle Fcrst∗≤12\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{{{\text{crst}}}}^{*} \le 12$$\end{document}. This model suggests that melting of the deep mantle and upward migration of the generated magma strongly affect the thermal history of the Moon. The model we developed here is a good starting point for constructing more realistic models of the thermal history of the Moon where the effects of heat and mass transport by mantle convection are also considered.


Introduction
The thermal history of the Moon has been a long-standing issue in studies of terrestrial planets (e.g., Toksöz et al. 1972;Solomon and Chaiken 1976;Cassen et al. 1979;Kirk and Stevenson 1989;Shearer et al. 2006;de Vries 2012;Evans et al. 2014). Geological and geodetic observations from spacecraft have revealed that the Moon expanded globally in the early stage of its history Sawada et al. 2016) and then contracted globally in the later stage (Watters et al. 2010;Matsuyama et al. 2021). The period when its radius reached the maximum coincides with that when the mare volcanism was active: the volcanism was weak before 4 Gyr ago, became active between 3.3 and 3.8 Gyr ago with the peak at ~ 3.6 Gyr ago, and then mildly continued until 1-2 Gyr ago (e.g., Hiesinger et al. 2003;Whitten and Head 2015). To understand these features of the history of the Moon, we developed a spherically symmetric model of the lunar mantle evolution where magma-generation and migration as well as transport of heat, mass, and incompatible heat-producing elements (HPEs) by migrating magma are considered.
In most of earlier one-dimensional numerical models of the thermal history, the planetary expansion/contraction of the Moon is assumed to occur thermally, and the early expansion is reproduced when the initial temperature in the deep mantle is 1100 K or less (e.g., Solomon and Chaiken 1976): the Moon monotonously contracts throughout its history when the mantle is initially hotter because of secular cooling of the mantle. In a model where the effect of volume change from compositional differentiation is also considered, the early expansion is reproduced even when the initial temperature in the deep mantle is as high as 1200 K (Kirk and Stevenson 1989). However, the initial temperature of 1200 K or less in the deep mantle is still considerably lower than that predicted from the giant impact and mantle overturn hypotheses of the Moon (e.g., Stevenson 1987;Alley and Parmentier 1998;Pritchard and Stevenson 2000;Canup 2004;Rufu et al. 2017;Boukaré et al. 2018;Lock et al. 2018).
The early expansion of the Moon is hard to account for by three-dimensional models of mantle evolution that have been developed in the literature to understand the mare volcanism, too. Some models succeeded in reproducing the peak of mare volcanism at around 3.6 Gyr ago (e.g., Zhong et al. 2000;Stegman et al. 2003;Zhang et al. 2013a, b). Some other models succeeded in reproducing volcanism that continues for billions of years as is observed for the lunar mare volcanism (Konrad and Spohn 1997;Spohn et al. 2001). By considering the blanket effect of the crust, Ziethe et al. (2009) suggest that the volcanism continues for an even longer period, up to 3.5 billion years. Asymmetric models where the uppermost mantle in the nearside is more enriched in HPEs than that of the farside account for the long-lasting mare volcanism in the Procellarum KREEP Terrane (Wieczorek and Phillips 2000;Laneuville et al. 2013Laneuville et al. , 2018. The early expansion calculated in these models is, however, too small to account for the observed global expansion of ~ 0.5-4.9 km until 3.8 Gyr ago Laneuville et al. 2013;Zhang et al. 2013a;Sawada et al. 2016), or the early expansion continues too long, 1 billion years or longer (Zhang et al. 2013b).
In these previous models, only the thermal expansion/contraction has been considered, and the planetary expansion by melting has been neglected. The occurrence of dykes and volcanism in the early Moon (e.g., Hiesinger et al. 2003;Andrews-Hanna et al. 2013;Needham and Kring 2017), however, suggests that the effects of volume change from melting are important. Magma migration and transport of heat, mass, and HPEs by migrating magma should also be considered (e.g., Shearer et al. 2006). In this study, we extend the earlier one-dimensional spherically symmetric model of the lunar mantle (e.g., Solomon and Toksöz 1973;Solomon and Chaiken 1976;Kirk and Stevenson 1989) by adding magma-generation and migration and evaluate their effects on the radial change of the planet and volcanism; we start our calculation from a mantle structure predicted from earlier models of the mantle overturn (e.g., Hess and Parmentier 1995;Alley and Parmentier 1998). The model presented here is the first step toward a more comprehensive model of the lunar mantle evolution where mantle convection is also considered.

Model description
To understand the history of the Moon, we numerically calculate a thermal and structural evolution of internally heated one-dimensional spherically symmetric mantle caused by generation and migration of magma, thermal conduction, and internal heating. We take account of transport of heat, mass, and incompatible HPEs by migrating magma, too. The outer radius of the spherical shell r p is 1740 km, the radius of the Moon, while the inner radius r c is 383 km, the radius of the lunar core (e.g., Williams et al. 2001;Weber et al. 2011;Yan et al. 2015), and the crust of 43 km in thickness enriched in HPEs is placed along the surface boundary . The mantle materials are modeled as a mixture of olivine-rich materials A and the ilmenite bearing cumulates (IBC) B (e.g., Ringwood and Kesson 1976); their composition is denoted as A ξ * B 1−ξ * . The binary system is assumed to constitute a eutectic system with the eutectic composition of A 0.1 B 0.9 , which corresponds to the basaltic composition enriched in the IBC component (see Ogawa (2021) for the detail of the thermodynamic formulation). The composition of magma calculated below is almost always eutectic. Magma is generated when the temperature exceeds the solidus temperature. Generated magma migrates upward as a permeable flow through the coexisting matrix driven by the buoyancy of the magma (McKenzie 1984;Kameyama et al. 1996;Ogawa 2018). The temperature is fixed at T sur = 270 K on the surface boundary, while the core is modeled as a heat bath of uniform temperature. The thermal diffusivity of the crust is about half that of the mantle and the crust serves as a blanketing layer (Ziethe et al. 2009). Initially, the deep mantle is more enriched in HPEs and IBC-rich component than the mid-mantle as earlier models of mantle differentiation by the magma ocean and the subsequent mantle overturn suggest (e.g., Ringwood and Kesson 1976;Hess and Parmentier 1995). Some earlier studies suggest that a part of IBC-rich cumulates enriched in HPEs remain just beneath the crust after mantle overturn (Yu et al. 2019;Zhao et al. 2019). Instead of explicitly simulating these remains of cumulates, we assume that the crust is uniformly enriched in HPEs.

The basic equations
The continuity equation is: Here, U is the velocity of the matrix, u that of magma, and φ * the melt-content; U and u have an only radial component because of the assumed spherical symmetry. The relative velocity u − U is proportional to the density difference between magma and the coexisting matrix �ρ as: where k φ is the permeability that depends on φ * as k φ0 φ * /φ * 0 3 (McKenzie 1984); the φ * -dependence is truncated at its value at φ * = 0.4 . Here, k φ0 is the reference permeability, the value of which is given in McKenzie (1984) and Miller et al. (2014); φ * 0 is the reference melt-content; µ the viscosity of magma; g the gravitational acceleration (the parameter values are given in Table 1). The density difference between solid-phase and melt-phase �ρ depends on the chemical composition of the phases as well as the depth as: Here, ξ * s and ξ * l are the content of the end-member A in the solid-phase and the melt-phase, respectively, while ρ 0 the reference density; �V l /V 0 expresses the amount of density reduction by melting and depends on the depth r p − r as: where r is the radial coordinate, and the values of D * 0 , D * 1 , ς are chosen so that the solidus temperature calculated from Eq. 4 by the Clausius-Clapeyron relationship (see Eq. 5 below) becomes close to the solidus temperature of mantle materials (Katz et al. 2003;Garcia et al. 2011a, b). We set β * = 0.135 based on the density of olivine-rich end-member A ( ρ ξ * =1 = 3300 kg m −3 ; Elkins-Tanton et al. 2011), that of the IBC end-member B ( ρ ξ * =0 = 3744 kg m −3 ; Snyder et al. 1992;Shearer et al. 2006), and that of the IBC-rich eutectic composition of A 0.1 B 0.9 ( ρ ξ * = 0.1 = 3700 kg m −3 ; Elkins- Tanton et al. 2011;Yu et al. 2019).
The solidus temperature T melt depends on the lithostatic pressure P = ρ 0 g r p − r as: T 0 is the solidus temperature at the surface (Katz et al. 2003), and h the latent heat of melting.
The energy equation is: where the enthalpy h is written as: , C P is the specific heat of mantle materials, T the temperature, k = ρ 0 C P κ the thermal conductivity, κ the thermal diffusivity, and q the internal heating rate. We calculated T and φ * from the value of h calculated from Eq. 7 under the assumption that thermodynamic equilibrium always holds (see again Ogawa (2021) for the detail of the thermodynamic formation of our α m Thermal expansivity in the mantle 3 × 10 −5 K −1 α c Thermal expansivity in the core 9 × 10 −5 K −1 model). We assumed that matrix disintegrates and that a strong turbulent diffusion occurs with the eddy diffusivity of κ edd = 50 κ in largely molten regions with φ * > 0.4 . κ edd is assumed to gradually increase with increasing φ * as φ * 3 to avoid numerical instability. (We confirmed that the radial change we calculated was affected by less than 7% even when κ edd was increased from 50 to 100 κ at φ * > 0.4 ). Migration of magma is calculated in the energy equation as a transport of the latent heat of melting by migrating magma (the first term on the right-hand side of Eq. 7). We also assumed that the thermal diffusivity in the crust κ crst is lower than that in the mantle κ mantle to take account of the blanketing effect of the crust and regolith layer (e.g., Ziethe et al. 2009), that is: Following the method described in Ziethe et al. (2009), we assumed that κ crst /κ mantle = 0.48.
The core is modeled as a heat bath that has a uniform temperature T c , which changes with time as: where C pc is the specific heat of the core, ρ core the core density, V the volume of the core, and S the surface area of the core-mantle boundary. The heat flux of the coremantle boundary f is estimated at the bottom of the mantle as: The internal heating rate changes with time as: where q 0 = 14.70 p W kg −1 is the average initial heating rate at 4.4 Ga calculated from the total amount of (9) κ = κ crst at r > r crst , κ mantle at r < r crst . (10) HPEs in the current Moon, and the decay constant (see Table 2). q * tr changes with time as: , D * is the partition coefficient of HPEs between melt-phase and solid-phase, and q * tr is normalized so that V q * tr dV = 1 holds, where the integration extends over the entire spherical shell. The value of D * is fixed at 0.01 unless otherwise mentioned (see Table 1 and below).
In the estimate of the radial expansion of the Moon, we consider the effects of thermal expansion (Solomon and Chaiken 1976) and expansion by melting: the radial change R is given by: where α c and α m are the thermal expansivity of the core and that of the mantle, respectively; T c , T , and �φ * are the deviation of the core temperature, temperature in the mantle, and melt-content from the values assumed in the initial condition. Figure 1 shows an example of the initial distributions of enthalpy, temperature, internal heating rate, and composition. This initial condition is motivated by earlier models of the mantle overturn that is expected to have occurred after solidification of the magma ocean (e.g., Snyder et al. 1992;Alley and Parmentier 1998;Elkins-Tanton et al. 2011;Boukaré et al. 2018). Crystal fractionation in magma ocean is expected to have formed a dense layer composed of IBC and urKREEP (materials enriched in K, rare-earth elements, and P) at the top of the mantle. The gravitationally unstable mantle is suggested to have overturned to form a layer enriched in these materials at the base of the mantle (hereinafter, termed the "overturned layer"; e.g., Ringwood and Kesson 1976;Hess and Parmentier 1995;Zhao et al. 2019). The materials in the overturned layer are enriched in HPEs; the overturned materials with the composition of A 0.1 B 0.9 are assumed to be 7.5 times more enriched in HPEs than the bulk Moon q 0 (Hess and Parmentier 1995). We estimated the (13) distribution of the overturned materials in the mantle after the mantle overturn based on the distributions of internal heat source q = q init calculated in the models of Yu et al. (2019) and Zhao et al. (2019); the initial distribution of q * tr in Eq. 13 is calculated as q init /q 0 . We first approximated the distribution of HPEs they calculated by:

Initial distribution of temperature, HPEs, and bulk composition
where q m = 2.77 pW kg −1 is the internal heating rate of the olivine-rich materials beneath the layer of IBC and urKREEP before mantle overturn (Yu et al. 2019), l the thickness of the basal layer enriched in the overturned materials, and q is calculated from the assumption that the excess internal heating rate in the overturned layer Q l defined by:

satisfies,
Here, q crst is the average initial heating rate of the crust which is treated as a free parameter (see Fig. 1a). Then, we calculated the initial distribution of chemical composition ξ * b in the mantle from q init as: and ξ * b = 1 in the crust. Here, ξ * e = 0.1 is the eutectic composition which is the composition of the magma. On the other hand, the initial distribution of enthalpy h is: where (We illustrate the assumed enthalpy distribution in terms of the equivalent temperature T eq = h/C p in Fig. 1c). Here, the initial temperature in the bulk of the mantle T M is a free parameter, E sur and r l are constants arbitrarily chosen to be 6 × 10 4 J kg −1 and 550 km (Hess and Parmentier 1995), respectively, and the initial temperature of the core is T c = 1900 K (e.g., Ziethe et al. 2009;Zhang et al. 2013a, b); we assumed that the core is hotter than the mantle 4.4 Gyr ago (e.g., Alley and Parmentier 1998;Boukaré et al. 2018).

Free parameters
There are four free parameters in our model: the initial temperature in the mid-mantle T M , the ratio of reference permeability to melt-viscosity k φ0 /µ (see Eq. 2), the ratio of the concentration of the HPEs in the crust to the mantle F * crst , and thickness of the overturned layer l . (The asterisk stands for non-dimensional quantities). Instead of k φ0 /µ and l , we will use below their normalized forms M * = k φ0 ρ 0 g(r p −r c) µκ and L * = l/ r p − r c , respectively (Table 3). We varied the normalized permeability M * in the range estimated from the parameter values of Table 1, F * crst is estimated from a value adapted in earlier models (Konrad and Spohn 1997;Spohn et al. 2001), and L * is assumed from a value of the "overturned layer" in earlier models (Fig. 6

Thermal and structural evolution of the mantle The reference model
We present the radial change of planet calculated in the reference case in Fig. 2a. The assumed parameter values are T M = 1550 K , M * = 28 , F * crst = 8 , and L * = 1/5.5 . The radius of the planet changes with time by two reasons, thermal expansion and melting (see Eq. 15), as indicated by the green and light blue lines in the figure, respectively. Thermal expansion consistently causes radial contraction of the planet over the calculated 4.4 Gyr history. On the other hand, the history of radial change by melting consists of three stages: in Stage I, the planet rapidly expands until around 0.6 Gyr; in Stage II, the planet slowly expands until around 1.3 Gyr; in Stage III, the planet gradually contracts. In total, the planet expands by 1.23 km because of mantle melting in Stage I but shrinks in Stage III because of thermal contraction and solidification of the mantle; the planet shrinks by 1.08 km for the last 800 million years of the calculated history (the purple line in Fig. 2a).
A closer look at Fig. 2c and e shows more clearly how mantle melting contributes to the calculated radial expansion/contraction history of the planet. The temperature of the deep mantle rapidly increases owing to strong internal heating by the HPEs-enriched overturned layer nearby the core-mantle boundary assumed in the initial condition, and a large amount of magma is generated there. The upward extension of the partially molten region in the deep mantle induced by internal heating and upward migration of the generated magma (Fig. 2c, e) causes the conspicuous radial expansion of Stage I (Fig. 2a). The maximum melt-content of the region is around 0.4 and occurs at the top of the region. The migrating magma concentrates HPEs and the end-member B enriched materials to the top of the region ( Fig. 2d and f ). When the magma enriched in HPEs ascends to the uppermost mantle (radius r = 1300-1400 km), however, the effect of conductive cooling from the cold surface boundary overcomes the internal heating by the HPEs in the magma, and the magma solidifies (1.0 Gyr in Fig. 2b); the radial expansion of the planet almost stops accordingly (Stage II in Fig. 2a). According to decay and depletion of the HPEs, the partially molten region then gradually shrinks (Stage III in Fig. 2a).
An interesting feature of the melt-distribution shown in Fig. 2c, e is that two regions with elevated melt-content develop in the mantle after around 1.6 Gyr, one at the radius r around 600 km and the other at r around 1300 km. The shallower region arises because magma migrates upward and causes strong internal heating there by concentrating HPEs. On the other hand, the deeper region arises because upward migrating magma has extracted the compositionally dense end-member B from the region by 1.6 Gyr (see Fig. 2f ). Because of the depletion in the end-member B, the density of matrix becomes almost equal to that of magma that is more enriched in the end-member B (see, Appendix 1 for more detail). Owing to the almost neutral buoyancy of magma in the deep mantle, the partially molten region persists there until today (Fig. 2c).

Radial expansion/contraction by magma
To more clearly show how the radial expansion in Stage I and contraction in Stage III of the reference model occur, we calculated three more models, Models 1-3. In Model 1, we switched off magma migration in the reference model by assuming M * = 0 and considered only the effects of thermal diffusion, internal heating, and The ratio of the concentration of the HPEs in the crust to that in the mantle ( HPEs crst /HPEs mantle ) 4-32 L * The thickness of the overturned layer after mantle overturn (see Section "Initial distribution of temperature, HPEs, and bulk composition") 1/5.5-1 U et al. Earth, Planets and Space (2022) 74:78 melt-generation. In Model 2, we switched on heat transport by magma migration in Model 1 by assigning the default value assumed in the reference model to M * . We still neglected HPE-and material-transport by assuming D * = 1 in Eq. 13 and fixing ξ * b at its initial value; we also assumed β * = 0 in Eq. 3 to turn off the effect of compositional density contrast on magma ascent. In Model 3, we added HPE-transport to Model 2 by assigning the default value to D * (i.e., D * = 0.01 ). The difference between the reference model and Model 3 is whether or not the transport of the end-member B by migrating magma (Eq. 14) is taken into account.
In Model 1, the planet first shrinks until around 0.4 Gyr and then expands owing to melting of the deep mantle by internal heating (Fig. 3). A totally molten region develops at the base of the mantle after 0.4 Gyr and becomes Fig. 2 The reference model. a The plot of radius changes against time, b temperature profile as a function of time, and the evolution of the distribution of c temperature, d internal heating rate, e melt-content, and f bulk composition ξ b . The assumed parameter values are T M = 1550 K, M * = 28, F * crst = 8, L * = 1/5.5 . In a the green and light blue lines indicate the contribution of thermal expansion and melting, respectively, to the total radius change (the purple line). In b the gray and light-gray areas are the temperature distributions in today's lunar mantle inferred by Khan et al. (2006) and Khan et al. (2014), respectively; the gray lines are those estimated by Karato (2013) under the assumption that the mantle consists of dry olivine (the solid line) and that containing or wet olivine with the H 2 O-content of 100 ppm (the dash line). In c, the contour lines show the distribution of melt-content with the contour interval of 0.05 starting from 0 (indicated by the dashed line); in e the contour lines show the distribution of temperature thicker with time because the effect of conductive cooling from the cold surface boundary is negligible at such a great depth; the totally molten region occupies the deeper part of the mantle with radius r less than around 1000 km at 4.4 Gyr.
The large-scale melting shown in Model 1 does not occur in Model 2 (Fig. 4) where heat transport by migrating magma is considered. Migrating magma extracts heat from the deep mantle and transports it to the uppermost mantle. The carried heat is further transported to the cold surface boundary by conduction. The melt-content in the deep mantle is less than 0.1, and there is no totally molten region as observed in Model 1. Consequently, the planet radially expands only during the period of ~ 0.3 to ~ 0.7 Gyr.
When HPE-transport by migrating magma is also considered (Model 3 in Fig. 5), the partially molten region persists in the mantle only until around 2 Gyr and the entire mantle becomes solid after that. Most of the HPEs in the region migrate upward to the uppermost mantle within the first 0.5 Gyr (Fig. 5d). The magma solidifies quickly at the top of the partially molten region due to the strong conductive cooling from the surface boundary, and the region shrinks  Fig. 3; transport of heat, and the end-members A and B by migrating magma is not considered. In b the contour lines show the distribution of melt with the contour intervals of 0.01 starting from 0 (the dashed lines) U et al. Earth, Planets and Space (2022) 74:78 ( Fig. 5b). Note that the early radial expansion by melting observed in Model 3 is 0.17 km, considerably smaller than that observed in the reference model shown in Fig. 2a. The smaller expansion in Model 3 is a consequence of a faster upward migration of magma calculated in Model 3 compared with that calculated in the reference model. Magma velocity is proportional to the density difference between magma and the coexisting matrix (see Eq. 3, Appendix 1). The compositional density difference (the first term on the right-hand side of Eq. 3) is, however, neglected and melt is always substantially less dense than matrix, leading to the faster upward migration of magma in Model 3. The changes of the chemical composition by magma migration also have important influences on the history of the lunar radial expansion/contraction history.

Radial expansion/contraction at various parameter values
We explored how the radial expansion in Stage I and contraction in Stage III depend on the four free parameters, T M , F * crst , L * , and M * by carrying out an extensive parameter search as shown in Fig. 6.
The initial mid-mantle temperature T M influences the timing and magnitude of radial expansion in Stage I (Fig. 6a). A higher T M leads to an earlier occurrence of expansion with a smaller amplitude and a longer period of contraction (Stage III). This tendency arises because a higher T M implies a higher melt-content in the initial mantle, which causes earlier expansion of the partially molten region and larger radial contraction due to solidification in the shallow mantle. F * crst and L * strongly affect the presence or absence of the expansion. Lower F * crst and L * imply more HPEs in the deep mantle, leading to more expansion in Stage I (Fig. 6b, c). The permeability M * strongly affects the magnitude of expansion/contraction when M * ≤ 4 as shown in Fig. 6d because a higher permeability induces quicker drainage of magma from the mantle. (Note that the case at M * = 0 in Fig. 6d is Model 1 in Fig. 3). At M * ≥ 4 , however, this effect of increasing M * becomes saturated and the history of planetary expansion/contraction becomes almost independent of M * .

Discussion
The spherically symmetric mantle calculated in our model radially shrinks for the first ~ 300 million years because the magma generated in the shallow mantle by the assumed high initial temperature solidifies. Then, the mantle radially expands by up to 1.49 km ( M * = 12   U et al. Earth, Planets and Space (2022) 74:78 in Fig. 6d) for several hundred million years as magma is generated in the deep mantle and migrates upward to the uppermost mantle. There, magma is cooled and solidifies by conduction from the surface boundary. Migrating magma extracts HPEs from the deep mantle. After 1-2 Gyr since the start of the calculated history, the mantle radially contracts (~ 1 km for the last 800 million years) owing to its solidification and cooling that result from the HPE-depletion at depth (see Fig. 7).

Comparison with earlier models
The most important feature of our model is that the mantle radially expands by up to 1.49 km for the first hundred million years in many of the cases where the initial temperature in the mid-mantle T M is as high as ~ 1600 K because of the expansion caused by melting. In earlier one-dimensional models where expansion/contraction of the planet occurs only thermally, the initial temperature in the deep mantle must have been less than 600-1200 K at most to account for the early expansion of the Moon (Solomon and Chaiken 1976;Kirk and Stevenson 1989). Otherwise, the earlier researchers found that the calculated mantle monotonously shrinks or only slightly expands at the beginning of its history at most. The same conclusion has been reached in three-dimensional models where heat transport by mantle convection is considered, too. Although the global expansion of the early Moon is as large as 2 km in some earlier models (Zhang et al. 2013a, b), the consumption of latent heat by melting is neglected and the effect of thermal expansion is overestimated in these models. Besides, the early expansion of the model (see H50E100MR2500 of Zhang et al. 2013b, for example) continues for so long as 1 Gyr as is the case for our model that starts with the initial mantle temperature T M of 1400 K (Fig. 6a); the period of early expansion is substantially longer than that for the Moon. In a model where HPEs are assumed to be enriched in the uppermost mantle of the nearside (Laneuville et al. 2013), the farside mantle shrinks by about 5 km although up to 7 km expansion is observed on the nearside, and this model does not account for the observed global expansion. In our model where planetary expansion by melting is considered, we found a planetary expansion by up to 1.49 km for the first several hundred million years of the calculated history of the mantle (Figs. 2, 6), suggesting that melting of the mantle has played a crucial role in the observed expansion/contraction history of the Moon (Watters et al. 2010;Andrews-Hanna et al. 2013;Sawada et al. 2016;Matsuyama et al. 2021). Our models also show that a careful modeling of magma-generation and migration is crucial to realistically reproduce the melting history of the lunar mantle. In earlier one-dimensional models where the initial temperature is more than 1000 K lower than the solidus temperature in the deep mantle (Wood 1972;Solomon and Toksöz 1973;Toksöz and Solomon 1973), a partially molten region is first induced by the high initial temperature in the uppermost mantle and then shifts toward the deep mantle as it is internally heated by HPEs; the mantle becomes totally molten at its base today in some models (see Fig. 6 of Solomon and Toksöz (1973), for example). This melting history is similar to that of Model 1 where the mantle thermally evolves only by thermal diffusion and internal heating (Fig. 3). The melting history, however, substantially depends on magma-generation and transport of heat, HPEs, and materials by migrating magma. When heat transport by migrating magma is added to Model 1, a partially molten region extends upward to the uppermost mantle (Model 2; Fig. 4). When HPE-transport by migrating magma is further considered, the deep mantle becomes depleted in HPEs within around 0.5 Gyr since the beginning of the calculation, and the deep mantle becomes solid within 2 Gyr (Model 3; Fig. 5). In the reference model where changes in the chemical composition by magma segregation are added to Model 3, depletion of matrix in the end-member B allows a partially molten region to persist for 4 billion years in the deep mantle (Fig. 2). To account for the melting history of the lunar mantle that exerts control over the history of planetary expansion/contraction and mare volcanism, a careful modeling of magma-generation and migration and the resulting transport of heat, HPEs, and materials by magma is indispensable.

Comparison with the observed features of the Moon
The evolution of partially molten region observed in the reference model that accounts for the radial expansion/ contraction history of the Moon (see Fig. 2) meshes with the history of mare volcanism of the Moon. After the initial deepening to the depth of around 400 km until ~ 0.4 Gyr, the top of a partially molten region in our model becomes shallower and reaches the depth of around 300 km at around 0.8-1.5 Gyr owing to magma migration, and then again becomes deeper with time (Fig. 2c,  e). The deepening of the top of partially molten region until around 0.4 Gyr is consistent with the period when mare volcanism on the Moon was not so active (e.g., Hiesinger et al. 2003;Whitten and Head 2015), while its shallow depth at around 0.8-1.5 Gyr accounts for the peak of mare volcanism, which occurred at ~ 3.6 Gyr ago (e.g., Hiesinger et al. 2000Hiesinger et al. , 2003. The shrinking partially molten region observed after ~ 1.5 Gyr in our model may correspond to the deepening of the source region of mare volcanism inferred from the variation of Ti concentration in mare basalts (e.g., Hess et al. 1978;Pieters et al. 1980;Taylor 1982;Morota et al. 2011). Although a thermal diffusion model where magma generation is not considered also shows the deepening of partially molten region (Wieczorek and Phillips 2000), this earlier model does not account for the weak volcanism of the Moon before about 4 Gyr ago. (Note that this argument is based on an implicit assumption that the activity of the mare volcanism reflects the depth and amount of magma generation in the uppermost mantle. Solomon (1975), however, suggests that the mare volcanism declined since around 3 Ga despite that magma was generated in the mantle because the stress-state in the lithosphere changed and impeded extrusion of the generated magma. Besides, the thickness of the crust and the density contrast between magma and the crust also affect the rate of extrusive magmatism (e.g., Morota et al. 2009;Head and Wilson 2017;Taguchi et al. 2017;Wilson and Head 2017). A more careful modeling of magma migration in the crust and the lithosphere is called for to ultimately resolve the issue of the volcanic history of the Moon).
The profile of horizontally average temperature calculated at 4.4 Gyr in Fig. 2b is consistent with that of today's Moon obtained from analyses of its seismic data and electrical conductivity as a whole (e.g., Karato 2013; Khan et al. 2006Khan et al. , 2014. The presence of a partially molten zone at 4.4 Gyr in our model is also consistent with that in today's mantle of the Moon inferred from its seismic observations (e.g., Latham et al. 1973;Nakamura et al. 1973;Weber et al. 2011). However, the calculated temperature profile is higher than the observed one in the depth range from 700 to 1200 km. Besides, the molten region in our model occurs in the mid-mantle and is ~ 500 km in thickness with the degree of melting of up to 4% (Fig. 2c, e), while geophysical evidence suggests that the partially molten region occurs at the base of the mantle with its thickness around 150 km and its degree of melting 5-30% (Weber et al. 2011). The differences in the temperature profile and the depth, melt content, and the thickness of partially molten region call for further improvements of the models. In particular, it is important to consider the effects of heat transport by mantle convection.

Possible effects of mantle convection
Constraints on the thermo-chemical state of the lunar mantle just after planetary formation from the thermal history of the Moon have been studied for a long time (e.g., Solomon and Chaiken 1976;Pritchard and Stevenson 2000). Among the models presented in Section "Radial expansion/contraction at various parameter values", those calculated at the initial temperature in the mid-mantle T M ≈ 1600 K , the initial crustal fraction of HPEs F * crst ≤ 12 , and the extent of the enriched layer at the base of the mantle L * ≤ 1/2 are consistent with the observed expansion/contraction history of the Moon (Watters et al. 2010;Andrews-Hanna et al. 2013;Sawada et al. 2016;Matsuyama et al. 2021). To constrain these quantities more reliably, it is important to add the effects of mantle convection to our model. Many threedimensional models suggest that evolution of the mantle caused by mantle convection also plays an important role in the thermal and volcanic history of the Moon (e.g., Konrad and Spohn 1997;Spohn et al. 2001;Ziethe et al. 2009;Laneuville et al. 2013Laneuville et al. , 2018. In particular, the thermal history of the deep mantle is strongly influenced by mantle convection (Stegman et al. 2003;Zhang et al. 2013a, b): migrating magma which is generated in the deep mantle could drive mantle convection of the Moon owing to the buoyancy of the magma; compositional heterogeneity of the mantle caused by magma segregation can also drive mantle convection (see Appendix 1). In the future, numerical models of mantle convection where magma generation and magma migration are considered will be needed for comprehending the thermal and volcanic history of the Moon.

Conclusions
We developed a one-dimensional spherically symmetric numerical model of the thermal history of the lunar mantle, taking magma-generation and migration into account. We found that transport of heat, incompatible heat-producing elements (HPEs), and materials by migrating magma play a crucial role in the thermal history. The planet radially expands by up to 1.49 km during the first around 600 million years in the calculated history as magma is generated in the deep mantle and migrates upward to the uppermost mantle. The migrating magma extracts HPEs from the deep mantle, and the planet then radially contracts due to solidification and cooling of the mantle. The radial expansion/contraction history in our model is consistent with the observed Fig. 8 The evolution of the distribution of density difference between magma and matrix �ρ (Eq. 3) history of the Moon, and the melting history of the mantle in our model accounts for the mare volcanism of the Moon. In future works, it is important to extend our model by taking account of heat and mass transport by mantle convection.