Lunar Regolith Temperature Variation in the Rümker Region Based on the Real-Time Illumination

Abstract

Chang’E-5 will be China’s first sample−return mission. The proposed landing site is at the late-Eratosthenian-aged Rümker region of the lunar nearside. During this mission, a driller will be sunk into the lunar regolith to collect samples from depths up to two meters. This mission provides an ideal opportunity to investigate the lunar regolith temperature variation, which is important to the drilling program. This study focuses on the temperature variation of lunar regolith, especially the subsurface temperature. Such temperature information is crucial to both the engineering needs of the drilling program and interpretation of future heat-flow measurements at the lunar landing site. Based on the real-time illumination, and particularly the terrain obscuration, a one-dimensional heat equation was applied to estimate the temperature variation over the whole landing region. Our results confirm that while solar illumination strongly affects the surface temperature, such effect becomes weak at increasing depths. The skin depth of diurnal temperature variations is restricted to the uppermost ~5 cm, and the temperature of regolith deeper than ~0.6 m is controlled by the interior heat flow. At such a depth, China’s future lunar exploration is adequate to measure the inner heat flow, considering the drilling depth will be close to 2 m.

1. Introduction

The Rümker region is located in the northern Oceanus Procellarum (Figure 1) and is the proposed candidate landing region for the Chang’E-5 mission [1,2,3,4], since this region has a long and complex volcanic history. In the forthcoming Chang’E-5 mission (CE-5), the lander will carry a driller to the lunar regolith in order to collect regolith samples up to 2 kg. The penetration depth is planned to reach 2 m [5,6]. This mission will not intend to probe subsurface heat flow, but will still permit exploration of the subsurface regolith thermal properties of the moon. The lunar subsurface temperature was directly measured during the Apollo missions [7]. Samples returned from these missions suggest depth-dependent density and thermal properties [8]. In addition to in situ measurements, thermal infrared data from lunar orbiters have also been used to investigate the thermal properties of the lunar regolith [9,10,11,12,13].

Utilizing data from the Lunar Reconnaissance Orbiter (LRO) Diviner Lunar Radiometer Experiment, lunar thermal state was globally mapped, specifically the first measurements of thermal emissions in the lunar polar areas [9]. Based on these data, the research [12] studied the lunar equatorial surface temperatures and regolith properties. Using the same data, the work [13] examined the thermophysical properties of the topmost 10 cm of the global lunar regolith, and these results were consistent with predictions based on the one-dimensional model. The one-dimensional heat conduction equation was also used to investigate the thermal stability of near-surface ground ice on Mars [14]. In order to fit the temperature measurements, the study [15] constructed a one-dimensional two-layer conducting model. This model has a top layer close to 2 cm and is consistent with the in situ measurements from the Apollo mission. Other research further verifies the reasonableness of the one-dimensional approximation [12,13,16,17,18].

Standard finite difference approximation has been employed to validate the one-dimensional governing equation. Its deficiency lies in the treatment of the boundary condition, specifically in the three-dimensional problem. In such a case, a finite element method (FEM) is, in general, a better choice compared with the finite difference method. However, regarding the lunar regolith, the FEM encounters difficulties as the FEM grid is hard to generate due to the extremely small aspect ratio between the depth and width over the regolith. Therefore, the one-dimensional heat equation is generally considered, and the finite difference method is widely employed. In practice, the boundary condition greatly affects results in the one-dimensional approximation, particularly the solar illumination [12,13,16,17,18].

The lunar illumination aspect was first studied by Bussey et al. [19] to discover permanently shaded regions at poles. Later studies mainly focused on polar illumination by considering terrain obscuration [20,21,22]. This research found potential landing sites for future explorations but did not further study subsurface temperature variation. The recent studies [13,23] explored the lunar regolith thermophysical properties. They used a theoretical formula of the lunar orbital elements to estimate the solar hour angle and the incident solar flux for one lunar day. This theoretical hour angle and solar flux may not correspond to real-time values due to precession and nutation of the Moon’s rotation and are difficultly associated with a specified UTC (Universal Time Coordinated) time. Besides, their formula did not consider the effect of terrain obscuration, which needs to be considered in the temperature estimation. In practical drilling, the mechanical properties of lunar regolith play an important role in the processes of cutting regolith and removing cuttings from the bottom of a borehole [5]. As to the mechanical properties, one of the dependents is the temperature distribution. Therefore, it is necessary to investigate the regolith temperature.

In this study, we estimate subsurface regolith temperature distribution to support the design of drill tool and future lunar in situ exploration. For this purpose, we used the SPICE system [24] to obtain the subsolar point and elevation angle. Combined with the high-resolution Digital Terrain Model (DTM) from Lunar Orbiter Laser Altimeter (LOLA) [25,26], we estimated the real-time illumination state of a target location and the effect of terrain obscuration. Further combined with the recent one-dimensional thermophysical model, we explored the localized subsurface temperature beneath the Rümker region. We arrange this paper as follows. We introduce the theory in Section 2; we present the results and discussion in Section 3; and finally, we summarize our work and draw conclusions in Section 4.

2. Methods

2.1. Heat Conduction Equation

For the temperature T of lunar regolith, it varies with the regolith depth z and solar flux. The flux changes with the lunar local time t. Assuming the conducting layers have a thermal conductivity k, and supposing ∇ denotes gradient operator, the governing equation for the temperature is a form of heat conduction as follows:

$$\rho c_p\frac{\partial T}{\partial t}=\nabla \cdot k\nabla T,$$

(1)

in which ρ and c_p are the density and specific heat of the layers, respectively.

In the one-dimensional case, Equation (1) can be simplified as:

$$\rho c_p\frac{\partial T}{\partial t}=k\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right),$$

(2)

Recent work, especially work based on Diviner [12], discovered the accurate depth-dependent density of ρ. According to the study [13], we define density as:

$$\rho \left(z\right)={\rho }_d-\left({\rho }_d-{\rho }_s\right)\exp \left(-\frac{z}{H}\right),$$

(3)

in which ρ_s is the surface layer density value 1100 kg⋅m⁻³ [13], and ρ_d is the density values of the bottom layer 1800 kg⋅m⁻³ taken from the study [27]. H is a parameter set to 0.06 m and is used to govern the increase of regolith density from surface to bottom following previous work [12,13].

According to previous research [15,28,29], the thermal conductivity is a function of density and temperature, defined as:

$$k\left(T,\rho \right)=\left[k_d-\left(k_d-k_s\right)\frac{{\rho }_d-\rho }{{\rho }_d-{\rho }_s}\right]\left[1+x{\left(\frac{T}{350}\right)}^3\right],$$

(4)

in which x denotes radiative conductivity parameters (x = 2.7 according to the study [13]), k_s and k_d are the conductivities at the surface and bottom layers (k_s = 7.4 × 10⁻⁴ Wm⁻¹K⁻¹, k_d = 3.4 × 10⁻³ Wm⁻¹K⁻¹ according to the study [13]), and ρ is the density of regolith from Equation (3). Not only conductivity but also heat capacity varies with temperature; the study [13] derived a polynomial fit based on the data from the researches [30,31]. We here use their formula, expressed as

$$c_p\left(T\right)=c_0+c_{1}T+c_2{T}^2+c_3{T}^3+c_4{T}^4,$$

(5)

in which c₀, c₁, c₂, c₃ and c₄ are expanded coefficients. According to the study [13], we take their values as c₀ = −3.6125 J⋅kg⁻¹⋅k⁻¹, c₁ = 2.7431 J⋅kg⁻¹⋅k⁻², c₂ = 2.3616 × 10⁻³ J⋅kg⁻¹⋅k⁻³, c₃ = −1.234 × 10⁻⁵ J⋅kg⁻¹⋅k⁻⁴, and c₄ = 8.9093 × 10⁻⁹ J⋅kg⁻¹⋅k⁻⁵.

2.2. Boundary Conditions

There are two boundary conditions that must be set in order to solve Equation (2). One boundary condition is at the surface, and the other one is at the bottom layer of the regolith. At the surface, conduction and absorbed insolation are balanced against infrared emission to space. This principle is expressed by the following formula:

$$k{\frac{\partial T}{\partial z}|}_{z=0}+Q_s={\epsilon }_0\sigma T_s^4,$$

(6)

where the variable Q_s represents the solar heating rate, ε₀ denotes infrared emissivity (here equal to 0.95 according to studies [12,13]), and σ signifies Stefan–Boltzmann constant (σ = 5.67 × 10⁻⁸ Wm⁻²K⁻⁴). Following the studies [12,32] and assuming incident solar flux F_sun, we can get the solar heating rate as follows:

$$Q_s=\left[1-A_0-a{\left(\frac{4\alpha }{\pi }\right)}^3-b{\left(\frac{2\alpha }{\pi }\right)}^8\right]F_{sun},$$

(7)

in which a and b are the best fit constants. They are updated with the values (a = 0.06, b = 0.25) from the study [13]. A₀ signifies Bond albedo at normal solar incidence, which is an overall average reflection coefficient of a celestial body. According to the study [12], we take the albedo A₀ = 0.12. The incident flux F_sun depends on solar incident angle α, and distance r between the Moon and the Sun, expressed as ${\mathrm{F}}_{\mathrm{sun}}=\frac{{\mathrm{Sr}}_0^2}{{\mathrm{r}}^2}\cos \alpha$ The symbol S indicates the solar constant, and we take 1361 W⋅m⁻² from the research [33] as its value. The constant r₀ represents the mean distance between the Sun and the Earth, with the value of 1.49598261 × 10¹¹ m corresponding to S.

As shown in Equations (6) and (7), the solar incidence angle α must be estimated to evaluate solar heating rate. This angle depends on local surface slope λ, aspect, zenith angle θ_z, the surface azimuth γ, and the solar azimuth γ_s. The research [34] gives a detailed discussion on this topic, expressing the solar incidence angle α as follows:

$$\mathit{cos}\alpha =\mathit{cos}{\theta }_z\mathit{cos}\lambda +\mathit{sin}{\theta }_z\mathit{sin}\lambda \mathit{cos}\left(\gamma -{\gamma }_z\right),$$

(8)

Given the hour angle h between the target position A and subsolar point, the zenith angle θ_z following the study [34] is expressed as:

$$\mathit{cos}{\theta }_z=\mathit{sin}\delta \mathit{sin}\phi +\mathit{cos}\delta \mathit{cos}\phi cosh,$$

(9)

where δ is the latitude of the subsolar point, and φ designates the latitude of the target position A. For solar azimuth angle γ_s, we adopted the geometric formula from the work [34].

At the bottom layer z = z_b (here, z_b = 2 m), the heat flow Q is determined by the temperature gradient as follows:

$${\frac{\partial T}{\partial z}|}_{z=z_b}=\frac{Q}{k},$$

(10)

where k denotes the thermal conductivity shown in Equation (4). The value of heat flow (Q = 18 mW m⁻²) was taken from the Apollo measurements [7].

2.3. Illumination

The incident flux in Equation (7) depends on solar illumination. To obtain the flux, we need to know the real-time subsolar point and elevation angle. These parameters are estimated using the SPICE system [24], developed by the NASA Navigation and Ancillary Information Facility (NAIF). The elevation angle is used to determine illumination aspect [192021,35], defined as the angle between the horizon and sunlight direction. As is shown in Figure 2, the dotted line AP signifies the horizon, and the sunlight direction by the line AD. An elevation angle less than zero means the sun is sinking below the horizon at nightfall. In such cases, the boundary conditions for heat conduction in Equation (2) will not include solar flux; however, in the daytime, the elevation angle is required to solve Equation (2).

The effect of terrain obscuration is non-negligible when the elevation angle is considered in the daytime. Such an effect can be determined according to the subsolar point and terrain maps. As shown in Figure 2, the sunlight direction is estimated according to subsolar point. Supposing A denotes the target position on the Moon, we can determine a candidate point B, whether blocking A or not. Figure 2 gives a schematic for the terrain obscuration. The dashed arc denotes an equivalent curved surface with the same radius of point A. B signifies the point used to determine the terrain obscuration, and its projection on the equivalent surface is represented as C. β is the central angle between the points A and B, and its value can be calculated through the two points’ longitudes and latitudes.

According to the study [35], we can get the length of line DC in Figure 2 as follows:

$$DC=\frac{AC\cdot \sin \left(\frac{\beta }{2}+\theta \right)}{\mathit{cos}\left(\beta +\theta \right)},$$

(11)

The length of line BC can be estimated according to DTM from LOLA. The LOLA team assembled the LOLA data as a global Digital Elevation Model (DEM) in cylindrical and polar stereographic projects. We employed the cylindrical version and applied the LDEM_512_00N_45N_000_090 and LDEM_512_45N_90N_000_090 (1/512° × 1/512°) in our calculations. Supposing sunlight can arrive at point A, the length of line DC needs to be larger than BC. Considering curvature of the lunar surface, we estimated several candidate points along the arc of AC. In practice, a value of β lower than 8 degrees is sufficient to determine the terrain obscuration effect [35]. In the case of sunlight even blocked during the local daytime, we will not consider solar heating in the boundary conditions of Equations (1) and (2).

3. Results and Discussion

The real-time sunlight direction and elevation angle are estimated with the SPICE system. They are used in our illumination program to test terrain obscuration. To validate this program, Figure 3 gives a comparison between the observed image (Figure 3a) and the synthetic image (Figure 3b) from LOLA.

The illumination aspect is displayed over the Aristarchus crater, centered at (23.73°N, 312.51°E). Figure 3a comes from the narrow-angle camera (NAC) image (http://lroc.sese.asu.edu/posts/291), displaying illumination observed during the morning with the Sun shining from the right. Figure 3b exhibits the corresponding illumination estimated by our program; its contour lines indicate that the minimum elevation of the Aristarchus crater is close to −3.3 km. The color bar demonstrates the range of relative intensity of illumination. The maximum intensity of illumination in Figure 3b arises at the inner wall of Aristarchus crater and is similar to the intensity of illumination in Figure 3a. The shadowed areas in the two figures resemble each other, which mirrors obscuration by terrain features. The resemblance of illumination in the two figures indicates the effectiveness of our program; we employed it to investigate subsurface temperature variation on the Moon.

To demonstrate the terrain obscuration affection on temperature, we also studied the surface temperature over Aristarchus with (Figure 4a) and without (Figure 4b) terrain obscuration in Figure 4 at the lunar local time t_m = 6:39:31.

Figure 4a,b denotes the surface temperature with and without consideration of terrain obscuration, and Figure 4c represents their temperature residual. Figure 4a indicates that the temperature on the right inner wall is quite cold at 90 K, due to the obscured sunlight, while Figure 4b demonstrates a relatively high temperature to 190 K at the same location. Figure 4c shows the great residual of temperature existing at the right inner wall, where the sunlight is obscured by the surrounding terrain. The right inner wall in Figure 4b is in fact shadowed due to no sunlight arriving; thus, the temperature showing there is unbelievable. It can be concluded that the ignorance of terrain obscuration will lead to untruthful results.

To further explore the temperature variation with depth, Figure 5 gives temperature variation with depth at the point (313°E, 24°N), which is shown in Figure 3 with a sign of white three-star.

Figure 5 denotew profiles of temperature–depth at the lunar local time tm = 06:39:31, tm = 12:00:00, and tm = 17:00:00. The solid black line with circle signifies temperature variation with terrain obscuration, while the dashed pink line with square represents the case without terrain obscuration. Figure 5a indicates that the terrain obscuration affects temperature greatly near to the surface, but the rest of the panels imply that these two cases with and without terrain obscuration have the same performance with the sunlight arriving. It is can be concluded that the terrain obscuration affects the near-surface temperature greatly, but this effect disappears with the depth increased and the sunlight not obscured.

To investigate temperature variation over Rümker region with time, not only a boundary condition, but also an initial condition is needed to solve the heat equation. According to the study [13], an appropriate initial temperature influences temperature variation at depth over time. To reduce errors, we continuously implemented estimation for several lunar days, each one corresponding to about one month on Earth. We selected the most recent lunar day to analyze. In order to take the lunar day associated with UTC on Earth, we took the analyzed lunar days expressed as UTC months. We also supplied the corresponding lunar local time. Although the CE-5 is scheduled to launch in 2020, the exact date of landing implementation is difficult to predict now. Considering previous Chinese lunar exploration projects such as Chang’E 3 (CE−3) and Chang’E 4 (CE−4), we took five months for the project implementation. The errors generated by initial conditions will be reduced by the last month of the mission, so we selected the last month as the research period and analyzed the temperature. The research period will last from UTC 2020-0-3T00:00:01 to 2020-06-01T23:45:01, corresponding to lunar local time from 04:23:01 of the first day to 04:22:31 of the next day. We sampled the days and give the corresponding surface temperature in Figure 6.

The temperature distribution at the time of UTC 2020-0-3T17:15:01 is displayed in Figure 6a, which shows a minimum temperature close to 80 K at night. As the sunlight in Figure 6b arrives in the morning, the temperature increases in the research area at the time of UTC 2020-05-05T04:30:01. As the rising of the Sun, it continues to increase at the time of UTC 2020-06-06T13:15:01 in Figure 6c. At noon, it even climbs to a maximum value of more than 400 K in Figure 6d, corresponding to the time of UTC 2020-05-13T06:30:01. After that, the temperature is declining as the sun is setting, which is exhibited in Figure 6e at the time of UTC 2020-05-18T17:00:01. Finally, the lunar surface becomes quite cold as the sun is sinking behind the horizon. This case is illustrated in Figure 6f, corresponding to the time of UTC 2020-05-19T20:00:01.

Adopting the same research period, we could analyze subsurface temperature variation. We tested various depths and found that the subsurface temperature does not change at increasing depths. A shallow depth varies from the surface temperature, but at deeper levels, the temperature is nearly constant. After several tests, we found a trade-off depth of 0.0656 m, which demonstrates varying but not significantly different temperatures. This depth illustrates subsurface temperature variation, as shown in Figure 7.

As can be seen from Figure 7a–c, the subsurface temperature changes a little, but remains in a range close to 200 K, unlike the variation seen in Figure 6a–c. The subsurface temperature increases gradually until noon, as shown in Figure 7d. It climbs to a maximum value in Figure 7e, close to 280 K, even when the sun is setting. It does not yet decline greatly in the evening of Figure 5f. This result confirms that the lunar regolith tends to slow down thermal conduction [13,16]. Such thermal insulation could be ascribed to the porous, pulverized, and granular attributes of the regolith, as investigated by the study [23]. Its improvement is the prediction of a lower thermal conductivity but a higher specific heat than previous research [9,10,11,12,13]. Therefore, the subsurface temperature varies more slowly than surface. The differences of temperature variation between Figure 6 and Figure 7 are a mirror of such insulation of lunar regolith.

To further probe the subsurface temperature, we investigated the temperature variations in the latitude and longitude directions in Figure 8 and Figure 9, respectively. They all pass through the CE-5 landing sites proposed by the research [2]. These landing sites are indicated with red star and red square, centered at (43.68°N, 49.85°W) and (41.41°N, 58.53°W) in Figure 1. Figure 8 and Figure 9 give temperature variation with time on cross sections and display their results at the same time as those in Figure 6 and Figure 7. Figure 8a–c shows the temperature variation in the longitude direction. The illumination heats the surface layers within 0.05 m as the sun rises. At noon, as seen in Figure 8d, the surface temperature reaches the maximum value, but the subsurface temperature varies little, especially those layers below 0.25 m. After midday, the high temperature at the surface conducts down to the lunar inner layers, which heat up as shown. As the sun sets in Figure 8e,f, the surface temperature declines greatly, but subsurface temperature falls slowly. These results demonstrate the thermal insulation properties of the lunar regolith.

Figure 8a–f shows several stripes of temperature. According to the study [13], a thermal inertia depends on local topographic features. It is used to describe the resistance of materials to changes in temperature. The temperature stripes in Figure 8a–f are a mirror of different topographic features. To display such correlations, Figure 8d gives a black profile line of the surface topography along latitude. This line varies in a range relative to the topography and is not to scale by distance. We also added such a profile of surface topography in Figure 9d in the latitude direction to demonstrate its correlation to the temperature stripes. Such stripes can be found in Figure 9, especially in Figure 9d–f. The stripes with significant variation by depth are on parts of Mons Rümker, which is located on the left part of Figure 9a–f. From the correlation between the topographic terrain and temperature stripes, it could be concluded that the temperature stripes are a mirror of anisotropic thermal inertia. As the sun rises in Figure 9a–c, the surface temperature increases, but inner temperature at layers under 0.25 m changes little. Similarly, in Figure 9d–f, the surface temperature declines gradually, but inner temperature varies slowly, especially those beneath layer of 0.25 m; these results indicate that the lunar regolith has thermal insulation properties. Thus, we provide results within the depth of 0.275 m in Figure 8 and Figure 9, although the temperature estimated ranges within the depth of 2 m.

To investigate the subsurface temperature variation beneath 0.25 m further, Figure 10 gives the variations at the proposed two landing sites [2]. As to the point of red square in Figure 1, Figure 10a indicates that the subsurface temperature at −0.2679 m varies between 235.97 K and 238.64 K in the whole research period, but alters a little below this depth. It even keeps in constant at the depth of −0.5691 m. Such constant subsurface temperature, as shown in Figure 10b, is consistent across the other landing sites, marked with a red star in Figure 1. Considering all the results shown in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10, it can be concluded that the subsurface temperature stays constant under the depth of −0.6 m. As more lunar heat flow measurements needed are to constrain lunar thermal evolution [36], remeasurement at different places is needed in future explorations. Based on the CE-5 drilling experience, China’s next exploration could consider a direct measurement of heat flow, and a drilling depth close to 2 m is adequate to explore the inner heat flow beneath −0.6 m.

The depth of constant temperature layer found above would depend on H parameter and lunar inner heat flow Q. We next take an analysis of these varied parameters. In the previous study [13], the H parameter ranged from 0.02 m to 0.09 m. In this paper, we used a mean value H = 0.06 m. As to the heat flow [17], the measured heat flow from the Apollo 15 was 21 ± 3 mWm⁻², but a low value of 15 ± 2 mWm⁻² was found for the Apollo 17. The later studies [37,38] found a global value low to 12 mWm⁻². We used a mean value of heat flow Q = 18 mWm⁻². As to the layer of constant temperature found above, it is necessary to discuss the uncertainties of the H parameter and heat flow. We here give a profile of temperature−depth for the CE-5 landing site centered at (41.41°N, 58.53°W). The profile is shown in Figure 11, which includes three panels (a–c) representing the lunar local time t_m = 05:00:00, t_m = 12:00:00, and t_m = 17:00:00.

The solid line with circle represents the result for the H parameter equal to 0.06 m, which is the value used in our calculation. The dashed line with the cross denotes the case for H = 0.02 m, while the dashed blue line with square is for H = 0.09 m. The dashed blue line with the upper triangle implies the case with a heat flow low to 12 mWm⁻², while the dashed pink line with the lower triangle indicates the case with a large heat flow Q = 21 mWm⁻². It can be found that the small and large heat flow changes the profiles little, compared with the mean value Q = 18 mWm⁻². The apparent deviation of the profile is mainly contributed by the small H parameter (H = 0.02 m) in the afternoon. However, all the profiles merge with each other with the depth increased, especially low to −0.3 m. This result suggests that the changes of H and Q have no influence on temperature below −0.6 m.

4. Conclusions

In this paper, we studied the lunar regolith temperature variation in the proposed Rümker landing region of CE-5. We investigated the real-time illumination based on the NASA’s SPICE system. Combined with local surface slope and aspect, zenith angle, surface azimuth, and solar azimuth estimations, the real-time solar heating rate was evaluated and used to calculate the relative intensity of illumination in the Aristarchus crater. Our synthetic figure closely resembles the observed NAC image, indicating the effectiveness of our calculations. We also used the real-time solar heating flux to study the temperature variation. The results demonstrate that solar illumination affects the surface temperature, becoming weaker with the increasing depth. In general, the subsurface temperature changes abruptly with the solar illumination around −5 cm. However, temperature anomalies like stripes occur in some places; these can be ascribed to anisotropic thermal inertia. The subsurface temperature fluctuates a little as the depth increases, but remains constant beneath a depth of −0.5691 m. The CE-5 drilling depth will approach 2 m, and thus is adequate to explore the heat flow beneath −0.6 m. This work will be helpful for future Chang’E landing missions to South Polar Region to investigate the subsurface interior temperature variation.