Modeling Thermal Emission under Lunar Surface Environmental Conditions

Thermal emission spectra can provide key insights into the composition and thermophysical properties of the regolith on the Moon and other airless bodies. However, under lunar surface environmental conditions, the uppermost millimeters of the regolith (from which thermal emission originates) cannot be characterized by a single temperature, leading to changes in spectral characteristics that should be accounted for in interpreting thermal emission measurements. Here, we develop and apply a Monte Carlo radiative transfer method to model thermal emission from particulate media with varying, nonisothermal subsurface temperature profiles. We model emission spectra for three major lunar mineral phases (pyroxene, olivine, and plagioclase), and investigate the effects of particle size and packing density. Modeled spectra are compared to lab measurements acquired under both ambient and simulated lunar conditions. We find that in some cases, the model provides useful constraints on the magnitude of the temperature profile established in a lab sample under lunar-like conditions, whereas in other cases, lab spectra are not well represented by the linear temperature profiles considered in this work. The model is generally successful at predicting changes in spectral contrast under lunar-like conditions, but less successful in accurately predicting shifts in the position of the Christiansen feature emissivity maximum; we illustrate and discuss the validity of the modeling approach for a range of different cases. Model results can also be used to quantify the depth within which observed thermal emission originates; this depth depends on composition and grain size, and ranges from ∼100 to 1000 μm for representative packing densities.


Introduction
Infrared emission from a planetary surface carries information about the composition and physical properties of the surface, as well as its temperature. Particulate regoliths on airless bodies, such as the Moon, are often highly insulating and, as a result, can sustain steep temperature gradients within the epiregolith (Mendell & Noble 2010), the near-surface boundary layer (∼1 mm deep) from which "surface" optical reflectance and thermal emission originate. Thus, remote sensing and laboratory measurements of particulate media in near-vacuum often measure thermal emission spectra that cannot be characterized by a single temperature, complicating the derivation of surface composition and thermophysical properties (such as particle size and porosity), which in turn are important for understanding the processes that have shaped a planetary surface. Epiregolith temperature gradients on airless bodies may also play a role in driving temperature-sensitive surface processes, such as volatile transport (Sarantos & Tsavachidis 2020. Several previous works have investigated the thermal emission properties of lunar regolith and regolith analogs under simulated lunar and asteroid conditions, through both laboratory measurements (e.g., Logan & Hunt 1970;Donaldson Hanna et al. 2017;Shirley & Glotch 2019;Breitenfeld et al. 2021;Donaldson Hanna et al. 2021) and numerical models (e.g., Henderson & Jakosky 1994;Millán et al. 2011). These studies have generally found a shift in characteristic spectral features (notably a shift of the Christiansen feature, the midinfrared emissivity maximum, to shorter wavelengths) and an increase in overall spectral contrast under lunar-like conditions. The degree to which these spectral features change is related to the magnitude of the thermal gradient within the sample and its albedo. However, there remains some disagreement between models as to the actual magnitude of near-surface thermal gradients on the Moon, Mercury, and other airless bodies (e.g., Henderson & Jakosky 1997;Hale & Hapke 2002).
In this work, we develop a Monte Carlo radiative transfer method to examine how epiregolith thermal gradients shape the thermal emission spectra of three major lunar mineral phases, pyroxene (enstatite), olivine (forsterite), and plagioclase (labradorite). Comparable previous models have focused on terrestrial mineral and rock phases, quartz (Henderson & Jakosky 1997;Millán et al. 2011) and basalt (Henderson & Jakosky 1997), that are not the most representative of the lunar crust. We apply our modeling approach to constrain the magnitude of the thermal gradients established under simulated (laboratory) lunar conditions, investigate the effects of particle size and packing density, and discuss implications for understanding the near-surface thermal environments of the Moon and other airless bodies.
The laboratory spectra that appear below were acquired in environment chambers capable of simulating lunar surface conditions, located at the Johns Hopkins Applied Physics Laboratory (Greenhagen et al. 2020), the University of Oxford (Thomas et al. 2012;Donaldson Hanna et al. 2021), and Stony Brook University (Shirley & Glotch 2019). Other laboratory facilities with these capabilities include Brown University (Donaldson Hanna et al. 2017) and the Planetary Spectroscopy Laboratory at DLR (Maturilli et al. 2008). Thermal gradients are typically created within samples by removing the atmospheric gases within the chamber (pressure <10 −4 mbar), heating the samples from below using cup heaters, and heating the samples from above using a solar-like lamp to brightness temperatures similar to those found on the Moon (∼290-400 K). To ensure that these thermal gradients are like those on the Moon, laboratory measurements of Apollo regolith samples under simulated lunar conditions have been compared to Diviner Lunar Radiometer observations of the sampling sites at which the Apollo soils were collected (Donaldson Hanna et al. 2017;Greenhagen et al. 2019). Laboratory measurements under ambient (isothermal) conditions are also possible, by filling chambers to pressures ∼1000 mbar and heating samples from below. To measure the thermal infrared emissivity spectra, chambers are connected to Fourier transform infrared (FTIR) spectrometers that have detectors and beamsplitters capable of measuring across the ∼5-50 μm spectral range.

Numerical Method
In this work, we simulate thermal emission spectra using a radiative transfer model in an approach that is similar to previous studies, with some differences. Henderson & Jakosky (1994, 1997) use a two-stream approximation that solves a pair of equations describing the intensities of upward and downward radiation in a plane-parallel medium. Similarly, Millán et al. (2011) use a discrete ordinate code that solves the radiative transfer equation for 4 to 16 discretized streams (i.e., directions of propagation). Alternatively, Ito et al. (2017) compute emissivity for isothermal media using analytical expressions derived from Conel (1969) and Hapke (1993Hapke ( , 1996. Our method, the Regolith Boundary Layer (ReBL) model adopts a Monte Carlo radiative transfer approach that models the scattering and attenuation of a large number of representative "energy bundles" propagating through the medium of interest. The main difference between this and the methods described above is that the Monte Carlo approach involves simulating the physical processes underlying radiative transfer (transport, scattering, attenuation) rather than explicitly solving the radiative transfer equation. In the limit of a large number of simulated "energy bundles," the Monte Carlo solution should yield similar or greater accuracy compared to analytical solutions.
Like other radiative transfer models, ReBL requires some prior knowledge of the scattering and absorption properties of the modeled medium. Here, we compute the scattering coefficient (μ s ), absorption coefficient (μ a ), and asymmetry parameter (g) over a range of wavelengths for a given grain composition, size, and packing density using a Mie scattering code (Mätzler 2002). We assume a Henyey-Greenstein phase function for each scattering event, such that the probability of a scattering angle θ is given by q p( ): where g is an asymmetry parameter (−1g1), such that g = 0 corresponds to isotropic scattering, while g > 0 indicates that forward scattering dominates and g < 0 indicates that backward scattering dominates.
Mie theory is strictly applicable only to isolated, spherical scatterers, which is seldom the case in a planetary regolith, but has nonetheless often been applied to model planetary surfaces (e.g., Henderson & Jakosky 1997;Millán et al. 2011) due to its relative computational simplicity. While we use the scattering and absorption coefficients (μ s , μ a ) here, it is also possible to work with the corresponding cross sections (σ s , σ a ) or dimensionless efficiencies (Q s , Q a ) instead. The scattering (or absorption) cross section, σ s(a) = Q s(a) A, and the corresponding coefficient μ s(a) = σ s(a) n, where A is the cross-sectional area of a scattering particle, and n is the number of particles per unit volume.
As shown schematically in Figure 1, ReBL models the epiregolith as a series of plane-parallel isothermal layers of equal thickness (though varying layer thicknesses are easily accommodated). The basic Monte Carlo radiative transfer algorithm (based on Modest 2003) is as follows: 1. The initial positions and directions of propagation of representative energy bundles are selected through random number draws, assuming uniform, isotropic emission within each isothermal layer of thickness Δz.

Each bundle is moved through a distance
where the extinction coefficient, μ e ≡ μ s + μ a , and R l is a random number between 0 and 1. (Note that although the simulated medium is onedimensional, the propagation of energy bundles is threedimensional.) The code then determines whether the bundle should be absorbed or scattered, or whether it has crossed a computational domain boundary. 3. If ω < R a (where ω ≡ μ s /μ e is the single scattering albedo and R a is a random number between 0 and 1), the bundle is absorbed. Otherwise, a scattering angle is sampled from the Henyey-Greenstein probability distribution function. 4. If the bundle has crossed the surface (the upper boundary of the computational domain), its contribution to the measured emission spectrum is recorded. If it crosses the lower boundary of the computational domain, propagation is terminated. (The domain is specified to be sufficiently deep that a negligible amount of energy from the deepest layer escapes the medium. Thus, we assume that if a bundle crosses the lower boundary, the probability that it will escape the medium is negligible.) 5. The output of the code is a set of "weights," w(λ, z) = N esc (λ, z)/N emit (λ, z), corresponding to the fraction of thermal energy at the wavelength emitted within the layer at depth z (i.e., between z and z + Δz) that escapes the medium. Here, N esc (λ, z) is the number of energy bundles that escape the medium, and N emit (λ, z) is the total number of energy bundles initialized (at wavelength within the layer at depth z).
ReBL can be applied in a number of different ways. Given subsurface temperature profiles, the computed weights w(λ, z) can be used to simulate emission spectra, as discussed in this work. Like other radiative transfer models, ReBL can also be combined with a one-dimensional conductive heat transfer model to compute planetary surface and subsurface temperatures (Prem et al. 2019).
The results presented here were obtained by tracking 10 5 energy bundles per wavelength per layer. (The accuracy of a Monte Carlo simulation increases, i.e., statistical noise decreases, as the number of simulated energy bundles increases. For this work, 10 5 bundles per wavelength per layer provided an acceptable balance between accuracy and computational cost.) The simulated medium is 1 mm deep and is discretized into 100 layers, each 10 μm thick (in order to resolve changes in temperature over this length scale). We note that it is also possible to adopt a nonlinear discretization scheme (e.g., if increased spatial resolution of the uppermost layers is desired); in this case, the number of "energy bundles" released within each layer should be scaled by layer thickness in order to maintain the same level of accuracy. Simulated spectra have a wavelength resolution of <0.1 μm (the resolution of available optical constant data).
Given a subsurface temperature profile T(z), the spectral intensity (W m −2 /μm) of thermal emission from the surface at wavelength is given by where n is the number of layers, w i (λ) is the weight associated with layer i at wavelength, Δz is layer thickness (assumed to be constant here), and B(λ, T i ) is the Planck function at wavelength and layer temperature T i . The summation term is the total flux emitted within the volume of a layer (Modest 2003) ) is the fraction of energy at wavelength isotropically emitted from a surface at the maximum modeled depth (z max ) that escapes the medium, computed according to the algorithm above, but considering energy bundles initialized at a surface (the lower boundary of the computational domain) rather than within a volume. The computational domain is sufficiently deep that l  w z , 0 max ( ) , but this term is included in calculations for completeness.
The brightness temperature spectrum thus obtained can be converted to an emission spectrum. Effective emissivity is defined as By this definition, ò eff = 1 at λ CF (the wavelength corresponding to the Christiansen feature). Under ambient (isothermal) conditions, "effective emissivity" is equivalent to "emissivity." Under anisothermal, lunar-like conditions, when the medium cannot be characterized by a single brightness temperature, the effective emissivity at a given wavelength is the ratio of energy emitted from the surface at that wavelength to energy emitted at the same wavelength by a blackbody with a brightness temperature equal to the brightness temperature of the Christiansen feature.

An Illustrative Isothermal Spectrum
To illustrate the use of ReBL, we simulate a medium composed of 28 μm diameter grains of enstatite, with a packing density f = 0.37 (corresponding to a regolith density of 1.1 g cm −3 , per Hayne et al. 2017). In order to calculate Mie scattering parameters (μ s , μ a , g), we used optical constants measured by Zeidler et al. (2015) for the three orthogonal crystallographic orientations of enstatite at 300 K. We do not presently account for the temperature dependence of optical constants, since the measurements of Zeidler et al. indicate that variations are modest over the range of temperatures considered here, also consistent with previous measurements of infrared reflectance by Hinrichs & Lucey (2002). Figure 2 shows isothermal (i.e., constant temperature versus depth) emission spectra constructed using the Mie scattering Figure 2. Isothermal spectra modeled using the Mie scattering properties of 28 μm diameter grains with a packing density of 0.37, and optical constants corresponding to three orthogonal crystallographic orientations of enstatite, dubbed x, y, and z (dotted/dashed lines). Solid lines indicate a simple average of x, y, and z emissivities, and an isothermal spectrum modeled using an average of x, y, and z scattering properties.
properties of each crystallographic orientation, and of a mixture of equal parts of each orientation (i.e., using averaged scattering properties). Since the relationship between scattering properties and emissivity is nonlinear, the latter differs from a simple average of the emissivities corresponding to the three crystallographic orientations. Applying this difference as a correction factor, we can then compute a best-fit mixture of orientations that approximates lab measurements acquired under ambient (isothermal) conditions such that e l e l e l e l where ε eff,lab (λ) is the measured effective emissivity, ε eff,model,x/y/z (λ) are the modeled emissivities for x/y/z crystallographic orientations, f x/y/z are the weights assigned to each orientation, and F is a correction factor defined as e l e l e e l º + + where ε eff,model,avg (λ) is the modeled emissivity using averaged scattering properties. Figure 3 compares model results to measurements acquired under ambient conditions in the Simulated Airless Body Emission Laboratory (SABEL; Greenhagen et al. 2020) for enstatite grains 25-32 μm in size, washed and dried to remove clinging fines. The modeled spectrum is based on the Mie scattering properties of 28 μm grains. In this case, the best-fit weights obtained using Equation (5) were f x = 0.21, f y = 0.43, and f z = 0.36. Weighted average scattering properties were obtained as follows: ( · · · · · · ) ( ) / Subscripts x/y/z are used to denote Mie scattering properties computed using optical constants corresponding to those respective crystallographic orientations.
It can be seen from Figure 3 that the modeled spectrum reproduces some aspects of the lab data well; the positions and relative contrast of spectral features are generally in good agreement. The most notable discrepancies are at shorter wavelengths, where modeled emissivity is significantly lower than measured. This unrealistically steep drop in emissivity at shorter wavelengths has also been noted in other models that apply Mie theory to model particulate media (e.g., Ito et al. 2017). This likely indicates that at these short wavelengths, approximation of the medium as composed of isolated, spherical particles becomes less valid. The model also considers only a single particle size rather than a size range, in addition to which optical constants drawn from the literature may not precisely match those of the laboratory sample. We also assume the same mixture of orthogonal crystallographic orientations for each modeled layer-this is another parameter that could be adjusted to obtain a better fit. In order to account for differences between modeled and measured spectra, we define an isothermal correction factor: where ε eff,lab and ε eff,model denote lab and modeled effective emissivity, respectively. F iso,λ is then applied as a scaling factor when modeling thermal emission under anisothermal conditions in subsequent sections.
It can be seen from Figure 3 that lab emissivities for enstatite at some longer wavelengths (>12 μm) are higher than at the Christiansen frequency, where emissivity is normalized to unity. Measurement conditions when these data were collected were optimized for the 5-12 μm range, and the high values at longer wavelengths may be a calibration artifact. Therefore, we focus on wavelengths <12 μm when discussing differences between ambient and anisothermal enstatite spectra below.

Isothermal Sensitivity Depths
The weights, w(λ, z), computed by ReBL can also be used to characterize a "sensitivity depth" for infrared emission measurements. For instance, Figure 4 shows the depth within which 50%, 90%, and 99% of measured thermal emission from a surface originates, as a function of wavelength, for the isothermal case described in Section 2.1. The sensitivity depth is qualitatively similar to the thermalization length (the mean distance traveled by a photon through a scattering medium before absorption), also shown in Figure 4, defined per Henderson & Jakosky (1997) as It should be noted that sensitivity depths depend on the subsurface temperature profile; if deeper layers are warmer than shallower layers, the sensitivity depth increases, particularly at more transparent (larger d t ) wavelengths. The depth resolution in Figure 4 is limited by the finite thickness of modeled planeparallel layers. However, layer thicknesses that are smaller than the grain size may not be physically meaningful, particularly at the uppermost surface, where regolith microstructure may be far from plane-parallel. In the remainder of the text, we use the term "isothermal sensitivity depth" to refer to the depth within which 99% of thermal emission from a medium originates under isothermal conditions (i.e., the "99% flux" line in Figure 4). The isothermal sensitivity depth is a function of wavelength.

Results and Discussion
In this section, we examine the relationship between epiregolith thermal gradients and spectral characteristics by modeling the infrared emissivity of three mineral phases of interest: pyroxene (enstatite), olivine (forsterite), and plagioclase (labradorite), with a variety of grain sizes and packing densities.
The simplest useful way to parameterize thermal gradients may be based on the magnitude of change in temperature (ΔT) and the depth (Δz) over which this occurs, assuming a linear increase in temperature from the surface to a constant temperature at depth. For each modeled medium, we consider several values of ΔT that span a broad range of possible temperature gradients. The range of Δz is informed by the maximum isothermal sensitivity depth (Section 2.2), which in turn is dependent on composition, particle size, and packing density. In the simulations below, the temperature at depth is set to be 400 K. However, it is worth noting that in anisothermal cases, the shape of the modeled emission spectrum is influenced not only by the shape of the temperature profile, but also by the magnitude of temperatures. This effect is also seen in lab data (Donaldson Hanna et al. 2017), but is not explored in detail here since it is a relatively small effect over the temperature ranges investigated.

Pyroxene (Enstatite)
We first consider the baseline case presented in Section 2 (28 μm enstatite grains, f = 0.37), for which the maximum isothermal sensitivity depth is 360 μm. We therefore model a set of nine temperature profiles with ΔT = 100 K, 75 K, and 50 K over depths of Δz = 360 μm, 180 μm, and 90 μm. Figure 5 shows a representative set of four modeled spectra, scaled by the correction factor defined in Equation (10) such that the modeled/measured isothermal spectra are identical. Similar to lab data, modeled spectra show an increase in overall spectral contrast as the epiregolith thermal gradient becomes steeper.
It can be seen from Figure 5 and Table 1 that there is a general trend of increasing shift of the Christiansen feature to shorter wavelengths with increasing thermal gradient, but the magnitude of this shift is generally smaller than seen in the lab data. While the modeled ΔCF for the 100 K/90 μm case agrees with lab data, the overall spectral contrast in this case is much higher than observed in the simulated lunar environment. The difficulty in matching the CF position is attributable to the fact that the discrepancies between modeled and measured spectra are most pronounced at shorter wavelengths (see Figure 3), and not fully compensated for by the application of a simple correction factor. To quantify changes in spectral contrast, we look at the local emissivity minima corresponding to two of the most clearly visible fundamental vibrational bands (also known as the reststrahlen bands), at wavelengths of ∼9 μm and ∼11.5 μm (designated RB1 and RB4, respectively, in Table 1, per Donaldson Hanna et al. 2012a). Based on this measure, we find that the lab simulated lunar environment (SLE) spectrum is best approximated by the case (shown in Figure 5) in which temperature varies by 100 K over 360 μm.
It should be noted that there are families of temperature profiles with different brightness temperatures but near-identical effective emissivity values, such as the 75 K/360 μm and 50 K/180 μm cases shown in Figure 6. In such situations, knowledge of brightness temperature can provide an additional constraint on the subsurface thermal gradient. . The depth within which 50%, 90%, and 99% of thermal emission from an isothermal medium originates, as a function of wavelength. Also shown, for comparison, is the analytically computed thermalization length (see Equation (11)).

Figure 5.
Modeled and measured thermal emission spectra for ∼28 μm enstatite grains. Modeled spectra are corrected by applying a scaling factor such that modeled and measured isothermal spectra match. Measured spectra were acquired under ambient (isothermal) and simulated lunar conditions. The model results shown correspond to four different anisothermal cases, in which temperature increases linearly by ΔT down to a depth Δz, below which temperature is constant.

Packing Density Effects
Variations in packing density (f, i.e., volume filling factor) can be modeled by linearly scaling scattering and absorption coefficients obtained using Mie theory. Since Mie theory assumes that each modeled grain scatters and absorbs radiation independently, the more densely packed the medium, the higher the extinction coefficient (i.e., the shorter the extinction length). In addition to the f = 0.37 case discussed above, we also modeled 28 μm enstatite grains with packing densities of f = 0.27 and 0.43, corresponding to regolith densities of 0.8 g cm −3 (per Henderson & Jakosky 1997) and 1.3 g cm −3 (per Hale & Hapke 2002), respectively. Since the packing density of the laboratory sample is not precisely known, the same ambient laboratory spectrum was used to obtain correction factors (per Equation (10)) in all cases.
Maximum isothermal sensitivity depths (Section 2.2) were found to be 500 μm for f = 0.27, and 300 μm for f = 0.43, compared to 360 μm for the baseline f = 0.37 case. Temperature gradients of ΔT = 100 K, 75 K, and 50 K over depths of Δz = 1.0, 0.5, and 0.25 times the maximum isothermal sensitivity depth resulted in near-identical sets of modeled emission spectra (compared to Figure 5 and Table 1) in all three cases. This finding reinforces the utility of the maximum isothermal sensitivity depth as a characteristic length Figure 6. Modeled brightness temperature and effective emissivity (uncorrected, and corrected per Equation (10)) for a medium of ∼28 μm enstatite grains, with two different temperature profiles. Despite different temperatures, both of these cases result in near-identical values of effective emissivity. Note. The change in spectral contrast under SLE conditions (indicated by the emissivity at reststrahlen bands RB1 and RB4) is best matched by a temperature change ΔT = 100 K over a depth Δz = 360 μm, although the modeled shift in the location of the Christiansen feature (ΔCF) is smaller than observed.   scale, and implies that epiregolith thickness increases with decreasing particle packing density. The general trend of shorter extinction lengths (and, thus, sensitivity depths) for more densely packed media should hold true even if more physically rigorous methods are used to compute scattering and absorption properties. However, we note that close-packing also tends to suppress forward scattering (e.g., Mishchenko 1994;Ito et al. 2018), which could lead to subtle differences in measured and modeled spectra that are not captured in this work.

Particle Size Effects
Figure 7 compares modeled and measured lab spectra for two additional size fractions of enstatite. Similar to the approach adopted in Section 2.1, the 32-45 μm size fraction is modeled using Mie scattering properties for 38 μm spheres, while the <25 μm size fraction is modeled using Mie scattering properties for 10 μm spheres.
The best-fit mixtures of crystallographic orientations were found to be f x = 0.22, f y = 0.36, and f z = 0.42 for the 38 μm case, and f x = 0.27, f y = 0.31, and f z = 0.42 for the 10 μm case. The resulting modeled isothermal spectra, before the application of a correction factor, are indicated by dashed lines in Figure 7. It can be seen that the measured and modeled (uncorrected) isothermal spectra agree more closely for larger grains (Figure 7(b)). This can be attributed to the tendency of Mie scattering models to overemphasize transparency features (such as that at ∼12.5 μm) that are more pronounced for smaller grain sizes (Ito et al. 2017). While the application of a correction factor to modeled anisothermal spectra can compensate for this to some extent, discrepancies in spectral shape between measured (light pink) and modeled (dark red) anisothermal spectra are also more noticeable for smaller grains (Figure 7(a)).
It is interesting to note that for all three enstatite grain sizes considered, measured RB1 (∼9 μm) and RB4 (∼11.5 μm) emissivity values are most closely matched by a modeled thermal gradient of ΔT = 100 K over a depth Δz = the maximum isothermal sensitivity depth = 480 μm for 38 μm grains, and 230 μm for 10 μm grains (see Figure 7). The measured location of the Christiansen feature remains challenging to model accurately.

Olivine (Forsterite)
We next apply the procedure described in Section 2.1 to model ambient and anisothermal olivine spectra, comparing to lab data for San Carlos olivine acquired at the University of Oxford. We first consider a sample size fraction of 45-75 μm, which we model using Mie scattering properties for 60 μm diameter spheres, with a packing density f = 0.37, and optical constants for San Carlos olivine at 300 K from Zeidler et al. (2015). It can be seen from Figure 8(a) that the resulting best-fit ambient spectrum shows less qualitative agreement with the measured spectrum than for enstatite (see Figure 3), particularly in the vicinity of the reststrahlen band region near 10.5 μm. The reason for this discrepancy appears to be related to the fact that the modeled emission spectra of the three orthogonal crystallographic orientations differ the most in this regionmaking it challenging to find to find a best-fit mixture of orientations that matches this and other regions of the spectrum equally well. Figure 8(b) shows the modeled isothermal sensitivity depth (see Section 2.2) and thermalization length for this medium. Due to the larger grain size, sensitivity depths and thermalization lengths are also larger for this medium compared to that in Figure 4. In both of these cases, the sensitivity depth (defined as the depth within which 99% of "surface" thermal emission originates) is ∼2.7-3 times the thermalization length.
Using the difference between modeled and measured isothermal spectra as a correction factor, we model a set of six anisothermal temperature profiles with ΔT = 150 K and 100 K, over depths of Δz = 680 μm, 340 μm, and 170 μm.  Table 2 shows the location of the Christiansen feature, and the emissivity at the reststrahlen bands near 10.5 μm and 19 μm (designated RB3 and RB5, per Donaldson Hanna et al. 2012a). Of these six profiles, the overall increase in spectral contrast observed in lab measurements under simulated lunar conditions is best approximated by a linear temperature gradient of ΔT = 150 K over a depth Δz = 680 μm (the maximum isothermal sensitivity depth), as shown in Figure 9(a).
In contrast with results for enstatite, there are some significant qualitative discrepancies in the trends seen in model results compared to lab measurements. In particular, model results show a shift of the Christiansen feature to longer wavelengths and a decrease in emissivity at ∼24 μm under anisothermal conditions, contrary to observed trends. This discrepancy in behavior near the emissivity peak is likely due to the poor agreement in the ∼10 μm region of model results obtained assuming Mie theory (Figure 8(a)). Discrepancies at longer wavelengths may be related to particle size. Lab measurements indicate that whereas smaller particles exhibit a decrease in emissivity in this region under lunar-like conditions, larger particles exhibit an increase in emissivity (Shirley & Glotch 2019). However, this behavior does not appear to be adequately represented by Mie scattering properties.

Particle Size Effects
Figure 10 compares modeled and measured spectra, as well as isothermal sensitivity depths, for smaller and larger size fractions of olivine: a <45 μm size fraction, modeled using Mie scattering properties for 20 μm spheres, and a 75-125 μm size fraction, modeled using Mie scattering properties for 100 μm spheres. Discrepancies between modeled and measured spectra are similar in nature to those observed for the Figure 10. Comparison of measured and modeled spectra for (a) smaller (∼20 μm) and (c) larger (∼100 μm) olivine grains; (b) and (d) indicate the corresponding isothermal sensitivity depths and thermalization lengths, based on Mie scattering properties for 20 μm and 100 μm grains, respectively. intermediate size fraction discussed above. It should also be noted that the 100 μm particles (the largest size modeled here) are significantly larger than the 10 μm thickness of modeled isothermal layers; in reality, sub-grain-scale temperature gradients are likely to be minimal due to solid state conduction. The maximum isothermal sensitivity depth in this case also approaches the depth of the modeled medium, such that in nonisothermal scenarios, thermal radiation originating from even deeper within the medium may be measurable. Figure 11 shows representative, modeled anisothermal spectra compared to measurements under ambient and lunarlike conditions for small (<45 μm) and large (75-125 μm) grains. The maximum isothermal sensitivity depths for the two size fractions are 260 μm and 990 μm, respectively. The model yields generally better results for the smaller grains, reproducing qualitatively the observed changes in spectral contrast and shift of the emissivity maximum to shorter wavelengths under simulated lunar conditions. The anisothermal spectrum for the larger size fraction is more difficult to model; the qualitative discrepancies seen in the intermediate case (Figure 9) at shorter and longer wavelengths become more severe, and none of the parameterized temperature profiles matches measurements particularly well.

Plagioclase (Labradorite)
Finally, we apply the modeling approach developed above to labradorite, a plagioclase mineral, leveraging lab measurements  The results shown in Figure 12(a) are comparable to those obtained for olivine and enstatite: the locations of modeled spectral features generally agree with measurements, although there are discrepancies in the magnitude of effective emissivity. One notable difference is the appearance of a double peak near the emissivity maximum in the modeled spectrum, attributable to prominent features in spectra calculated from the three principal complex refractive indices. This double peak is not as apparent for larger or smaller size fractions. The maximum isothermal sensitivity depth is found to be 300 μm (Figure 12(b)). Figure 13 illustrates the effects of modeled anisothermality on spectral characteristics. Notably, none of several modeled linear temperature gradients reproduced the emission spectrum measured under simulated lunar conditions. In addition to the lack of agreement in both the position and direction of shift of the Christiansen feature, there is a marked difference in spectral slope between 10 and 25 μm, with the modeled anisothermal spectra showing relatively flat, rather than increasing, emissivity in this wavelength range.
There are a few potential reasons for the discrepancy in spectral slope. One may be the presence of clinging fines; Shirley & Glotch's (2019) measurements of different size fractions of labradorite show that the "upward" slope observed is more pronounced for smaller grains. Alternatively, it may be the case that the temperature profile established for this labradorite sample in the lab is significantly nonlinear, and not well represented by a modeled linear gradient. This possibility is also suggested by measurements presented by Donaldson Hanna et al. (2012b), which indicate that changing how a sample is heated/cooled under vacuum conditions (leading to a different thermal gradient) results in differences in measured emission spectra. Measurements of a range of fine particulate plagioclase minerals by Donaldson Hanna et al. (2012b) show a rise in emissivity between 10-25 μm only when a sample is heated from below in a cooled vacuum chamber, but not when the chamber is maintained at room temperature. The measurements shown in Figure 13 were acquired in a simulated lunar environment in which samples were heated from above (to simulate solar illumination) and below in a cooled chamber (Shirley & Glotch 2019). Figure 14 compares modeled and measured spectra, as well as isothermal sensitivity depths, for smaller and larger size fractions of labradorite: a <32 μm size fraction, modeled using Mie scattering properties for 16 μm spheres, and a 63-90 μm size fraction, modeled using Mie scattering properties for 76.5 μm spheres. Both modeled and measured spectra for the 63-90 μm grains (Figure 14(c)) are similar to the 32-63 μm grains (Figure 12(a)); the modeled best-fit mixture of orientations is also near-identical. Notably, the best-fit mixture of orientations for the <32 μm case includes no contribution from one of the principal refractive indices (unlikely to be the case in reality). Maximum isothermal sensitivity depths approximately scale with grain size. Figure 15 shows representative, modeled anisothermal spectra compared to measurements under ambient and lunarlike conditions for small (<32 μm) and large (63-90 μm) grains. In both cases, the model still has difficulty capturing the position and direction of shift of the Christiansen feature. However, model results for the smaller size fraction (Figure 15(a)) show a better agreement in the overall spectral slope between ∼10-25 μm, suggesting that the clinging fines thought to be present may contribute to the similar upward slope seen in measurements of larger size fractions.

Summary and Conclusions
This work presents a Monte Carlo radiative transfer approach to modeling thermal emission under conditions similar to those found on the Moon and other airless bodies, where low thermal conductivity can give rise to significant thermal gradients within the epiregolith-the uppermost layer of the surface to which remote sensing observations are sensitive.
We applied this approach to model lab analog regoliths composed of three different minerals (corresponding to major lunar compositional phases) and different grain size fractions, studied under both ambient and simulated lunar conditions. The availability of ambient (isothermal) spectra for each of these cases allows us to correct to some extent for the effect of model simplifications, such as the use of a single intermediate grain size to approximate the behavior of a size fraction, and the use of Mie theory to approximate the scattering properties of a medium of irregular, closely packed grains.
We model a series of anisothermal spectra, corresponding to different, linear temperature gradients, and find that this approach can provide constraints on the temperature profiles established in the lab setting, and help to predict (with varying accuracy) how ambient spectra change under simulated lunar conditions. In general, this modeling approach is able to capture the changes in spectral contrast, but is less successful at modeling shifts in the position of the Christiansen feature emissivity maximum. The temperature profiles that most closely match measured spectra under lunar-like conditions are similar in magnitude to those calculated by Henderson & Jakosky (1997) and Millán et al. (2011) using a forward-modeling approach. The model results presented here also illustrate how emission spectra may change as thermal gradients vary with latitude, and over the course of a lunar day (e.g., Hale & Hapke 2002).
In addition to emission spectra, Monte Carlo radiative transfer calculations also provide an estimate of the depth of the epiregolith thermal gradient. We define the isothermal sensitivity depth (a function of wavelength) as the depth within which 99% of thermal emission from a medium (at a given wavelength) originates under isothermal conditions. The maximum isothermal sensitivity depth is characteristic of the depth of the epiregolith thermal gradient. This parameter varies with composition and particle size as shown in Table 3, and increases with decreasing packing density. The sensitivity depth would increase in the presence of a positive epiregolith thermal gradient (i.e., increasing temperature with depth within the uppermost subsurface) as is thought to be the case under lunar-like conditions.
The modeling approach presented here involves several simplifications. In addition to the assumption of monodisperse Mie scatterers, we neglect any variation of optical constants with temperature as well as any variation of physical properties (e.g., packing density, relative proportions of orthogonal crystallographic orientations) with depth. The results presented in Section 3 demonstrate that the net effect of these simplifications can be accounted for to some extent by leveraging available lab measurements of emission spectra under ambient conditions. However, the agreement between modeled and measured spectra could potentially be improved by revisiting these simplifications. One advantage of the Monte Carlo approach, and a target for future research, is that the methods described here can readily be adapted to accommodate nonuniform regolith properties, or   Note. Isothermal sensitivity depth is defined as the depth within which 99% of thermal emission from a medium (at a given wavelength) originates under isothermal conditions. The maximum isothermal sensitivity depth is a measure of the depth of the epiregolith thermal gradient.