A Microphysical Thermal Model for the Lunar Regolith: Investigating the Latitudinal Dependence of Regolith Properties

The microphysical structure of the lunar regolith provides information on the geologic history of the Moon. We used remote sensing measurements of thermal emission and a thermophysical model to determine the microphysical properties of the lunar regolith. We expand upon previous investigations by developing a microphysical thermal model, which more directly simulates regolith properties, such as grain size and volume filling factor. The modeled temperatures are matched with surface temperatures measured by the Diviner Lunar Radiometer Experiment on board the Lunar Reconnaissance Orbiter. The maria and highlands are investigated separately and characterized in the model by a difference in albedo and grain density. We find similar regolith temperatures for both terrains, which can be well described by similar volume filling factor profiles and mean grain sizes obtained from returned Apollo samples. We also investigate a significantly lower thermal conductivity for highlands, which formally also gives a very good solution, but in a parameter range that is well outside the Apollo data. We then study the latitudinal dependence of regolith properties up to ±80° latitude. When assuming constant regolith properties, we find that a variation of the solar incidence‐dependent albedo can reduce the initially observed latitudinal gradient between model and Diviner measurements significantly. A better match between measurements and model can be achieved by a variation in intrinsic regolith properties with a decrease in bulk density with increasing latitude. We find that a variation in grain size alone cannot explain the Diviner measurements at higher latitudes.


Introduction
The surface of the Moon is covered in regolith (Langevin & Arnold, 1977;McKay et al., 1991), loose material formed through the weathering of the local rock by meteoroid bombardment, space weathering (Pieters & Noble, 2016) and thermal erosion (Delbo et al., 2014).The lunar surface is composed of two different terrains, the highlands and the maria.The maria are dark, geologically young and flat areas formed out of low-lying basins filled with basaltic lava from the interior after massive impacts had penetrated the lunar crust (Hiesinger & Jaumann, 2014).The highlands are older areas, more heavily cratered, mountainous and consist of anorthositic and plagioclase-rich magnesian-suite rocks (Hiesinger & Jaumann, 2014).
During the Apollo era, the properties of the lunar regolith were studied by several in situ experiments and laboratory measurements on returned samples.Core tube samples generally showed an increase in bulk density with increasing depth (Carrier et al., 1991) with drill core bulk densities ranging from 1,500 kg m 3 to 2,000 kg m 3 in the first 2 m (Carrier, 1974).Sieving of lunar regolith particles showed that the particles follow a log-normal grain size distribution with grain sizes ranging from a few micrometers to several millimeters (McKay et al., 1991).The mean grain sizes (particle diameter) of investigated samples collected at the Apollo landing sites range from 40 to 280 μm (Carrier, 1973;Heiken et al., 1973;McKay et al., 1974), with most means falling between 45 and 100 μm (McKay et al., 1991).In the Apollo heat flow experiments, lunar surface and subsurface temperatures were measured (Keihm & Langseth, 1973;Keihm et al., 1973), with mean surface temperatures reaching almost 400 K during the day and dropping below 100 K during the night near the lunar equator, while temperatures at greater depths (∼1 m) remained constant, suggesting the presence of an insulating surficial layer.The thermal conductivity of the regolith was measured in situ (Langseth et al., 1972(Langseth et al., , 1976) ) and in the laboratory (Cremers, 1972(Cremers, , 1975;;Cremers & Birkebak, 1971;Cremers & Hsia, 1973, 1974;Cremers et al., 1970) and found to be on the order of 10 2 W m 1 K 1 and 10 3 W m 1 K 1 , respectively.However, after the Apollo era, the global distribution of regolith properties was still not well understood.Besides in situ measurements, remote sensing measurements combined with models provide a powerful tool for characterizing lunar regolith in unsampled areas.In particular, radiometric measurements of the thermal emission of the (sub-)surface allow the microphysical structure (e.g., porosity and grain size distribution) and the thermophysical properties of the regolith to be determined by comparison with synthetic emission profiles derived from thermophysical models, since the subsurface temperature profiles and surface temperatures are sensitive to these properties.The lunar surface thermal emission has been measured by the Diviner Lunar Radiometer Experiment (Paige, Foote, et al., 2010) on board the Lunar Reconnaissance Orbiter (LRO) at high resolution.These measurements have been interpreted by a number of thermophysical models with different complexity levels to study the equatorial regolith properties (Vasavada et al., 2012), the global surficial thermal conductivity (Yu & Fa, 2016) and the global bulk density and thermal conductivity variation with depth (Hayne et al., 2017, their heat1d-model).The heat1d-model has been modified by Martinez and Siegler (2021) for application to the lowest temperatures at high latitudes and has been coupled to a 3D surface model including shadowing and scattering effects by King et al. (2020, their O3DTM model).Recently, Schorghofer et al. (2024) applied the heat1d-model in combination with a grain-size dependent expression for the radiative thermal conductivity (Gundlach & Blum, 2013;Sakatani et al., 2017) to explore the inference of regolith grain size.All these models are based on the description of the variation of bulk density and thermal conductivity with depth using the empirical H-parameter (Hayne et al., 2017;Vasavada et al., 2012).In general, the models used to describe the regolith properties have a critical influence on the modeled temperatures.
In this work, we expand upon the previous investigations, by developing a one-dimensional microphysical thermal model ("1DMTM"), which more directly simulates regolith properties, such as volume filling factor and grain size.The stratification is modeled as a function of regolith grain size and lithostatic pressure (Schräpler et al., 2015) and the thermal conductivity is a function of grain size, volume filling factor, and temperature (Gundlach & Blum, 2012).We compare the simulated surface temperatures with Diviner measurements and derive regolith grain size and volume filling factor as a function of latitude, with separate investigation of the lunar highlands and maria.
The thermal emission measured by Diviner is not only composed of contributions from the regolith, but also from rocks sitting on top of the regolith or from rocks embedded into it.Since rocks possess a higher thermal conductivity and bulk density and thus a higher thermal inertia than regolith, they heat up slower during the daytime and cool slower during the night.Hence, rocks stay much warmer during nighttime and increase the thermal emission.In this study, we use the regolith temperatures derived by Bandfield et al. (2011Bandfield et al. ( , 2017)), who separated the effect of surface rocks from regolith using the Diviner measurements.By comparing the Planck radiance as a function of wavelength for multiple Diviner channels and different surface rock area fractions within each Diviner pixel, they derived rock abundances and regolith temperatures for latitudes between ±80°and local times between 7:30 p.m. and 5:30 a.m.These restrictions arise because at higher latitudes and during the daytime unresolved local slopes appear as an anisothermality that is incorrectly interpreted by the algorithm as increased rock abundance.Another limit is the sensitivity to smaller rocks.The derived rock abundance corresponds to the fractional areal coverage of ∼0.5 m rocks and the technique is not sensitive to smaller rocks.Thus, the regolith temperatures still include the contribution of small rocks, and may also be affected by buried rocks within ∼10 cm of the surface (Elder et al., 2017).The regolith temperature data product used in this study is publicly available at the Geosciences Node of the Planetary Data System, covers the period from June 2009 to September 2016 and has the highest available resolution of 128 pixel per degree.
We further down-selected the Diviner nighttime regolith temperatures to achieve optimal data quality by applying filters based on (a) albedo, (b) rock abundance, and (c) local slope.First, the data were filtered by albedo to reduce the uncertainty due to albedo variation on the regolith temperatures by choosing only those pixels that are within (0.995, 1.005) × (1 A n ), with A n being the mean normal albedo measured by the Lunar Orbiter Laser Altimeter (LOLA; Lucey et al., 2014;Smith et al., 2010) on board LRO at zero phase angle and 1,064 nm wavelength.In addition, this filter allows to distinguish between highlands and maria, as both regions show a distinct difference in albedo with mean values of A n,M = 0.17 and A n,H = 0.31, respectively.Second, the data were filtered by rock abundance.Although this study already uses the Diviner regolith temperatures with the effects of rock removed, the technique used by Bandfield et al. (2011) involves a number of simplifications, including the modeling of rocks, and selecting pixels with low rock abundance reduces uncertainties.In addition, many small rocks are present at pixels with increased rock abundance, which cannot be detected by the technique, but may affect regolith temperatures (Byron et al., 2023).Here, we used the filter to select only pixels with a rock abundance ra ≤0.3%, which is below the global average rock frequency of 0.4% derived from Bandfield et al. (2011) and still gives us a sufficient number of pixels to evaluate.Third, a filter based on local slope based on the 128 pixel-perdegree gridded LROC global digital terrain model (Scholten et al., 2012) was applied to reduce the effect of resolved topography on the measurements with a selected maximum slope angle of α ≤ 5°.Finally, only the regolith temperatures at which all three filter criteria are met were selected.We find that generally, the maria have a higher rock abundance and smaller slope angles in comparison to the highlands.
The Diviner regolith temperatures were collected for every latitude ±0.5°at the pixels selected by the above described filter.Assuming that the regolith properties are the same in the northern and southern hemispheres, the data for negative and positive latitudes were combined for the maria and highlands pixels separately.Figure 1 shows the binned mean nighttime regolith temperatures as a function of latitude.The statistical error due to the standard error of the mean regolith temperatures is on the order of ∼10 2 K and therefore negligible (for more Figure 1.The filtered and binned mean Diviner regolith temperatures versus local time expressed in hours past noon as a function of latitude for the lunar maria (left) and highlands (right).The regolith temperatures were selected at the corresponding latitude (see legend) ±0.5°and are only available for nighttime and latitudes up to ±80°.The lunar maria do not cover latitudes higher than ∼50°.details see Appendix A).The regolith temperatures are found to be relatively similar (i.e., within ±0.5 K) in the highlands and maria.

Heat Transfer Model
The evolution of temperature T with depth x and time t is described by the one-dimensional heat transfer equation, which takes as input the physical and thermophysical properties of the modeled material, namely the depthdependent bulk density ρ, the temperature-dependent specific heat capacity c p and the thermal conductivity λ that may depend on depth and temperature.The models used in this study to describe these properties are presented in Sections 3.2, 3.3, and 3.4.To solve the heat transfer equation numerically, a forward time, centered space finite difference approximation was applied and the simulated subsurface was divided into several layers of thickness Δx, where the change in temperature was calculated for each time step Δt.The numerical solution is presented in more detail in Appendix B1.
At the surface, the change in temperature depends on the absorbed solar radiation, the heat transported in or out of the surface layer from the subsurface and the thermal emission to space.Thus, the upper boundary condition reads with the bolometric Bond albedo A θ (the index indicates that its value is a function of the solar incidence angle θ), the incoming solar flux F ⊙ , the infrared emissivity ϵ, the Stefan-Boltzmann constant σ and the surface temperature T 0 .The bolometric Bond albedo A θ depends on the solar incidence angle θ and is described in more detail in Section 3.5.The incoming solar flux equals zero during nighttime and is during the day a function of the solar constant I E , the heliocentric distance of the simulated body in AU r h , its rotation period T period and the considered latitude β.The rotation period of the Moon corresponds to one synodic month (∼29.5306days) and its distance to the Sun is approximated to be 1 AU at all times.We neglect here both the eccentricity of Earth's and Moon's orbit and the inclination of the Moon's equator against the ecliptic (∼1.5°).The values of all model parameters are summarized in Appendix B2 in Table B1.
Details about the grid implementation, equilibration times and numerical stability are presented in Appendix B1.

Bulk Density and Volume Filling Factor Model
The bulk density and volume filling factor are depth-dependent and described by the stratification model for planetary regolith of Schräpler et al. (2015).This empirical model builds upon the pressure dependence of the volume filling factor of cohesive-particle layers found by Güttler et al. (2009), with the volume filling factor ϕ, the pressure p, the minimum volume filling factor at the surface ϕ s and the maximum volume filling factor at deeper layers ϕ d , the turnover pressure p m and the logarithmic transition width from low to high volume filling factor Δ. The turnover pressure p m (r) is a function of grain radius r p m (r) ∝ r 4 3 (5) and is the formal inflection point of the ϕ-(log p)-curve because The basic assumption of the stratification model in Schräpler et al. (2015) is that the pressure in the regolith is determined by the lithostatic stress and therefore the regolith is compacted by the weight of the overlying regolith layers.Thus, by calculating the lithostatic pressure as a function of depth, p(x), and inserting this relation into Equation 4, the volume filling factor is obtained as a function of depth ϕ(x).However, only the inverse, the depth as a function of the volume filling factor x(ϕ), can be calculated and reads Δln( 10) with the gravitational acceleration g, the material density ρ grain of a single regolith grain and the Gaussian hypergeometric function H2F1 (Schräpler et al., 2015).In practice, the volume filling factor at the desired depths is determined by interpolating the result of Equation 7.
We note that Equations 17 and 18 in Schräpler et al. (2015) were not published in the correct form, so we provide here the correct expressions (Equations 7 and 8, respectively).The pressure as a function of volume filling factor p (ϕ) reads in the correct form An alternative to Equation 4 was recently provided by Tatsuuma et al. (2023) who derived A comparison between Equations 8 and 9 shows that for ϕ ≫ ϕ s their functional behavior for ϕ → ϕ d is very similar.However, the exponent in Equation 8, Δln(10) = 1.3 for Δ = 0.58 (Güttler et al., 2009;Schräpler et al., 2015), is smaller than the exponent 3 in Equation 9 so that we treat Δ as a free parameter.
Finally, the regolith bulk density and volume filling factor are related by the material density, which is evaluated separately for the highland and maria regions from measurements of Kiefer et al. (2012).The bulk density at the surface ρ s is determined from the respective volume filling factor.At the surface, a fixed value for the volume filling factor of ϕ s = 0.15 is chosen, which on the one hand represents random ballistic deposition (Blum & Schräpler, 2004;Jullien & Meakin, 1988), and on the other hand was empirically derived by Szabo et al. (2022) from modeling energetic neutral atom emission of the upper lunar regolith layers.We calculated the mean free path of thermal emission and absorption within a packing with such a low volume filling factor of ϕ s = 0.15 and for regolith grains with r = 50 μm considering the number density and effective cross-section to l ∼ 1 mm.This means that the location of absorption and emission has a certain uncertainty of around 1 mm.However, in order to view and resolve this more precisely, the model would have to be calculated in 3D.The volume filling factor and bulk density in the deeper layers are treated as a free parameter in the simulations.The Apollo drill cores showed a variation in bulk density not only with the landing site but also within one landing site, which can be explained by differences in grain density, grain shape, and grain size distribution (Carrier et al., 1991).We note, that high stresses and/or dynamical processes can lead to compaction so that values of ϕ d → 1 are also feasible, given that the particle-size distribution is sufficiently wide.
Figure 2 shows the resulting stratification profiles for two different regolith grain radii and three logarithmic transition widths, but fixed deep layer volume filling factor.The grain radius uniquely determines the depth of the turning point (Equation 6) and thus, how easily a regolith layer can be compressed under its own weight.According to Equation 5, the depth at which the transition from low to high volume filling factors occurs decreases with increasing grain radius.On the other hand, the transition width Δ determines the steepness of the transition from the loose packing at the surface to the maximum volume filling factor in the deeper layers.The steepness decreases with increasing transition width Δ.We note that for small grain radii and especially for large transition widths Δ much larger depths are needed to reach the maximum volume filling factor and bulk density and therefore the corresponding value at ∼1-2 m depth does not always equal the maximum volume filling factor ϕ d and bulk density ρ d .

Thermal Conductivity Model
There are three modes of heat transfer in a granular medium: transport through the solid network of grains, radiation inside the void spaces of the packing, and gas diffusion.In the case of lunar regolith, the latter is not important, due to the low gas density on the Moon.Gundlach and Blum (2012) describe the heat conductivity λ of a granular packing of equal-sized spheres in vacuum as a function of particle radius r, volume filling factor ϕ and temperature T through (11) All parameter values can be found in Appendix B2 in Table B1.We note that in their original paper, Gundlach and Blum (2012) assumed that the surface free energy γ was linearly temperature dependent.However, newer research suggests that this assumption was wrong.Pinon et al. (2016) experimentally and theoretically showed that in the temperature range relevant for the lunar surface (∼100 400 K), the adhesion force in Equation 8 in Gundlach and Blum (2012) decreases by about 20%.As the surface free energy enters our Equation 11 only as the cubic root, the reduction effect shrinks to about 7%.Pinon et al. (2016) also state that for contact areas greater than 10 4 nm 2 , the temperature dependency basically vanishes.Thus, we disregard any temperature dependence of γ in this work.
The radiative term in Equation 11 is strongly temperature-dependent and is a function of grain radius and volume filling factor, which determine the mean free path of the photons.In a granular medium, the thermal conductivity is decreased in comparison to the solid material (with a conductivity of λ solid ), due to the reduced number of contacts in the packing, the small contact area, the packing structure and an empirical scaling factor χ. The latter was introduced by Gundlach and Blum (2013) to fit the thermal conductivity model to laboratory measurements of lunar regolith samples and to account for the irregular shape and non-monodisperse size-distribution of the grains.
In this work, we re-perform this calibration by introducing, first, a temperature dependence of the thermal conductivity of the solid material λ solid and, second, a separate fit for material from the highlands and the maria.
The introduction of a temperature dependence of λ solid is important, because the thermal conductivity should approach zero as the temperature approaches 0 K (Woods-Robinson et al., 2019) and during nighttime, regolith temperatures below 70 K at 80°latitude (see Figure 1) are observed.The details of the derived temperature dependence can be found in Appendix B3.
Following the approach in Gundlach and Blum (2013), the thermal conductivity model is then fitted to laboratory thermal conductivity measurements of returned Apollo regolith samples.These thermal conductivity measurements as a function of temperature were conducted for returned samples from each Apollo mission, from Apollo 11 to Apollo 16 (Cremers, 1972(Cremers, , 1975;;Cremers & Birkebak, 1971;Cremers & Hsia, 1973, 1974;Cremers et al., 1970).However, clear association with the maria or highland regions can only be assigned to the Apollo 11 (maria), Apollo 12 (maria) and Apollo 16 (highlands) landing sites.The Apollo 14 landing site is located at large impact ejecta ridges containing highland debris and maria fragments, while the Apollo 15 landing site is located at a highland-mare boundary (Vaniman et al., 1991).For the calibration of the model, the volume filling factor and grain radius of the measured regolith samples are needed.Cremers et al. (1970), Cremers and Birkebak (1971), and Cremers and Hsia (1974) provide information only on the bulk density of the measured sample, which is converted to a volume filling factor via the mean grain densities for maria and highland material derived from Kiefer et al. (2012).Information on the grain size is taken from the Lunar Sample Compendium (Meyer, 2005).Details on the exact values and the fitting procedure are listed in Appendix B3.The variation of the scaling parameter χ leads to best-fit values of χ = 0.69 for samples from the maria (Apollo 11 and 12), hereafter "maria model," and χ = 0.14 for the highland sample (Apollo 16), hereafter "low conductivity highlands model." Figure 3 illustrates the model fits along with the laboratory thermal conductivity measurements.We note already at this point that the thermal conductivity of the highland sample is much lower than the measured values for the samples from the maria.This difference is even more pronounced due to the different bulk densities of the samples.The data from Apollo 11 and 12 were measured at a bulk density of 1,300 kg m 3 , while the sample from Apollo 16 was measured at a bulk density of 1,500 kg m 3 .Therefore, we added another scenario, the "high conductivity highlands model," in which the highlands are given the same thermal conductivity as the maria with χ = 0.69.As a result, the "maria model" and "high conductivity highlands model" assume the same thermal conductivity, while the "low conductivity highlands model" assumes a lower thermal conductivity for the highlands.Making this case distinction allows us to compare the effects of a distinct difference in thermal conductivity between maria and highlands with a globally uniform thermal conductivity.(Cremers & Birkebak, 1971;Cremers & Hsia, 1974;Cremers et al., 1970).The Apollo 11 and 12 data were measured at a bulk density of 1,300 kg m 3 and are associated with material from maria, while the Apollo 16 sample was measured at a bulk density of 1,500 kg m 3 and is associated with highland material.The corresponding maria model fit and low conductivity highlands model fit are illustrated.In this study, we add a third scenario, the high conductivity highlands model, where the highlands are given the same thermal conductivity as the maria with χ = 0.69.

Specific Heat Capacity Model
The specific heat capacity is temperature-dependent and fitted to measurements of regolith samples from the Apollo missions (Hemingway et al., 1973;Robie et al., 1970).It is modeled as a fourth-order polynomial with coefficients taken from Hayne et al. (2017), which are listed in Table B1.The regolith in the maria and the highlands is described by the same fit as the measurements vary only insignificantly.

Albedo Model
The bolometric Bond albedo A θ is a function of solar incidence angle θ and describes the ratio of reflected solar radiation to incident solar radiation.It was derived by Feng et al. (2020) from Diviner noontime surface temperatures as a function of latitude.
A θ = A 0 + 0.6126 1.597 cosθ + 2.15 cos 2 θ 1.636 cos 3 θ + 0.4704 cos 4 θ) (13) .The difference between the maria model and the high conductivity highlands model is due to an assigned difference in albedo and grain density, whereas the low conductivity highlands model is assigned an additional difference in thermal conductivity.The temperature profile as well as the specific heat capacity and thermal conductivity profile are illustrated for two different local times, noon and the time just before sunrise.It is important to note two aspects of thermal conductivity.First, it is a function of volume filling factor (see Section 3.3) and therefore the maria and high conductivity model only show similar thermal conductivity profiles for similar volume filling factor profiles.Second, the thermal conductivity is an interface property and therefore evaluated at the layer boundaries (see Appendix B1).
The bolometric Bond albedo at zero solar incidence angle was found to be A 0,highlands = 0.12 for the highlands and A 0,maria = 0.07 for the maria.

Highlands and Maria
In this study, the regolith in the highlands and maria is characterized by a difference in grain density (see Section 3.2) and albedo (see Section 3.5) (maria and high conductivity highlands model) as well as a difference in thermal conductivity (see Section 3.3) (low conductivity highlands model).The result is an initial difference, without any variation of the free parameters, in the modeled stratification, thermophysical properties and temperature, which is illustrated in Figure 4.The difference in grain density results in a difference in bulk density at fixed volume filling factor (lower panels) and vice versa (upper panels).We note that the volume filling factor profiles for maria and highlands for fixed surface and deep layer volume filling factor are not completely identical, because the volume filling factor as a function of depth also depends on grain density (see Equation 7).

Parameter Sensitivity
In Appendix B4 we investigate the sensitivity of the three free parameters in the thermal model-the regolith grain radius r, the logarithmic transition width Δ and the deep layer density ρ d -to the modeled nighttime temperatures along with other model parameters.

Derivation of Grain Size and Volume Filling Factor for Highlands and Maria at the Lunar Equator
We start at the lunar equator and find the best fitting parameter set by determining the minimum root-mean-square deviation (RMSD), with the measured regolith temperature T n , modeled regolith temperature Tn and number of compared temperature pairs N. Figure 5 shows the filtered Diviner regolith temperatures (see Section 2 for details) at the lunar equator for maria and highlands together with the modeled regolith temperatures that result in the minimum RMSD.The minimum RMSD was found for the following parameter combinations: Maria r = 30 μm, Δ = 0.4, ρ d = 2,200 kg m 3 with RMSD = 0.31 K, high conductivity highlands r = 40 μm, Δ = 0.1, ρ d = 1,700 kg m 3 with RMSD = 0.17 K and low conductivity highlands r = 70 μm, Δ = 0.2, ρ d = 2,500 kg m 3 with RMSD = 0.21 K.For comparison, as a measure of the lowest possible RMSD, that is, the intrinsic noise of the averaged Diviner regolith temperatures, a third-order polynomial was fitted to the Diviner regolith temperatures, yielding RMSD M = 0.25 K and RMSD H = 0.17 K for maria and highlands, respectively.The lower panels of Figure 5 show the temperature difference between the modeled and measured regolith temperatures, which is very similar to the temperature difference between the polynomial fit and the measurements.Therefore, the developed microphysical model is able to describe the equatorial nighttime regolith temperatures with an RMSD very close or similar to the lowest possible value.
However, as discussed in Appendix B4.2, there is a degeneracy among the three free parameters: the grain radius, the transition width, and the deep layer density.Numerically, the degeneracy can be determined by the rank of the Jacobian matrix for the least squares fit, and it indicates that the available data does not constrain the parameters alone.To characterize these degenerate solutions, we created heat maps, which are presented in Figure 6, to visualize the resulting RMSD for various combinations of these parameters.For the maria and high conductivity highlands model, the transition width Δ tends to increase as the deep layer density increases, while the grain radius slightly decreases.This means that, for similar surface temperatures, the transition from surface to deep layer density becomes more gradual at higher deep layer densities.Conversely, at lower deep layer densities, the transition becomes sharper, and the point at which this transition occurs (determined by the grain radius) is found at shallower depths.In general, it is difficult to describe the degeneracy explicitly, since the grain radius directly affects not only the stratification but also the thermal conductivity, and furthermore, the thermal conductivity is also a function of stratification (volume filling factor).For the low conductivity highlands model, very high bulk densities and a sharp transition in the stratification profile (indicated by a very small Δ) are preferred.The degenerate solutions are limited by an upper constraint imposed by the deep layer density, which is influenced by the bulk density of the highland crust (Han et al., 2014;Wieczorek et al., 2013).
In the following subsections, we will discuss the derived grain size values and their meaning as well as the volume filling factors and bulk densities suggested by the model.To better constrain the non-unique solutions, we compare the derived values with grain sizes and bulk densities inferred from Apollo samples.In addition, we will discuss the applied difference in thermal conductivity between the low conductivity and high conductivity highlands (and maria) and the resulting effect on the derived regolith properties.

Derived Bulk Densities and Volume Filling Factors
All thermophysical model solutions for the lunar equator are illustrated in Figure 6 and the resulting bulk densities and volume filling factors will be discussed in the following.For the maria, deep layer bulk densities in the range ρ d = 1,800-2,500 kg m 3 lead to a RMSD <1 K between model and measurement, which translates into a deep layer volume filling factor of ϕ d = 0.54-0.74.In the case of the high conductivity highlands, solutions with RMSD <1 K can be found for deep layer bulk densities in the range ρ d = 1,700-2,300 kg m 3 and at ρ d = 2,500 kg m 3 , which translates into a volume filling factor of ϕ d = 0.59-0.87.In contrast, the low conductivity highlands case suggests only very high deep layer bulk densities in the range ρ d = 2,400-2,500 kg m 3 , leading to a deep layer volume filling factor of approximately ϕ d = 0.83-0.87.The difference in derived bulk density between low and high conductivity highlands is a consequence of their difference in thermal conductivity.
In the case of the lower thermal conductivity, similar surface temperatures are achieved with higher bulk densities.This will be discussed in more detail in Section 4.1.3.
A lower limit of the deep layer volume filling factor is given by ϕ = 0.56 for random loose packing (Onoda & Liniger, 1990).As an upper limit, random close packing of monodisperse spheres gives ϕ = 0.64, but the maximum volume filling factor increases with increasing polydispersity.For log-normal size distributions, spheres can pack to a maximum filling factor of ϕ = 0.87 (maximum ϕ d value from above) when they possess a logarithmic width of ∼1 dex (Anzivino et al., 2023), which is in agreement with the measured size distribution   , 1991).For the highlands, we also add the mean grain radius of 40-140 μm from Apollo 16 samples with 40 μm as lower limit (dashed line), as determined by Heiken et al. (1973).As a constraint on the derived bulk density, we draw boundaries between ρ bulk = 1,500 kg m 3 and ρ bulk = 2,000 kg m 3 (center column), as obtained from the Apollo drill core experiments (Carrier, 1974).For illustration purposes, since the bulk density is a function of the grain radius, the transition width, the deep layer density and the depth below the surface, the bulk density is evaluated at 1 m depth and for the upper and lower limit of the grain radius, namely 20 μm (solid line) and 50 μm (dotted line).
( McKay et al., 1991).However, such a jammed packing can only be achieved under compression, vibration or thermal cycling so that the constant meteorite bombardment and the day-night temperature variations have to be considered.
The returned Apollo regolith drill cores, in contrast, showed a maximum regolith bulk density of ≈2,000 kg m 3 (Carrier, 1974).We used the bulk density inferred from the drill cores as a constraint on the derived bulk density and drew boundaries between ρ bulk = 1,500 kg m 3 and ρ bulk = 2,000 kg m 3 in Figure 6 (center column).Since the bulk density model is a function of the grain radius, the transition width, the deep layer density and the depth below the surface, we chose to evaluate the bulk density at 1 m depth and for an upper and lower limit of the grain radius, namely 20 and 50 μm (see Section 4.1.2for more details).For the high conductivity highlands, the initial best-fit solution is within these boundaries and only some solutions are excluded.However, the upper constraint on the bulk density eliminates all solutions for the low conductivity highlands model and the majority of the solutions for the maria model.The best solution for the maria would then be r = 40 μm, Δ = 0.3, ρ d = 2,000 kg m 3 with RMSD = 0.42 K.This deep layer density translates into a deep layer volume filling factor of ϕ d = 0.60, which is very close to the deep layer volume filling factor for the best-fit of the high conductivity highlands model ϕ d = 0.59.In fact, assuming that the grain size-frequency distribution is fairly similar in the maria and highlands, the regolith should have similar packing in both regions and the limit of the densest packing at great depth should be the same.The predicted difference in bulk density between maria and high conductivity highlands is a result of the difference in assigned grain density (see Section 3.2).This difference would be less pronounced for more similar grain densities.Recent data from regolith samples from the Chang'E 5 mission indicate ρ grain ≈ 3,200 kg m 3 for material from the maria (Li et al., 2022), which is slightly less than the value of ρ grain ≈ 3,360 kg m 3 used in this work.
Figure 7 shows the derived volume filling factors and bulk densities as a function of depth below the surface for maria, high conductivity and low conductivity highlands and their three best fits (no constraints applied).The profiles are compared with the mean bootprint depth and best estimates of regolith bulk density derived by Mitchell et al. (1974).From the analysis of more than 700 astronaut footprints, the average bootprint depth in intercrater areas is ∼0.7 cm (Carrier et al., 1991;Houston et al., 1972;Mitchell et al., 1973).This means that the regolith can only have reached its maximum compaction below this depth and this gives us another constraint on the derived stratification profile.Furthermore, Mitchell et al. (1974) derived best estimate values for the bulk density from the Apollo core tube samples in distinct depth intervals and suggested a bulk density of 1,580 ± 50 kg m 3 for the depth interval 0 30 cm and a bulk density of 1,740 ± 50 kg m 3 for the depth interval 30 60 cm.The profiles for maria and high conductivity highlands show similar volume filling factors, but different bulk densities due to their assigned difference in grain density (see Section 3.2).As a result, only the high conductivity highlands show a reasonable agreement with the best estimate values of bulk density from Mitchell et al. (1974), whereas the maria show higher bulk densities.

Derived Grain Radius
The developed microphysical thermal model describes the regolith stratification and thermal conductivity in terms of grain radius.Although the thermophysical model solutions are degenerate (see Figure 6), we find that we can explain the Diviner regolith temperatures at the lunar equator well with mean grain radius values obtained from returned Apollo samples, which range from 20 to 140 μm (Carrier, 1973;Heiken et al., 1973;McKay et al., 1974), with most means falling between ∼20 and 50 μm (McKay et al., 1991).For the maria and a RMSD <1 K, we find a best-fitting grain radius of r = 20-70 μm and for the low conductivity highlands r = 50-100 μm.
In the case of the high conductivity highlands, we find a solution for the entire range of investigated values r = 10-100 μm, but the best three solutions are between 20 and 50 μm, which corresponds to the range in which most of the mean grain radii of the Apollo samples fall.To better constrain the non-unique solutions, this range of the most likely grain radius from returned Apollo samples is illustrated in Figure 6.For the highlands, we also add the mean grain radius of 40-140 μm determined by Heiken et al. (1973) from Apollo 16 samples with 40 μm as the lower limit.Carrier (1973) lists a mean grain radius for the Apollo 11 and 12 samples (associated with maria) of r ≈ 20-50 μm, which equals the range in which most of the mean values of the Apollo samples fall.However, accurate statistical comparisons between maria and highlands are difficult due to different experimental techniques.When analyzing the grain radius of the Apollo 11, 12, 14 and 15 samples, Carrier (1973) found little statistical difference and argues that the particle size distribution for most of the lunar surface has reached a steady state in which comminution is balanced by agglutination and particle size is therefore relatively insensitive to exposure age.
The non-unique solutions can be further constrained by the bulk density from Apollo drill core experiments, which is discussed in more detail in Section 4.1.1.Since the bulk density model is a function of all three parameters, the grain radius, the transition width and the deep layer density, it is not straightforward to draw exact boundaries in Figure 6, but roughly speaking, the upper limit of ρ bulk = 2,000 kg m 3 excludes solutions with RMSD <1 K at the smallest grain size values of r = 20 μm for maria and r = 10-20 μm for the high conductivity highlands as well as all solutions for the low conductivity highlands.This narrows down the best-fitting grain radius for all solutions with RMSD <1 K to r = 30-70 μm for the maria and to r = 30-100 μm for the high conductivity highlands.In both cases, the best fit is achieved with a grain radius of r = 40 μm.
The recent work of Schorghofer et al. (2024) analyzed a small region in the maria at the lunar equator and derived a most likely grain radius of 90 μm, which is slightly above the range determined in this work for equatorial maria regions.However, it should be noted that their approach only considers the dependence of the radiative thermal conductivity on the grain size, while this work additionally employs a network conductivity and stratification dependent on the grain size.B4) at 0°latitude and higher latitudes (50°maria, 80°highlands).The best fit shown in Figure 5 is highlighted with a dashed line.The maria and high conductivity highlands model show similar volume filling factor profiles, but a difference in bulk density due to their different grain density.The bulk densities of the Apollo core tubes (Mitchell et al., 1974) are shown by the vertical bars and the average bootprint depth (Carrier et al., 1991;Houston et al., 1972;Mitchell et al., 1973) is marked by the hatched box at the top panels.
In reality the regolith particles follow a log-normal grain size distribution with grain sizes ranging from a few micrometers to several millimeters (McKay et al., 1991) and thus, the assumption of a single mean grain radius presents a simplification.The smaller grains can fill the voids between larger grains modifying the volume filling factor, the mean free path and the number of contacts between the grains.In this study, the network and radiative thermal conductivity as well as the stratification is a function of grain size, which does not allow a direct physical interpretation of the derived mean grain radius.Ryan et al. (2020) describe the radiative thermal conductivity of polydisperse particles with a mean grain size corresponding to the Sauter diameter, which represents the surface area-weighted mean.Sakatani et al. (2018) model the thermal conductivity of a polydisperse granular packing as monodisperse grains identified as the volumetric median grain size.Another study by Xiao et al. (2022) identifies the derived mean grain size as the number-weighted grain size.

Caveats: Applied Thermal Conductivity Values
In the following, we will discuss the thermal conductivity values applied in this work and the hypotheses of a difference in thermal conductivity between maria and highlands, which was realized through the low conductivity highlands model.In general, the thermal conductivity is a critical input parameter for thermal models.Laboratory data of returned lunar samples and regolith simulants suggest low thermal conductivities (λ < 0.005 W m 1 K 1 for T < 400 K) and many theoretical thermal conductivity models are calibrated to these laboratory data (e.g., Reiss, 2018;Sakatani et al., 2018;Vasavada et al., 2012;Woods-Robinson et al., 2019, this work).The heat1d model (Hayne et al., 2017) is not calibrated to the lab data and is based on a different density dependence for thermal conductivity and thus shows significantly higher thermal conductivities for higher densities (λ > 0.005 W m 1 K 1 for T < 400 K).In general, there is a discrepancy between thermal conductivity measured in the laboratory and in situ, with in situ measurements being an order of magnitude higher (Langseth et al., 1972(Langseth et al., , 1976)).The incorporation and discussion of alternative thermal conductivity values, for example, the measured in situ thermal conductivity data will be part of a future work.
The exact dependence of thermal conductivity on the volume filling factor (bulk density) and temperature as well as the calibration to absolute values critically influences the derived regolith properties.In this work, the maria and high conductivity highlands model are assigned the same thermal conductivity as a function of volume filling factor and the derived regolith properties suggest a stratification with similar volume filling factors (see Section 4.1.1 and Figure 7 for more details).In contrast, the low conductivity highlands model, which is based on the thermal conductivity measurement of the Apollo 16 sample (Cremers & Hsia, 1974), which is at least a factor of 2-3 lower than that of the maria samples (see Figure 3), leads to very high and possibly unrealistic bulk densities.In this case, the lower thermal conductivity is thus compensated for by an increase in bulk density, in particular by increasing the deep layer density and decreasing the transition width, but also to some extent by increasing the regolith grain radius.
In the following, the hypotheses of a lower thermal conductivity for the highlands in contrast to a uniform conductivity for maria and highlands is discussed in more detail.First, we would like to motivate additional thermal conductivity measurements for returned samples from the lunar highlands to verify the measurement by Cremers and Hsia (1974).Qualitatively, the difference in thermal conductivity is consistent with differences in the composition of the two terrains, as the maria are regions with an elevated abundance of FeO and TiO 2 , which have high thermal conductivity.The main components of both terrains are basalt and anorthosite, respectively, and the thermal conductivity of anorthosite (≈1.7 W m 1 K 1 at T ≈ 275 K; Birch & Clark, 1940) is lower than the one of basalt (≈2.6 W m 1 K 1 at T ≈ 275 K; Halbert & Parnell, 2022) (≈2.5 W m 1 K 1 at T ≈ 275 K; Clauser & Huenges, 1995).However, Fujii and Osako (1973) measured the thermal conductivity in vacuum over a wide temperature range T = 100-600 K of an Apollo 11 and 14 rock with similar porosity, but difference in material density and Ti and Fe concentration, and found a fairly similar thermal conductivity.Another property influencing thermal conductivity of granular materials is particle shape, with more irregular or rough particles leading to a lower network thermal conductivity (Gundlach & Blum, 2013).There are few studies investigating the shape of regolith particles using X-ray computed tomography (e.g., Baidya et al., 2022;Katagiri et al., 2010) and the number of particles investigated is quite low and comparisons are difficult.Finally, we would like to mention the discovery of a cold spot at the Apollo 16 landing site by Powell, Williams, et al. (2023), close to where the sample of the highland thermal conductivity measurement was taken from.Cold spots show reduced nighttime temperatures and can be explained by a decompaction of the upper regolith layer.
In summary, the available data suggest a lower thermal conductivity for the main component of the highlands material, anorthosite, in comparison to basalt, the main component of the maria.However, a factor of 2-3 or even more between highland and maria thermal conductivity seems unrealistic and additional and comparable measurements are needed to investigate this further.If the difference in thermal conductivity were real, this would lead to much higher bulk densities in the highlands than suggested by the Apollo drill core experiments, and to a significant difference between maria and highlands in terms of volume filling factor as well as bulk density.We therefore conclude that at least such a significant difference in thermal conductivity between highlands and maria is very unlikely.

Latitudinal Dependence of Lunar Regolith Properties
We investigate the latitudinal dependence of intrinsic regolith properties by fitting the thermophysical model to Diviner regolith temperatures as a function of latitude up to ±80°(see Section 2 for details on the data set).Initially, a fixed parameter set is assumed, using the best-fit regolith properties determined at the lunar equator for maria, high conductivity highlands and low conductivity highlands, respectively (see Section 4.1).Although the solutions are non-unique at the lunar equator (Figure 6), the observed trend with latitude is the same for the different possible parameter combinations.The modeled temperatures assuming constant regolith properties compared with the Diviner data set are shown in Figure 8.We find that the Diviner regolith temperatures can be well described for low latitudes, but significant systematic deviations arise at a latitude >30°and increase with increasing latitude.At the highest latitudes, the modeled regolith temperatures are up to ∼3.5 K higher than the measured ones.We note that using Diviner nighttime bolometric temperatures, which represent the wavelengthintegrated flux from all Diviner channels (Paige, Siegler, et al., 2010), cannot eliminate or even reduce the observed latitudinal gradient.Figure 6 in Powell, Horvath, et al. (2023) compares the median regolith and bolomeric temperatures at midnight as a function of latitude and shows a nearly constant offset (slightly decreasing with increasing latitude) with the bolometric temperatures being ∼1 K higher than the regolith temperatures.In order to reduce the latitudinal gradient, the difference between nighttime regolith temperatures and bolometric temperatures would need to increase with increasing latitude.
Earlier investigations involving Diviner data and thermal modeling showed contrasting results.The H-parameter maps created by Hayne et al. (2017) did not show a latitudinal gradient.In contrast, Yu and Fa (2016) found a strong latitudinal dependence of the derived surficial thermal conductivity, which they attributed to a variation of intrinsic regolith properties.First, in Section 4.2.1, we will discuss processes that may contribute to regolith formation that can be associated with a latitudinal dependence and whether we expect this to change regolith properties.Second, in Section 4.2.2 we will investigate whether the temperature deviations between model and observations at higher latitudes shown in Figure 8 can be removed by applying an alternative dependence of the albedo on solar incidence angle.This is motivated by the fact that one difference between the previously presented works is the albedo model used in each case.Finally, in Section 4.2.3 we use the microphysical thermal model to predict and discuss a variation of intrinsic regolith properties with latitude.

Latitudinal Dependence of Regolith Formation and Evolution Processes
One of the processes involved in regolith formation is space weathering, due to the effects of solar wind and micrometeoroid impacts.Space weathering effects are reduced toward the lunar poles, because the maximum flux of solar wind particles and micrometeorites occurs near the lunar equator and decreases strongly with increasing latitude (Cremonese et al., 2013;Kallio et al., 2019).However, larger impactors show the opposite effect.The relative flux of near-earth objects increases from the lunar equator to the poles by a factor ∼1.6 (Robertson et al., 2021).This is in qualitative agreement with crater counts on the Moon (Robertson et al., 2021).However, the latitudinal dependence of the incoming fluxes only affect the regolith gardening rate, but not the physical properties of the regolith.This might change if the impact velocity is latitude-dependent.The largest dependency arises from the apex/antapex dichotomy and, thus, on the source of the impactors (Zuluaga & Sucerquia, 2018).Averaged over all possible impactor sources, this effect should be rather small.On the other hand, the average impact angle seems to be slightly dependent on the lunar latitude, with impacts on the poles being slightly steeper (∼60°) than at the equator (∼40°; Robertson et al., 2021).We conclude that we would not expect a large variation in regolith grain size or bulk density over the lunar latitudes, because at any lunar latitude, the regolith should be in a steady state after millions of years of constant exposure and bombardment.This is in line with Carrier (1973), who reviewed the measured grain size distributions of returned Apollo samples and concluded that the regolith has reached a steady state.
A second formation process for regolith is thermal fatigue, which has been proposed to be a very effective mechanism on small asteroids (Delbo et al., 2014), but should be less efficient on the Moon because of the much slower rotation period.The diurnal temperature variation on the lunar surface decreases from the equator toward ±80°latitude by approximately a factor of two.Experiments by Chen et al. (2006) showed that thermal cycling leads to compaction in granular material and thermal cycling is the mentioned process that can affect the regolith bulk density as a function of latitude.Metzger et al. (2018) performed thermal cycling experiments with a lunar regolith simulant and found an increase in bulk density of 1.3% for a temperature amplitude of ΔT = 187 K and of 0.26% for the control group with ΔT = 0 K. Divoux (2010) reviews thermal compaction processes and shows that even after thousands of thermal cycles no steady state is reached.Although this is a large number of cycles, it is unclear how to extrapolate these results to the lunar environment and millions to billions of years of evolution.Future work and in situ data at high latitudes are needed to investigate this further.Due to the lower day-night temperature variations at higher latitudes, it is conceivable that the deep layer densities may be somewhat smaller toward the poles than in the equatorial regions.We will get back to this point in Section 4.2.3.

Variation of Solar Incidence Dependent Albedo
The bolometric Bond albedo is a function of solar incidence angle θ, but the exact dependence is subject to a large uncertainty.Because the reflectance measurements of the Diviner solar channels (Paige, Foote, et al., 2010) and the Lunar Orbiter Laser Altimeter (LOLA; Smith et al., 2010) on board LRO are both acquired normal to the surface, they cannot be used to infer the angular dependence.Yu and Fa (2016) and Hayne et al. (2017) both adopted the empirical expression from Keihm (1984) to describe the angular dependence of the bolometric Bond albedo, but used different values for the coefficients a and b.In this work, we applied the description by Feng et al. (2020), who used Diviner surface temperatures acquired at noon as a function of latitude to determine the angular dependence (see Section 3.5 and Equation 13).However, at high latitudes shadows are present due to roughness and therefore the albedo inferred by Feng et al. (2020) must include these shadowing effects and contains some degree of uncertainty.Figure 9 gives an overview on the functional behavior of these albedos.
We now intend to investigate whether the observed increase in RMSD with latitude (see Figure 8) can be removed by applying an alternative dependence of the albedo on the solar incidence angle.Qualitatively, we find that a steeper increase in bolometric Bond albedo than predicted by Equation 13for the highest incidence angles is required to reduce the modeled surface temperatures.We apply Equation 15and allow a variation of the coefficients a and b as well as the albedo at zero incidence angle A 0 .The best fit values are obtained by minimizing the total RMSD for the modeled and measured regolith temperatures (taking all latitudes into account) and are found to be a = 0.10, b = 0.02, and A 0,H = 0.09 for the highlands high conductivity model case.We also tried to start with the best thermophysical model fit at 80°latitude and varied the albedo function for the lower incidence angles, but this resulted in a bolometric Bond albedo that is nearly zero at low latitudes and therefore very unrealistic.
For the highlands low conductivity model as well as for the maria, we adopt the same values for the coefficients a and b as found for the highlands high conductivity model.However, in the case of the maria, the value of the albedo at zero incidence angle A 0 needs to be modified according to the LOLA normal albedo measurements.As already shown in Feng et al. (2020), the LOLA normal albedo measurements show two peaks, one for the mean albedo of the highlands with A n,H = 0.31, and one for maria with A n,M = 0.17, respectively.Because LOLA measures the normal albedo at a wavelength of 1,046 nm, the values have to be scaled down to convert to the broadband solar albedo at zero incidence angle.Previous works found scaling factors of f = 0.39 (Feng et al., 2020) and f = 0.49 (Hayne et al., 2017).With the best-fitting albedo at zero incidence angle for highlands found in this work A 0,H = 0.09, the scaling factor is f = 0.29.Thus, the albedo at zero incidence angle for the maria is calculated to be A 0,M = 0.05.The solid curve in Figure 9 shows the resulting bolometric Bond albedo as a function of solar incidence angle compared to the albedo models used in previous works (Feng et al., 2020;Hayne et al., 2017;Vasavada et al., 2012;Yu & Fa, 2016).
We also include the albedo function derived from laboratory measurements of Apollo samples by Foote et al. (2020), which in contrast does not show a steep increase in albedo at high incidence angles.In fact, there is a gap between the albedo measured in the laboratory and the albedo used in thermal models to fit thermal emission measurements.One reason for this could be that the extreme surface roughness of the Moon is difficult to reproduce in the laboratory.A simple energy balance argument suggests roughness is the source of the apparent discrepancy: at high solar incidence angles, rough surfaces will produce a wider temperature distribution, due to shadowing, which intrinsically emits less infrared radiation than a smoother surface.Therefore, to balance absorbed insolation with emitted infrared radiation over the complete diurnal cycle, the rough surface must have higher albedo; this effect should increase drastically at higher latitudes.
If the albedo were to increase less with solar incidence angle, the discrepancy between modeled and measured regolith temperatures would be even more pronounced at higher latitudes, assuming constant regolith properties.
Figure 9.The bolometric Bond albedo as a function of incidence angle, as applied in previous works (Feng et al., 2020;Hayne et al., 2017;Vasavada et al., 2012;Yu & Fa, 2016) and measured in the laboratory (Foote et al., 2020).The solid curve shows the albedo function used in this work.The steeper increase of albedo with incidence angle is able to significantly reduce the observed deviation between modeled and measured regolith temperatures at higher latitudes.
Apart from the increase with solar incidence angle, there is also a difference in the albedo value at zero incidence angle between previous works and this work.Previous works found A 0,M = 0.07 for the maria and A 0,H = 0.12-0.16for the highlands.Foote et al. (2020) compared modeled temperatures applying their albedo model with Diviner bolometric temperatures at the corresponding landing sites and found a good agreement for the maria, but a deviation for the highlands, which would require an albedo value of A 0,H = 0.093, closer to the value found in this work.
Figure 10 shows the modeled regolith temperatures as well as the RMSD as a function of latitude after the albedo variation.The steeper increase of albedo with solar incidence angle almost completely removes the previously observed systematic offset between model and measurement.The RMSD remains below 1 K for all latitudes.However, there remains a small systematic deviation in the slope of the modeled nighttime temperatures at higher latitudes, which is manifested in a small increase of the RMSD values with increasing latitude.
In summary, varying the bolometric Bond albedo as a function of incidence angle allows us to significantly reduce, but not eliminate, the observed latitudinal gradient in RMSD between model and measurement.The albedo correction is supported by the fact that the onset of the initial latitudinal gradient (∼30°) coincides with the angle of repose of granular media and therefore likely is a roughness effect.With respect to the strong latitudinal gradient found in the work of Yu and Fa (2016), this study supports the suggestion of Hayne et al. (2017) that the gradient they observed could be a solar incidence-dependent albedo effect.

Variation of Intrinsic Regolith Properties With Latitude
With the microphysical model, we can make a prediction about the intrinsic variation of regolith properties with increasing latitude if the albedo correction were not applied.These predictions are valuable for planning future missions targeting higher latitudes and can be compared with future in situ measurements and returned samples.We investigate the best-fitting regolith properties at the highest latitudes (50°latitude for maria, 80°latitude for highlands) by varying the grain radius, the transition width and the deep layer density in the microphysical model.Figure B4 in Appendix B5 shows the results in the form of heat maps.The comparison of these heat maps with the ones obtained for the lunar equator (Figure 6) shows for the higher latitudes solutions with lower deep layer densities, larger transition widths and an expansion of the solution space to larger grain radii at higher latitudes.
The resulting bulk density and volume filling factor profiles for the three best fits for the maria, low and high Figure 10.The Diviner regolith temperatures as a function of latitude are plotted along with the modeled regolith temperatures (TPM) assuming constant regolith properties and applying the solar incidence-dependent albedo variation.The steeper increase of albedo with solar incidence angle reduces the observed systematic difference between model and measurement significantly, but a small deviation in the slope at higher latitudes remains.conductivity highlands model are illustrated in Figure 7 together with the ones obtained at the lunar equator.Due to the rather small variation in latitudes in the case of maria, all density profiles are fairly similar.In contrast, the best solutions for the highlands at ± 80°clearly show that the deep layer density is smaller and the transition from low to high densities is somewhat smoother than at the equator.

Journal of Geophysical
Due to the larger extent of latitudes, we will focus in the following on the highlands.In Figure 11, we show for the high conductivity highlands model the nighttime temperatures (top row) and the resulting RMSD values (bottom row) as a function of latitude when allowing a variation of only one parameter with latitude (middle row) and holding the other two parameters constant.We find that a variation of the regolith grain radius with latitude cannot remove the systematic increase of the RMSD with increasing latitude (left column).However, varying the regolith deep layer density (right column) or the logarithmic transition width (middle column) results in a better fit to the Diviner data with latitude than varying the bolometric Bond albedo as a function of solar incidence angle.
While we find that the albedo variation provides a poorer fit, we also note that the fixed form of the albedo function (Equation 15) does not allow us to choose any curve, in contrast to the free variation of the logarithmic transition width and deep layer density for each latitude.Figure 12 shows the resulting RMSD from the variation of solar incidence-dependent albedo and the variation of intrinsic regolith properties as a function of latitude.For comparison the RMSD of a third-order polynomial fit to the Diviner data is also calculated and the difference between the respective parameter variation and the RMSD of the polynomial fit is shown in the lower panel of Figure 12.The variation of solar incidence dependent albedo results in a better fit to the Diviner data than the Nighttime temperatures for the highland high conductivity model and latitudes between ±0°and ±80°under variation of the regolith properties, using the standard albedo (Section 3.5).In the thermal model, one of the three free parameters is varied, while the other two are kept constant at the best-fit parameter values of r = 40 μm, Δ = 0.1, ρ d = 1,700 kg m 3 .The grain radius (left column) is varied in steps of 2 μm, the transition width (center column) in steps of 0.02 and the deep layer bulk density (right column) in steps of 10 kg m 3 .The resulting RMSD between model and measurement as a function of latitude is shown in the three bottom graphs.
variation in regolith grain radius.Up to a latitude of ∼60°the variation of solar incidence dependent albedo and transition width show similar RMSD, but for the highest latitudes, the variation in transition width gives a better result.The predicted increase in transition width with increasing latitude leads to a slower increase of bulk density with depth (see Figure 7).Overall, the variation of the deep layer density yields the best result and removes the systematic deviations between model and observations almost completely.We performed the same parameter variation for the low conductivity highlands model and find the same behavior and predictions (see Figure B5).
Thus, this version of the model predicts a decrease of the deep layer density with increasing latitude resulting in an overall lower bulk density.As pointed out in Section 4.2.1, thermal cycling is more effective at higher day-night temperature differences (i.e., toward the lunar equator).At a depth of ∼2 cm, at which the local bulk density has almost saturated to the deep layer density, the temperature difference in the high conductivity highlands between noon and sunrise is ∼160 K at the equator and ∼75 K at ±80°latitude, respectively.Thus, the diurnal thermal expansion is roughly a factor of two lower at the high latitudes so that the decrease in bulk density and volume filling factor with increasing latitude is at least qualitatively in agreement with this effect.

Conclusion
In this work, we used a microphysical thermal model for the lunar regolith, which describes (a) the increase of bulk density and volume filling factor with increasing depth (i.e., increasing lithostatic pressure) as a function of grain size and (b) the thermal conductivity as a function of grain size, temperature and volume filling factor.The regolith in the lunar maria and highlands was investigated separately and characterized by the material differences in albedo and grain density (maria model and high conductivity highlands model) and an additional difference in thermal conductivity (low conductivity highlands model).The modeled thermal conductivity as a function of temperature was fitted to laboratory measurements of returned Apollo samples.The comparison of modeled temperatures with regolith temperatures measured by Diviner allowed us to determine the best values for the regolith grain radius, the logarithmic transition width from surface to deep layer bulk density and the deep layer density.However, these three model parameters are not uniquely constrained by the Diviner data.Each parameter set is considered to be a representative of the solution space describing all possible parameters combining model approach and data set.This prevented us from predicting exact values.
Our key findings are summarized in the following: 1.The regolith temperatures measured by Diviner are very similar for maria and highlands (Figure 1).The regolith temperatures in the highlands and maria can be well described by the mean grain size values obtained from returned Apollo samples.Assuming similar thermal conductivity for highlands and maria (Figure 3, upper curve), the derived volume filling factor profiles are fairly similar, but the resulting bulk density is higher for maria due to the higher grain density (Figure 7).When investigating a significantly lower thermal conductivity for the highlands (Figure 3, lower curve), this also formally results in a very good fit to the Diviner data, but with a bulk density that lies well outside the range of the Apollo data.2. We investigated the regolith properties up to ±80°latitude.First, we showed that for constant regolith properties, namely those derived at the equator, a variation of the solar incidence-dependent albedo can significantly reduce the initially observed deviation between model and Diviner measurements at higher latitudes (Figure 10).The onset of the initial latitudinal gradient coincides with the angle of repose of granular media and therefore likely is a roughness effect.However, a small systematic deviation between model and observation remains at high latitudes and the resulting albedo is slightly lower at low latitudes and slightly higher at high latitudes than any assumption made before (Figure 9).Whether this is realistic or not needs to be investigated in future.3.With the microphysical model at hand, we make a prediction on the intrinsic latitudinal variation of the regolith properties, under the assumption that the applied solar incidence-dependent albedo from the literature is correct (Figure 11).We predict a decrease in bulk density and volume filling factor with increasing latitude.In more detail, the best solution is identified to be a decrease in deep layer density above >30°latitude.The second best solution is a slight increase in logarithmic transition width with increasing latitude and both variations can fit the data better than the albedo correction.In addition, we showed that a variation in grain radius alone cannot explain the Diviner measurements at higher latitudes.These predictions are valuable for future mission planning and can be compared with in situ measurements and returned samples from the upcoming missions to higher latitude regions.

Outlook
This section summarizes the gaps in our knowledge of the lunar regolith that are relevant to this work and provides suggestions for future measurements and modeling.
1. Since thermal conductivity is a critical input parameter for thermal models, accurate calibration of the thermal conductivity of lunar regolith is needed, in particular studies of the possible difference in thermal conductivity between material from the maria and the highlands.However, laboratory measurements may not accurately represent lunar conditions, while in situ measurements are an order of magnitude higher (Langseth et al., 1972(Langseth et al., , 1976)), but may suffer from effects such as compaction (Grott et al., 2010).2. The stratification model applied in this study is based on laboratory measurements of cohesive-particle layers, which are described by an analytical function (Equation 4) including the empirical transition width Δ (Güttler et al., 2009).A fully physical model for the volume filling factor as a function of pressure is needed.The model should be complemented by a non-invasive in situ measurement of the stratification of lunar regolith in the first 10 cm.Here, measurements from the Chang'E-2 microwave radiometer may be helpful, which are more sensitive to deeper regolith layers, and therefore can help to constrain the stratification and thermophysical properties of the lunar regolith (Feng et al., 2020).The data could also uniquely constrain the parameter set investigated in this study and solve the degeneracy problem.However, for the interpretation of the microwave data, a radiative transfer model is needed and therefore the optical properties of the regolith must be known (Bürger, Glißmann, et al., 2023;Feng et al., 2020).3. The interpretation of the Diviner measurements is based on the assumption that only the uppermost regolith particle emits the measured radiation.However, due to the fairy castle structure in the uppermost regolith layers, the uppermost particle is not precisely defined, and since the Diviner wavelengths are not necessarily smaller than the regolith particle size, the thermal emission of the underlying particles could also play a role.4. The bolometric Bond albedo as a function of the solar incidence angle has a critical influence on the modeled regolith temperatures, but is subject to a large uncertainty.Diviner solar channel measurements during emission phase function campaigns could be analyzed to constrain the bistatic phase function of the lunar regolith over the solar spectrum.
For the approximation of the spatial derivatives, the grid set-up will be introduced first.
The layout for the general case of a non-uniform grid is illustrated in Figure B1 and the approach follows the method presented in Patankar (1980).The grid points 0, 1, …, i 1, i, i + 1, …, N are located at respective depths of x 0 , x 1 , …, x i 1 , x i , x i+1 , …, x N .Focusing on grid point i and its two neighboring grid points i 1 and i + 1, the following approximation for the spatial derivatives can be made Note.For the free parameters, the deep layer bulk density, the logarithmic transition width and the regolith grain radius, the range of values studied is listed.
Journal of Geophysical Research: Planets The first derivatives at the interfaces (equal the layer boundaries) i + 1 2 and i 1 2 are with the temperatures T i 1 , T i , T i+1 at the respective grid points i 1, i, i + 1 and the corresponding depths x i 1 , x i , x i+1 .The grid space Δx i between the two interfaces for grid point i can be expressed as and the discretization for the spatial derivatives therefore reads Plugging Equations B2 and B6 back into Equation B1 and rearranging the terms then leads to the following discretization of the heat transfer equation, which allows to derive the change in temperature ΔT i = T i (t + Δt) T i (t) in layer i with time step Δt.
The thermal conductivity is an interface property and gets evaluated at the respective layer boundary.Since the thermal conductivity is a function of temperature and volume filling factor, these properties are derived at the interface location by linear interpolation over the temperatures and volume filling factors of the adjacent layers, and the resulting thermal conductivity is calculated.
At the upper boundary, the first grid point i = 0 is located directly at the surface x 0 = 0 and the change in surface temperature ΔT 0 is calculated by with the resulting energy loss or gain ΔE in the first layer depending on the absorbed solar radiation, the heat conducted in or out of the surface layer from the subsurface and the thermal emission to space.The discretization of the heat conduction term at the upper boundary reads To reduce the runtime of the simulation, a non-uniform grid with an equal grid spacing in logarithmic space is chosen.This grid is finer at the surface, where temperature variations with depth are steep, and coarser in deeper layers, where variations are slower.For the grid implementation in this study, the grid spacing between the first two grid points is set to Δx 0 = 2 mm, the total depth of the simulation is fixed to x N = 2 m, and the number of grid points is set to N = 50.An investigation of the surface grid spacing shows that the selected value of Δx 0 = 2 mm and a two-and four-times better resolution lead to the same result within 0.1 K for local times between 7:30 p.m. and 5:30 a.m.
At the start of the simulation, the temperature in each layer T i is initialized following the approach presented in Hayne et al. ( 2017) with T N and T 0 being the temperature at the lower and upper boundary, respectively.At the lower boundary, where the depth is greater than the depth affected by the diurnal temperature variations, the temperature is set constant and is chosen to be represented by the equilibrium mean annual temperature calculated from the balance of incoming solar insolation and thermal emission.In this study, the equilibrium mean annual temperature is slightly scaled down by to match the Apollo 15 and 17 sub-surface temperature measurements of 252 and 256 K, respectively (Keihm & Langseth, 1973;Langseth et al., 1973).The temperature at the upper boundary is represented by instantaneous radiative equilibrium at noon The equilibration time of the model is determined by calculating the mean difference in surface temperature between two successive runs over a lunar rotation period (one synodic month).The average variation in surface temperature is found to be <10 2 K after 21 rotation periods and we chose to evaluate the final temperatures after a model run of 24 lunar rotation periods (∼2 years).In this study, we focus on the analysis of surface temperatures.However, for the analysis of simulated temperatures in deep layers (∼1 m), longer equilibration times are required to remove the effect of temperature initialization.
A measure for the numerical stability in an explicit scheme is the Fourier number Fo, which must be less than 0.5 in order for the simulation to be stable.It can be derived for the general case of a non-uniform grid and a temperature-dependent thermal conductivity as follows Since the thermal conductivity λ, density ρ, heat capacity c p and layer thickness vary over depth and/or time, the Fourier number is not a constant through one simulation run.We chose a fixed time step of Δt = 25 s, which results in acceptable Fourier numbers for all simulation runs.

B2. Thermal Model: Parameters
The values of all input parameters for the thermophysical model are summarized in Table B1.

B3.1. Temperature Dependent Thermal Conductivity of the Solid Material
For the introduction of a temperature-dependent thermal conductivity of the solid material, we looked at thermal conductivity measurements of meteorites at low temperatures conducted by Opeil et al. (2012).Each measurements was fitted with a fifth-order polynomial and extrapolated for temperatures T ≥ 300 K, taking into account the local slope at T = 300 K.The resulting curves were then normalized to λ = 1 W m 1 K 1 at T = 200 K and a master curve was fitted to the individual curves λ solid (T) = a + bT + cT 2 + dT 3 + eT 4 + f T 5 (B15) with fit-coefficients a = 0, b = 1.8866 × 10 2 W m 1 K 2 , c = 1.3863 × 10 4 W m 1 K 3 , d = 4.8674 × 10 7 W m 1 K 4 , e = 8.0380 × 10 10 W m 1 K 5 , f = 5.0906 × 10 13 W m 1 K 6 .Figure B2 shows the measurement curves together with the resulting master curve.It is expected that the thermal conductivity approaches zero for T → 0 K (Woods-Robinson et al., 2019) and indeed the Opeil et al. (2012) measurements show a sharp drop in thermal conductivity at low temperatures.To avoid nonphysical results at the lowest temperatures, the zero intercept of the master curve was enforced by setting a = 0. Note that because of the somewhat arbitrary normalization of the master curve, when fitting the model to the measurements, the scaling factor χ (see Equation 11) now accounts not only for the irregularity and non-monodisperse size-distribution of the regolith particles, but also for the scaling of λ solid (T ).

B3.2. Calibration to Laboratory Measurements
The thermal conductivity model is fitted to the laboratory measurements (see Figure 3) by varying the scaling factor χ (see Equation 11).Initially, two different fits are investigated, the maria model and the low conductivity highlands model.For the calibration, the volume filling factor and the grain radius of the corresponding samples are required.These are taken from the literature and are summarized in Table B2.For the Apollo 16 sample, we only give an upper limit of the mean grain radius, because the measured sample listed in Graf (1993) contains particles up to a size of 16 mm, however, the sample for the thermal conductivity measurements contains only a few grams and therefore the average grain size is likely lower.In this case, we tested different grain sizes in the low conductivity highland model fits and found the best fit for r = 20 μm and a resulting χ = 0.14.For the maria model and the Apollo 11 sample, we obtain χ = 0.73 and for Apollo 12, we get χ = 0.65, which results in a mean value for the maria of χ = 0.69.The additional high conductivity highlands model assumes the same thermal conductivity for maria and highlands and therefore is based on the fit of the maria model with χ = 0.69.The fit results are shown in Figure 3.

B4.1. Daytime and Nighttime Surface Temperatures
As was already shown by Vasavada et al. (2012) in great detail, some thermal model parameters have separate influences on daytime and nighttime surface temperatures.During the day, thermal conduction is negligible compared to radiation, since almost all of the incoming energy is re-radiated.The absorbed solar energy is a function of the albedo, while the radiated energy depends on the emissivity.Therefore, both of these parameters  Opeil et al. (2012).The measurements were fitted with a fifth-order polynomial, extrapolated for temperatures T ≥ 300 K and the resulting curves were then normalized to λ = 1 W m 1 K 1 at T = 200 K. Note.
have an effect on daytime surface temperatures.During the night, the incoming solar energy is zero and the amount of energy conducted in and radiated from the surface layer is very similar.Hence, modeled nighttime surface temperatures are sensitive to model parameters, such as thermal conductivity and bulk density, while these parameters have only a small effect on daytime surface temperatures.In consequence, if we want to learn about regolith properties, like grain radius, that influence the stratification and thermal conductivity, we need to focus on nighttime surface temperatures.

B4.2. Model Parameters
In this subsection, we show the sensitivity of the three free parameters in the thermophysical model-the grain radius, the logarithmic transition width, and the bulk density of the deep layer-to the modeled nighttime regolith temperatures.In addition, we investigate the influence of other model parameters that were assumed to be constant in this study.
Figure B3 shows the resulting nighttime surface temperatures for a variation of the respective parameter values by ± 25%.An increase in grain radius and deep layer density or decrease in transition width leads to increased nighttime surface temperatures, as the bulk density increases more rapidly in the upper layers (see Figure 2).In contrast, the variation of the surface bulk density has a negligible effect on the modeled nighttime surface temperatures.We also tested an explicit variation of the network thermal conductivity and radiative conductivity and found that network conductivity dominates during nighttime.As mentioned earlier, the albedo at zero incidence angle has little effect on nighttime temperatures, but plays an important role for daytime surface temperatures.The nighttime surface temperatures are sensitive to emissivity, which affects the whole level of temperature, but does not change the slope.
We note that the three free parameters predominantly governing the shape of the bulk density profile-the grain radius, the logarithmic transition width, and the bulk density of the deep layer-all affect the whole level of nighttime surface temperature and the slope of the nighttime cooling curve to some degree.The resulting degeneracy when fitting the modeled temperatures to Diviner surface temperatures is shown and discussed in Section 4.1.The investigated range in grain radius is r = 10-100 μm, motivated by the range in mean grain radius from the returned Apollo samples.For the transition width Δ, values between Δ = 0.1 and Δ = 1.2 were tested, including the value Δ = 0.58, experimentally found by Güttler et al. (2009).For the deep layer bulk density, we chose as an upper boundary the bulk density of highland crust with ρ d ≲ 2,500 kg m 3 (Han et al., 2014;Wieczorek et al., 2013).

Figure 2 .
Figure 2. The volume filling factor and bulk density as a function of depth for regolith layers consisting of highland particles with a radius of 20 μm (left) and 80 μm (right).The particle radius determines the depth of the turning point and thus, how easily a regolith layer can be compressed under its own weight.The transition width Δ determines the steepness of the transition from the surface volume filling factor ϕ s to the deep layer volume filling factor ϕ d .

Figure 3 .
Figure 3. Measured thermal conductivity of returned Apollo samples as a function of temperature(Cremers & Birkebak, 1971;Cremers & Hsia, 1974;Cremers et al., 1970).The Apollo 11 and 12 data were measured at a bulk density of 1,300 kg m 3 and are associated with material from maria, while the Apollo 16 sample was measured at a bulk density of 1,500 kg m 3 and is associated with highland material.The corresponding maria model fit and low conductivity highlands model fit are illustrated.In this study, we add a third scenario, the high conductivity highlands model, where the highlands are given the same thermal conductivity as the maria with χ = 0.69.

Figure 4 .
Figure 4. Modeled regolith properties and resulting temperatures for the maria, high conductivity and low conductivity highlands as a function of depth below the surface.The thermophysical model was run for a fixed parameter set of r = 50 μm, Δ = 0.4, and ρ d = 1,800 kg m 3 (upper panels) as well as for a fixed deep layer volume filling factor of Φ d = 0.64 (lower panels) at 0°latitude.The difference between the maria model and the high conductivity highlands model is due to an assigned difference in albedo and grain density, whereas the low conductivity highlands model is assigned an additional difference in thermal conductivity.The temperature profile as well as the specific heat capacity and thermal conductivity profile are illustrated for two different local times, noon and the time just before sunrise.It is important to note two aspects of thermal conductivity.First, it is a function of volume filling factor (see Section 3.3) and therefore the maria and high conductivity model only show similar thermal conductivity profiles for similar volume filling factor profiles.Second, the thermal conductivity is an interface property and therefore evaluated at the layer boundaries (see Appendix B1).

Figure 5 .
Figure 5.The filtered Diviner regolith temperatures and their mean values (see Section 2 for details) at the lunar equator for maria and highlands are illustrated together with the modeled regolith temperatures that result in the minimum RMSD.The minimum RMSD was found for the following parameter combinations: maria r = 30 μm, Δ = 0.4, ρ d = 2,200 kg m 3 with RMSD = 0.31 K, high conductivity highlands r = 40 μm, Δ = 0.1, ρ d = 1,700 kg m 3 with RMSD = 0.17 K and low conductivity highlands r = 70 μm, Δ = 0.2, ρ d = 2,500 kg m 3 with RMSD = 0.21 K.However, it is important to note that the solutions show a degeneracy, which is illustrated in Figure6.The lower panels show the temperature difference between model and measurement ΔT = T TPM T Diviner (solid curves).For comparison, as a measure of the theoretical best fit, a third-order polynomial was fitted to the Diviner regolith temperatures and the resulting temperature differences are plotted as the dashed curves.

Figure 6 .
Figure6.Derived regolith properties at 0°latitude for maria (upper panel), high conductivity highlands (center panel) and low conductivity highlands (lower panel).The three free parameters in the thermophysical model, the grain radius, the transition width and the deep layer density, and the comparison with Diviner regolith temperatures results in non-unique solutions.The heat maps show the minimum RMSD for the combinations of two of each of the three free parameters, with the third parameter being variable.For better readability, all RMSD ≥1 K are colored in black.To better constrain the non-unique solutions, the range of the most likely grain radius from returned Apollo samples is illustrated (left and right columns), ranging from 20 μm (solid line) to 50 μm (dotted line)(McKay et al., 1991).For the highlands, we also add the mean grain radius of 40-140 μm from Apollo 16 samples with 40 μm as lower limit (dashed line), as determined byHeiken et al. (1973).As a constraint on the derived bulk density, we draw boundaries between ρ bulk = 1,500 kg m 3 and ρ bulk = 2,000 kg m 3 (center column), as obtained from the Apollo drill core experiments(Carrier, 1974).For illustration purposes, since the bulk density is a function of the grain radius, the transition width, the deep layer density and the depth below the surface, the bulk density is evaluated at 1 m depth and for the upper and lower limit of the grain radius, namely 20 μm (solid line) and 50 μm (dotted line).

Figure 7 .
Figure7.The derived volume filling factors and bulk densities as a function of depth below the surface for maria, high conductivity and low conductivity highlands and their three best fits (see Figures6 and B4) at 0°latitude and higher latitudes (50°maria, 80°highlands).The best fit shown in Figure5is highlighted with a dashed line.The maria and high conductivity highlands model show similar volume filling factor profiles, but a difference in bulk density due to their different grain density.The bulk densities of the Apollo core tubes(Mitchell et al., 1974) are shown by the vertical bars and the average bootprint depth(Carrier et al., 1991;Houston et al., 1972;Mitchell et al., 1973) is marked by the hatched box at the top panels.

Figure 8 .
Figure 8.Comparison between the modeled regolith temperatures (curves labeled TPM) and the Diviner regolith temperatures as a function of latitude (symbols, see legend), assuming constant regolith properties as derived for the lunar equator (see Section 4.1).The lower panels show the resulting root-mean-square deviation between model and measurement as a function of latitude.

Figure 11 .
Figure11.Nighttime temperatures for the highland high conductivity model and latitudes between ±0°and ±80°under variation of the regolith properties, using the standard albedo (Section 3.5).In the thermal model, one of the three free parameters is varied, while the other two are kept constant at the best-fit parameter values of r = 40 μm, Δ = 0.1, ρ d = 1,700 kg m 3 .The grain radius (left column) is varied in steps of 2 μm, the transition width (center column) in steps of 0.02 and the deep layer bulk density (right column) in steps of 10 kg m 3 .The resulting RMSD between model and measurement as a function of latitude is shown in the three bottom graphs.

Figure 12 .
Figure12.The upper panel compares the resulting RMSD from the variation of solar incidence dependent albedo (Figures 9 and 10) and the variation of intrinsic regolith properties (Figure11) as a function of latitude.The lower panel shows the difference in RMSD between the respective variation and the RMSD of a third-order polynomial fit to the Diviner data.A variation in deep layer density (predicts decrease with latitude, see Figure11) yields the best solution.

Figure B2 .
Figure B2.The master curve describing the temperature dependence of the thermal conductivity of the solid material, derived from measurements ofOpeil et al. (2012).The measurements were fitted with a fifth-order polynomial, extrapolated for temperatures T ≥ 300 K and the resulting curves were then normalized to λ = 1 W m 1 K 1 at T = 200 K.

Figure B3 .
Figure B3.The sensitivity of the modeled nighttime regolith temperatures to both free and fixed model parameters.For the analysis, a reference temperature profile for maria, with r = 50 μm, Δ = 0.4, ρ d = 1,800 kg m 3 at 0°latitude and other model parameters according to the values in TableB1, was created and the respective parameter values were varied by ±25%.An exception to the variation is the emissivity, where only a maximum value of ϵ = 1 is possible.

Figure B4 .
Figure B4.Derived regolith properties at higher latitudes for maria (upper panel), high conductivity highlands (center panel) and low conductivity highlands (lower panel).The lunar maria do not extend over latitudes above ∼50°, while the highlands in this work are studied up to a latitude of 80°.The heat maps show the minimum RMSD for the combinations of two of each of the three free parameters, with the third parameter being variable.For better readability, all RMSD ≥1 K are colored in black.

Table B2
Summary of the Literature Values for the Calibration of the Thermal Conductivity Model