Issue 
A&A
Volume 627, July 2019



Article Number  A129  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201935514  
Published online  12 July 2019 
Modeling nearsurface temperatures of airless bodies with application to the Moon
^{1}
Technische Universität Berlin, Institute of Planetary Geodesy, 10623 Berlin, Germany
email: philipp.glaeser@tuberlin.de
^{2}
University of Stuttgart, Dept. of Hydromechanics and Modelling of Hydrosystems, 70569 Stuttgart, Germany
email: dennis.glaeser@iws.unistuttgart.de
Received:
19
March
2019
Accepted:
5
June
2019
In this study we present a model to determine surface and subsurface temperatures of airless bodies in the solar system. To precisely model direct sunlight we incorporated the solar limb darkening effect of the solar disk. Scattered sunlight and thermal reradiation from nearby planets is also considered in our model. We further consider multiple scattering of reflected sunlight and thermal reradiation on the modeled object itself. The finite volume method is applied to solve the model for which we present full derivations for the governing equations that control scattering and heat diffusion into the subsurface. We assessed errors stemming from the chosen discretization of the depth profile, the window size from which scattering is considered, as well as from the chosen integration stepsize and the spatial resolution of the Digital Terrain Model (DTM). Exemplarily, we determine surface and subsurface (2 m depth) temperatures for the lunar polar areas. Topography of the lunar poles is modeled by measurements of the Lunar Orbiter Laser Altimeter (LOLA). We integrated temperatures over a 18.6year time frame using 180 m pixel^{−1} LOLA DTMs of the poles, a 60 × 60 km window, and a 12 h integration timestep. The resulting preliminary temperature maps for the lunar poles are presented. Further, we show that our model agrees with temperatures obtained by the Diviner lunar radiometer experiment.
Key words: Moon / methods: numerical / conduction / radiation mechanisms: thermal
© 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 Earthdays 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 surfacebound 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 subsolar 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 nearpolar 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 temperaturedependent 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 insitu 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 cratershapes while also considering thermal reradiation. Later, analytical and numerical thermal models were solved for idealized cratershapes near the poles of Mercury and the Moon to investigate surface temperatures and nearsurface ice taking into account multiple scattering effects (Paige et al. 1992; Ingersoll et al. 1992; Salvail & Fanale 1994). A more sophisticated study using a twolayer model with depth and temperaturedependent thermophysical properties taking into account multiple scattering was carried out by Vasavada et al. (1999) for flat surfaces and surfaces within bowlshaped and flatfloored polar impact craters. Thermal models were also used to study the effect of roughness on surface temperatures, coldtrap 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 dwarfplanets) 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 modelcalculated temperatures and compared them to measured Diviner bolometric temperatures using a 500 mscale 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 nearsurface regolith. Bauch et al. (2014) used 100 m pixel^{−1} 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 nonenergyconserving 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 Earthshine or no radiogenic heat source in the subsurface has been considered. Furthermore, some of the presented models require timesteps 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 highresolution lunar topography derived from the Lunar Orbiter Laser Altimeter (LOLA) data (Smith et al. 2017) in order to determine surface and subsurface 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 reradiation from nearby terrain, heat conduction in the subsurface, radiation into space, a constant radiogenic heat source in the subsurface, 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 nearsurface 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 wellknown heat equation:
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 threedimensional, 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:
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 solarsystem 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 solarsystem objects. In order to allow for a compact notation, let us introduce the set 𝒞 = {ϕ_{1}, ϕ_{2}, …, ϕ_{n}} of solarsystem 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 solarsystem object ϕ ∈ 𝒞, at a location x ∈ ∂ϕ on its surface ∂ϕ, on the basis of the Stefan–Boltzmann law:
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
where S_{0} 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 timedependent distance of ϕ to the Sun in astronomical units (AU). Making use of (6) and (7) we now decompose the heat flux q such that
where the terms and refer to scattered sunlight and thermal radiation that reaches ℳ from the considered solarsystem objects ϕ_{k}. The former is modeled by
Here, δ = δ(x,y) ∈ {0, 1} describes whether or not the terrain at location y ∈ ∂ϕ on the solarsystem 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
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:
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
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 finitevolume 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 finitevolume 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. Finitevolume discretization
In this section, we outline the derivation of the discrete formulation of the model problem. We employ a finitevolume 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
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,
At the core of finite volume schemes lies the construction of the discrete flux approximation,
for which we are using a classical twopoint 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)
We note that we consider control volumewise 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 cellwise 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
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:
Finally, approximating the time derivative in Eq. (15) by a finite difference, the discrete formulation reads
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 finitevolume 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:
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. Onedimensional approximation
We assume that the nearsurface 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, ddimensional setting (with d ∈ {1, 2, 3}). Starting from a 3D domain Ω ⊂ ℝ^{3} (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 ith subdomain 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 ith subdomain. Equations (2)–(5) are then solved separately in each subdomain , 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 subdomains Ω_{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 subdomains Ω_{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 subdomains within the area of influence. We note that in this discrete setting Eq. (11), evaluated for subdomain i, reads
where the subscripts i and j refer to quantities evaluated on and , respectively, and where we have introduced the discrete analogon (around subdomain i) of the area of influence, ℛ_{i} (see Sect. 2.2). denotes the temperature at the surface of subdomain j, that is, on .
3.4. Approximation of flux density from solarsystem objects
In order to determine reflected sunlight coming from nearby solarsystem objects, for example from Earth, we divide the visible celestial disk in 18 tendegree 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 solarsystem object to the Sun d_{ϕ, sun} (we note that if the segment is not lit or not visible):
Fig. 5. 18 tendegree zonal segments of the disk of the solarsystem 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 tendegree zonal segments to model heat flux from solar system objects:
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 solarsystem 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 R_{visible}, 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):
The emission angle ψ can be calculated from the angles θ and γ:
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_{2} = 0.93 and ν_{2} = −0.23 (Cox 2000). Setting w_{0} = 1 − u_{2} − ν_{2}, w_{1} = u_{2} and w_{2} = ν_{2} leads to
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:
with W_{1} = −0.47 and W_{2} = −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,
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 R_{j} and P_{j}, respectively. The points C_{j} are placed halfway between two successive radii except for the C_{1} 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:
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,
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,
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,
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 finitevolume 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 subsurface 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^{−2}, 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^{−8} W m^{−2} K^{−4} and in our study we set the solar irradiance to S_{0} = 1366 W m^{−2} (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 depthdependent as described in Vasavada et al. (2012),
where ρ_{top} = 1300 kg m^{−3} and ρ_{bottom} = 1800 kg m^{−3}. The discovery of the nonlinear dependency of density with depth when evaluating Diviner temperature maps led to the introduction of the Hparameter (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 Hparameter 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 micrometeor impacting and the resulting overturn within the top layer and compaction of the underlying layers (Speyerer et al. 2016).
Fig. 8. Depthdependent 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 temperaturedependent 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:
where Λ_{top} = 0.0006 W m^{−1}K^{1}, Λ_{bottom} = 0.007 W m^{−1}K^{1}and 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 Hparameter 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 temperaturedependent,
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 longwave 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:
Fig. 9. Left: average OLR value assigned to each of the 18 tendegree 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 tendegree 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 longwave radiation of segment k on Earth, , in W m^{−2}.
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^{−1} 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 farfield 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^{−1} 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, colorcoded 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. LRODiviner
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 coboresighted telescopes with a total of nine 21element 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):
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:
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 subdetector 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 colorcoded 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 heatflow 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^{5} (≈1%). We use 29 cells for the discretization of the subsurface, 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 tradeoff 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 600cell profile as our benchmark profile, the relative error of our 29cell profile amounts to 0.19% and is considered negligible. Consequently, we adopted the 29cell irregularly spaced grid for our study due to vast computational benefits.
Fig. 13. Temperatures evaluated over a tenday 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. 
The three different profiles that were tested: two profiles with cells of constant size and one with an irregularly spaced grid.
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^{−1} 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 colorcoded 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 selfheating 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^{−1} (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^{−1} 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^{−1} 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 Ydirection) 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 Ydirection) 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 (L1norm) 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 colorcoded 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^{−1}; for example, 50 m pixel^{−1}, 100 m pixel^{−1}, 200 m pixel^{−1}, and 400 m pixel^{−1}. In order to find a suitable timestep for this study we similarly tested four different timesteps: 3 h, 6 h, 12 h, and 24 h. As a benchmark, the most highly resolved test case of 50 m pixel^{−1} DTM resolution combined with a 3 h integration timestep 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 reradiation, 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 timesteps 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 timestep (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 timestep of 12 h). Consequently temperatures are far more sensitive to changing DTM resolutions (Figs. 17a,b). For our chosen resolution of 180 m pixel^{−1} an average error of 3.7 K pixel^{−1} is found (Fig. 17b) and in some rare cases differences can reach up to about ±100 K (evaluated at 200 m pixel^{−1} and 12 h). Further increasing DTM resolution would significantly improve the solution, especially at resolutions higher than 100 m pixel^{−1} 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 timestep of 12 h and a DTM resolution of 180 m pixel^{−1}. 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 (L1norm) 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 timestep 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 colorcoded. (b) Profile showing the mean absolute error with increasing resolutions at a fixed timestep 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.6year period (lunar precessional cycle) with a 12 h timestep while deriving temperatures for the surface and subsurface. As described before we used a LOLA DTM with 180 m pixel^{−1} 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 subsurface has a temperature of 50 K. After integrating over the 18.6year 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 subsurface 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 subsurface, for example depthtoice 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 nearside (negative Ydirection 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 reradiation of warmer terrain and subsequent heating. At the south pole Shackleton’s nearside 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 colorcoded by temperature. 
To quantify the influence of direct and scattered sunlight as well as reradiation 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^{−2} which is roughly only ≈36% of the solar constant which is received at the subsolar 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 reradiation. 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 doubleshadowed 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 colorcoded by height and maps (b–f) are colorcoded 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^{−1}. After additionally accounting for an offset to the footprints of the Diviner measurement in the xydirection 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^{−1} (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 colorcoded 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 subsurface 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 highresolution 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. S4A and S4B 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^{−1}. Future studies need to address this issue and further increase the DTM resolution. As was also demonstrated in Sect. 6.4 the subsurface temperatures are not yet equilibrated and need further iterations for studies interested in, for example, depthtoice investigations. The incorporation of variable albedos (Lemelin et al. 2016) and thermophysical parameters, for example through a variable Hparameter (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 stepsize 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/21). 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
 Acton, C. H. 1996, Planet. Space Sci., 44, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Aharonson, O., & Schorghofer, N. 2006, J. Geophys. Res., 111,E11007 [NASA ADS] [CrossRef] [Google Scholar]
 Araki, H., Tazawa, S., Noda, H., et al. 2009, Science, 323, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Bauch, K. E., Hiesinger, H., Helbert, J., Robinson, M. S., & Scholten, F. 2014, Planet. Space Sci., 101, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Buhl, D., Welch, W. J., & Rea, D. G. 1968, J. Geophys. Res. (1896–1977), 73, 7593 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, A. N. 2000, Allen’s Astrophysical Quantities [Google Scholar]
 Davidsson, B. J. R., & Rickman, H. 2014, Icarus, 243, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Delbo, M., Mueller, M., Emery, J. P., Rozitis, B., & Capria, M. T. 2015, Asteroids IV, 107, 205 [Google Scholar]
 Droniou, J. 2014, Math. Models Methods Appl. Sci., 24, 1575 [CrossRef] [MathSciNet] [Google Scholar]
 Folkner, W. M., Williams, J. G., & Boggs, D. H. 2009, Interplanet. Netw. Prog. Rep., 178, 1 [Google Scholar]
 Gaster, B., Howes, L., Kaeli, D. R., Mistry, P., & Schaa, D. 2011, Heterogeneous Computing with OpenCL, 1st edn. (San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.) [Google Scholar]
 Gläser, P., Scholten, F., De Rosa, D., et al. 2014, Icarus, 243, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Gläser, P., Oberst, J., Neumann, G. A., et al. 2018, Planet. Space Sci., 162, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Harris, A. W., & Lagerros, J. S. V. 2002, in Asteroids in the Thermal Infrared, eds. W. F. Bottke, Jr., A. Cellino, P. Paolicchi, & R. P. Binzel, 205 [Google Scholar]
 Hayne, P. O., & Aharonson, O. 2015, J. Geophys. Res. (Planets), 120, 1567 [NASA ADS] [CrossRef] [Google Scholar]
 Hayne, P. O., Greenhagen, B. T., Foote, M. C., et al. 2010, Science, 330, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Hayne, P. O., Ghent, R., Bandfield, J. L., et al. 2013, Lunar Planet. Sci. Conf., 44, 3003 [NASA ADS] [Google Scholar]
 Hayne, P. O., Bandfield, J. L., Vasavada, A. R., et al. 2016, in New Views of the Moon 2, LPI Contrib., 1911, 6065 [NASA ADS] [Google Scholar]
 Hayne, P. O., Bandfield, J. L., Siegler, M. A., et al. 2017, J. Geophys. Res. (Planets), 122, 2371 [Google Scholar]
 Hemingway, B. S., Robie, R. A., & Wilson, W. H. 1973, Lunar Planet. Sci. Conf. Proc., 4, 2481 [NASA ADS] [Google Scholar]
 Hu, X., Shi, X., Sierks, H., et al. 2017, MNRAS, 469, S295 [CrossRef] [Google Scholar]
 Ingersoll, A. P., Svitek, T., & Murray, B. C. 1992, Icarus, 100, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Jaeger, J. C., & Harper, A. 1950, Nature, 166, 1026 [NASA ADS] [CrossRef] [Google Scholar]
 Jamsa, S., Peltoniemi, J. I., & Lumme, K. 1993, A&A, 271, 319 [NASA ADS] [Google Scholar]
 Keihm, S. J., Peters, K., Langseth, M. G., & Chute, J. L. 1973, Earth Planet. Sci. Lett., 19, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Kuhrt, E., & Giese, B. 1989, Icarus, 81, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Lagerros, J. S. V. 1997, A&A, 325, 1226 [NASA ADS] [Google Scholar]
 Langseth, M. G., Keihm, S. J., & Peters, K. 1976, in Lunar and Planetary Science Conference Proccedings, ed. D. C. Kinsler, 7, 3143 [NASA ADS] [Google Scholar]
 Lemelin, M., Lucey, P. G., Neumann, G. A., et al. 2016, Icarus, 273, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Linsky, J. L. 1966, Icarus, 5, 606 [NASA ADS] [CrossRef] [Google Scholar]
 Mazarico, E., Neumann, G. A., Smith, D. E., Zuber, M. T., & Torrence, M. H. 2011, Icarus, 211, 1066 [NASA ADS] [CrossRef] [Google Scholar]
 Mitchell, D. L., & de Pater, I. 1994, Icarus, 110, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Paige, D. A., Wood, S. E., & Vasavada, A. R. 1992, Science, 258, 643 [NASA ADS] [CrossRef] [Google Scholar]
 Paige, D. A., Foote, M. C., Greenhagen, B. T., et al. 2010a, Space Sci. Rev., 150, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Paige, D. A., Siegler, M. A., Zhang, J. A., et al. 2010b, Science, 330, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Pettit, E., & Nicholson, S. B. 1930, ApJ, 71, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Piqueux, S., Hayne, P. O., Elder, C. M., et al. 2016, Lunar Planet. Sci. Conf., 7, 1762 [Google Scholar]
 Rubanenko, L., & Aharonson, O. 2017, Icarus, 296, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Salvail, J. R., & Fanale, F. P. 1994, Icarus, 111, 441 [NASA ADS] [CrossRef] [Google Scholar]
 Scholten, F., Oberst, J., Matz, K. D., et al. 2012, J. Geophys. Res. (Planets), 117, E00H17 [Google Scholar]
 Siegler, M., Paige, D., Williams, J.P., & Bills, B. 2015, Icarus, 255, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, D. E., Zuber, M. T., Neumann, G. A., et al. 2017, Icarus, 283, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Speyerer, E. J., Povilaitis, R. Z., Robinson, M. S., Thomas, P. C., & Wagner, R. V. 2016, Nature, 538, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Stubbs, T. J., & Wang, Y. 2012, Icarus, 217, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Titus, T. N., & Cushing, G. E. 2012, LPI Contributions, 1673, 95 [NASA ADS] [Google Scholar]
 Trenberth, K. E., & Stepaniak, D. P. 2003, J. Clim., 16, 3706 [NASA ADS] [CrossRef] [Google Scholar]
 Vasavada, A. R., Paige, D. A., & Wood, S. E. 1999, Icarus, 141, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Vasavada, A. R., Bandfield, J. L., Greenhagen, B. T., et al. 2012, J. Geophys. Res. (Planets), 117, E00H18 [CrossRef] [Google Scholar]
 Williams, J.P., SeftonNash, E., & Paige, D. A. 2016, Icarus, 273, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, J.P., Paige, D. A., Greenhagen, B. T., & SeftonNash, E. 2017, Icarus, 283, 300 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
The three different profiles that were tested: two profiles with cells of constant size and one with an irregularly spaced grid.
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 subdomains Ω_{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 tendegree zonal segments of the disk of the solarsystem 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 R_{visible}, 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 R_{j} and P_{j}, respectively. The points C_{j} are placed halfway between two successive radii except for the C_{1} 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. Depthdependent 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 tendegree 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 tendegree 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, colorcoded 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 colorcoded by temperature. 

In the text 
Fig. 13. Temperatures evaluated over a tenday 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 colorcoded by temperature. 

In the text 
Fig. 15. Mean absolute error (L1norm) 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 colorcoded 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 (L1norm) 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 timestep 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 colorcoded. (b) Profile showing the mean absolute error with increasing resolutions at a fixed timestep 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 colorcoded 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 colorcoded by height and maps (b–f) are colorcoded 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 colorcoded 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 stepsize 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.