© ESO 2019

1. Introduction

The thermal environment of the Moon can be divided up into three main regimes: day, night, and polar (Paige et al. 2010a). Due to the slow rotation of the Moon, the lunar day and the lunar night each last about 14 Earth-days which leads to quite extreme surface temperatures. Contrary to Earth, on the Moon no efficient transport of heat on a global scale exists due to its only tenuous surface-bound exosphere and the lack of an ocean. Furthermore, the surface is highly insulating and hence the daytime surface temperature is almost in equilibrium with the solar flux and reaches ≈400 K near the sub-solar region and cools down to ≈100 K during the night (Williams et al. 2017). The third regime, the polar temperature regime, is again quite different to the day and night regimes. Due to the small lunar obliquity of 1.54° the polar areas experience extreme illumination conditions. Many near-polar crater floors are deep enough to be in permanent shadow, and are referred to as permanently shadowed regions (PSRs). Here, the coldest lunar temperatures near 20 K can be found (Siegler et al. 2015). In addition, areas receiving almost continuous illumination exist which are found on the highest topographic locations near the poles where sunlight is not obstructed by nearby terrain (Gläser et al. 2014). Average maximum and minimum temperatures are found to be 202 K and 50 K, respectively (Williams et al. 2017).

Early measurements of the lunar thermal environment date back as early as 1930, when Pettit & Nicholson (1930) measured the diurnal surface temperatures and temperature profiles during lunar eclipses. Jaeger & Harper (1950) concluded from the temperature profiles during eclipse measurements that this could be explained by a temperature-dependent conductivity parameter or a nonhomogeneous surface, both of which being assumed to be constant until then. They introduced the idea that a 2 mm surface layer of dust with a relatively low thermal inertia lies on top of a layer with higher thermal inertia. Remote temperature measurements from Earth and spacecrafts as well as in-situ measurements and sampling from Apollo and Surveyor missions eventually lead to a much better understanding of the thermophysical properties of the lunar surface (Linsky 1966; Hemingway et al. 1973; Keihm et al. 1973; Langseth et al. 1976).

Early efforts were made by Buhl et al. (1968) to model lunar surface temperatures using idealized crater-shapes while also considering thermal re-radiation. Later, analytical and numerical thermal models were solved for idealized crater-shapes near the poles of Mercury and the Moon to investigate surface temperatures and near-surface ice taking into account multiple scattering effects (Paige et al. 1992; Ingersoll et al. 1992; Salvail & Fanale 1994). A more sophisticated study using a two-layer model with depth- and temperature-dependent thermophysical properties taking into account multiple scattering was carried out by Vasavada et al. (1999) for flat surfaces and surfaces within bowl-shaped and flat-floored polar impact craters. Thermal models were also used to study the effect of roughness on surface temperatures, cold-trap formation, and ice stability (Jamsa et al. 1993; Aharonson & Schorghofer 2006; Davidsson & Rickman 2014; Rubanenko & Aharonson 2017). Further, the thermal environment of small bodies (e.g., asteroids, comets, small planetary satellites, and dwarf-planets) were investigated in various studies applying thermal modeling (Kuhrt & Giese 1989; Lagerros 1997; Delbo et al. 2015; Hayne & Aharonson 2015; Hu et al. 2017). Most recently the Lunar Reconnaissance Orbiter (LRO) Diviner Lunar Radiometer, referred to as Diviner, mapped the global thermal state of the entire lunar surface for the first time (Paige et al. 2010a). Paige et al. (2010b) derived model-calculated temperatures and compared them to measured Diviner bolometric temperatures using a 500 m-scale triangular mesh based on topography derived from Kaguya LALT laser altimeter data (Araki et al. 2009). Vasavada et al. (2012) used the Diviner measurements to derive new estimates for the thermophysical properties of the lunar near-surface regolith. Bauch et al. (2014) used 100 m pixel⁻¹ topography derived from stereo observations made by the LRO Wide Angle Camera (WAC; Scholten et al. 2012) to model temperatures at the Apollo 17 landing site. Finally, Hayne et al. (2017) used a sophisticated thermal model to derive thermophysical properties of the lunar regolith from Diviner measurements.

The above- mentioned studies differ from the one presented here in at least one of the following ways: No (or simplified) terrain scattering or at least no multiple scattering approach has been considered (e.g., Aharonson & Schorghofer 2006; Kuhrt & Giese 1989; Hayne & Aharonson 2015; Hayne et al. 2017; Hu et al. 2017); only a constant, linear, or discontinuous model for the thermophysical parameters has been introduced (e.g., Paige et al. 1992, 2010b; Delbo et al. 2015; Rubanenko & Aharonson 2017); no heat conduction has been considered (e.g., Jamsa et al. 1993); additionally, in several of these studies idealized topography has been used, locally nonenergy-conserving finite difference schemes have been applied to solve the heat equation, Sun has been introduced as a point source, and no solar limb darkening effect, no Earth-shine or no radiogenic heat source in the sub-surface has been considered. Furthermore, some of the presented models require time-steps in the order of seconds and are unusable for the present study (e.g., Davidsson & Rickman 2014). In this study we introduce a thermal model which is applied to high-resolution lunar topography derived from the Lunar Orbiter Laser Altimeter (LOLA) data (Smith et al. 2017) in order to determine surface and sub-surface temperatures. Our model incorporates direct solar insolation while considering the solar limb darkening effect and obstruction by the local horizon, multiple scattering of sunlight and thermal re-radiation from nearby terrain, heat conduction in the sub-surface, radiation into space, a constant radiogenic heat source in the sub-surface, as well as scattered sunlight and thermal radiation coming from Earth. Furthermore, our method allows nonlinear dependences of the thermophysical parameters with depth and temperature to be incorporated, while being both globally and locally energy conserving.

2. Mathematical model

In this work we present a model for the temporal evolution of the near-surface temperatures on rocky solar system objects that do not have a considerable atmosphere. We base the presented model on the fundamental axiom of conservation of energy, which for the considered system can be formulated by means of the well-known heat equation:

(1)

where ρ is density, c is specific heat capacity, Λ is thermal conductivity, and T is temperature.

2.1. Governing equations

Let us use ℳ to denote the rocky object to be modeled and Ω ⊂ ℳ the model domain on which computations are carried out. We note that we allow Ω to be one-, two-, or three-dimensional, that is, Ω ⊂ ℝ^d, d ∈ {1, 2, 3}. We use Γ to denote the boundary of Ω, which is further decomposed into disjoint partitions Γ = Γ_s ∪ Γ_b ∪ Γ_l (see Fig. 1). Furthermore, n is the unit outer normal vector on Γ. Using this notation, the system of governing equations with boundary conditions reads as follows:

(2)

(3)

(4)

(5)

Fig. 1. Illustration of the model domain Ω (for d = 3) and boundary Γ, which is further decomposed into the part coinciding with the surface of the rocky object, Γs, the lateral boundary Γl, and the bottom boundary Γb. While there are generally no restrictions on the choice of the model domain, in this work we construct Ω by vertical extrusion of parts of the surface of the rocky object with a constant depth. Therefore, the surface topography is also illustrated on the bottom boundary Γb.

Equations (2)–(5) state that the energy conservation equation must be fulfilled within the domain Ω, which interacts with its surroundings via the heat fluxes q and g specified on the boundary segments Γ_s and Γ_b, respectively. These boundary conditions are explained in detail in the following section.

2.2. Boundary conditions

While we consider homogeneous Neumann boundary conditions on the lateral boundary Γ_l, the heat exchange with deeper parts of the rocky object is modeled via the function g on Γ_b (see Eq. (4)). We note that in this work we use a constant heat flux, that is, g = const. (see Sect. 4). We incorporate the interaction with space, Sun, and other solar-system objects in the function q on the surface boundary Γ_s (see Eq. (3)). These interactions comprise radiation into space, solar insolation, and contributions stemming from reflected sunlight and thermal radiation from further solar-system objects. In order to allow for a compact notation, let us introduce the set 𝒞 = {ϕ₁, ϕ₂, …, ϕ_n} of solar-system objects ϕ_k, k ∈ {1, …, n} from which reflected sunlight and thermal radiation is taken into consideration (we note that this includes the modeled object itself, i.e., ℳ ∈ 𝒞).

Using this notation we define the radiant emittance of a solar-system object ϕ ∈ 𝒞, at a location x ∈ ∂ϕ on its surface ∂ϕ, on the basis of the Stefan–Boltzmann law:

(6)

Here, ϵ^ϕ = ϵ^ϕ(x) is the hemispherical emissivity, σ_SB is the Stefan–Boltzmann constant, and refers to the bolometric surface temperature at position x and time t. Furthermore, we define the direct solar insolation

(7)

where S₀ is the solar irradiance and f_sun = f_sun(x, t)∈[0, 1] is the normalized visible fraction of the solar disk taking into account solar limb darkening (see Sect. 3.5) and obstruction by the horizon (Gläser et al. 2014) on ϕ as seen from the respective location x at time t. Here, ι_⊙ = ι_⊙(x, t) and A^ϕ = A^ϕ(x) are the solar incidence angle and the solar spectrum average albedo of the surface, respectively, and d_ϕ, sun = d_ϕ, sun(t) is the time-dependent distance of ϕ to the Sun in astronomical units (AU). Making use of (6) and (7) we now decompose the heat flux q such that

(8)

where the terms and refer to scattered sunlight and thermal radiation that reaches ℳ from the considered solar-system objects ϕ_k. The former is modeled by

(9)

Here, δ = δ(x,y) ∈ {0, 1} describes whether or not the terrain at location y ∈ ∂ϕ on the solar-system object is visible from location x ∈ Γ_s on the modeled object where η = η(x,y) is the respective emission angle (see Fig. 2 and view factors, e.g., Vasavada et al. 1999). Accordingly, ι = ι(x,y) denotes the incidence angle at position x for scattered light and radiation coming from location y. The distance between the two locations is denoted by d_ϕ, ℳ and the reflected light from location y is defined as

(10)

Fig. 2. Geometric relations between infinitesimal areas dA around the locations x and y. The normal vectors are denoted n(x) and n(y), respectively. The graphic shows how the surface at position y is illuminated by the Sun at an incidence angle ι⊙ and reflects light at an emission angle η which is received at location x at a distance d and incidence angle ι.

Similarly, we can now define the irradiance stemming from thermal radiation:

(11)

where the absorptivity α^ℳ = α^ℳ(x) describes how much of the incident radiation is absorbed at location x and where is defined after Eq. (6). We note that by Kirchhoff’s law the emissivity ϵ and absorptivity α must be equal, that is, ϵ^ℳ(x) = α^ℳ(x).

As mentioned above, ℳ ∈ 𝒞, and thus scattered light and radiation, is also taken into account from nearby terrain on ℳ itself. However, a large portion of the surface ∂ℳ is not visible from a location x ∈ ∂ℳ; for example due to the curvature of the modeled object there is a physical limit on how far terrain features can be seen. Therefore, we introduce an area of influence ℛ = ℛ(x) ⊂ Γ_s ⊂ ∂ℳ, which we use as a region of integration in Eqs. (9) and (11) instead of the entire surface ∂ℳ. Thus, the expressions for scattered light and radiation from nearby terrain on ℳ itself read

(12)

(13)

The distance d is the distance between location x and y, that is, d = d(x, y)=‖x−y‖.

3. Numerical model

We use this chapter to describe the finite-volume scheme that is applied to solve Eqs. (2)–(5) as well as the discrete approximations of the heat fluxes presented in Sect. 2.2. An extensive discussion of various finite-volume schemes can be found for example in Droniou (2014), however we provide a compact presentation of the scheme used in this work together with the resulting discrete form of our governing Eq. (2).

3.1. Finite-volume discretization

In this section, we outline the derivation of the discrete formulation of the model problem. We employ a finite-volume scheme for the spatial discretization, for which we now introduce a set of control volumes 𝒯 such that Ω = ∪_{K ∈ 𝒯}K. Integration of Eq. (2) over a control volume K ∈ 𝒯 yields

(14)

We note that the divergence theorem has been used here. Let us now introduce the set of faces 𝒮 of the discretization and for each cell K ∈ 𝒯 the subset 𝒮_K such that ∂K = ∪_{σ ∈ 𝒮K}σ. The second integral of the left side of Eq. (14) can then be split into the integrals over the faces of the cell, that is,

(15)

At the core of finite volume schemes lies the construction of the discrete flux approximation,

(16)

for which we are using a classical two-point flux approximation. Therefore, we introduce intermediate face unknowns on the face centers x_σ, which are subsequently eliminated by enforcing both continuity of temperature and fluxes across the faces σ. For an interior face σ between the two cells K, L, that is, , these conditions read (we denote by 𝒮_int the set of interior faces)

(17)

We note that we consider control volume-wise constant temperatures T_K, densities ρ_K, heat capacities c_K, and thermal conductivities Λ_K. The last three of the above parameters are possibly functions of both T and the position x. In this work, we evaluate the cell-wise values for these three parameters as functions of the cell centers x_K and the cell temperatures T_K. Using this notation we write the fluxes

(18)

with |σ| being the measure, i.e., the area of the face σ and where we have defined as the outer unit normal vector on σ as seen from cell K, d_K, σ = x_σ − x_K, and (see Fig. 3).

Fig. 3. Illustration of the geometrical setting between two cells K, L ∈ 𝒯. Depicted are the geometrical quantities required to describe the flux given in Eq. (18) across the face .

Assembling F_L, σ accordingly, insertion into Eq. (17) and elimination of T_σ yields the final expression for the discrete flux across a face σ between the two cells K and L:

(19)

Finally, approximating the time derivative in Eq. (15) by a finite difference, the discrete formulation reads

(20)

which must hold for each cell K ∈ 𝒯. We note that the superscripts t and t + Δt indicate time level on which the respective terms are evaluated; that is, in this work we are using an implicit time discretization scheme.

3.2. Surface temperature reconstruction

Equations (6) and (11) require the evaluation of temperatures at the surface Γ_s. Since the chosen finite-volume scheme does not readily provide the surface temperatures, we have to reconstruct them by locally balancing the participating heat fluxes, that is, on the basis of a local energy balance. We consider a control volume K of the computational grid that has a face σ on Γ_s, i.e., σ ⊂ Γ_s. Using the definition of the flux F_K, σ from Eq. (18) we can write the local balance of heat fluxes:

(21)

which essentially states that the discrete flux across σ must fulfill the boundary condition (3). Both F_K, σ and q depend on the sought surface temperature , while the latter even shows a nonlinear dependency. Furthermore, q depends on the surface temperatures of all boundary faces σ ∈ 𝒮 : σ ⊂ ℛ(x_σ), that is, all faces within the area of influence.

3.3. One-dimensional approximation

We assume that the near-surface temperature gradients on ℳ in the vertical direction are expected to be much larger than in the horizontal directions. Therefore, we adopt a 1D description neglecting horizontal heat conduction (Bauch et al. 2014; Kuhrt & Giese 1989). So far we have presented the mathematical and numerical models for a general, d-dimensional setting (with d ∈ {1, 2, 3}). Starting from a 3D domain Ω ⊂ ℝ³ (see Fig. 1), we now approximate this by a set of vertical columns, Ω_i ⊂ Ω, such that Ω ≈ ∪_iΩ_i (see Fig. 4). We use to denote the surface segment of the i-th sub-domain and introduce the measure . We then define the 1D approximations as the center lines of the columns Ω_i with associated measure a_i such that . We also note that is equal to the depth of the i-th sub-domain. Equations (2)–(5) are then solved separately in each sub-domain , which is advantageous from a computational point of view, as it makes the system perfectly suitable for parallelization.

Fig. 4. Illustration of the computational grid and the splitting of the domain Ω into a set of vertical sub-domains Ωi, for which two exemplary columns Ωi and Ωj with discretization are depicted. We note that Ωi are still 3D objects, while the actual computational domains are strictly 1D, extruded by the measure . Furthermore, the discrete area of influence ℛi ⊂ Γs around Ωi is illustrated, within which the scattering of sunlight and thermal radiation of neighboring surface segments , j ≠ i, is considered.

However, the individual sub-domains Ω_i or are still coupled through thermal radiation, that is, the evaluation of the integral in Eq. (13) requires knowledge of the surface temperatures in adjacent sub-domains within the area of influence. We note that in this discrete setting Eq. (11), evaluated for sub-domain i, reads

(22)

where the subscripts i and j refer to quantities evaluated on and , respectively, and where we have introduced the discrete analogon (around sub-domain i) of the area of influence, ℛ_i (see Sect. 2.2). denotes the temperature at the surface of sub-domain j, that is, on .

3.4. Approximation of flux density from solar-system objects

In order to determine reflected sunlight coming from nearby solar-system objects, for example from Earth, we divide the visible celestial disk in 18 ten-degree zonal segments. We then calculate the actual terminator on the disk which is visible from the observer at location at time t and determine the surface area of the lit part of each visible segment , k ∈ {1, …18} (see Fig. 5). The amount of reflected sunlight is then determined based on the incidence angle ι_⊙k, albedo and the distance of the solar-system object to the Sun d_ϕ, sun (we note that if the segment is not lit or not visible):

(23)

Fig. 5. 18 ten-degree zonal segments of the disk of the solar-system object along with the terminator and horizon as seen from the observer location . In this example no sunlit segments below 20°S are visible to the observer due to obstruction by the horizon.

Similarly, we adopt the 18 ten-degree zonal segments to model heat flux from solar system objects:

(24)

We note that here is the visible radiating area of segment k with surface temperature , independent of whether or not it is lit. Furthermore, η_ki and ι_ik in Eqs. (23) and (24) denote the emission angle from the scattering or radiating segment k and the incident angle of the incident scattered light and radiation from segment k at . The distance between the modeled and solar-system object is d_ik.

3.5. Approximation of f_sun

As mentioned previously in Sect. 2.2 we compiled f_sun from the fraction of the visible solar disk and solar limb darkening. The determination of the fraction of the visible solar disk is described in Gläser et al. (2014). Solar limb darkening describes the gradual decrease in brightness of the solar disk from its center towards the limb (see Fig. 6). The darkening is due to different observed depths of the solar atmosphere depending on the emission angle at which a point on the disk is observed (angle ψ and point C in Fig. 6). At the center, the deepest, warmest, and hence brightest layers can be observed.

Fig. 6. Solar limb darkening. This is the phenomenon by which the Sun darkens from its center to its limb. The actual visible solar disk is defined by Rvisible, the apparent radius of the Sun depending on the distance of the observer. All points closer than the plane defined by the apparent radius are visible to the observer and are observed under emission angles ψ, ψ ∈ [0, 90].

The ratio between the brightness at the center to any point on the visible solar disk is incorporated in the parameter f_sun using the formula given by Cox (2000):

(25)

The emission angle ψ can be calculated from the angles θ and γ:

(26)

The brightness of the Sun is dependent on the actual wavelength, which in our study is fixed to λ = 550 nm at the center of the visible spectrum. Subsequently the parameters at this wavelength are u₂ = 0.93 and ν₂ = −0.23 (Cox 2000). Setting w₀ = 1 − u₂ − ν₂, w₁ = u₂ and w₂ = ν₂ leads to

(27)

At an angle ψ = 0 the ratio is defined to be 1 and so the formula can be further reduced to consist of only two parameters:

(28)

with W₁ = −0.47 and W₂ = −0.23.

To derive a value for f_sun in our model (see Eq. (7)) we discretized the solar disk in five concentric rings (Fig. 7). The radii are all multiples of the radius of the innermost ring, that is,

(29)

Fig. 7. Solar disk discretized in five concentric rings from which we determine the actual visible solar disk ratio while considering solar limb darkening and obstruction by the horizon. The radius and area of each ring are denoted Rj and Pj, respectively. The points Cj are placed halfway between two successive radii except for the C1 which is placed right in the center of the disk. In this example only the brightest, most inner ring of the Sun is completely visible. All other rings are only partially visible due to obstruction by the horizon.

For each ring we then calculated an average brightness ratio assumed to coincide with points C_j using Eq. (28). Furthermore, we determined the visible areas P_j of each ring, where f_j ∈ [0, 1] is the fraction of the visible ring area:

(30)

To compile an area weighted average brightness value I_total for the entire visible solar disk we sum over the product of the average brightness per ring with its area,

(31)

The retrieved brightness value needs to be scaled by I_max which is the maximum brightness in case the entire solar disk is visible, that is,

(32)

The final value for f_sun is the ratio of I_total to I_max, which now incorporates solar limb darkening while taking into account the part of the solar disk that is seen and to what extent,

(33)

When ≈20% or ≈80% of the solar disk is visible the limb darkening effect has its biggest impact. Here a difference in irradiance of ≈10% is reached compared to a solar disk of constant brightness which roughly also corresponds to a 10 K difference in surface temperature.

3.6. Solution strategy and multiple scattering

We solve the governing Eqs. (2)–(5) on a set of 1D domains (see Sect. 3.3) using a finite-volume scheme for the spatial and a backward Euler scheme for the temporal discretization. Due to the nonlinear dependences of the parameters involved (see Sect. 4.2), Newton’s method is employed as the nonlinear solver in each time step together with the Gauss–Seidel method as the iterative linear solver in each Newton iteration. As mentioned in Sect. 3.3, the individual domains are partly coupled due to the reflected sunlight and thermal radiation from neighboring terrain taken into account within an area of influence. Here, we iterate until from one iteration to the next the relative change of the norm of all incoming reflected sunlight and thermal radiation from nearby terrain for the entire scene is less than 0.001% (abort criterion). Regardless of the abort criterion however we run at least three iterations to ensure several propagation cycles of energy transport by scattering and radiation. This way we ensure propagation of thermal radiation and scattered light into shadowed terrain which is especially important to derive meaningful temperature estimates inside PSRs. Here, scattering and thermal radiation are the only heat source and single scattering would not sufficiently describe the energy transport and surface temperatures.

The software is implemented for General Purpose Computation on Graphics Processing Unit (GPGPU) using OpenCL (Gaster et al. 2011) which allows for fast computation of large areas and long timescales.

4. Model parameters – the lunar case

In this study we aim to evaluate surface and sub-surface temperatures near the lunar poles. All assumptions, parameters, and constants used to achieve this task are introduced in this section.

4.1. Constants

On the lower domain boundary Γ_b we set g = 0.016 W m⁻², representing a constant radiogenic heat source stemming from radioactive decay in the lunar interior (Langseth et al. 1976; Vasavada et al. 2012). We introduce the solar spectrum average albedo to be constant A^moon(x)=A^moon = 0.2 (Paige et al. 2010b). We assume lunar emissivity and absorptivity to be constants, that is, independent of changing material properties and therefore location, ϵ^moon(x)=ϵ^moon = α^moon(x)=α^moon = 0.95 (Hayne et al. 2010). The Stefan–Boltzmann constant is σ_SB = 5.670367 × 10⁻⁸ W m⁻² K⁻⁴ and in our study we set the solar irradiance to S₀ = 1366 W m⁻² (Stubbs & Wang 2012).

4.2. Thermophysical parameters

Fits of the thermophysical parameters to Apollo core samples as well as to recent Diviner temperature measurements have shown a nonlinear dependency on either depth z or temperature T or both of them (Vasavada et al. 2012; Hemingway et al. 1973). In our approach we adopted density to be modeled depth-dependent as described in Vasavada et al. (2012),

(34)

where ρ_top = 1300 kg m⁻³ and ρ_bottom = 1800 kg m⁻³. The discovery of the nonlinear dependency of density with depth when evaluating Diviner temperature maps led to the introduction of the H-parameter (Hayne et al. 2013, 2016; Piqueux et al. 2016). This parameter describes the vertical variations of the thermal properties of the regolith within the upper few centimeters and supersedes the previously used model where only two discrete layers were assumed (Vasavada et al. 1999; Mitchell & de Pater 1994). Typical values for the H-parameter lie between 0 and 10 cm, which mainly depends on the age of the surface, with the background value of 6 cm used in this model (Hayne et al. 2013). Depending on the depth z, all densities between ρ_top and ρ_bot are possible, see Fig. 8. The strong increase of density in the upper layers of regolith can be explained by the constant micro-meteor impacting and the resulting overturn within the top layer and compaction of the underlying layers (Speyerer et al. 2016).

Fig. 8. Depth-dependent density profile as in Vasavada et al. (2012). The strong nonlinearity with depth can be seen in the top few centimeters of regolith where density increases rapidly. The profile was evaluated at 29 irregularly spaced cells (cell centers in red).

The thermal conductivity is modeled to be depth- and temperature-dependent and is also adopted from Vasavada et al. (2012). They used the same exponential gradient as for the density parameter since density and conductivity are possibly physically related, that is, solid conduction increases with density:

(35)

where Λ_top = 0.0006 W m⁻¹K^-1, Λ_bottom = 0.007 W m⁻¹K^-1and X = 2.7, which describes the ratio between radiative and solid (phonon) conduction. Solid conduction is dominant in the lower layer, whereas thermal radiation dominates in the upper layer at temperatures above 350 K (Vasavada et al. 1999). The H-parameter is again chosen to be the background value of 6 cm.

The specific heat capacity is adopted from Hemingway et al. (1973) which was derived through a fit to Apollo core sample data and is modeled to be temperature-dependent,

(36)

4.3. Thermal radiation from Earth

The determination of thermal radiation from Earth to the Moon would entail a second thermal model running simultaneously to the one for the Moon. To avoid the computational overhead we assume Earth has a constant temperature at each zonal segment independent of the time of day. Following the results published in Trenberth & Stepaniak (2003), where the authors determined the average outgoing long-wave radiation (OLR) of Earth based on the years 1985–1989, we use the average OLR from March 2000 through March 2018 (personal communication Kevin Trenberth and Yongxin Zhang, both at UCAR: University Corporation for Atmospheric Research). For each of the 18 segments we calculate the mean OLR and assign it to the entire segment (see Fig. 9). In our lunar case study, Eq. (24) now reads:

(37)

Fig. 9. Left: average OLR value assigned to each of the 18 ten-degree zonal segments. We note the different average OLR near the north and south poles of Earth and the maximum radiation at the tropics. In this example no thermal radiation from segments south of 70°S would be received at the lunar location due to obstruction by the horizon. Right: average zonal OLR (black line, courtesy of Kevin Trenberth and Yongxin Zhang) and the average for each ten-degree zonal segment used in this study (red crosses).

We note that the distance d between the Moon and Earth is given in meters and the outgoing long-wave radiation of segment k on Earth, , in W m⁻².

5. Data sets

5.1. LRO–LOLA

To precisely model the illumination states near the poles we derived north and south polar LOLA DTMs each spanning 400 × 400 km and with 20 m pixel⁻¹ resolution (see Fig. 10 and Gläser et al. 2018). We then derived horizon maps and illumination states for the central 50 × 50 km subsets of the DTMs, defined as our region of interest (ROI; Figs. 10c,d), while taking into account possible obstructions of sunlight by the near- and far-field topography (Gläser et al. 2014). Here we included information from the DE421 ephemeris (Folkner et al. 2009) and applied SPICE routines (Acton 1996) to derive the exact position and size of the solar disk with respect to the lunar surface (Gläser et al. 2014). To demonstrate the accuracy of the LOLA DTM and the implemented illumination method we simulated illumination for our ROI at the lunar south pole (Fig. 11b). The time was chosen to be identical to the time of image acquisition of the LRO Wide Angle Camera (WAC) image M172692029ME shown in Fig. 11a. Besides the differences in brightness due to variable albedos which are not accounted for in our simulation, we find that the correspondence of illuminated and shadowed pixels between the two images is excellent (visual inspection). We note that the resolution of the simulated image is five times higher than the WAC image. Comparisons with more highly resolved LRO Narrow Angle Camera (NAC) images of 50 cm pixel⁻¹ can be found in (Gläser et al. 2014).

Fig. 10. (a) 400 × 400 km north polar DTM with outlines of ROI. (b) 400 × 400 km south polar DTM with outlines of ROI. (c) 50 × 50 km ROI at the north pole. (d) 50 × 50 km ROI at the south pole. Maps are displayed in gnomonic map projection, color-coded by heights and are centered at the lunar poles which is marked as a black cross. The lunar polar LOLA DTMs have a resolution of 20 m pixel−1.

Fig. 11. 50 × 50 km region of interest centered at the lunar south pole. (a) Cutout of WAC image M172692029ME displayed at 100 m pixel−1 resolution. (b) Synthetically illuminated LOLA DTM at 20 m pixel−1 resolution using the image acquisition time of the WAC image shown in (a) as the input time for simulation. Maps are displayed in gnomonic map projection.

5.2. LRO-Diviner

We chose temperature maps recorded by Diviner (Paige et al. 2010a) as our benchmark for comparison with our simulated surface temperatures. The Diviner instrument consists of two co-boresighted telescopes with a total of nine 21-element thermopile detector arrays at the focal planes. Each array has a separate spectral filter which allows measurement of the reflected solar and emitted infrared radiation in nine spectral channels in the range 0.3 − 400 μm. Channels 1 and 2 measure the reflected solar radiation from the lunar surface and are not used in this study. Channels 3–5 are located near 8 μm to precisely identify the spectral location of the Christiansen feature. Channels 6–9 cover the thermally emitted radiation of the lunar surface from approximately 12.5–50 μm. From a 50 km orbit altitude the footprint of each pixel is 320 m along track and 160 m across track with a swath width of 3.4 km. To infer bolometric brightness temperatures of the lunar surface we integrated the measured brightness temperatures in Diviner’s infrared channels 3–9 as described in the supporting online material of Paige et al. (2010b):

(38)

where T_bri, i is the measured brightness temperature in channel i and l is a weight stemming from the ratio of the Planck function B(λ, T_bri) evaluated for the spectral limits of channel i to B(λ, T_bri) over all wavelengths:

(39)

For illustrative purposes, we calculated six orbits of Diviner data over the north pole from February 2, 2010 and made a mosaic (see Fig. 12). In general, the temperature measurements of Diviner decrease in accuracy at low temperatures due to offset uncertainty, while the decrease in accuracy at high temperatures is due to gain uncertainty (Paige et al. 2010a). Further, when sub-detector scale anisothermality occurs, the derived brightness temperatures differ in each channel (Vasavada et al. 2012). This nonhomogeneous temperature distribution within the field of view of each detector increases when large illumination angles enhance the influence of roughness, shadowing, and topography. In our study however we concentrate on the polar areas where the lowest lunar temperatures are found on crater floors next to almost constantly illuminated and hot crater rims. Also, the illumination angles are very shallow such that roughness and topography are the main drivers of illumination. Although the polar areas are challenging for measuring surface temperature due to the aforementioned reasons, we chose the Diviner data set as our benchmark since it is the best available data source.

Fig. 12. Mosaic of six orbits of Diviner data from February 2, 2010, near the lunar north pole. Map is displayed in gnomonic map projection and color-coded by temperature.

6. Preliminary results

In this section we present initial results gained while implementing, testing, and validating our model.

6.1. Optimal choice of depth and discretization

The skin depth , where κ = Λ/ρc is the thermal diffusivity and p is the lunar synodic rotation period (≈29.53 d), describes the depth at which the intensity of the thermal wave of the diurnal surface temperature cycle decays to 1/e (≈37%) (Titus & Cushing 2012). The diurnal skin depth of the Moon is on the order of ≈4 − 10 cm (Hayne et al. 2017). The annual heat wave penetrates about three to four times deeper than the diurnal heat wave, for example, ≈40 cm. Kuhrt & Giese (1989) suggest that the lower boundary be placed at 2.5 skin depths (≈100 cm) in order to include the effects of the annual heat wave. In fact, Apollo heat-flow experiments have shown annual temperature variations at depths up to 120 cm (Langseth et al. 1976), and for this reason we place the lower boundary in our study at a depth of 2 m to comfortably include all diurnal and annual effects. Here, the intensity of the annual heat wave will have roughly decayed to 1/e⁵ (≈1%). We use 29 cells for the discretization of the sub-surface, where the first ten cells have a depth of 5 mm, followed by five cells with a depth of 1 cm, another five cells with a depth of 2 cm, two cells with a depth of 5 cm, and another two cells with a depth of 10 cm. The last five cells have a depth of 30 cm. The fine and irregular spacing is the result of a trade-off between sufficiently sampling the nonlinear progression of the lunar thermophysical parameters, which vary drastically in the first few centimeters of the regolith, and computation time, which scales with the amount of cells used (see also Fig. 8 where the irregular spaced grid is shown in red). For illustration purposes, we calculated temperatures over a period of ten days sampled every two hours for the upper 200 cm of regolith with an initial temperature distribution of 100 K using different sampling intervals (Fig. 13). If we consider thirty cells with regular sampling (≈6.6 cm cell depths) the thermophysical parameters are not resolved properly and heavy smoothing occurs. The temperature profile using 600 cells (≈0.33 cm cell depths) properly samples the nonlinear region of the thermophysical parameters in the upper layers. In the lower layers however there is little variation in the parameters and temperature and the sampling rate is excessive. The last profile is the result of the irregular sampling with 29 cells, as described above. We used the composite Simpson’s rule to determine the area beneath our three temperature profiles shown in Fig. 13 (see Table 1). Adopting the 600-cell profile as our benchmark profile, the relative error of our 29-cell profile amounts to 0.19% and is considered negligible. Consequently, we adopted the 29-cell irregularly spaced grid for our study due to vast computational benefits.

Fig. 13. Temperatures evaluated over a ten-day period starting from a uniform 100 K temperature distribution using different discretizations. Temperature profiles with constant cell sizes are shown for the upper 5 cm for 30 (red) and 600 (blue) cells and one profile is shown using an irregularly spaced grid with 29 cells (black). We note that for all profiles we additionally show the corresponding surface temperatures, that is, at 0 cm depth. The differences between the irregularly spaced grid with 29 cells and 30 cells of constant size are significant, whereas compared to the results obtained with a grid using 600 cells of constant size they are negligible. The irregularly spaced grid was adopted in this study.

6.2. Impact of area of influence

Due to memory and computational limitations, scattering and radiation from nearby terrain cannot be taken into account from all visible terrain patches; for example, from a sunlit crater rim as far as 150 km away. To asses the size of our chosen area of influence, that is, the extent to which distance scattering and radiation from nearby terrain is to be considered, we carried out a test using a lower DTM resolution of 250 m pixel⁻¹ to ensure that we remain within our computational boundaries. As a benchmark we defined our area of influence to be within a window of 80 × 80 km and compared it to smaller window sizes of 2 × 2 km, 10 × 10 km, 20 × 20 km, 30 × 30 km, 40 × 40 km, and 60 × 60 km. For each window size we integrated temperatures over 2 years with a 12 h step size and compared the final temperature map (January 1, 2010, at midnight (UTC)) to our benchmark map which was compiled using the 80 × 80 km window (see Fig. 14).

Fig. 14. Differences in surface temperature for the lunar north pole on January 1, 2010, using different window sizes for which scattering and radiation are considered. A window size of 80  ×  80 km serves as a benchmark. Differences are presented for window sizes of (a) 2 × 2 km, (b) 10 × 10 km, (c) 20 × 20 km, (d) 30 × 30 km, (e) 40 × 40 km, and (f) 60 × 60 km. Maps are displayed in gnomonic map projection and are color-coded by temperature.

It can be seen that considering scattering and radiation only within a 2 × 2 km window may not be sufficient for some regions (Fig. 14a). About 40% of the surface temperatures differ more than 1 K and 15% differ by more than 10 K compared to our benchmark map. On average each pixel differs by ≈4.3 K (Fig. 15). The difference maps for window sizes 10 × 10 km and 20 × 20 km do show a significant improvement at terrain features related in size, for example, inside the two larger craters; Hinshelwood and Whipple (Figs. 14b,c). Here, self heating due to radiation and reflectance from one crater wall to the opposite site is now more appropriately accounted for. For smaller terrain features, for example craters below ≈5 − 10 km radius, for which self heating from inside the crater is fully accounted for at these window sizes, the remaining differences must all come from features beyond the chosen window size. This effect is prominent in the crater field on Peary’s crater floor shown in the lower right quadrant in Figs. 14b,c. The average temperature differences per pixel are 2.5 K and 0.8 K for the 10 × 10 km and 20 × 20 km window sizes, respectively (Fig. 15). Increasing the window size further to 30 × 30 km and 40 × 40 km covers all self-heating processes within the craters shown in the scene and only 1.13% of the surface temperature differences in the 40 × 40 km difference map exceed 2 K (Fig. 14e). On average, the temperature differences per pixel decrease to 0.3 K and 0.1 K for the 30 × 30 km and 40 × 40 km windows, respectively (Fig. 15). However, some of the small craters located on Peary’s crater floor still show large temperature differences to our benchmark map. The final map (Fig. 14f) shows surface temperature differences to our benchmark considering a 60 × 60 km window. Although the area covered by the chosen window is 1.8 times smaller than our benchmark, 99.9% of it can be considered identical (< 1 K) and only two values exceed the 2 K difference with the maximum difference being 2.2 K. The average difference value per pixel is 0.016 K (Fig. 15). Based on our benchmark test and limitations in memory we prefer a 60 × 60 km window for our study. However, considering scattering and radiation from terrain within a 60 × 60 km window and using a DTM at a resolution of 20 m pixel⁻¹ (see Sect. 5) to derive temperatures for a ROI of 50 × 50 km we would need several terabytes of memory to handle the visibility maps. Such memory requirements are far beyond our current capabilities. We decided to keep the size of the region of interest at 50 × 50 km and were therefore obliged to lower the DTM resolution to 180 m pixel⁻¹ in order to stay within our computational capabilities (24 GB GPU memory). For illustrative purposes, we calculated scattering and radiation within our chosen 60 × 60 km window and 180 m pixel⁻¹ resolution for three locations, which are one pixel apart from each other (see Fig. 16). It can be seen that besides current illumination and temperature states, irradiance through scattered sunlight and thermal radiation also varies greatly with location. The first location is on a crater wall superposed on Peary’s crater rim. Here reflected sunlight and thermal radiation from the opposite crater wall but also from Peary’s crater rim have a substantial influence on the temperature at the respective location. The second location is just one pixel away (180 m in negative Y-direction) and is situated precisely on the crater rim where relatively little reflected light and thermal radiation from within the crater and from Peary’s rim can be received. In total however, more terrain is visible than in the first location. The third location (180 m in negative Y-direction) is just outside the small crater and receives radiation from an even larger set of features than the first two locations. Nevertheless, the first location can receive two times more radiation from nearby terrain than the second and seven times the amount received by the third location.

Fig. 15. Mean absolute error (L1-norm) for each map presented in Fig. 14 over the size of the window areas. This reveals how fast the errors are decreasing with increasing window sizes. Here the errors of the six window sizes used in our benchmark and the benchmark itself are shown (black crosses). Data points are evaluated at window sizes (left to right) 2 × 2 km, 10 × 10 km, 20 × 20 km, 30 × 30 km, 40 × 40 km, 60 × 60 km, and 80 × 80 km.

Fig. 16. Visible terrain calculated for three different locations within our chosen 60 × 60 km window. The maps are displayed in gnomonic map projection and color-coded by a dimensionless value describing the magnitude of influence that a pixel has on the center of the window (blue cross). Although the locations are only one pixel apart from each other there is a vast difference in magnitudes and visible terrain.

6.3. Convergence

The convergence of the model was tested in both space and time. We chose to test four different DTM resolutions placed around our adopted DTM resolution of 180 m pixel⁻¹; for example, 50 m pixel⁻¹, 100 m pixel⁻¹, 200 m pixel⁻¹, and 400 m pixel⁻¹. In order to find a suitable time-step for this study we similarly tested four different time-steps: 3 h, 6 h, 12 h, and 24 h. As a benchmark, the most highly resolved test case of 50 m pixel⁻¹ DTM resolution combined with a 3 h integration time-step was chosen, to which we compared all other cases. The test period was set to 1 year and the test area to a 10 × 10 km region centered at the south pole. For simplification, we only consider direct illumination and neglect scattering, thermal re-radiation, and radiation from Earth. All test cases started from the same initial temperature distribution and had the same depth profiles as described in Sect. 6.1. In total, 16 test cases were evaluated and compared to our benchmark; see Fig. 17a. It can be noted that the resulting temperatures only differ slightly with changing time resolution since the Moon is a slow rotator and all presented time-steps model illumination conditions sufficiently. For instance, there is only a 10% difference in the mean absolute errors between the results for a 24 h and 3 h time-step (at a fixed DTM resolution of 200 m) but there is over 1000% difference in mean absolute errors between a resolution of 50 m and 400 m (at a fixed time-step of 12 h). Consequently temperatures are far more sensitive to changing DTM resolutions (Figs. 17a,b). For our chosen resolution of 180 m pixel⁻¹ an average error of 3.7 K pixel⁻¹ is found (Fig. 17b) and in some rare cases differences can reach up to about ±100 K (evaluated at 200 m pixel⁻¹ and 12 h). Further increasing DTM resolution would significantly improve the solution, especially at resolutions higher than 100 m pixel⁻¹ as can be seen in Fig. 17b. Whereas the window size from which scattering is considered influences mostly temperatures within craters and PSRs (compare Fig. 14), the DTM resolution has a distinct and relatively large effect on temperatures at terrain edges; at crater rims for example. This is due to the smoothing effect of edges with decreasing resolution which results in changing incidence angles and therefore temperature. The discretization of time and space can be selected based on Figs. 15 and 17 and available computational capabilities. For this study however, which aims at presenting the fundamentals of the model independent of a scientific question at hand, we adopted, as a compromise between computational effort and accuracy, a time-step of 12 h and a DTM resolution of 180 m pixel⁻¹. Based on the scientific goal, it is crucial for future studies to assess the window size and DTM resolution with respect to the introduced errors presented in Sects. 6.2 and 6.3.

Fig. 17. (a) Mean absolute errors (L1-norm) of 16 test cases (black crosses) calculated using different resolution in time and space. The data point at 50 m DTM resolution and 3 h integration time-step serves as a benchmark; cf. lower left cross. The black dashed line shows the location of the profile displayed in (b) and the black circle indicates our chosen DTM resolution of 180 m pixel−1. The mean absolute errors are color-coded. (b) Profile showing the mean absolute error with increasing resolutions at a fixed time-step of 12 h. Our adopted DTM resolution is indicated as a black circle.

6.4. Equilibration

In order to derive temperatures for the 50 × 50 km ROIs centered at each pole, we illuminated our LOLA DTMs for a 18.6-year period (lunar precessional cycle) with a 12 h time-step while deriving temperatures for the surface and sub-surface. As described before we used a LOLA DTM with 180 m pixel⁻¹ resolution and a 60 × 60 km window to model scattering and radiation from nearby terrain. We started from a homogeneous temperature distribution for each profile from 80 K at the surface with a 1 K decrease in temperature for each node; for example, the last node (29th node) in the sub-surface has a temperature of 50 K. After integrating over the 18.6-year period the resulting temperature distribution was then used as an initial temperature and the process was repeated. By comparing corresponding profiles of the two resulting maps for each pole we could conclude that surface temperatures had equilibrated. The lowest sub-surface nodes (at 1.85 m depth) however were not yet equilibrated and on average still changed by 4 K. For our study, where we only anticipate presenting surface temperatures, we can consider the solution as equilibrated and for future studies involving temperatures in the sub-surface, for example depth-to-ice investigations, further iterations are needed.

6.5. Temperature maps for the lunar poles

Figure 18 shows preliminary temperature maps for the ROIs at the lunar poles for January 1, 2010, at midnight (UTC). We note that both poles are illuminated at this time and the Sun is shining on the lunar near-side (negative Y-direction at the north pole and vice versa at the south pole). While at the north pole Peary’s crater floor is mainly in shadow, its crater rim is fully illuminated and reaches temperatures of ≈300 K (Figs. 18a,c). Although many smaller craters are not illuminated at this particular time or are even never illuminated at all, for example they are located in PSRs like the craters on Hinshelwood’s crater floor, one can still see their outlines in the temperature map due to reflected sunlight and re-radiation of warmer terrain and subsequent heating. At the south pole Shackleton’s near-side crater rim and the ridge towards the de Gerlache crater are strongly illuminated and also reach temperatures of ≈300 K (Figs. 18b,d). The relatively high temperatures of ≈70 K within Shackleton itself are remarkable since previous studies have identified Shackleton as a PSR (Mazarico et al. 2011; Gläser et al. 2014). Heat transport towards the inside of Shackleton seems to be quite effective.

Fig. 18. Simulated illumination and temperature for the ROIs centered at the lunar poles on January 1, 2010. (a) Illumination of the north polar ROI. (b) Illumination of the south polar ROI. (c) Temperature of the north polar ROI. (d) Temperature of the south polar ROI. Maps (a–d) are displayed in gnomonic map projection. Maps (c) and (d) are color-coded by temperature.

To quantify the influence of direct and scattered sunlight as well as re-radiation from nearby terrain and Earth shine on the total irradiance we present, for illustrative purposes, the components coming from each single source on January 1, 2010, at midnight (UTC) for the south polar ROI (Fig. 19). At this particular time, direct insolation reaches ≈500 W m⁻² which is roughly only ≈36% of the solar constant which is received at the sub-solar point, for example at the lunar equatorial regions (Fig. 19b). This is due to the large incidence angles at the poles which scale the solar irradiance with the cosine of it; see Eq. (7). Scattered sunlight and thermal radiation coming from nearby terrain are mostly acting on terrain that is not illuminated (Figs. 19c,d). Here mainly opposite sites of illuminated and hence warm crater rims experience irradiance through (multiple) scattering and re-radiation. Although only a very thin part of Shackleton’s inner crater rim is illuminated at this time, there is still a nonzero heat flux in the interior, showing the strength of our multiple scattering approach; see Sect. 3.6. All but very few pixels (≈0.1%) experience some nonzero heat flux stemming from scattering and radiation from nearby terrain. Some of those locations that do not receive any additional heat flux coincide with the most illuminated locations near the lunar south pole which are located on Shackleton rim and on the ridge towards the de Gerlache crater; cf. Gläser et al. (2014). These areas are all relatively flat and are situated at high altitudes where no additional radiation from surrounding terrain can be received. The other few locations with zero heat flux are double-shadowed craters within PSRs for which no heat flux was found even after multiple scattering. We note that in our example thermal radiation from nearby terrain is approximately four times higher than irradiance from scattering. Scattering of sunlight and thermal radiation from Earth play a minor role in terms of magnitude but they usually irradiate surface patches completely independently of those illuminated directly or indirectly by the Sun, for instance when Earth is in opposition to the Sun as seen from the Moon. Furthermore, Earth can reach much higher elevations (≈6.68°) than the Sun (≈1.54°) meaning that it can potentially directly illuminate areas within PSRs. In our example however, Earth and the Sun appear very close to each other in the lunar sky (just one day after Earth is in conjunction with the Sun) and hence the directions of incoming radiations are similar and so are the irradiated surface patches (Figs. 19b and e,f). We note that in this example a significantly larger part of the lunar surface is irradiated by Earth shine than by direct sunlight.

Fig. 19. Single components of the total irradiance (b–f); (a) is shown for context. (a) Topography of the north polar ROI. (b) Direct sunlight. (c) Scattered sunlight from terrain. (d) Thermal radiation from terrain. (e) Scattered sunlight from Earth. (f) Thermal radiation from Earth. Map (a) is color-coded by height and maps (b–f) are color-coded by irradiance. Maps (a–f) are displayed in gnomonic map projection.

6.6. Comparison to Diviner

In order to validate our temperatures derived with the presented model we compare them to surface temperatures measured by Diviner. For the comparison we chose the Diviner measurement taken at February 2, 2010, at 01:20 (UTC) near the lunar north pole during orbit 2763 (also part of the mosaic shown in Fig. 12). At the south pole we chose the Diviner measurement taken on February 2, 2010, at 00:30 (UTC) from the previous orbit 2762. We then ran our simulation to the exact same times as when the Diviner measurements were taken. The resulting maps were cropped and resampled to match Diviner’s footprint and resolution of 250 m pixel⁻¹. After additionally accounting for an offset to the footprints of the Diviner measurement in the xy-direction of −0.3 pixel, −0.9 pixel for the north and −1.9 pixel, −0.3 pixel for the south polar DTM we derived difference maps (Figs. 20d,h). For the most part, the difference maps are flat when coinciding with areas where temperatures do not change drastically. However, at light–shadow boundaries the differences are more pronounced. In Sect. 6.3 we found that the used DTM resolution has the biggest effect on temperatures modeled at terrain edges. Since light–shadow boundaries often coincide with abrupt changes in topography such as crater rims, this is an indication that at least some of the remaining temperature differences can be attributed to the used DTM resolution of 180 m pixel⁻¹ (compare Figs. 20a,c,d where the largest temperature differences correlate with light–shadow boundaries at changing topography). Furthermore, the few larger differences of about ±100 K can be explained by the used DTM resolution as was found in Sect. 6.3. Overall there are less light–shadow boundaries at the south polar scene compared to the rougher surface morphology at the north pole resulting in smaller differences between the model and Diviner. Nevertheless, a correlation between the two maps at the north pole of 94.4% is achieved (Fig. 21a). At the south pole the correlation reaches 98.0% which emphasizes the great similarity between our modeled temperatures and the actual temperature measurements of Diviner (Fig. 21b).

Fig. 20. Our model temperatures compared to temperatures derived from Diviner data. (a–d) North pole and (e–h) South pole. (a),(e) DTM for which temperatures were modeled. (b),(f) Model temperatures at exact times of Diviner measurements. (c),(g) Diviner measurements taken on February 2, 2010, in (a) from orbit 2763 at 01:20 am and in (b) from orbit 2762 at 00:30 am. (d),(h) The difference between the model temperatures and the Diviner measurements. We note that no interpolation was used. Maps are displayed in polar gnomonic map projection and are color-coded by height (a),(e) and temperature (b–d),(f–h). For context the footprints of the Diviner measurements are indicated as black lines in (a),(b) and (e),(f).

7. Discussion

We present a method to calculate surface and sub-surface temperatures and have successfully derived temperature maps of the lunar polar areas. We analyzed the impact that the chosen window size from which scattering is considered, the disretization of the depth profile, and the used DTM and time resolution have on the resulting temperature and assessed the respective errors. Furthermore, the high-resolution LOLA DTMs used in this study allowed us to determine temperatures for terrain at a significantly increased resolution compared to a similar study carried out by Paige et al. (2010b). We also used updated thermophysical parameters of the lunar soil from Vasavada et al. (2012). Diviner represents the benchmark for lunar temperatures and hence our model temperatures could be validated in an absolute sense. The differences between our model and Diviner are shown in Fig. 20 and reveal a very good correspondence. Indeed, the two examples given correlate by ≈94.4% (north pole case) and ≈98.0% (south pole case; Figs. 21a,b). When contrasting Figs. 21a,b with similar Figs. S4-A and S4-B given in the supporting online material of the study by Paige et al. (2010b), we find great similarity. Nevertheless, we note that our model seems to simulate lower temperatures more accurately. Similar to Paige et al. (2010b), there is some evidence at temperatures between 50 and 100 K for a cluster of warmer Diviner measurements than what our model predicts. However, contrary to the findings of these latter authors, this is only observed for a faint cluster and only at the north pole case where the majority of the modeled temperatures in this range still correlate at a ratio of ≈1 : 1. This is a strong indication that the applied model can reproduce the actual lunar thermal behavior over a large temperature range. Despite the promising results, there are several remaining error sources for both the model and the derivation of Diviner measurements, which should be addressed in future studies. In this study, for instance, we used a window size of 60 × 60 km from which scattering was considered. Here we were able to show that remaining errors can be neglected (Fig. 15). However, our chosen DTM resolution, which was constrained by computational capabilities, still introduces large mean absolute errors of 3.7 K pixel⁻¹. Future studies need to address this issue and further increase the DTM resolution. As was also demonstrated in Sect. 6.4 the sub-surface temperatures are not yet equilibrated and need further iterations for studies interested in, for example, depth-to-ice investigations. The incorporation of variable albedos (Lemelin et al. 2016) and thermophysical parameters, for example through a variable H-parameter (Hayne et al. 2017), could further improve the model. Additionally, the typical integration time for a Diviner measurement is 128 ms during which LRO travels ≈200 m (velocity 1.66 km). That smearing effect (Williams et al. 2016) is also not taken into consideration in our study and might add to some of the resulting temperature differences seen in Fig. a,b. Additionally, we only used ≈1% of the data points for our plot as compared to their study (Paige et al. 2010b) and more data points need to be evaluated and compared in the future.

Fig. 21. Scatter density plot with 2 K step-size of Figs. 20c,g over the corresponding terrain in Figs. 20b,f to visualize the good overall agreement of the two maps. The derived Pearson’s correlation coefficients reach (a) 94.4% for the comparison of the model and Diviner measurement at the north pole and (b) 98.0% at the south pole.

8. Conclusion

Our model produces temperatures that are in accordance with Diviner measurements. All major effects affecting temperature were modeled and only minor improvements to the model itself can be made in the future which are expected to alter the results only slightly. The model parameters used, however, can have a significant effect on the resulting temperatures. For instance, in future studies we intend to improve our model by using variable albedo and thermophysical parameters as well as DTMs at higher resolutions. Also the introduction of a beaming parameter (thermal equivalent of the opposition surge in the visible) will have an additional impact on the modeled temperatures (Harris et al. 2002).

Acknowledgments

P. Gläser was funded by a Grant of the German Research Foundation (GL 865/2-1). We also gratefully acknowledge the support of NVIDIA Corporation with the donation of the Quadro P6000 GPU used for this research. The data used in Eq. (37) and Fig. 9 was kindly provided by Kevin Trenberth and Yongxin Zhang from the University Corporation for Atmospheric Research. Last, we wish to thank the LOLA, LROC and Diviner Science Teams for releasing such wonderful data products.

References

All Tables

All Figures

Fig. 1. Illustration of the model domain Ω (for d = 3) and boundary Γ, which is further decomposed into the part coinciding with the surface of the rocky object, Γs, the lateral boundary Γl, and the bottom boundary Γb. While there are generally no restrictions on the choice of the model domain, in this work we construct Ω by vertical extrusion of parts of the surface of the rocky object with a constant depth. Therefore, the surface topography is also illustrated on the bottom boundary Γb.

In the text

	Fig. 2. Geometric relations between infinitesimal areas dA around the locations x and y. The normal vectors are denoted n(x) and n(y), respectively. The graphic shows how the surface at position y is illuminated by the Sun at an incidence angle ι⊙ and reflects light at an emission angle η which is received at location x at a distance d and incidence angle ι.
In the text

	Fig. 3. Illustration of the geometrical setting between two cells K, L ∈ 𝒯. Depicted are the geometrical quantities required to describe the flux given in Eq. (18) across the face .
In the text

Fig. 4. Illustration of the computational grid and the splitting of the domain Ω into a set of vertical sub-domains Ωi, for which two exemplary columns Ωi and Ωj with discretization are depicted. We note that Ωi are still 3D objects, while the actual computational domains are strictly 1D, extruded by the measure . Furthermore, the discrete area of influence ℛi ⊂ Γs around Ωi is illustrated, within which the scattering of sunlight and thermal radiation of neighboring surface segments , j ≠ i, is considered.

In the text

	Fig. 5. 18 ten-degree zonal segments of the disk of the solar-system object along with the terminator and horizon as seen from the observer location . In this example no sunlit segments below 20°S are visible to the observer due to obstruction by the horizon.
In the text

	Fig. 6. Solar limb darkening. This is the phenomenon by which the Sun darkens from its center to its limb. The actual visible solar disk is defined by Rvisible, the apparent radius of the Sun depending on the distance of the observer. All points closer than the plane defined by the apparent radius are visible to the observer and are observed under emission angles ψ, ψ ∈ [0, 90].
In the text

Fig. 7. Solar disk discretized in five concentric rings from which we determine the actual visible solar disk ratio while considering solar limb darkening and obstruction by the horizon. The radius and area of each ring are denoted Rj and Pj, respectively. The points Cj are placed halfway between two successive radii except for the C1 which is placed right in the center of the disk. In this example only the brightest, most inner ring of the Sun is completely visible. All other rings are only partially visible due to obstruction by the horizon.

In the text

	Fig. 8. Depth-dependent density profile as in Vasavada et al. (2012). The strong nonlinearity with depth can be seen in the top few centimeters of regolith where density increases rapidly. The profile was evaluated at 29 irregularly spaced cells (cell centers in red).
In the text

Fig. 9. Left: average OLR value assigned to each of the 18 ten-degree zonal segments. We note the different average OLR near the north and south poles of Earth and the maximum radiation at the tropics. In this example no thermal radiation from segments south of 70°S would be received at the lunar location due to obstruction by the horizon. Right: average zonal OLR (black line, courtesy of Kevin Trenberth and Yongxin Zhang) and the average for each ten-degree zonal segment used in this study (red crosses).

In the text

	Fig. 10. (a) 400 × 400 km north polar DTM with outlines of ROI. (b) 400 × 400 km south polar DTM with outlines of ROI. (c) 50 × 50 km ROI at the north pole. (d) 50 × 50 km ROI at the south pole. Maps are displayed in gnomonic map projection, color-coded by heights and are centered at the lunar poles which is marked as a black cross. The lunar polar LOLA DTMs have a resolution of 20 m pixel−1.
In the text

	Fig. 11. 50 × 50 km region of interest centered at the lunar south pole. (a) Cutout of WAC image M172692029ME displayed at 100 m pixel−1 resolution. (b) Synthetically illuminated LOLA DTM at 20 m pixel−1 resolution using the image acquisition time of the WAC image shown in (a) as the input time for simulation. Maps are displayed in gnomonic map projection.
In the text

	Fig. 12. Mosaic of six orbits of Diviner data from February 2, 2010, near the lunar north pole. Map is displayed in gnomonic map projection and color-coded by temperature.
In the text

Fig. 13. Temperatures evaluated over a ten-day period starting from a uniform 100 K temperature distribution using different discretizations. Temperature profiles with constant cell sizes are shown for the upper 5 cm for 30 (red) and 600 (blue) cells and one profile is shown using an irregularly spaced grid with 29 cells (black). We note that for all profiles we additionally show the corresponding surface temperatures, that is, at 0 cm depth. The differences between the irregularly spaced grid with 29 cells and 30 cells of constant size are significant, whereas compared to the results obtained with a grid using 600 cells of constant size they are negligible. The irregularly spaced grid was adopted in this study.

In the text

Fig. 14. Differences in surface temperature for the lunar north pole on January 1, 2010, using different window sizes for which scattering and radiation are considered. A window size of 80  ×  80 km serves as a benchmark. Differences are presented for window sizes of (a) 2 × 2 km, (b) 10 × 10 km, (c) 20 × 20 km, (d) 30 × 30 km, (e) 40 × 40 km, and (f) 60 × 60 km. Maps are displayed in gnomonic map projection and are color-coded by temperature.

In the text

Fig. 15. Mean absolute error (L1-norm) for each map presented in Fig. 14 over the size of the window areas. This reveals how fast the errors are decreasing with increasing window sizes. Here the errors of the six window sizes used in our benchmark and the benchmark itself are shown (black crosses). Data points are evaluated at window sizes (left to right) 2 × 2 km, 10 × 10 km, 20 × 20 km, 30 × 30 km, 40 × 40 km, 60 × 60 km, and 80 × 80 km.

In the text

Fig. 16. Visible terrain calculated for three different locations within our chosen 60 × 60 km window. The maps are displayed in gnomonic map projection and color-coded by a dimensionless value describing the magnitude of influence that a pixel has on the center of the window (blue cross). Although the locations are only one pixel apart from each other there is a vast difference in magnitudes and visible terrain.

In the text

Fig. 17. (a) Mean absolute errors (L1-norm) of 16 test cases (black crosses) calculated using different resolution in time and space. The data point at 50 m DTM resolution and 3 h integration time-step serves as a benchmark; cf. lower left cross. The black dashed line shows the location of the profile displayed in (b) and the black circle indicates our chosen DTM resolution of 180 m pixel−1. The mean absolute errors are color-coded. (b) Profile showing the mean absolute error with increasing resolutions at a fixed time-step of 12 h. Our adopted DTM resolution is indicated as a black circle.

In the text

	Fig. 18. Simulated illumination and temperature for the ROIs centered at the lunar poles on January 1, 2010. (a) Illumination of the north polar ROI. (b) Illumination of the south polar ROI. (c) Temperature of the north polar ROI. (d) Temperature of the south polar ROI. Maps (a–d) are displayed in gnomonic map projection. Maps (c) and (d) are color-coded by temperature.
In the text

Fig. 19. Single components of the total irradiance (b–f); (a) is shown for context. (a) Topography of the north polar ROI. (b) Direct sunlight. (c) Scattered sunlight from terrain. (d) Thermal radiation from terrain. (e) Scattered sunlight from Earth. (f) Thermal radiation from Earth. Map (a) is color-coded by height and maps (b–f) are color-coded by irradiance. Maps (a–f) are displayed in gnomonic map projection.

In the text

Fig. 20. Our model temperatures compared to temperatures derived from Diviner data. (a–d) North pole and (e–h) South pole. (a),(e) DTM for which temperatures were modeled. (b),(f) Model temperatures at exact times of Diviner measurements. (c),(g) Diviner measurements taken on February 2, 2010, in (a) from orbit 2763 at 01:20 am and in (b) from orbit 2762 at 00:30 am. (d),(h) The difference between the model temperatures and the Diviner measurements. We note that no interpolation was used. Maps are displayed in polar gnomonic map projection and are color-coded by height (a),(e) and temperature (b–d),(f–h). For context the footprints of the Diviner measurements are indicated as black lines in (a),(b) and (e),(f).

In the text

	Fig. 21. Scatter density plot with 2 K step-size of Figs. 20c,g over the corresponding terrain in Figs. 20b,f to visualize the good overall agreement of the two maps. The derived Pearson’s correlation coefficients reach (a) 94.4% for the comparison of the model and Diviner measurement at the north pole and (b) 98.0% at the south pole.
In the text