Application of Finite Difference Method for determining lunar regolith diurnal temperature distribution

. This study was performed in order to verify viability of using finite difference method and proposed simple astrometrical model for modelling heat transfer in lunar regolith. The concept was examined by developing FD model of heat flow for upper 0,9 m of lunar regolith, and comparing obtained results with in situ measurements provided by Apollo 15 and 17 heat flow experiments. The model was based on FDM approximation of Fourier’s law for one dimensional transient heat flow. Both constant and temperature-dependent thermophysical properties of lunar regolith were obtained from in situ measurements. Thermal boundary conditions were assumed on in situ measurements and on remote sensing based analytical model. In order to approximate Sun's position at lunar sky, simple analytical astrometric model of lunar rotation was developed. Matlab 2012a was used to conduct the calculations. Stable solutions were obtained for latitudes between 0 and 80°. Satisfactory agreement between Apollo 15 and 17 in situ measurements and FDM modelling was observed. A conclusion was reached, that both FDM and proposed astrometrical model are to be successfully applied for modelling heat transfer in lunar regolith.


New perspectives for space exploration
An increased global interest in space exploration is to be observed for several last years. Large national space agencies as well as privately owned companies conduct research and develop technologies to be used for both exploratory and commercial space applications. Noticeable part of those efforts are focused on re-establishing human presence on the Moon [1-6].

Finite difference method application in Building Physics
Finite difference method (FDM) is a numerical method for solving differential equations by approximating them with difference equations [13,14]. An example of this method's numerous applications are to be found in Building Physics where it is used for solving conductive heat transport problems. In his book, Pogorzelski presents detailed derivation of FDM equation from Fourier's Law of Heat Conduction, where numerous examples of heat flow problems in civil engineering are elaborated [13]. Kisilewicz presents successful FDM modelling of conductive and radiative heat transfer combination in transparent building barrier [15]. A vast comparative analysis of external walls' thermal properties was conducted by Kaczmarzyk, where thermal responses for short temperature changes for multiple combinations of structural and thermal insulating materials were determined by FDM modelling [16].
A general, one-dimensional, transient heat flow FDM formula for temperature at iinterface for time step t is given by following equation [15]: 2 ≤ m ≤ 6 (4) Equation (4). FDM numerical stability condition. This condition has to be met for every modelled material layer, otherwise a numerical instability will occur [13,14]. This restriction. Is a noticeable inconvenience in FDM modelling.

Selected physical properties of lunar regolith
According to data provided by manned Apollo lunar missions, lunar soil has very low thermal conductivity, that changes significantly with temperature, bulk density and applied pressure [12,[17][18][19][20][21]23]. Specific heat of lunar regolith is also highly temperature dependent [22].
Bulk density of lunar regolith is relatively low at the surface, but increases logarithmically with depth from 1,300 kg/m 3 to about 1,900 kg/m 3 [12,21].
Since early 70's several models of lunar regolith thermal conductivity temperature dependence have been proposed [12,[17][18][19][20][21]23]. There is a general agreement between thermal conductivities calculated using all these models for dense regolith, i.e. for depths greater several centimetre's. However, there are significant differences between thermal conductivities calculated for loose regolith, so called fluff, that exists down to two centimetres beneath Moon's surface [12,23].
Nonetheless, it is to be observed, that the thin, near-surface fluff layer has very low, and highly temperature dependent thermal conductivity.
It is to be pointed out, that the very process of boring into lunar regolith disturbed and compacted the material that surrounded later inserted bore stems [12,23]. Moreover, there were only four boreholes made for temperature measurements during the Apollo programme [24], thus it is very difficult, if possible, to provide accurate thermal model of lunar soil basing only on these data.
For the purpose of this study, we chosen regolith thermal conductivity model introduced by Vasavada et al. as it reproduces near-equatorial lunar surface temperatures with high accuracy [21].

Thermal boundary conditions
FDM requires input data, on which to start and conduct further calculation's. In case of this model these data include surface temperature, fixed temperature at a specified boundary depth, and initial temperature distribution in so called zeroth time step. As the latter is usually unknown, almost any values may be assumed. As the calculations proceeds, heat flows through analysed strata and eventually final results should be obtained, regardless of initially assumed temperature distribution [13][14][15].
Basing on in situ measurements, constant temperature 252 K and 254 K were assumed at the boundary depths of 0.9 m and 0.67 m for Apollo 15 and Apollo 17 landing site respectively [12].
Temperature of lunar surface was calculated according to Hurley et al. as a function of time, selenographic latitude and altitude of the Sun [25]. In order to determine the latter, without using ephemerides or additional software, a simplified analytical model of the Mon's rotation was devised.

Analytical model of lunar rotation
The Earth and its natural satellite -the Moon constitute a complex, dynamical system. Hence, the determination of the accurate ephemerides for selected celestial bodies on the Moon's sky is a non-trivial problem, and complicated mathematical formulae are needed to be used. However, it is possible to build a simple model, that will allow to estimate Sun's position in the lunar sky with precision of 1-2°, what is fully acceptable for engineering purposes.
The plane of the lunar orbit is inclined to the ecliptic by about 5.15°, and the rotational axis of the Moon is inclined to its orbital axis by 6.69° [26]. In 1722 Jacques Cassini discovered, that the rotational axis of the Moon precesses with the same rate as its orbital plane, but is shifted by 180° in phase [27]. Therefore, the angle between the ecliptic and the lunar equator remains almost constant (1.54°). For the simplicity of the model it may be assumed, that the equatorial plane of the Moon is coplanar with the ecliptic plane, i.e. the rotational axis of the Moon is perpendicular to the ecliptic plane. Moreover, it may also be assumed, that the Moon is located in the barycenter of the Earth-Moon system. The radius of lunar orbit is significantly smaller than the radius of the Earth's orbit, hence all effects caused by the trigonometric parallax of the Sun on the Moon's firmament may be neglected in this model. To summarise, proposed model of the Earth-Moon system could be reduced to the Moon located in the Earth's orbit with the rotational axis perpendicular to the ecliptic plane. The following paragraph describes how it is possible to calculate the azimuth and the altitude of the Sun on the Moon sky using the proposed model.
Let it be assumed, that the observer is located on 0° selenographic latitude and on the selenographic longitude φ>0 o (for φ>0 o the problem is symmetrical with respect to the Moon's equator). In this simplified model the rotational axis of the Moon is perpendicular to its orbital plane, hence for a given selenographic longitude the ecliptic great circle intersects the horizon at an angle 90 o -φ.
At the moment t 0 the Sun is located in the maximum elevation equal to 90 o -φ (the arc PS). After the time Δt=t'-t 0 the Sun will be located in P', and will draw the arc PP'=ε. At the time t' solar altitude h is equal to the arc P'G. The elevation arc is a part of the great circle ZP'G, and can be used for the azimuth (A) estimation, which is equal to arc SG (here it is assumed that the azimuth is counted clockwise from the south direction). Also, the arc P'W is equal to 90 o -φ as a section of the arc PW. Now using the law of sines for the parallactic triangle P'GW we have: This directly leads to the estimation of the altitude , as: h= arcsin(cosφ cos ε) (6) From the law of sines for the parallactic triangle ZPP' it is: Using the formula (3) the azimuth A may be described as: A= arcsin( sin ε cos h ) With formulae (2) & (4) available, it is possible to describe the Sun's position in the horizontal coordinate system for a given arc ε.
For engineering purposes it is not necessary to calculate the Sun's position for the real time t, as only general description of the altitude as the function of time and the selenographic latitude is of interest.
It is also possible to perform a simple estimation of the length of the arc ε. The Moon's rotation period T K ≈27.322 Eearth's days, and lunar year in our model is identical with Earth's year P K ≈365.256 days. If Δt is given in days we can simply estimate ε using following formula: ε= Δt (1/T K − 1/ P K ) (9) as the combination of orbital and rotational motion of the Moon. Equation (9) may be also simply described as: where T S is so called the Moon synodic month i.e. the period of lunar phases, of mean value 29.53 days.
In order to verify accuracy of described analytical model, its results were compared with precisely calculated Sun's positions, determined by Stellarium astronomical software. As seen in table 1. there are just minor differences between calculated positions of the Sun at Moon's sky and actual values. We consider such accuracy as fully satisfactory for lunar engineering purposes. This simple astrometric AL model may be easily incorporated into any software, and provide reliable astrometric AL data for energy calculations at the Moon's surface.

FDM numerical stability
At any, relatively low depth beneath the Moon's surface, regolith temperature results from thermal balance between periodically changing surface solar radiation and internal heat production by radioactive decay. These two heat sources reach an equilibrium at some depth, and produce constant temperature, that slowly increases with depth [12].
At any time, subsurface regolith temperatures cannot exceed some temperature range, determined by two extreme, steady state temperature distributions. This temperature range will be the widest at the surface, narrowing down with depth. In order to determine this temperature range as a function of depth, two steady state temperature distributions were calculated, each for the highest and for the lowest surface temperature, i.e. for lunar noon and pre-dawn respectively.
As lunar regolith thermal properties change with temperature, an iterative approach was used here.
Having those temperature ranges, it was necessary to calculate extreme values of temperature dependent, regolith thermal diffusivity for each discrete layer. Due to time step thermal diffusivity dependence and FDM stability condition (figures 2. to 4.), FDM will not yield stable solution if thermal diffusivity of any single discrete layer changes more than three times during a diurnal cycle. It was calculated, that maximum possible thermal diffusivity changes for single discrete layer is least at the equator (1.90 times) and greatest at 80° (2.67 times). As regolith thermal diffusivity changes less than 3.0 times, it was concluded, that FDM may be successfully applied for heat flow modelling in lunar regolith.
On the other hand, possible thermal diffusivity changes between all the discrete layers (7.42 times at the equator and up to 9.24 times at 80°) do not allow to apply singlethickness discrete layers. Thus, variable layer thicknesses were used, between 1 mm and 2 mm.
Such theoretically possible temperature ranges were determined for latitudes between 0 and 80°, with 10° intervals. Lunar regolith steady state temperature distributions for lunar noon and pre-dawn at 20° latitude. Theoretically possible regolith temperatures at given depths are limited by these two, extreme distributions. As actual regolith temperatures are contained in a field enclosed by two graph lines, meeting FDM stability condition for that limited region guarantees stable solution for FDM model at specified selenographic latitude.

Discretization
According to FDM numerical stability condition (Fig. 2) analysed upper 0.8 m of lunar regolith was divided into elementary layers.
For analysed Apollo 15 and 17 landing sites (Rima Hadley 25.0°N, 3.0°E and Tauruslittrow 20.0°N, 31.0°E respectively) elementary layer thickness was 0.001 m for uppermost 25 mm of the regolith and 0.0015 m for the remaining strata, resulting in exactly two minute time step. Such high density mesh was necessary for accurate determination of high temperature gradients, expected to appear in near-surface depths.

Model details
In order to conduct considerable amount of necessary calculations and to automate the process, a suitable program was created using MATLAB 2012A.
Thermal properties of discrete layers were updated every time step, basing on temperatures from previous time steps. Very short duration of assumed time step (2 minutes) guaranteed high accuracy and near-continuous character of the computation. The computation proceeded until periodically repetitive temperature distribution was reached.

Results and discussion
Results of conducted computation were saved as a table representation of regolith temperature as a function of depth and time, that passed since sunrise. Selected data were depicted in the charts below.
Using the same thermal conductivity model [21] for both Apollo 15 and 17 sites resulted in obtaining very similar overall character of temperature distribution for both analysed sites. There are some differences resulting from ~5° difference in the sites' latitude, but the shape of temperature curves remains the same. That is why above charts present only results for Apollo 15, while results for Apollo 17 are only elaborated According to Apollo in situ measurements [20] surface temperatures for Apollo 15 site varies between 84 K and 380 K, and between 88 K and 385 K for Apollo 17 site. Results of conducted FDM simulation give extreme temperature range 81-379 K for Apollo 15 and 83-384 K for Apollo 17. This small difference fits within lunar applied global surface temperature model. It is to be observed, that thermal inertia of the fluff is so low that the daytime surface temperature remains almost in thermal equilibrium with absorbed solar radiation, It implies, that surface temperature model results agree with observations, as long as Sun's position is calculated accurately. This fact proofs reliability of the astrometric model introduced in this paper. Nonetheless, lunar surface temperature model [25] used in this analysis does not consider any physical properties of local surface, nor its slope. This should be changed, for more detailed, future analyses.   3. to 6. presents, that large temperature amplitudes observed at the surface, rapidly attenuates within the first few centimetres of the regolith, remaining relatively stable at depth larger than 50 cm. This observation is confirmed by Malla and Brown, who obtained similar results using fourth-order Runge-Kutta technique of numerical integration [12].
Unfortunately, there were no precise in situ temperature measurements taken down to 5cm, due to relatively large thermocouples span along used bore stems [20,23,28]. Apollo 15 in situ measurements shown, that nearly 300 K surface temperature amplitude is almost completely attenuated below 50 cm [28]. In comparison, at the same location, our FDM model yields temperature amplitude 299 K at the surface and just 0.15 K at the depth of 50 cm.
On the other hand, not every set of compared values matches so well within few K tolerance. For example, between surface and a depth of 35 cm about 45 K and 40 K increase of mean temperature was registered at Rima Hadley and Taurus-Littrow respectively [20]. Our FDM model yielded lower values of This increase: 35 K and 37 K respectively. There is satisfactory agreement between our FDM model and Apollo 17 heat flow experiment results, what cannot be stated for Apollo 15 results. However, analysis conducted by [28] show, that temperature distribution at Rima Hadley was disturbed by thermal perturbations caused by nearby topographic features, and cannot be treated as globally representative data. All temperature distributions and thermal conductivity profiles presented above indicates, that lunar regolith, especially its uppermost, "fluff" layer exhibits extraordinary thermal insulating properties. Even peak values of regolith thermal conductivity (0.0076 Wm -1 K -1 ) are comparable with state-of-the-art man-made thermal insulators [29]. This exceptionally advantageous property of lunar regolith makes the material a perfect resource for thermal shields for future lunar buildings.
As noted by other authors [12,23] some differences between assumed lunar soil thermophysical properties and resulting modelled temperature distributions are totally acceptable, as lunar regolith remains relatively not well examined material to date.

Conclusions
This paper demonstrates, that FDM may be successfully applied for lunar regolith thermal analysis at least up to 80° latitude. Introduced simple analytical astrometric model approximates position of the Sun in lunar sky with precision that is fully satisfactory for lunar engineering purposes.
Results of performed FDM heat flow analysis are in general agreement with in situ experiments and simulations conducted using other methods. Presented FDM model may be Improved by applying internal planetary heat flow instead of using fixed temperature at specified depth.
It has been also concluded, that lunar surface temperature should be calculated using surface physical properties' and actual solar irradiance, instead of global, remote sensing based model.
Due to its excellent thermal insulating properties, lunar regolith may be used as thermal shield for future lunar engineering structures.
More extensive in situ experiments are required to properly investigate thermal properties of lunar regolith and to create more representative models of the material.