Inversion of the Lunar Subsurface Rock Abundance Using CE-2 Microwave Brightness Temperature Data

Abstract

The rock strongly affects the surface and subsurface temperature due to its different thermophysical properties compared to the lunar regolith. The brightness temperature (TB) data observed by Chang’E-1 (CE-1) and Chang’E-2 (CE-2) microwave radiometers (MRM) give us a chance to retrieve the lunar subsurface rock abundance (RA). In this paper, a thermal conductivity model with an undetermined parameter $\beta$ of the mixture has been employed to estimate the physical temperature profile of the mixed layer (rock and regolith). Parameter $\beta$ and the physical temperature profile of the mixed layer are constrained by the Diviner Channel 7 observations. Then, the subsurface RA on the 16 large (Diameter > 20 km) Copernican-age craters of the Moon is extracted from the average nighttime TB of the CE-2 37 GHz channel based on our previous rocky TB model. Two conclusions can be derived from the results: (1) the subsurface RA values are usually greater than the surface RA values retrieved from Diviner observations of the studied craters; (2) the spatial distribution of subsurface RA extracted from CE-2 MRM data is not necessarily consistent with the surface RA detected by Diviner data. For example, there are similar RA spatial distributions on both the surface and subsurface in Giordano Bruno, Necho, and Aristarchus craters. However, the distribution of subsurface RA is obviously different from that of surface RA for Copernicus, Ohm, Sharonov, and Tycho craters.

1. Introduction

Rocks are one of the predominant microrelief forms of the lunar surface. Understanding the distribution of rock is significant for the following reasons: (1) The rocks on the lunar surface are potential hazards to the lander or rover [1,2]. Therefore, the knowledge of the distribution of rock abundance (RA) on the lunar surface can help us to determine a suitable landing site in the future. (2) The relationship between rocks and craters’ impact melt and ejecta and crater morphology implies that the distribution of the rocks on the lunar surface may include the potential formation and evolution information of the lunar crust layer [3,4,5]. (3) The correlation between RA on the lunar surface and crater model ages provides a new method for dating young lunar craters [6,7]. (4) The rocks on the lunar surface and subsurface can provide a good understanding of the lunar thermal evolution [8,9].

To date, there have been many estimations of lunar rocks using both qualitative and quantitative methods. A lot of data are used to estimate the lunar rocks in a qualitative way, including thermal infrared data [5,8,10,11], optical data [3,4,5,12,13], and radar data [14,15]. Meanwhile, thermal infrared data [8] and optical data [12] can also be used to quantitatively constrain the lunar rocks. As far as we know, the first global map of surface RA, with the meter and larger size, have been made using Diviner infrared thanks to the large thermal emission differences between lunar regolith and rocks. At present, we know little about the subsurface RA of the Moon. The microwave TBs of the global Moon, detected by the CE-1 and CE-2 MRMs at four channels: 3, 7.8, 19.35, and 37 GHz [1], give us a good chance to constrain the subsurface RA due to their longer penetration compared to thermal infrared data. Using the differences in both the dielectric loss and the thermal properties of the lunar rocks and regolith, the subsurface RA could be quantitatively constrained by the microwave TBs [16,17]. Wei et al. [18] recently retrieved the subsurface RA at several craters from the CE-2 TB data by combing 7.8 GHz and 19.35 GHz channels. In their work, the same mixture model as that of Hu et al. [17] is used to calculate the mixture’s thermal physical properties, including density, thermal conductivity, and heat capacity at all depths.

The rocks are proven to have a significant influence on the lunar TBs [9,16,17,19,20,21]. Our previous rock model, proposed by Hu et al. [17], is employed in this paper, including the rock effects on both the physical temperature and the dielectric loss of the lunar subsurface layers. In addition, a thermal conductivity model with an undetermined parameter $\beta$ of the mixture was employed to estimate the physical temperature profile of the mixed layer (rock and regolith). The data applied in this paper are the level 2C TB data after pre-calibrations collected from CE-2. Their calibration and a discussion of their accuracy can be found in [22,23,24,25,26,27,28]. Recent recalibration works [26,27] suggest that the pre-calibrated TB data at 37 GHz can well fit both the infrared and TB models and the correction is smaller than other channels.

The inversion algorithm is another factor that affects the inversion results. Meng [29] inverted the thickness of the lunar regolith using the Backpropagation Neural Network method. However, the training data selection is a difficult problem due as there are few samples for network training. The Root Mean Square (RMS) algorithm [8,30] is usually used to solve the retrieved problem. The first map of global surface RA distribution was obtained using this algorithm. Fa et al. [30] also used this algorithm to estimate the thickness of the lunar regolith. Therefore, the RMS method is employed in this work.

In this paper, we first analyzed the effects of the parameter $\beta$ on the surface temperature and constrained the best fitting value of $\beta$ using the observed data of the seventh channel of Diviner. Then, we discussed the effect of RA on surface temperature and microwave TB based on this best-fitting value of $\beta$. Finally, we extracted the subsurface RA values on the 16 large (Diameter > 20 km) Copernican-age craters of the Moon from the CE-2 37 GHz midnight TB. We also conducted a simple uncertainty analysis.

2. Data Sets

The Diviner experiment is a nadir-pointed, push broom scanning radiometer with two spectral channels for reflected solar radiation, each from 0.35 to 2.8 μm, and seven channels for infrared emission, spanning from 7.55 to 400 μm [31]. It is designed to globally map surface albedo and temperature over the lunar day and year. The rock-free regolith surface temperature was estimated using the algorithm described by Bandfield et al. [8]. The data were collected from 5 July 2009 to 10 October 2016. Data were restricted to local times of from 1930 to 0530, with solar incidence angles greater than 90 degrees, with a spatial resolution of 128 pixels per degree (ppd). The Diviner Lunar Radiometer Experiment Reduced Data Records (RDR) data at a single channel (T₇, 25–50 μm) were used to constrain the surface boundary condition and determine the parameter $\beta$ due to its high signal-to-noise ratio over the full range of equatorial surface temperatures [32].

The Diviner’s measurements were also used to infer surface RA. The global map of surface RA of meter size and larger was drawn by Bandfield et al. [8] with the Diviner measurements. The global surface RA data have a spatial resolution of 128 ppd bin size (~125–250 m/pixel). Most of the average surface RA values within the $\pm 60°$ latitudes are smaller than 0.01. In some fresh craters, the average surface RA values are greater than 0.01, such as the Giordano Bruno crater, which has an average surface RA of up to 0.1.

The iron (FeO) and titanium (TiO₂) contents on the lunar surface affect the microwave TB by changing the dielectric loss tangent of the materials in the mixed layers [33]. They can be obtained from returned samples or remote sensing analysis. They have also been derived from the gamma-ray wavelength range of the lunar prospector data with a spatial resolution of 150 km × 150 km [34]. Recently, Sato et al. [35] derived a new TiO₂ map based on the near-global Wide Angle Camera (WAC) aboard the LRO multispectral mosaic. In this work, the used lunar surface global iron (FeO) and titanium (TiO₂) contents map (100 m/pixel) were derived by Lucey et al. [36,37,38,39] from Clementine ultraviolet-visible data.

The level 2C MRM data after pre-calibrations, collected from the CE-1 and CE-2, are published in the Lunar and Planetary Data Release System. A detailed description of data accuracy and data calibration can be found in [25,26,27,28,40,41]. The average normalized CE-2’s 37 GHz TBs at midnight (removing the lunar phase effect) is used to constrain the subsurface RA, which has a spatial resolution of 5 pixels per degree (ppd).

3. Models

Here, we treat each location of the Moon as a discrete vertical column of a mixture layer consisting of rock and regolith. The microwave TB of each vertical column strongly depends on its physical temperature profile and dielectric loss characteristics. Physical temperature is significantly affected by thermal inertia at all depths. Thermal inertia is a key quantity that is widely used in planetary science because it is related to the intrinsic physical properties of surface and subsurface materials. In general, the rock has a larger thermal inertia and loss tangent than the lunar soil (regolith and dust). In this work, we use “subsurface RA” to express the average volumetric fraction of rocks or rock-like materials within the column of each location detected by CE-2 MRM data. The subsurface RA can be extracted from CE-2 MRM data by minimizing the differences between simulated microwave TB and the observed TB of CE-2 37 GHz. Here, the scattering effect of the meter size and larger rock fragments could be neglected due to the wavelength (millimeter scale) of the 37 GHz channel being small [16].

The one-dimensional heat conduction model employed to calculate the physical temperature profile is described in this section. The surface energy flux and albedo affected the surface temperature are presented. The mixed models of physical properties are also described. As we do not have precise estimations of the lunar regolith thickness for Copernican-age craters, the maximum penetration depth $d_0$ of 37 GHz into the mixed layer is about 0.5 m. This can be estimated from the permittivity of the mixed layer with the formula $d_0=\lambda \sqrt{{\epsilon }_{mix}^{\prime }}/{\epsilon ″}_{mix}$, where $\lambda$ is the wavelength (m), ${\epsilon }_{mix}^{\prime }$ and ${\epsilon ″}_{mix}$ denote the real and imaginary part of the permittivity in the mixed layers, respectively. Therefore, we can simply assume the thickness of the mixed layer is about 1 m, which is deep enough to cover the penetration depth.

The rock TB model proposed by Hu et al. [17] is used to calculate the microwave TB of each column. In order to avoid repetitive descriptions, the multilayer microwave radiation transfer model is not described. Here, we only present the loss tangent models of the lunar regolith and rocks because they are the main model influence on microwave TB calculation.

3.1. Heat Conduction Model

The temperature profile can be obtained using the finite difference method to solve the 1-D time-dependent heat conduction equation [42,43,44,45]:

$$\rho \left(x,T\right)C\left(x,T\right)\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left[K\left(x,T\right)\frac{\partial T}{\partial x}\right]$$

(1)

where x (m) and $t$ (h) stand for depth and time; $\rho (x,T)$, $C(x,T)$, and $K(x,T)$ stand for the density ($\mathrm{g}·\mathrm{c}{\mathrm{m}}^{-3}$), heat capacity ($\mathrm{W}·{\mathrm{g}}^{-1}·{\mathrm{K}}^{-1}$), and thermal conductivity ($\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1}$) of the material, respectively. They are the functions of the lunar physical temperature $T\left(\mathrm{K}\right)$ and vertical variable of depth $x$ (m). The upper boundary condition depends on the energy on the lunar surface, which can be written as:

$$K_s\frac{\partial T}{{\partial }_x}{\left.\right|}_{x=0}+Q_s=\epsilon {\sigma }_B{T}_s^4$$

(2)

where $K_s(\partial T)/(\partial x)$ is the heat energy transmitted from the surface to the subsurface, $T_s$, $\epsilon =0.95$, and ${\sigma }_B(\mathrm{W}·{\mathrm{m}}^{-2}·{\mathrm{K}}^{-4})$ are the surface temperature, the infrared emissivity of lunar regolith, and the Stefan–Boltzmann constant, $Q_s=\left(1-A\right)F_{\odot }$ is the surface energy flux, which is equal to the solar heating rate, and $F_{\odot } (\mathrm{W}·{\mathrm{m}}^{-2})$ is the incident solar flux, which can be written as [9,17,44,46,47]

$$F_{\odot }\left(t\right)=\frac{S_{sun}}{4\pi ·R^2}\left(sin\varphi sin\delta +cos\varphi cos\delta \cos h\right)$$

(3)

where $S_{sun}=3.846\times {10}^{26}$ W is the solar luminosity [48], $\varphi$ is latitude, $\delta$ is the solar declination angle, which is usually set as 0 [49], and the hour angle is $h=2\pi t/P$, with $P$ as the length of the synodic day. R is the distance between Sun and Moon, which varies over time and has the following form

$$R=r_1·\cos \left(\frac{2\pi t}{29.53\times 12}+\pi \right)+r_2$$

(4)

where $r_1=3.85\times {10}^8$ m, $r_2=1.496\times {10}^{11}$ m, they are the distance from the Moon and Sun to Earth, respectively. The dependence between albedo and solar incidence angle $\theta$ is given by [9,19,32]

$$A\left(\theta \right)=A_0+a{\left(\frac{\theta }{\pi /4}\right)}^3+b{\left(\frac{\theta }{\pi /2}\right)}^8$$

(5)

where $A_0=0.12$ is the normal bolometric Bond albedo; coefficients $a$ and $b$ have a Hayne’s value of $a=0.06$ and $b=0.25$. At the lower boundary $d,$ the condition can be written as:

$$\frac{\partial T}{\partial x}{\left.\right|}_{x=d}=-\frac{Q}{K_d}\ll 1$$

(6)

where $Q \left(\mathrm{W}·{\mathrm{m}}^{-2}\right)$ is the internal lunar heat source; $K_d$ and $\partial T/\partial x{\left.\right|}_{x=d}$ are the thermal conductivity at $x=d$, which is sufficiently deep below the subsurface, and the temperature gradient, respectively.

3.2. Physical Parameters of the Mixed Layer

The relationship between the density $\rho$, subsurface $RA,$ and depth $x$ can be written as [17]:

$${\rho }_{mix}\left(x\right)=\left[{\rho }_d-\left({\rho }_d-{\rho }_s\right)e^{\frac{-x}{H}} \right]\left(1-RA\right)+RA·{\rho }_r$$

(7)

where ${\rho }_{mix}$ represents the mixed density, (${\rho }_d-\left({\rho }_d-{\rho }_s\right)e^{\frac{-x}{H}}$) is the density of regolith at a depth $x$ m, $RA$ represents the volumetric fraction of rocks at the penetration depth of CE-2 37 GHz, and $H=0.06$ (m) is the average H-parameter. ${\rho }_s=1300 \left(\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-3}\right), {\rho }_d=1800 (\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-3}),$ and ${\rho }_r=2940 (\mathrm{k}\mathrm{g}·{\mathrm{m}}^{-3})$ stand for the density of the lunar surface, balance depth, and rocks, respectively, which are obtained from Horai et al. [50]. From their investigation, the heat capacity of rocks can be written as:

$$C_{rock}={\alpha }_1+{\alpha }_{2}T+{\alpha }_3{T}^2+{\alpha }_4{T}^3$$

(8)

where ${\alpha }_1=-154.9, {\alpha }_2=4.983, {\alpha }_3=-8.207\times {10}^{-3}, {\alpha }_4=5.192\times {10}^{-6}$, and $C_{rock}$ is the average heat capacity of lunar rocks. The depth-dependent heat capacity of the mixed layer has the following forms:

$$C_{mix}=\left(c_1+c_{2}T+c_3{T}^2+c_4{T}^3+c_5{T}^4\right)\times \left(1-RA\right)+C_{rock}\times RA$$

(9)

where $C_{mix}$ represents the mixed heat capacity of the mixed layer, $\left(c_1+c_{2}T+c_3{T}^2+c_4{T}^3+c_5{T}^4\right)$ is the heat capacity of lunar regolith, and $c_1$, $c_2$, $c_3$, $c_4$, and $c_5$ are equal to −3.6125, 2.7431, $2.3616\times {10}^{-3}$, $-1.234\times {10}^{-5}$, and $8.9093\times {10}^{-9}$ [9], respectively.

The correlation between the lunar regolith thermal conductivity $K_{reg}$, depth, $x$, and temperature, $T$ can be written as:

$$K_{reg}\left(x\right)=K_d-\left(K_d-K_s\right)e^{\frac{-x}{H}}+2.7K_s·({\frac{T\left(x\right)}{350})}^3$$

(10)

where $K_s=0.00074$ and ${ K}_d=0.0034$ are the thermal conductivity of the regolith on the surface and the bottom, respectively [9]. In the work of [17,18], the same linear mixture model as used for density and heat capacity was also used to determine the thermal conductivity of the mixed layer. However, Horai et al. [51] suggested that the mixed thermal conductivity of two-phase materials is nonlinearly mixed. Here, we simplify the mixed layer consisting of two-phase material rock and regolith, and the mixing forms can be expressed as:

$$K_u=K_r\left(1+\frac{1-RA}{\frac{RA}{3}-1/(1-r)}\right)$$

(11)

$$K_l=K_r\left(r+\frac{r\times RA}{(1-RA)/3+\frac{r}{1-r}}\right)$$

(12)

$$r=\frac{K_{reg}}{K_r}$$

(13)

where $K_u$ and $K_l$ are the theoretical upper and lower bounds, normalized by $K_r$, where $K_{reg }$ is the thermal conductivity of regolith, $K_r=1.491 (\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1})$ is the thermal conductivity of rocks [50]. According to their results, if the two phases are incompletely randomly admixed, the mixed thermal conductivity depends on the degree of random admixing. The mixed thermal conductivity can be expressed as follows:

$$K_{mix}=K_u\times \beta +K_l\times \left(1-\beta \right)$$

(14)

where $\beta$ is the mixture weight of rocks in a mixed layer.

3.3. Loss Tangent Model

The loss tangent quantifies the transparency of material to electromagnetic radiation. Low loss tangents indicate very transparent media, while high loss tangents indicate relatively opaque material. Mathematically, the loss tangent is simply a shorthand combination of real, $\epsilon \mathrm{’}$, and imaginary, $\epsilon \mathrm{”}$, dielectric constants, where $tan\delta =\epsilon \mathrm{”}/\epsilon \mathrm{’}$. The real dielectric constant ($\epsilon \mathrm{’}$) controls the speed of the electromagnetic waves through the material, while the imaginary dielectric constant ($\epsilon \mathrm{”}$) dominates energy attenuation.

The loss tangent of materials has been discussed many times [20,33,47,52]. The dielectric loss tangent models of both the mare and highland lunar regolith employed here are derived from Siegler et al. [20]:

$${\overline{tan\delta }}_{mare}={10}^{[0.312\rho -2.65+Ti\times (f^{-0.0025}-0.958\left)\right]}$$

(15)

$${\overline{tan\delta }}_{highlands}={10}^{[0.312\rho -3.79+f^{-0.069}]}$$

(16)

where the value of $\mathrm{w}\mathrm{t}.\%\mathrm{T}\mathrm{i}{\mathrm{O}}_2$ content, $Ti$, is taken from Lucey et al. [39]. According to Siegler, Equation (15)’s fit is only valid at TiO₂ values over 2 wt.%; therefore, Equation (16) is used to calculate the loss tangent at the location with a TiO₂ smaller than <2 wt.%. To distinguish the regolith and rocks, the real parts of the permittivity for regolith and rocks are ${\epsilon }^{\prime }={1.871}^{\rho }$, ${\epsilon }^{\prime }={1.919}^{\rho }$, respectively. The dielectric loss tangent model of the rock used here is derived from [33]:

$${\overline{tan\delta }}_{rock}={10}^{0.312\rho -3.26+0.038\times \left(Ti+FeO\right)}$$

(17)

where the value of $\mathrm{w}\mathrm{t}.\%\mathrm{F}\mathrm{e}\mathrm{O}$ content, $FeO$ is also taken from Lucey et al. [39]. According to Hashin et al. [53], the dielectric loss tangent of the mixed layer can be modeled similarly with Equations (11)–(14), where the thermal conductivity is replaced with the corresponding loss tangent.

4. Method

In order to extract the subsurface RA from the CE-2 microwave TB data, the effect of subsurface RA on physical temperature and microwave TB has to be estimated. Before estimating these effects, we have to confirm the most satisfactory mixture weight $\beta$. Two regions with a lower and higher surface RA are selected to estimate the best-fitting $\beta$ value. Then, based on this value, the influence of subsurface RA on physical temperature and microwave TB is estimated. Finally, the subsurface RA can be obtained by minimizing the differences between average simulated TB and the observed TB of CE-2.

4.1. Determination of the Parameter $\beta$

To estimate the effect of parameter $\beta$ on the physical temperature, a vertical uniform RA distribution is assumed at the lunar center ($0°\mathrm{N},0°\mathrm{E}$) due to its relatively stable geological background compared to other regions. Here, we suppose a uniform value that is equal to the average surface RA of 0.0027 obtained from Diviner. Figure 1 shows the influence of the parameter $\beta$ on the modeled surface temperature. We can see that the nighttime surface temperature increases with increasing parameter $\beta$, and the surface temperature variation can reach up to 13 K at midnight (local time 0:00) when $\beta$ increases from 0 to 0.7. To obtain a satisfactory parameter $\beta$, the surface temperatures observed by a single channel ($T_7, 25-41 \mathsf{\mu }\mathrm{m}$) are selected to constrain due to their high signal-to-noise ratio [32]. We used observations of $T_7$ at two different regions to determine a general parameter $\beta$. The first region is the lunar center ($0\pm 0.1°\mathrm{N},0\pm 0.1°\mathrm{E}$), which has a relatively smaller average surface RA value (0.0027). The other region is the CE-3 landing site ($44.12\pm 0.1°\mathrm{N},340.5\pm 0.1°\mathrm{E}$), which has a relatively larger surface RA value (0.005).

The RMS method was employed to determine the best-fitting parameter $\beta$. We can calculate the RMS of surface temperature between the average observed and the simulated by changing parameter $\beta$. Figure 2 compares the nighttime surface temperature between the observed data and the simulated results. We can find that the modeled results matched very well with the observed data. Parameter $\beta$ is estimated to be small and contain the same number of $6.6\times {10}^{-5}$ at two different regions. That means that this value could be used to simulate the physical temperature at other regions of the Moon. Furthermore, our simulation shows that the variation in the TB at 37 GHz is far less than 1 K when the $\beta$ varies from 0 to 0.7. Therefore, we can take the value $\beta =6.6\times {10}^{-5}$ to model both the thermal conductivity and loss tangent of the mixed layer in our research.

4.2. Effects of Subsurface RA on the Surface Temperature and Microwave TB

Figure 3 shows the effect of the subsurface RA on the surface temperature and 37 GHz microwave TB. Figure 3a shows that the nighttime surface temperature is sensitive to the subsurface rocks, while the daytime temperature is not. In Figure 3a, we can see that the surface temperature increases with increasing subsurface RA, and the variation in nighttime surface temperature can reach up to 10 K at midnight when the subsurface RA increases from 0 to 0.2. However, the subsurface RA has a significant influence on both daytime and nighttime microwave TB. Figure 3b shows the sensitivity of microwave TB at 37 GHz to the subsurface rocks. We can see that the influence of the subsurface rocks on 37 GHz microwave TB changes with time. The nighttime TB value decreases with increasing subsurface RA and the daytime TB value shows the opposite pattern. In addition, the TB variation can reach up to 38 K at midnight when the subsurface RA increases from 0 to 0.2, which is larger than the surface temperature variation of 10 K. This means that the microwave TB is more sensitive to the subsurface rocks compared to the surface temperature.

4.3. Subsurface RA Inversion

Since the influence of the subsurface rocks on TB varies over time, we assume that the midnight effect is equal to the average effect. The microwave TB simulated by the surface RA at the lunar equator ($0°\mathrm{N}, 120°\mathrm{E}$) shows that the simulated midnight TB value is greater than the average observed TB of CE-2 37 GHz. The TB differences between the simulated and observed CE-2 at midnight can be up to 8 K. From Figure 3b, we can speculate that the RA value detected by CE-2 37 GHz microwave should be greater than that detected by Diviner infrared. Therefore, the average midnight TB observed by CE-2 37 GHz is selected to constrain the best-fit subsurface RA values. The RMS method is employed to estimate the most appropriate subsurface RA, which has the following forms:

$$\Delta TB=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left( T{B}_i^s\left(RA\right)-T{B}_i^o\right)}^2}$$

(18)

where $i$ indicates the $i$th channel, and N is the total channel number. Here, $N=1$ means that only the TB of the 37 GHz channel is chosen for comparison between theoretical $T{B}^s$ and CE-2 observation $T{B}^o$. The $T{B}^o$ data used here are the average normalized midnight microwave TB of CE-2. By calculating the RMS TB with changing subsurface $RA$ values, the best-fit $RA$ can be obtained until the RMS of $\Delta TB$ reaches the minimum. The best-fit subsurface $RA$ value at the lunar equator ($0°\mathrm{N},120°\mathrm{E}$) is 0.0676, and the RMS is $3.3\times {10}^{-4}$ K.

5. Results

The subsurface RA of the 16 large Copernican-aged craters was extracted from the average midnight microwave TBs of the CE-2 37 GHz. Figure 4 presents the locations of all the selected craters. The background mosaic map is the surface RA, retrieved from the Diviner.

Table 1. The characteristics of the selected craters.

Crater Names	Latitude (°N)	Longitude (°E)	Diameter (km)	Subsurface RA (Mean)	Surface RA (Mean)	RMS (K)
Aristarchus (A)	23.7	312.6	40	0.0929	0.032	0.001
Burg (B)	45	28.2	40	0.0556	0.0085	0.003
Byrgius A (BA)	−24.5	296.3	19	0.0636	0.0307	0.027
Copernicus (CP)	9.62	340	93	0.0747	0.0113	0.013
Crookes (CK)	−10.4	194.9	49	0.0668	0.018	0.002
Giordano Bruno (GB)	36	102.89	22	0.2636	0.0959	0.001
Glushko (G)	8.11	282.33	40	0.0835	0.0231	0.004
Jackson (J)	22.05	196.68	71	0.1119	0.0257	0.005
King (K)	4.96	120.49	76	0.1423	0.0151	0.025
Moore F (M)	37.4	185	24	0.1192	0.0327	0.001
Necho (N)	−5.25	123.24	30	0.1885	0.034	0.001
Ohm (O)	18.4	246.5	64	0.0398	0.0171	0.001
Sharonov (SH)	12.4	173.3	74	0.0274	0.0057	0.008
Stevinus (ST)	−32.5	54.2	75	0.0532	0.0067	0.007
Tycho (T)	−43.3	348.62	85	0.0757	0.0489	0.006
Vavilov (V)	−0.8	222.1	98	0.0446	0.0092	0.012

presents the latitude and longitude, diameter, mean surface RA retrieved from Diviner, the inverted average subsurface RA, and the mean RMS values of 16 craters. It can be seen from Table 1 that the maximum average RMS value between simulated and observed TB is 0.027 on the Byrgius A crater, which is small enough compared to the radiometric sensitivity of 0.5 K of CE-2 MRMs. Therefore, we think that the estimated results are reliable. From Table 1, we can also see that the average RA value of the subsurface layer of the crater is greater than the average RA value of the surface at the corresponding location. The Giordano Bruno (GB) crater has a maximum average subsurface RA of 0.2636, and the Sharonov (SH) crater has a minimum average subsurface RA of 0.0274; and the average RMS values at these two craters are 0.001 and 0.008, respectively.

Figure 5 shows the distribution of the subsurface RA of the selected craters. All of the images are the inverse distance weight (IDW) interpolation results produced by ArcMap. The background mosaic images are from the LRO wide-angle camera (WAC). The red color represents the higher subsurface RA, and the green color is the smaller subsurface RA. It can be seen from Figure 5 that the elevated subsurface RA is usually present over the crater’s rims, walls, and floors. In comparison with the distribution of surface RA and subsurface RA between the studied craters, we can see that the abundance distribution of surface rocks is not necessarily the same as that of subsurface rocks. For example, at the Giordano Bruno, Necho, and Aristarchus craters, there are similar RA distributions on both the surface and subsurface. However, there are obviously different distributions of RA at both the surface and subsurface at the Copernicus, Ohm, Sharonov, and Tycho craters. The surface RA at the central peak and the wall corner of Copernicus is larger than that at other places, but the subsurface RA at the south wall is larger than that at other places. The surface RA at the central peak, the north and east wall corner of Ohm is larger than that at other places, but the subsurface RA at the west wall is larger than that at other places. The surface RA at the northeast wall of Sharonov is larger than that at other places, but the subsurface RA at the south and west wall is larger than that at other places; the subsurface RA at the central peak and east wall of Tycho is larger than that at other places, but the subsurface RA at the north wall is larger than that at other places.

According to the analysis of the influence of rocks on the variation trends of surface temperature (Figure 3a) and microwave TB (Figure 3b), theoretically, the surface temperature and microwave TB at night show change conversely with the increase in subsurface RA. This is the reason that the abnormal “hot spots” [54] during a total eclipse on the lunar surface, revealed by infrared scanning images, turned out to be the “cold spots” on the nighttime microwave TB distribution observed by CE-1 and CE-2 MRM [11,55]. Upon inspection of the data distribution of Diviner infrared and CE-2 37 GHz TB of the studied craters, we found that when the “hot spots” in infrared coincide with the “cold spots” in the microwave, and the distribution of RA both on the surface and subsurface is consistent; otherwise, the distribution is inconsistent.

Figure 6 presents the relative distribution of nighttime regolith temperature obtained from the Diviner and nighttime microwave TB of CE-2 37 GHz at the craters where the RA distribution on the surface and subsurface is inconsistent. From Figure 6, we can see that the “hot spots” in the infrared are located at the central peak and wall corner of the Copernicus crater, at the north and east wall corners of the Ohm crater, at the northeast wall corner of the Sharonov crarer, and at the central peak and east wall corners of the Tycho crater. However, the “cold spots” in the microwave are located at the south, southwest, southwest, and north of the Copernicus, Ohm, Sharonov, and Tycho craters, respectively.

The surface RA detected by Diviner is defined by the areal percentage of rock slab in each pixel. According to Wei et al. [18], the sensitivity of the MRM data to the size of the rocks is much smaller than that of the Diviner infrared data. That is, the subsurface RA in our research not only includes the exposed surface rocks at a meter or larger scale detected by Diviner data but also the exposed surface rocks at centimeter scale and buried rocks at a centimeter or larger scale. This indicates that the differences between surface and subsurface RA are mostly caused by the different detection abilities of Diviner and CE MRMs. This conclusion is consistent with that of Ghent et al. [56], that the Diviner RA signature and radar signature are not necessarily consistent. For example, the Diviner RA signature correlates with the radar signature for the Giordano Bruno crater, while it is not correlated with the Copernicus crater.

From Figure 5, we can see that the subsurface RA detected by CE-2 37 GHz MRM data is usually clustered on one side of the impact craters. This distribution trend can also be found in Diviner RA for certain craters, including Giordano Bruno, Necho, and Jackson. In fact, there are two main sources of subsurface rocks detected by CE MRM data in the impact crater: the bedrock that is exposed or buried in situ, and the ejecta deposits sourced from other places. Therefore, this distinct asymmetry spatial distribution of subsurface RA can provide information about rocky ejecta and impact direction [56]. The ejecta reflected by subsurface RA suggest that the impact direction is from the northwest, southwest, and northeast for the Giordano Bruno, Necho, and Jackson craters, respectively. However, the Diviner RA is not consistent with this distribution trend in some craters. For example, the Diviner RA is almost evenly distributed around the floor bounds of Copernicus. In addition, the higher subsurface RA values are located in the south of the crater where the surface RA is relatively smaller. This discrepancy indicates that most of the rocks that lead to a decline in the low microwave TB values in the subsurface are covered by a sufficient depth of insulating regolith to hide them from Diviner’s view [8,56].

6. Discussion

Quantitative maps of subsurface RA at 16 large Copernican-aged craters were generated. Although the RMS between the theoretical and observed TB is small enough, there are still some parameters that affected the inversion results. For example, the measurement uncertainties, CE-2 MRM data, Diviner infrared data, TiO₂ and FeO contents, and Apollo sample’s thermal properties, constrain the inversion results. In this section, we will analyze the uncertainties generated by two factors: the H-parameter and the rock’s loss tangent.

6.1. Uncertainty Generated by H-Parameter

The average H-parameter of 0.06 m was used in our simulation, which was obtained from Hayne et al. [9]. According to their research, the H-parameter has a significant influence on the physical temperature. In addition, the H-parameter between nine Copernican-aged impact craters (King, Copernicus, Aristarchus, Jackson, Tycho, Necho, Byrigus A, Moore F, and Giordano Bruno) varies from 0.036 m to 0.063 m, as calculated by $H=H_0{t}_{Ma}^{0.1}$, with best-fit parameters $H_0=0.032$ m, where the crater age $t_{Ma}$ is given in millions of years. Here, we selected the King crater (5$°$N, 120.5$°$E) as a case study to discuss the uncertainty generated by the H-parameter. To highlight the effect of the H-parameter on microwave TB and subsurface RA, we calculated the microwave TB and subsurface RA within $\mathrm{H}=0.06\pm 0.024$ m.

Figure 7 shows the effect of three different H-parameters on a nighttime microwave TB of 37 GHz, and the inverted subsurface RA varies with an increase in H-parameter. It can be seen from Figure 7 that the midnight microwave TB variation values can reach up to 3 K when the H-parameter varies from 0.036 to 0.084 m, and the nighttime microwave TB increases with increasing H-parameter. In addition, the estimated subsurface RA also increases with increasing H-parameter. Moreover, we can find that when H-parameter varies from 0.036 to 0.084 m, and the inversion subsurface RA has an uncertainty of −0.02/ +0.01.

6.2. Uncertainty Generated by Loss Tangent of Rocks

Parameters such as the density, heat capacity, thermal conductivity, and loss tangent of rocks used in our study were obtained from the Apollo samples experiments. However, according to the spectral analysis of the lunar impact basins [57], the plagioclase and orthopyroxene abundance in some regions are different from those at the Apollo landing site regions. Although the relative dielectric permittivity was proven to be a function of density but not of chemistry, the loss tangent is a function of density and chemistry [33]. In addition, in Carrier III et al.’s [33] research, the loss tangent shows three theoretical fits with the data: the equivalent BHS formula, the two-dimensional regression to the power-law formula, and the three-dimensional regression of variance as a function of density and selected values of % ($Ti+FeO$). The loss tangent values calculated by the three fittings at the King crater (5$°$N, 120.5$°$E) are about 0.0035, 0.008, and 0.0224, respectively.

Figure 8 shows a comparison of the nighttime microwave TB of 37 GHz at three different loss tangent values and presents the effect of the rock’s loss tangent on inverted subsurface RA. It is worth noting that when we calculate the microwave TB for different rock loss tangents, the subsurface RA is set to 0.15. In Figure 8, the microwave TB variation values can reach up to −9 K at midnight when the rock loss tangent varies from 0.0035 to 0.0224, and the estimated subsurface RA decreases with increasing rock loss tangent. The uncertainty of the inversion subsurface RA generated by the different loss tangents of different rocks is about 0.04/−0.05.

7. Conclusions

Rocks are one of the basic materials that make up the Moon, which usually has more thermal inertia and a greater loss tangent than the lunar regolith. The rock TB model proposed by Hu et al. [17] was used to calculate the microwave TB of each vertical column. In this work, a new mixed model of two-phase thermal conductivity, consisting of rock and regolith, is used to replace the mixed thermal conductivity model of Hu et al. [17]. Then, by calculating the RMS TB with changing subsurface RA values, the best-fit subsurface RA can be obtained until the RMS of TB reaches the minimum. Finally, we inverted the subsurface RA on 16 large Copernican-aged craters of the Moon. The maximum average RMS value of these craters was about 0.027. This demonstrated that our inversion results are reliable.

In a comparison of the inversion subsurface RA values and surface RA values retrieved from Diviner, the results show that the average RA values in the subsurface are usually larger than those on the surface at the corresponding locations. Also, the spatial distribution of subsurface RA is not necessarily the same as that of surface RA. In addition, the subsurface rocks are usually clustered on one side of the impact craters. This distribution trend can also be found in the surface RA for certain craters, including Giordano Bruno, Necho, and Jackson. This distinct asymmetry spatial distribution of subsurface RA can provide information about rocky ejecta and impact direction [56]. However, in some craters, such as Tycho, the rock ejecta direction revealed by subsurface RA differs from that revealed by surface RA. This could indicate the different formation and evolution information of surface rocks and subsurface rocks; an interpretation will be presented in future work.