Modeling the Angular Dependence of Emissivity of Randomly Rough Surfaces

Directional emissivity (DE) describes how the emissivity of an isothermal surface changes with viewing angle across thermal infrared wavelengths. The Oxford Space Environment Goniometer (OSEG) is a novel instrument that has been specifically designed to measure the DE of regolith materials derived from planetary surfaces. The DE of Nextel high‐emissivity black paint was previously measured by the OSEG and showed that the measured emissivity decreases with increasing emission angle, from an emissivity of 0.97 ± 0.01 at 0° emission angle to an emissivity of 0.89 ± 0.01 at 71° emission angle. The Nextel target measured was isothermal (<0.1 K surface temperature variation), and the observed change in emissivity was due to Fresnel‐related effects and was not due to nonisothermal effects. Here we apply several increasingly complex modeling techniques to model the measured DE of Nextel black paint. The modeling techniques used here include the Hapke DE model, the Fresnel equations, a multiple slope Fresnel model, and a Monte Carlo ray tracing model. It was found that only the Monte Carlo ray tracing model could accurately fit the OSEG measured Nextel data. We show that this is because the Monte Carlo ray tracing model is the only model that fully accounts for the surface roughness of the Nextel surface by including multiple scattering effects.


Introduction
Three-dimensional (3-D) thermophysical models of airless bodies in the solar system are used to calculate the surface and subsurface temperature distribution across a planetary body. By comparing surface temperature measurements from an orbiting thermal infrared (TIR) radiometer instrument to modeled 3-D thermophysical surface temperatures, a number of properties can be inferred about the surface including composition, rock abundance, emissivity, albedo, density, and thermal conductivity (Bandfield, 2009;Bandfield et al., 2015). These properties are in turn useful for investigating and understanding numerous areas of planetary scientific research and for the characterization of future landing and in situ resource utilization sites. For example, surface particle size and roughness are important for understanding regolith development. Three-dimensional thermophysical models of the lunar surface are also vital to help constrain the distribution of volatile-rich materials in the polar regions.
Three-dimensional thermophysical models use topographical information to model the surface of an airless body as a mesh of surface elements. Each element is split into layer depths, and the one-dimensional thermal diffusion equation is then solved to obtain the temperature depth profiles for each element Rozitis & Green, 2011;Vasavada et al., 2012). The surface temperature is calculated using the incoming solar ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. flux, reflected visible flux, and reradiated TIR radiation from other sunlit elements. If the surface is flat, then there is no contribution from reflected visible or reradiated TIR radiation; however, if the surface is curved like a crater, then the reflected visible and reradiated TIR radiation emitted from an element can be absorbed by another element. For the lunar surface scattering and reradiation have little effect near the equator, where the main driver for the surface temperature is the high solar incidence angle; however, at the poles scattering and reradiation have a significant impact on the surface temperature . Due to the low obliquity of the Moon, the solar incidence angle remains low throughout the day at the lunar poles and large shadows are cast across the terrain. Thus, the way in which illuminated surface elements scatter visible and reradiate TIR radiation into other elements that are in shadow is important in modeling the surface temperature. By extension, the surface temperature then controls the distribution of volatile-rich materials in the polar regions.
Typically, 3-D thermophysical models assume that sunlit elements scatter visible and reemit TIR radiation in an isotropic fashion. This assumed isotropic scattering has previously been shown to be incorrect in the visible by a number of laboratory measurements Johnson et al., 2009;Pommerol et al., 2011). There is also strong evidence from remote sensing measurements (Bandfield et al., 2015) and a very limited number of laboratory measurements (Warren et al., 2017) that the isotropic assumption is also incorrect in the TIR. The assumption of isotropic scattering in both the TIR and visible is believed to be the reason for differences of up to tens of Kelvin in previous 3-D thermophysical models of the lunar polar surface when compared to infrared remote sensing measurements .
How TIR is reemitted from a surface is also significant for better constraining processes acting on the surfaces of asteroids and the orbital and spin evolution of asteroids. Due to an asteroid's spin, TIR radiation is emitted asymmetrically from an asteroid. This asymmetrical emission causes a net photon force and torque. The net force causes an asteroid's orbit to drift and is known as the Yarkovsky effect. The net torque changes an asteroid's spin rate and is known as the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect. An accurate 3-D thermophysical model is required to model the surface temperature distribution of an asteroid to allow the YORP and Yarkovsky effects to be accurately estimated. Such a thermal model must therefore include realistic scattering and reemission from the asteroid's surface. For modeling the YORP and Yarkovsky effects the direction of TIR emission and reemission from the asteroid surface is also important, since this dictates the direction of the photon force and torque that causes the YORP and Yarkovsky effects. In the asteroid community the directionality of TIR reemission or emission from a surface is often referred to as "thermal beaming" (Harris, 1998).
At visible and near-infrared wavelengths (0.1-1 μm), several goniometers including the Bloomsburg University Goniometer and the Physikalisches Institut Radiometer Experiments have measured the bidirectional reflectance distribution function (BRDF) of Apollo samples and of the lunar analogue material JSC-1AF Pommerol et al., 2011). The BRDF describes how visible and near-infrared light is scattered from a particulate surface for all viewing angles (Hapke, 2014). These measurements have shown that Apollo samples and JSC1-AF do not scatter light in an isotropic fashion Pommerol et al., 2011;Shepard, 2002). Instead, the intensity of scattered radiation generally decreases with increasing reflection angle. In the visible, numerous models exist to theoretically calculate the nonisotropic scattering BRDF of a surface and the most widely used model is that of Hapke (Hapke, 2014). Hapke defines the BRDF as the ratio of the scattered radiance at a detector to the collimated incident irradiance illuminating a surface.
where r(θ r , θ i , θ p ) is the BRDF and θ r , θ i , and θ p are the angles of reflection, incidence, and azimuth shown in Figure 1. L i is the irradiance incident on the surface at incidence angle θ i , and E r is the radiance reflected from the surface at reflection angle θ r .
The DE modeling described in this paper is for isothermal surfaces and does not account for more realistic rough surfaces that will likely be nonisothermal. For rough thermally heterogeneous surfaces variations in emission with observation angle would be caused by changes in the fraction of hot and cold surfaces that are visible as the observation angle is changed and by the Fresnel-driven effects described in this paper. When analyzing remotely sensed emissivity measurements of planetary surfaces, the DE modeling described in this paper will need to be accounted for as part of a multilayered problem.
Although there have been scattering function measurements of Apollo and JSC1-AF samples made in reflection across visible and near-infrared wavelengths, there have been no equivalent laboratory measurements made in emission across near and TIR wavelengths (1-25 μm). The BRDF is strictly defined only for reflection in the visible and near infrared and is not applicable to near, thermal, and far infrared (1-400 μm) light that is emitted or reemitted from a surface. Emitted radiation is emitted from a hot surface that is heated by a nondirectional source; this could be, for example, an electric heater placed under the surface. This is in contrast to reemitted radiation, which is emitted from a hot surface heated by a directional light source. In the case of airless bodies in the solar system this would be solar illumination, which is absorbed and reemitted. For infrared wavelengths, the directional emissivity (DE) describes how radiation is reflected, emitted, and reradiated from a surface for all viewing angles (Warren et al., 2017). In some previous works, the DE has also been referred to as the emission phase function (Bandfield et al., 2015).
This paper describes efforts to numerically model the DE measurements made of a surface heated from below. When a surface is heated from below, rather than from above with a light source, there is no incidence angle. Instead, the DE describes how the emission of the hot surface (>0 K) varies with emission and azimuth angles, where the azimuth angle is defined from an arbitrary plane, since there is no incidence plane. The DE in this case is the radiance emitted from the sample divided by the radiance emitted by an ideal blackbody that is under the same conditions such as wavelength range, surface temperature, and viewing angle (equation (2)).
where ε d is the DE, U b is the total radiance emitted by a blackbody at the same temperature as the surface, E is the radiance emitted by the surface, and θ e and θ p are the emission and azimuth angles, respectively, as shown in Figure 2. T is the temperature of the surface, and λ is the observed wavelength range.
Although there have been no DE laboratory measurements of lunar soils made in the infrared, as part of its third extended science mission, the Diviner Lunar Radiometer Experiment (Diviner) has been making offnadir measurements of the lunar surface (Bandfield et al., 2015). Diviner is an infrared radiometer onboard NASA's Lunar Reconnaissance Orbiter (LRO) that has been orbiting the Moon since its launch in June 2009 . Diviner is a nine-channel radiometer with two solar channels (0.35-2.8 μm), three narrow band channels near 8 μm, and four thermal/far infrared channels ranging from 12.5 to 400 μm. As well has making global nadir multispectral measurements of the lunar surface, Diviner has made multispectral off-nadir TIR measurements at several different surface locations at multiple local times of day.
The initial results from the Diviner off-nadir measurements have been reported in Bandfield et al. (2015) and show that the lunar surface is not an isotropic emitter in the thermal and far infrared spectrum. Diviner is able to make off-nadir measurements by using its two-axis pointing capability to target a surface location from multiple angles as it passes over the surface. Thus, the off-nadir measurements made by Diviner are similar to the type of observations made by a laboratory goniometer and can be compared to TIR DE goniometer measurements made in the laboratory.  By investigating in the laboratory how surface properties such as composition, particle size, thermal conductivity, and surface roughness affect the measured emissivity at different viewing angles, one could infer these properties from the off-nadir global measurements acquired by Diviner. Comparisons between these types of laboratory measurements and remote sensing observations will become increasingly important with the large amounts of off-nadir data expected from missions including NASA's Origins Spectral Interpretation Resource Identification Security Regolith Explorer and the Japanese Aerospace Exploration Agency's Hayabusa2. Both missions have TIR radiometers, spectrometers, and/or cameras onboard that will acquire measurements of their target asteroids at many different times of day and at many different viewing angles.  (2) the Fresnel equations. When neither of these models were able to fit the measured DE, the complexity of the modeling was increased to include surface roughness with (3) a multiple slope Fresnel model and finally multiple scattering and shadowing with (4) a Monte Carlo ray tracing model. A description of each model that was able to fit the OSEG DE measurements is described in section 2.

Hapke's Directional Emissivity Model
The first model used to model the DE of Nextel black paint was the Hapke (2014) DE model. This is different from the more commonly known Hapke BRDF model that attempts to describe how visible light is scattered from a particulate surface (Hapke, 2014). Hapke's BRDF model is complex with a number of additional parameters (11 parameters in the latest version of the model) that attempt to capture how properties, including the composition and porosity of a surface, affect the scattering behavior in the visible and near infrared.
Hapke developed a relatively simple model for the DE in the TIR with only one parameter, the single scattering albedo. The Hapke DE model predicts the angular dependence of the TIR radiance emitted by a surface of an optically thick particulate medium at uniform temperature T (equation (3); Hapke, 2014).
, ω is the single scattering albedo, μ = cos (θ e ), and θ e is the emission angle (defined in Figure 2). The Hapke modeled DE is plotted for several values of the single scattering albedo in Figure 3a. Hapke's model predicts that for materials with a small (<0.1) single scattering albedo the DE is approximately unity for all emission angles and hence the surface emits like a blackbody. Therefore, Hapke's emission model states that for high-emissivity materials the emissivity can be approximated independently of the emission angle. However, Hapke (2014) also states that only a few measurements of the DE have been published, and further measurements are required to test the model.
Hapke's DE model is valid for any particulate medium that is at a uniform temperature; thus, it should be valid for the Nextel target in Warren et al. (2017) as the sample's surface was at a uniform temperature (±0.1 K), and Nextel black paint is primarily made up of spherical polyurethane balls that range from 20 to 100 μm in diameter (described in more detail later). An unconstrained Matlab least squares fitting algorithm was used to fit the Hapke DE model to the measured DE of Nextel black paint with the single scattering albedo (w) as an unconstrained free parameter. However, Hapke's DE model was found not to fit the measured DE of Nextel black paint ( Figure 3b). Hapke's emission model predicts a more isotropic DE and does not capture the large (Δ0.1) drop in emissivity observed with increasing emission angle, particularly between 50°and 70°.

Fresnel Equations
The DE of an optically smooth surface can also be modeled using the Fresnel equations. The Fresnel equations for reflection and transmission use the optical constants of a medium (n and k) to describe the behavior of the parallel and perpendicular components of an electromagnetic wave incident on a dielectric or conducting surface. The optical constants of a medium are known as the real (n) and imaginary (k) parts of the refractive index and are related to the electromagnetic properties of the material by equations (4) and (5).
where μ m0 , ε m0 , and c 0 are the magnetic permeability, electrical permittivity, and the speed of light in vacuum, respectively. The μ m , ε e , and σ e are the relative magnetic permeability, relative electric permittivity, and electrical conductivity, respectively, of the medium, and ω f is the angular frequency of the electromagnetic wave. A material's complex refractive index is defined as n ¼ n þ ik. The imaginary part of the refractive index is known as the extinction coefficient as it characterizes the attenuation of electromagnetic waves through a medium. A large extinction coefficient means that the beam is quickly attenuated while passing through the medium, and a small extinction coefficient means that the medium is relatively transparent to the electromagnetic radiation. From equation (5), a perfect dielectric such as Nextel paint (σ e = 0) has an imaginary part of the refractive index equal to zero and hence its complex refractive index only contains a real part, n. In this case the Fresnel equations simplify to those shown in equations (6) and (7).

Journal of Geophysical Research: Planets
where n is the ratio of the refractive index between the two media, r p is the reflection of the parallel component, and r s is the reflection of the perpendicular component. Equations (6) and (7) are known as the Fresnel equations of reflection for an ideal dielectric medium. For unpolarized light, where there is an equal intensity from the parallel and perpendicular components, the reflection for an ideal medium can be calculated by taking the average of equations (6) and (7) to give (equation (8)) where θ r is the angle of reflection ( Figure 1) and r t is the reflected electromagnetic beam for unpolarised radiation. The DE of a perfect dielectric can then be calculated using Kirchoff's law (equation (9)).
where ε d is the DE at emission angle θ e , r t is the radiation reflected at angle θ r given in equation (8), θ r is the angle of reflection, and θ e is the angle of emission and in this case they are equal (i.e., the angle of reflection is the angle of emission).
According to the Fresnel equations, the DE of a perfect dielectric is dependent on only one free parameter, the ratio of the two media's refractive indices, n. The DE of a dielectric is isotropic for all values of emission angle when n is unity. For larger values of n the DE is approximately isotropic for θ e less than about 60°, past which the DE decreases with increasing emission angle. In addition, as n increases from unity the normal emissivity (DE at θ e = 0) decreases (Figure 4a).
An unconstrained Matlab (version 2016a) least squares fitting algorithm (Trust-Region) was used to fit the Fresnel equations to the DE measurement of Nextel black paint with the refractive index as the only free parameter (the refractive index of air was assumed to be unity). The best fit (Figure 4b) was found for a refractive index of 1.37, but the fit was poor with an R 2 value of 0.5. The Fresnel equations fit the measured DE well at emission angles <50°, but at larger emission angles the Fresnel equations predict a greater decrease with emission angle than is observed in the measurement.
The poor fit between the measured DE values of Nextel paint and the Fresnel predicted values has previously been observed by De Silva and Jones (1987), who suggested two possible reasons for this difference: (1) The Nextel paint was sufficiently transparent so that some emission from the aluminum substrate below the Nextel paint influenced the measurement or (2) the surface of Nextel black paint is sufficiently rough at the micron scale that an increase in DE with emission angle compared to an ideal dielectric is possible. As Nextel black paint is a dielectric with high emissivity, it is unlikely that any thermal radiation from the underlying metal surface affected the measurement. We assert that the surface roughness was the reason for the difference between the Fresnel-predicted DE response and the measured DE. A new more complex and realistic model (a multiple slope Fresnel model) that accounted for surface roughness was then used to try to better understand the effect of surface roughness on the DE.

Multiple Slope Fresnel Model
The Fresnel model was adapted to include micron-scale roughness effects by modeling the surface as a set of slopes with a Gaussian distribution. The Gaussian slope roughness model was taken from Bandfield et al. (2015) and used to model the DE of the Nextel target. The model calculates the radiance emitted from slopes of 0°to 90°at 1°intervals. The radiance of each slope is calculated using the Fresnel equations for a perfect dielectric, and its contribution to the total modeled radiance is weighted by the Gaussian statistical probability of its occurrence. This reduces surface slopes to a single parameter, the root-mean-square (RMS) slope, which is independent of length scale.
The probability distribution P for a given slope angle is given in Bandfield et al. (2015) as where θ 0 is the RMS slope angle and θ s is the angle of the slope from the surface normal. This describes the adirectional distribution of slopes that closely approximates a Gaussian distribution of unidirectional slopes for a RMS slope angle of θ 0 (Bandfield et al., 2015).
Using the probability distribution in equation (10), the model calculates the mixture of Fresnel DE functions for each slope angle in proportion to their contribution inside the measurement field of view based on the observational geometry. Therefore, slopes that face away from the observer will contribute proportionally less to the observed radiance than slopes that face the observer. This factor is accounted for in the model by the dot product of the observation vector and the vector normal to the local surface in equation (10) (Bandfield et al., 2015).
The best fit of the multiple slope Fresnel model to the measured DE of Nextel was found to have a better R 2 value (0.84; Figure 5) compared to the simple Fresnel model (R 2 value of 0.5) indicating that the poor fit of the simple Fresnel method is at least in part due to surface roughness effects. Using a constrained least squares fit, the best fit was found to be a RMS slope value of 23°and a refractive index of 1.35. The RMS slope was constrained to be between 0°and 90°, and the refractive index was constrained to be between 0 and 10. However, the fit of the multiple surface Fresnel model is still poor, which could be due to the model not fully accounting for surface roughness. Although the multiple slope model does account for roughness by accounting for a Gaussian distribution of slopes, it does not account for multiple scattering between surfaces or view shadowing. A new Monte Carlo ray tracing model that does account for multiple scattering and view shadowing was then used to further explore this effects on the DE.

Monte Carlo Ray Tracing
Both multiple scattering and view shadowing are inherently accounted for in ray tracing models where a packet of rays is traced backward from the detector through their interactions with the visible parts of the surface. View shadowing occurs when certain sections of a rough surface are blocked from the observer's view at a particular emission angle, and multiple scattering occurs when radiation interacts more than once against a rough surface ( Figure 6). View shadowing will affect the DE since at high emission angles; some slopes can be hidden from the observer's view (Figure 6), and large slopes that face away from the observer will be preferentially hidden from the observer's view.
From the point of view of the previously discussed multiple slope Fresnel model, the effects of view shadowing will modify the Gaussian slope distribution such that slopes facing away from the observer will not contribute to the emissivity. It is possible to account for view shadowing in a multiple slope Fresnel model by using a shadowing approximation methodology as outlined in Bandfield et al. (2015). The multiple slope Fresnel model used in the previous section (3.3) does not account for view shadowing. When the view shadowing approximation methodology was included in the multiple slope Fresnel model (see Bandfield et al., 2015, for model details), it was still found not to adequately fit the measured DE (R 2 value of 0.89 as seen in Figure 7).
The Monte Carlo ray tracing (MCRT) model works by tracing a bundle of rays backward through the system, that is, in this case rays are emitted from the detector at a particular emission angle onto the rough surface. When a ray hits the surface, the direction it reflected is calculated taking into account the local slope and using Snell's law. If the reflected ray hits the surface again, its new reflected direction is once again calculated using Snell's law. However, if the reflected ray does not interact with the surface again, then it is defined as "reflected from the surface" and the energy associated with that ray is defined as reflected from the surface. Totaling the amount of energy that is reflected from the surface and dividing by the total amount of energy the ray bundle had when emitted from the detector gives the directional reflectance of the surface at a particular emission angle. Kirchhoff's law (equation (9)) is then used to calculate the emissivity of the surface as one minus the reflectance at a particular emission angle.
The ray tracing algorithm and source code used to model the DE of the Nextel black painted surface is publicly available to download from the MySimLabs website (Bergström, 2012), and a full description of the

10.1029/2018JE005840
Journal of Geophysical Research: Planets algorithm is given in Bergström et al. (2007). Because the algorithm was designed to model the BRDF of samples, we have slightly modified the algorithm to calculate the DE from the angular reflection using Kirchhoff's law (equation (9)).
In MCRT models, random rough surfaces are often characterized in statistical terms using two distribution functions, the height distribution (P(ζ(x, y))) in the vertical direction and the autocovariance function (C τ ð Þ) in the lateral directions. The height distribution describes the height deviation from a certain mean reference level, and the autocovariance function describes the typical distance between two similar features (e.g., hills or valleys) (Bergström, 2012;Bergström et al., 2007Bergström et al., , 2008. If a random rough surface is described by the function z = ζ(x, y), then the height probability function is given in equation (11).
where H rms is the RMS height. The autocovariance function is normally defined using a Guassian function as where x 1 and x 2 are two difference points along the surface and τ is the correlation length (the typical lateral distance between two features on a surface). See Figure 8 for an example of surfaces with the same H rms but different correlation lengths.
In light scattering theory, ray tracing is an energy approximation where bundles of radiation rays are traced throughout their interactions with a rough surface until they leave the surface. At each interaction the Figure 7. Comparison between the best fit modeled directional emissivity using a multiple slope Fresnel model with and without a view shadowing approximation to the measured directional emissivity of Nextel paint. RMS = root-meansquare.

10.1029/2018JE005840
Journal of Geophysical Research: Planets surface is treated as locally smooth, so that each scattering event on the surface can be treated as a specular reflection and modeled using the Fresnel equation (described in section 2.2), which is known as the Fresnel approximation (Bergström et al., 2007). Tang et al. (1996) showed that the Fresnel approximation is valid for where H rms and τ are the RMS height and the correlation length of the surface, respectively, λ is the wavelength of the radiation, and θ i is the angle of incidence.
For modeling the OSEG experiments the angle of incidence is equal to the angle of emission (since the rays are traced back from the detector) and the OSEG measures over an emission angular range of 0-72°. The peak of the Planck function at the surface temperature (423 K) used to measure the DE of the Nextel paint in the OSEG is~7 μm. Putting these numbers into equation (13), the Fresnel approximation is valid at 72°e mission angle (the maximum measurable emission angle of OSEG) when the H rms is greater than 3.85 μm and the correlation length is greater than 1.9 μm. The correlation length and the RMS height of the Nextel target were measured using a noncontact profiler (described in section 2.4.1 below) and shown to be well above these values (13.5 and 95 μm, respectively), so the use of the MCRT model is valid over the angular range measured by the OSEG.

Nextel Surface Roughness Measurements
To model the DE of the Nextel target using the MCRT technique, the surface roughness profile was measured under a visible light microscope (Figure 9a). Surface roughness features were observed on scales approximately <100 μm in length. When illuminated with blue (488 nm) light, the polyurethane balls that make up Nextel fluoresce green (Figure 9b). The image in Figure 9b is 1.1 mm across and as such can be analyzed to show that the largest polyurethane ball is approximately 25 μm in diameter. The roughness profile of the Nextel target was further measured using a Polytech TopSens configurable Chromatic Confocal Sensor (Polytec Limited, 2015) that can measure surface roughness at 3-nm sampling (Figures 9c and 9d). Using the TopSens, the RMS height distribution of the Nextel target was found to be 13.5 μm and the correlation length was 95 μm. The measured surface roughness profile in the form of an array of heights (X, Y, Z) of the Nextel target was used in the MCRT model to produce a modeled DE. Using the measured surface roughness profile, the ray tracing model has only one free parameter, the real part of the refractive index (n). At an emission angle of 0°, the ray tracing algorithm showed that no multiple scattering or view shadowing occurs and the emissivity is defined only by the Fresnel formula, which implies that the real part of the refractive index of the material must be 1.43.
Using a refractive index value of 1.43, the MCRT model was run to calculate the DE shown in Figure 10c. The DE predicted by the MCRT model fits the measured DE of the Nextel black paint well with an R 2 value of 0.995. The MCRT model also records the average number of scattering events at each emission angle (Figure 10a) and the proportion of surface elements that are view shadow hidden (Figure 10b). The average number of scattering events between rough surface elements begins to increase with increasing emission angle from unity at~15°to 1.278 at 70°. The percent of surface elements that are view shadowed is zero for emission angles between 0°and 50°. At emission angles >50°, the percent of surface elements that are view shadowed increases with increasing emission angle from 0% at 50°to 1.27% at 70°. From Figures 10a  and 10b it is not clear which effect, view shadowing or multiple scattering, has a bigger impact on the DE. The following section investigates this in further detail.

Comparing the Fresnel Multiple Slope Model and the Monte Carlo Ray Tracing Model
In this section we investigate why the MCRT model could fit the measured DE of Nextel paint and why the multiple slope Fresnel model was not able to fit the measured DE. The MCRT model accounts for both view shadowing and multiple scattering, but the simple multiple slope Fresnel model used in section 2.3 does not account for either. The multiple slope Fresnel model was adapted to include view shadows (Figure 7) but was still found not to fit the measured DE of Nextel paint. This would strongly suggest that multiple scattering is critical for modeling the DE of a surface. However, the view shadowing included in the multiple slope Fresnel model was for a generic slope distribution rather than measured surface roughness of the Nextel target. In this section we show conclusively multiple scattering is required to accurately model the DE of a surface and we explore how large an effect multiple scattering and view shadowing have on the predicted DE.
To compare the predicted DE from the two models, we ran the models using the same input surface properties (e.g., Nextel surface profile) and optical constants. Since the multiple slope Fresnel model uses a Gaussian slope distribution as an input, the best fit Gaussian slope distribution was calculated from the measured Nextel surface profile ( Figure 11). Using this slope distribution and a real optical constant value of 1.43 (and an imaginary optical constant value of 0), the calculated DE from the two models is shown in Figure 12. Figure 12 also compares the predicted DE from the two models to the DE predicted by the Fresnel equations for the interested reader.
The MCRT and multiple slope Fresnel model predict different DE for the same input parameters (Figure 12). The MCRT model predicts the emissivity decreases at a slower rate with increasing emission angle when compared to the multiple slope Fresnel model. To investigate if this difference is caused by multiple scattering or view shadowing or a combination of both, the multiple scattering inside the MCRT model was switched off and the model was run again without multiple scattering-such that each ray was limited to only one scattering event. The DE calculated from MCRT model with and without including multiple scattering is shown in Figure 13 and is also compared to the multiple slope Fresnel model.
As seen in Figure 13, multiple scattering becomes increasingly important with increasing emission angle. Multiple scattering has no effect on the DE until~30°emission angle, at which point the DE decreases with increasing emission angle at a faster rate than when multiple scattering is not accounted for. Also, in Figure 13 we observe that multiple scattering accounts for the majority of the difference between the MCRT model and the multiple slope Fresnel model. The MCRT model without scattering and the multiple slope Fresnel model now agree within 0.001 emissivity at emission angles <60°. Between 60°and 70°there is a small difference between the two models, which increases with increasing emission angle starting with a difference of 0 at~60°increasing to a difference of 0.006 at an emission angle 70°. Since the only difference between these two models is view shadowing, it is reasonable to suggest that this difference must be due to view shadowing.

Journal of Geophysical Research: Planets
We can demonstrate the effects of view shadowing by modifying the MCRT models such that it records the local slope each ray hits at different emission angles. This shows how the slope distribution calculated by the MCRT model changes at different emission angles ( Figure 14). From Figure 14, the slope distribution is predicted to change at different emission angles because some slopes are preferentially hidden from view at certain emission angles (this is the process of view shadowing as described earlier in section 2.3). At lower emission angles slopes facing away from the observer are less likely to be seen, which skews the slope  distribution at lower emission angles ( Figure 14). When this changing slope distribution is accounted for in the multiple slope Fresnel model, the calculated DE is the same as calculated DE from the MCRT with no multiple scattering model ( Figure 15). We can therefore conclude that multiple scattering has the larger effect on the predicted DE for a rough surface than view shadowing. For the Nextel surface profile, both multiple scattering and view shadowing had the largest effect on the DE at an emission angle of 70°. At this emission angle multiple scattering changes the DE by 0.02 and view shadowing changes the DE by 0.006. Figure 14. Slope distributions at every emission angle. As the emission angle increases above~50°the slope distribution changes such that large slopes facing away from the observer (i.e., negative slopes) are less likely to be seen. Figure 15. A multiple slope Fresnel model that accounts for the changing slope distribution agrees without error with a MCRT model that does not account for scattering. MCRT = Monte Carlo ray tracing.

When Is a Surface Optically Smooth?
The DE of a theoretical perfectly smooth surface (i.e., with an H rms of 0) will follow a Fresnel function. The DE will be modified by adding roughness to the surface (i.e., with an H rms of >0) and can be accurately modeled by accounting for multiple scattering using an MCRT model. The difference between the DE of a smooth surface predicted by the Fresnel equations and the DE of a rough surface predicted by a MCRT model can be seen in Figure 16a. The emissivity of a smooth surface decreases at a faster rate with increasing emission angle. The "smoother" a surface is, the closer to a perfect Fresnel response the MCRT model predicts the DE to be. As the modeled surface approaches perfectly smooth, the MCRT model will predict that the DE of the surface will follow a Fresnel response. In this section we answer the question: how "smooth" does a surface need to be before the DE is predicted to follow a perfect Fresnel response? Meaning, when can a surface be thought of as optically smooth?
To quantify how smooth a surface is, we calculated the ratio between the correlation length (see equation (12) for definition) and the H rms , where a low ratio value indicates a smooth surface and high ratio value indicates a "rough" surface. For the surfaces shown in Figure 8, one can see the surface in Figure 8a is more smooth than a surface in Figure 8b. Figure 8a has a H rms of 1 μm and a correlation of 10 μm giving a ratio of 0.1, and Figure 8b has a correlation length of 10 μm and a H rms of 10 μm giving a ratio of 1.
A Matlab function that generated a random rough surface profile with different correlation and H rms lengths was used to generate hundreds of rough surfaces with varying ratios of smoothness. These surface roughness profiles were then used as an input to the MCRT model with a real part of the refractive index between 1.1 and 4 and imaginary part of 0. The predicted DE for a selection of the surfaces with different ratios of smoothness is shown in Figure 16a to demonstrate the effect of increasing surface roughness on the DE.
The total difference over the emission angular range 0-80°between the MCRT model predicted DE of the randomly generated rough surface profiles and the equivalent Fresnel response of the surface was calculated and plotted against the smoothness ratio (Figure 16b). In Figure 16b the difference is the total difference between 0°and 80°between each of the solid lines and the red dashed Fresnel line. The total difference between the Fresnel response and the modeled MCRT DE should be thought of as the difference between the DE of a rough surface and a perfectly smooth surface (i.e., when the total difference is close to 0, the surface can be described as optically smooth and will follow a Fresnel function). From Figure 16b any surface with a ratio of >40 can be thought of as optically smooth and the DE of the surface over the angular range 0-80°will follow the Fresnel equations.

Conclusions
Previous DE measurements of Nextel black paint have shown that the DE is not isotropic, rather the emissivity decreases with increasing emission angle. A decrease with increasing emission angle is predicted by the Fresnel equations, but the modeling work presented in this paper has shown that the measured DE of Nextel paint does not follow a Fresnel response as would be expected for an optically smooth surface. We have shown that DE is dependent on the surface roughness, and to accurately model DE, surface roughness must be appropriately accounted for in the model. We also found that modifying the Fresnel model to account for surface roughness by calculating the average Fresnel response of a Gaussian distribution of multiple slopes does not accurately fit the measured DE. The measured DE could only be fit accurately by using a MCRT model that accounted for multiple scattering and view shadowing. Further analysis shows that multiple scattering has a greater effect on the modeled DE than view shadowing. This work has important implications for remote sensing observations of airless bodies. In general, previous works have assumed the surfaces of airless bodies emit in an isotropic fashion with emission angle. This means that regardless of emission angle, the emissivity of a surface will always appear to be the same as its nadir (normal) emissivity. Our results show that the emissivity of a surface decreases with emission angle and as such off nadir and hemispherical emissivity will generally be lower than nadir orientated emissivity. This affects several types of TIR remote sensing investigations of airless bodies. For example, 1. Thermophysical studies. When calculating the total emitted radiance from a surface, if the surface is assumed to be an isotropic emitter and the nadir emissivity value is assumed to be equal to the hemispherical emissivity, then the total emitted radiance from the surface will be overestimated. For the case of the Moon, this could mean assuming an emissivity of 0.95 instead of 0.98 , which could make differences in modeled surface and subsurface temperatures in the region of 10-50 K. 2. Yarkovsky and YORP studies. Compared to an isotropic emitter, a nonisotropic surface emitter will have a different angular distribution of the emitted radiance. This will potentially affect the radiative pressure, which can lead to changes in the spin rate of small bodies. 3. TIR spectral compositional studies. If the emissivity of a surface is lower at high emission angles, this could also potentially change the shape of the surfaces spectral emissivity when viewed from large emission angles. In general, one would expect the change in emissivity with emission angle to result in changes in the depth of the spectral features. However, the scale of sensitivity would be expected to change with wavelength, which may add effects such as slopes into the spectra.
For the application of TIR remote sensing of airless bodies, the Fresnel-driven DE described in this paper is part of a multilayered problem. The surface of an airless body will not be isothermal as the modeled Nextel surface was in this paper but will be nonisothermal. This will cause the DE of the surface to appear to change with viewing angle as different surface temperatures are viewed at different viewing angles. The modeling technique described in this paper addresses the isothermal Fresnel-related DE effects but does not address the nonisothermal DE effects.
As well as having important implications for TIR remote sensing observations of airless bodies, measuring and modeling the DE as described in this paper could also be used as a technique for measuring the optical constants of materials. The measurement and subsequent modeling of the DE of surfaces have shown that the DE is dependent on the optical constants of the surface material and the surface roughness. It was shown that a surface with a ratio between its correlation length and its H rms of <40 can be thought of as optically smooth. If the DE of such a surface was measured in the OSEG, it will follow a Fresnel function and the optical constants of the surface material could be derived by finding the best fit Fresnel function to the measurement.