An alternative to Hapke’s macroscopic roughness correction

In this paper, we present an alternative to Hapke’s correction for the photometric effect of macroscopic roughness of planetary surfaces (Hapke, 1984). Our model includes the effects of both single-facet scattering and multi-facet scattering. The single-facet scattering model is derived directly from the definition of bidirectional reflectance of a rough surface. We model projected shadowing by adopting an approximation published elsewhere in the literature. Monte Carlo simulations of single-facet scattering demonstrate very good agreement with the single-facet scattering component of our model. Using 3D printed molds with known roughness statistics, we prepared rough mineral surfaces consisting of quartz and olivine in the laboratory. We recorded extensive measurements of the rough surface bidirectional reflectance over a wide range of illumination and view geometries at 1751 wavelength bands ranging from 500–2250 nm. By fitting the residual between the measured bidirectional reflectance and our single-facet scattering model, we developed an empirical approximation for multi-facet scattering. Our results show that the multi-facet scattering is proportional to the surface’s diffusive reflectance, RMS slope, and the cosine of the illumination angle. Furthermore, the multi-facet scattering is approximately Lambertian for phase angles less than 90 ◦ but becomes forward scattering at higher phase angles. Our empirical multi-facet scattering model incorporates all of these features. As a surface becomes more macroscopically rough, our data show that the overall reflectance distribution becomes progressively more backscattering in a way that is distinct from the shadow hiding opposition effect. We refer to this effect as the macroscopic roughness backscattering bias (MRBB), which affects the entire hemisphere rather than being localized to small phase angles. Our proposed model is more accurate than Hapke’s for phase angles less than 90 ◦ , and the accuracy of the two models is comparable at higher phase angles. However, researchers should be aware of an issue regarding the roughness parameter defined in Hapke’s model, namely the fact that ‘‘ ϑ ’’ is not actually equal to the mean slope angle of the surface for the probability density function used in the model. This is an issue which apparently has gone unnoticed in the literature, and has potentially caused the model to be applied incorrectly in past studies.


Introduction
A major challenge in the field of planetary remote sensing is to adequately model the effects of surface roughness on bidirectional reflectance measurements. As discussed in Shepard and Campbell (1998), the parameterization of surface roughness is a difficult problem due to the fact that surface statistics, such as root-mean-square (RMS) heights or slopes, vary depending on the length scale over which they are measured. For this reason, it is useful to draw a distinction between the wavelength scale roughness, particle scale roughness, and macroscopic roughness of a particulate medium. The wavelength scale roughness of individual particles and its effect on the single scattering phase function have been described extensively (Hapke, 1986(Hapke, , 2012b. At the particle scale, surface particles may shadow deeper particles within the medium, resulting in the well known shadow hiding opposition employ well defined structures, such as crater shaped depressions or holes (Van Diggelen, 1958;Hämin-Anttila et al., 1965;Veverka and Wasserman, 1972;Lumme and Bowell, 1981a,b) whereas others employ statistical quantities such as RMS heights or slopes (Hapke, 1984;Despan et al., 1998Despan et al., , 1999 or a fractal measure (Shepard and Campbell, 1998). An alternative approach is to inverse model the roughness based on empirical observations without utilizing a detailed forward model of rough surface reflectance. For example, Mushkin and Gillespie (2006) proposes a ''two look'' method in which an area is imaged under different illumination and/or emergence directions and the difference in shadow fraction of each pixel is interpreted as a proxy for surface roughness. Alternatively, Cuzzi et al. (2017) presents an approach wherein the shape of the observed phase curve of an object is used to infer a shadowing parameter. This parameter is then used to correct for the effect of shadowing on the object's apparent spherical albedo. Further research is needed to investigate the relative performance of these inverse modeling strategies with respect to forward modeling. Perhaps they may be viewed as complementary.
The most widely used model for the photometric effect of macroscopic roughness of planetary bodies was introduced in Hapke (1984). This model forms one component of Hapke's overall photometric model of soil reflectance (Hapke, 2012b). In his approach, the slope of an arbitrary point on the surface ϑ is modeled as a random variable with probability density function (ϑ). His roughness correction treats the bidirectional reflectance of the macroscopically smooth surface as known. It is assumed that the relevant facet size far exceeds the particle size, such that the bidirectional reflectance of each facet is equal to that of the macroscopically smooth surface. He then estimates the rough surface bidirectional reflectance as the product of the macroscopically smooth surface reflectance, at an effective incidence angle and effective emergence angle, and a shadowing function. The model is parameterized by the quantity ϑ. Several attempts have been made to implement and validate Hapke's roughness model in the laboratory, in observations of terrestrial and extraterrestrial terrains, and using computer simulations. In Helfenstein (1988), computer generated imagery of cratered surfaces were used to compare the photometrically determined value of ϑ to its true value. The results were found to be generally consistent with the predictions of Hapke's roughness correction, but with lower accuracy reported at high illumination angles and high phase angles. Roughness statistics of the lunar surface were assessed in Helfenstein and Shepard (1999) using stereophotogrammetry and using Hapke's roughness correction. It was found that the majority of photometrically detectable roughness occurred at length scales between 0.1 mm and 8 cm for these surfaces. The parameters of Hapke's photometric model, including the roughness parameter ϑ, were inverted from surface images of Mars obtained by the Spirit and Opportunity rovers (Johnson et al., 2006a,b), from the High Resolution Stereo Camera (HRSC) (Jehl et al., 2008;Pinet et al., 2014;Gwinner et al., 2016), and from the Compact Reconnaissance Imaging Spectrometer (CRISM) (Fernando et al., 2013(Fernando et al., , 2015. A high spatial resolution inversion of Hapke's photometric model was also performed on the Lavoisier lunar crater using the Advanced Moon Micro-Imager Experiment (AMIE) camera in Souchon et al. (2013). Most recently, Hapke's photometric model was used to measure ϑ of dwarf planets and large Kuiper belt objects using the New Horizons spacecraft (Verbiscer et al., 2022). However, no in-situ or stereophotogrammetric measurement of roughness was performed to quantitatively validate the inversion results in any of these studies. The ϑ parameter has also been inverted from laboratory measurements of bidirectional reflectance of prepared surfaces. In Souchon et al. (2011), granular volcanic sand samples were prepared by gently shaking trays of material to create level surfaces. The samples were not pressed or scraped, however, such that the surfaces remained macroscopically rough. Hapke's full photometric model was then inverted using a genetic algorithm to obtain ϑ. However, independent measurements of ϑ were not performed in this study, so the inversion results could not be quantitatively validated. In Badura and Bachmann (2019), macroscopically rough beach sand surfaces were prepared in the laboratory, and ϑ was obtained both via inversion of Hapke's photometric model and via stereophotogrammetry. Even for the surfaces that met the underlying model assumptions, the model exhibited large errors (in excess of 40%) at high illumination zenith angles and high phase angles. In Labarre et al. (2017Labarre et al. ( , 2019, the macroscopic roughness model was applied to airborne and satellite imagery of terrestrial volcanic terrains. The model's performance was generally better at moderately high albedos and relatively smooth surfaces. However, a systematic underestimation of roughness was observed at higher values of ϑ. As noted by the authors, many of the natural surfaces that they examined, such as slabs, hollows, and fractures in lava fields, had complex structures that are not well modeled by the assumed azimuthally symmetric form of the probability density function (ϑ) in Hapke's model.
One of the challenges with assessing Hapke's roughness model is that the length scale and physical meaning of ϑ is ambiguous. In Helfenstein (1988), it is suggested that ϑ should be viewed as an integral quantity for roughness scales up to the resolution of the imaging system. Alternatively, in Shepard and Campbell (1998), it is suggested that the most important roughness scale is the smallest scale over which shadows still exist, i.e. over which shadows are not filled in by multifacet scattering from neighboring portions of the surface. Laboratory measurements conducted in Cord et al. (2003Cord et al. ( , 2005, analysis of the lunar surface in Helfenstein and Shepard (1999), and analysis of terrestrial volcanic terrains in Labarre et al. (2017Labarre et al. ( , 2019 suggest that the photometrically dominant roughness scale ranges anywhere from the sub-millimeter to centimeter scale. Due to the constraints of our laboratory setup, the surfaces measured in this paper include roughness at only a single macroscopic length scale: 4 mm correlation length. Therefore, we are not in a position to comment on whether ϑ should be interpreted at any particular length scale. We consider the proper length scale interpretation of ϑ of actual planetary bodies to be an open question. Rather, our goal is to compare our proposed model against Hapke's for a fixed, well controlled, macroscopic roughness scale to assess their relative performance. For roughness to be considered macroscopic, the roughness length scale must significantly exceed the particle size. This is the length scale for which Hapke originally derived his macroscopic roughness correction, and is therefore the scale at which it is most appropriate to test his model. Several authors have suggested that the primary shortcoming in Hapke's macroscopic roughness correction is the lack of multi-facet scattering, which is particularly important for rough, high albedo surfaces (Buratti and Veverka, 1985;Shepard and Campbell, 1998;Cord et al., 2003Cord et al., , 2004. Hapke himself noted this, and proposed a modification to the model in Hapke (2012b) wherein the effective photometric ϑ is decreased in proportion to the material's diffusive reflectance.
The assumption is that multi-facet scattering partially ''fills in'' the shadowed regions, creating the appearance of a smoother surface. However, to our knowledge, this modification was not based on any measurements and has never been tested. Our measurements indicate that the ''filling in'' of shadows is only part of the story. As will be shown, multi-facet scattering increases the scattered radiance from all facets (whether they are in shadow or not) in a way that is roughly Lambertian over a substantial portion of the hemisphere. This effect can be a significant contribution to the total scattered radiance even at nadir illumination, when there are no shadows at all. An additional drawback of Hapke's roughness correction is an analytical issue that seems to have gone unnoticed in the literature: the parameter identified as ''ϑ'' is not actually equal to the mean value of ϑ for the probability density function (ϑ) specified in his model. The two quantities differ by a factor of approximately π∕2 (see Sections 3 and 4.2). This issue may partially account for the disparity between the model and measurements in past studies.
In this paper, we perform a rigorous assessment of Hapke's macroscopic roughness correction (including his later modification for multifacet scattering) and propose an alternative that includes a more exact D.J. Shiltz and C.M. Bachmann  treatment of single-facet scattering and an empirical approximation of multi-facet scattering. We begin in Section 2 with the definition of bidirectional reflectance for a rough surface based on fundamental radiometry. In Section 3, we discuss Hapke's macroscopic roughness correction in detail, including the analytical issues associated with his choice of probability density function for the facet slope. In Section 4, we derive the correct form of the probability density function and use it to derive an expression for the single-facet scattering component of the rough surface bidirectional reflectance. We validate the singlefacet scattering component by performing Monte Carlo simulations. In Section 5, we present the details of our experimental measurements of rough mineral surface reflectance using 3D printed molds with known roughness statistics and a robotic spectro-gonio-radiometer. We derive an empirical model for multi-facet scattering from these measurements. In Section 6, we compare our proposed model with Hapke's and show in detail how the macroscopic roughness of the surface affects the bidirectional reflectance across the hemisphere. Section 7 provides concluding remarks.

Bidirectional reflectance of a rough surface
Consider a surface facet whose geometry is depicted in Fig. 1. Incoming collimated irradiance strikes the surface and radiance is scattered from the surface to a detector. Angles and are the zenith angles of the incidence and emergence directions with respect to the vertical̂, while angles ι and ε are the zenith angles with respect to the facet normal̂. The azimuth angle between the incidence and emergence directions is ψ. The zenith and azimuth angles of the facet normal are ϑ and φ respectively. These angles are related by the law of cosines for spherical triangles: cos ι = cos cos ϑ + sin sin ϑ cos φ (1a) cos ε = cos cos ϑ + sin sin ϑ cos(ψ − φ) The phase angle is the angle between the incidence and emergence directions, which is the same in both the untilted and tilted coordinate systems. It may be expressed as cos = cos cos + sin sin cos ψ The bidirectional reflectance of a surface, denoted by , is defined as the ratio of scattered radiance in the emergence direction̂to the incident collimated irradiance in the plane orthogonal to the direction of incidencê (Hapke, 2012b). For a tilted facet of differential area , we may write the bidirectional reflectance in the tilted coordinate system as where is the differential flux scattered by the facet and is the differential solid angle subtended by the detector. The tilted area may be written in terms of the underlying horizontal area in the -plane as Therefore, the total radiant intensity scattered by single-facet scattering events is where is the portion of the surface that is both illuminated and visible. The projected area subtended by the rough surface from the perspective of the detector is cos , and so the total radiance scattered by single-facet scattering events is In Hapke's original model for macroscopic roughness, only the radiance due to single-facet scattering events is accounted for. However, our measurements show that multi-facet scattering is a significant contributor to the total radiance, especially at high albedo and when the illumination direction is close to nadir. Let̃s ingle and̃m ulti denote the rough surface bidirectional reflectance due to single-facet scattering and multi-facet scattering, respectively. The complete bidirectional reflectance of the rough surface can therefore be defined as ( , , ) =̃s ingle +̃m ulti = 1 ∬ (ι, ε, ) cos ε cos sec ϑ +̃m ulti (7)

Hapke's macroscopic roughness correction
In his derivation, Hapke begins by introducing the following assumptions: 1. The roughness scale is large enough such that geometric optics is valid. 2. The rough surface is comprised of facets whose physical size far exceeds the mean particle size. 3. The slope distribution function (ϑ), which describes the probability density function of the facet normal zenith angle ϑ, is independent of azimuth φ. Surfaces with a preferred orientation are not considered. 4. The slope distribution function is assumed to have the form where the roughness parameter ϑ is defined by and the probability density function is normalized with 5. The parameter ϑ is assumed to be small. In particular, any Taylor expansion terms involving order ϑ 3 and higher are neglected.
Using the above assumptions, Hapke derives an approximate analytical expression for the bidirectional reflectance of a rough surface in terms of ϑ. In particular, the rough surface bidirectional reflectance is expressed as the smooth surface reflectance at an effective incidence angle and effective emergence angle multiplied by a shadowing term. His expression, neglecting multi-facet scattering, is given bỹ where effective angles and are functions of , , ψ, and ϑ as outlined in the Appendix. In an attempt to model the effects of multi-facet scattering, Hapke proposed the following modification: where 0 is the diffusive reflectance of the medium (Hapke, 2012b). His reasoning was that multi-facet scattering would effectively reduce the photometric effects of roughness, and that this reduction would depend non-linearly on the material's albedo. However, to our knowledge this modification was not based on any measurements and has never been tested until now. A link to our Python implementation of the original and modified model is included in Section ''Data Availability''. As noted in the introduction, there appear to be analytical issues with Hapke's roughness correction that seem to have gone unnoticed in the literature. As defined in Eq. (9), the parameter identified as ''ϑ'' is not actually equal to the mean value of a random variable ϑ whose probability density function is (ϑ). This can be easily shown by taking the limit of the distribution as the surface becomes arbitrarily smooth. As ϑ → 0, the distribution of ϑ also tends to zero, and so we have sin ϑ ≈ ϑ and sec ϑ ≈ 1. This gives which is a Rayleigh distribution with mean (π∕2) tan ϑ. This means that the true mean value of the distribution, which we denote as ⟨ϑ⟩, differs from ''ϑ'' by a factor of approximately π∕2. Therefore, research efforts that have tried to measure ϑ directly have potentially misapplied the model by assuming erroneously that ϑ in Hapke's model is equal to the mean value of the surface's slope angle. Furthermore, the normalization condition given by Eq. (10) is only approximately correct for small values of ϑ. As shown in Fig. 2, even for moderately rough surfaces, the area under (ϑ) is significantly less than unity and therefore does not meet the requirement of a probability density function.
In addition to the issues with the slope distribution function (ϑ), it is not clear why the roughness correction should take the form given by Eq. (11), i.e. that the rough surface reflectance should be the product of the smooth surface reflectance at an effective orientation and a shadowing function. As described in Section 2, the roughness of a surface causes facets at a variety of incidence and emergence angles to be presented to the illumination source and detector simultaneously. To find the radiance reaching the detector, the contribution of each facet needs to be accounted for at each tilt angle, weighted by how likely a facet is to have a particular tilt orientation, how much projected area that facet has with respect to the detector, and how likely a facet is to be shadowed or obscured by other facets. The derivation of such an expression is presented in the next section.

Proposed model derivation
In this section, a complete derivation of our proposed macroscopic roughness correction is presented. We begin by listing all of the simplifying assumptions. It should be noted that these assumptions are no more restrictive than the assumptions made in Hapke's correction. However, since our statistical model is described in terms of autocorrelation functions, some additional discussion of the assumptions and their implications is warranted. We make the following assumptions: 1. The surface can be described statistically by a wide sense stationary random process, i.e. the autocorrelation function is not spatially variant. Of course, in any remote sensing imagery, all aspects of the scene -such as albedo, surface roughness, chemical composition, grain size, etc. -are spatially variant. Therefore, when such a scene is imaged, it is understood that the radiance detected in each pixel represents the average contribution across the corresponding area on the surface. Therefore, we do not assume that the surface autocorrelation function is invariant across the entire scene, since this is never the case. However, we do assume that it is invariant within the area corresponding to a single pixel of the imaging system. 2. The surface autocorrelation function is not directionally dependent. This is a more restrictive assumption and is less likely to be satisfied in practice. Many natural surfaces have a ridged structure with a preferred orientation. Such surfaces cannot be described by an autocorrelation function that is the same in every direction. However, it should be noted that this assumption is equivalent to the one made by Hapke that (ϑ) is independent of azimuth angle φ. 3. The height of an arbitrary point on the surface is normally distributed. The reason for this assumption is that jointly normal random variables exhibit certain properties that simplify their analysis, as will be shown in this section. Virtually all of the roughness corrections and shadowing functions in the literature make this assumption for this reason. Studies of the lunar surface using stereophotogrammetry have indicated that the surface elevation is approximately normally distributed (Lumme et al., 1985;Helfenstein and Shepard, 1999). 4. The autocorrelation function is twice differentiable at the origin. As will be shown, this is equivalent to requiring that an RMS slope exists. It should be noted that the autocorrelation function need not take any particular form, so long as it is not directionally dependent. In practice, the autocorrelation function need not be measured directly anyway. Our proposed roughness correction will be shown to depend only on the RMS slope along a linear transect, which is a much easier quantity to measure. In general, the RMS slope of a surface depends on the length scale over which it is measured (Shepard and Campbell, 1998). In this paper, the correlation length of the rough surfaces is 4 mm, which greatly exceeds the particle size for the materials used in this study (see Section 5), and therefore meets the assumptions of the facet scale as originally described by Hapke.

Slope distribution function
Given this set of assumptions, our first task is to rigorously derive the probability density function of the facet slope. Consider a set of three points on a rough surface viewed from above as depicted by Fig. 3. Point 0 lies at the origin, point 1 lies along the -axis, and point 2 lies along the -axis (the positive -axis points out of the page). The azimuth angles of the illumination source and detector projected onto the -plane are denoted by ϕ and ϕ , respectively.
We model the heights of the three points ⃗ = ( 0 , 1 , 2 ) as jointly normally distributed random variables with zero mean, variance σ 2 , and autocorrelation function ρ( ). The horizontal distance from point 0 to point 1, and from point 0 to point 2, is . The horizontal distance from point 1 to point 2 is √ 2 . Therefore, for autocorrelation function ρ( ), the covariance matrix describing ⃗ is are the correlations between the three heights. The joint probability density function of ⃗ can therefore be written as Define the random variables = ( 1 − 0 )∕ and = ( 2 − 0 )∕ , which are the slopes of the surface along the positive -axis and -axis respectively. The joint probability density function of ⃗ = ( , ) can be found as This integral may be evaluated by computing the determinant and inverse of and employing the following identity: which yields the fact that and are jointly normally distributed with zero mean and covariance Next, we take the limit of ⃗ ( , ) as → 0. To do this, we apply L'Hopital's rule twice and utilize the following properties of the autocorrelation function ρ( ): where is the RMS slope along any transect of the surface (Smith, 1967). Using these substitutions, the joint distribution of and can be written as This indicates that the slopes along the and axes are independent normally distributed random variables with zero mean and variance 2 . The normal vector̂can be expressed in terms of and aŝ Taking the dot product of the normal vector and the vertical gives Next, define the random variable = 2 + 2 . To find the probability density function of , we first definê= ∕ and̂= ∕ and definê=̂2 +̂2. Sincêand̂are standard independent normally distributed random variables (with means of zero and variances of unity), the probability density function of̂is a chi-squared distribution with two degrees of freedom: We note that = 2̂, so the distribution ( ) can be found via the method of transformations and is expressed as Next, define the random variable = 1 / √ + 1 . Its probability density function can be found by employing the method of transformations a second time: We can express the facet normal zenith angle as ϑ = cos −1 . By employing the method of transformations one last time, we arrive at the desired probability density function:

Comparison with (ϑ)
Next, we compare the exact distribution ϑ (ϑ) derived above with the approximate distribution (ϑ) used by Hapke. To do this, we take the limit of both distributions as the surface becomes arbitrarily smooth. This gives which is a Rayleigh distribution with mean √ π∕2 . As was shown in Eq. (13), (ϑ) also approaches a Rayleigh distribution, with mean π∕2 tan ϑ. Since both (ϑ) and ϑ (ϑ) approach Rayleigh distributions in the limiting case of an arbitrarily smooth surface, we can find a relationship between Hapke's roughness parameter ϑ and the RMS slope by ensuring the distributions approach one another. Setting Eq.'s (13) and (28) equal to each other gives This relationship allows a direct comparison between the roughness correction described in this paper (given in terms of ) and Hapke's roughness correction (given in terms of his parameter ϑ). It should again be noted that the parameter identified by Hapke as ϑ is not equal to the mean value of ϑ as defined by (ϑ). Indeed, for small ϑ, we have ⟨ϑ⟩ ≈ (π∕2) tan ϑ as shown by Eq. (13). Researchers should be aware of this when comparing experimentally measured slopes with the predictions of Hapke's roughness correction. For this reason, we advocate using the RMS slope as the roughness parameter to avoid ambiguity in notation.
The difference between the exact distribution ϑ (ϑ) and Hapke's approximation (ϑ) is shown in Fig. 4. As seen in the figure, the approximation (ϑ) performs well for small values of but is less accurate for larger values of . As increases, and therefore as ϑ increases by Eq. (29), the area under (ϑ) decreases from unity, which is consistent with Fig. 2.

Proposed roughness correction
Now that the correct slope distribution function has been identified, we may utilize it to evaluate the rough surface bidirectional reflectance, which was defined by Eq. (7). However, the function ϑ (ϑ) will not be used directly. Rather, we will start with the joint probability density function for and (the slopes along the positive and axes) given by Eq. (21), and then express angles ϑ, ι, and ε as random variables defined in terms of and . As a first step, let us define and , the slope of the surface in the direction of the illumination source and detector, projected onto the -plane. These slopes are related to and by The three tilt angles can then be expressed as Since ϑ, ι, and ε are now modeled as random variables, we replace the deterministic area-averaged integral in Eq. (7) by the expected value: where (¬ ) is the probability that a facet with slopes and is not in shadow. Following Hapke, we consider two types of shadows: tilt shadows and projected shadows. A facet is in tilt shadow if its normal vector̂makes an angle of more than 90 • with the direction of incidencêor with the direction of emergencê. We may equivalently express this condition in terms of the slopes in the incidence and emergence directions: if > cot or > cot , then the facet is in tilt shadow. A facet is in projected shadow if some portion of the surface obstructs the view from the facet to the source or to the detector. Let us rewrite (¬ ) explicitly in terms of tilt and projected shadows: where ¬ denotes the absence of projected shadow and ¬ denotes the absence of tilt shadow. Following Hapke and several other authors, we assume that facets that are not in tilt shadow have a probability of being in projected shadow that is independent of their slopes (Saunders, 1967;Smith, 1967;Hapke, 1984;Heitz et al., 2013). This means that (¬ | ¬ ) can be taken outside the integral in Eq. (32). However, the probability of being in tilt shadow depends directly on the facet slope, so (¬ ) cannot be taken outside the integral. Instead, we can directly account for tilt shadowing by using step functions. Let ( ) denote the unit step function, i.e. Shiltz and C.M. Bachmann We can then define the probability of not being in tilt shadow as An estimate of (¬ | ¬ ), denoted bŷ(¬ | ¬ ), was presented in Heitz et al. (2013) and is adopted in this paper. Unlike other models in the literature, this model is fully bistatic and accounts for the correlation effect: namely the fact that when ψ is close to 0 • , and are strongly correlated. It should be noted that the original publication of the model contained a few typesetting errors. These errors have been corrected and the equations have been re-written below in terms of the notation used in this paper. Their model defines an estimate of projected shadowing via the term̂(¬ | ¬ ) according to: In Heitz et al. (2013), validation of Eq. (37) was performed using Monte Carlo simulation which showed very good agreement with their model. We have performed additional Monte Carlo simulations (presented in Section 5.3) which corroborate these findings. As a final simplification, we may utilize Eq. (31) to rewrite the trigonometric terms in the integral in Eq. (32) as: This expression has a clear physical intuition: it is the ratio between the projected area subtended by a tilted facet and the projected area that would be subtended if the facet were flat. In this paper we refer to it as the projected area factor. We can now write the single-facet scattering portion of our proposed roughness correction as witĥ(¬ | ¬ ) given by Eq. (37), angles ϑ, ι, and ε given in terms of and in Eqs. (30) and (31), (¬ ) given by Eq. (36), and ⃗ ( , ) given by Eq. (21). The multi-facet scattering term is discussed in the next section. Since our roughness correction involves computing an integral, some additional discussion on computational considerations is warranted. We begin by noting that all of the terms in the integrand in Eq. (39) are well behaved over all and , with no asymptotes or discontinuities (other than those introduced by the step functions, which are easily integrable). Since ⃗ ( , ) is a 2D Gaussian distribution centered at the origin, the non-zero extent of the integrand is limited by the standard deviation . Therefore, the lower and upper limits of integration may be set to some small multiple of while retaining very high accuracy. In all of the calculations presented in this paper, the integral in Eq. (39) was computed using limits of integration of , ∈ [−5 , 5 ], a uniformly spaced grid of length 100 in both and , and trapezoidal rule. With any modern computer and a precompiled trapezoidal rule implementation (such as that offered by Python's Scipy library (Virtanen et al., 2020)) the roughness correction may be quickly computed without issue. A link to our Python implementation of the model is included in Section ''Data Availability''.

Monte Carlo simulation
In this section, a numerical Monte Carlo simulation of rough surface scattering is presented. The purpose of the simulation is to validate the single-facet scattering portion of our proposed model,̃P roposed, single . For a rough surface that satisfies the assumptions presented in the beginning of this section (i.e. wide sense stationary, no preferred orientation, and normally distributed heights), the only approximation made in the single-facet scattering portion of the model is the treatment of projected shadows via the term̂(¬ | ¬ ). Everything else iñP roposed, single follows directly from the definition of the true single-facet reflectancẽs ingle derived from fundamental radiometry in Eq. (7).
Multi-facet scattering is not included in this Monte Carlo simulation; we address multi-facet scattering later in Section 5.4 when we consider the residual between our single-facet scattering model and measurements. To numerically model multi-facet scattering would require a fully three dimensional simulation with iterative ray tracing to simulate multiple orders of scattering, such as the model presented in Shkuratov et al. (2005). Such a model would certainly be useful to study multifacet scattering more closely, but its development is beyond the scope of this paper. Instead, we will rely on direct measurements of the rough surface reflectance to obtain the multi-facet scattering model. To construct the model, consider the set of discrete points depicted in Fig. 5. We denote as a randomly generated vector of length 2 + 2. Points 0 , 1 , and 2 +1 comprise a set of heights associated with the simulated facet. Points 0 , … , define a set of heights for a transect toward the projected direction of the illumination source, while points 0 , +1 , … , 2 represent a set of heights for a transect toward the projected direction of the detector. Let the distance between two adjacent points along a transect in the -plane be . The matrix of the horizontal distances in the -plane between all of these points may be written as Eq. (40) (Box I), where the matrix is given by For a Gaussian autocorrelation function with correlation length , the covariance may be written as The random vector is generated as a set of jointly normally distributed random variables with zero mean and covariance Σ. Let be a matrix in which each row is a particular realization of the random vector . Each particular realization of the surface is indexed by ∈ {0, 1, … , − 1}. The realizations of random variables ϑ, ι, and ε for surface realization are denoted by ϑ , ι , and ε . These angles are given by: In these expressions, the first index is the surface realization index , while the second index refers to the position of the point (0 through 2 + 1) within that surface. For surface , the facet is in tilt shadow if its local illumination or emergence direction is greater than 90 • . We may write this condition for surface as The facet is in projected shadow if either of the two rays from the origin to the illumination source or from the origin to the detector intersects the surface at some point. We may write this condition as Therefore, we may write the shadowing condition of surface realization as The exact expression for rough surface single-facet scattering̃s ingle was derived in Eq. (7). We simulate this result by replacing the integral over the surface with a discrete summation over the illuminated and visible facets. The Monte Carlo model (which includes single-facet scattering only) for the rough surface bidirectional reflectance is given bỹ In this paper, the Monte Carlo simulations were performed using correlation length = 1. The simulated transects had length 10 and discretization scale = ∕20. This gives a value of = 200. For each view geometry, the number of independent realizations of the surface was = 100, 000. A link to our Python implementation of the model is included in Section ''Data Availability''.

Sample preparation
We performed experimental validation of our proposed roughness correction using two minerals, quartz and olivine. Geotechnical properties of the two minerals are listed in Table 1. To study the radiometric effects of surface roughness, we 3D printed 30 cm diameter plastic molds with randomly generated surfaces. We generated each surface using correlation length = 4 mm and a Gaussian autocorrelation . For the Gaussian autocorrelation function, the RMS slope can be determined using Eq. (20) as We produced four different molds, each with a different RMS height σ : 0 mm (smooth), 0.5 mm, 0.75 mm, and 1.0 mm. These σ values correspond with values of 0, 0.177, 0.265, and 0.354 respectively. After 3D printing, we coated each mold with a two part epoxy resin to smooth out the striations left by the printing process. All mineral samples, including the smooth samples, were prepared using an identical process, as depicted in Fig. 6. The container is first filled upside down, with the 3D printed mold forming the base. Next, the container is sealed and placed on a shaker table. The purpose of this step is to ensure that all of the small gaps in the mold surfaces are completely filled with material. After experimenting with different time intervals, we determined that 5 min was long enough to ensure a well-prepared surface without stratifying the material grain size. Next, the upper half of the container is removed and excess material is scraped away. A base plate is then installed and the container is flipped vertically. Finally, the 3D printed mold is carefully released, revealing the prepared mineral surface.

Bidirectional reflectance measurement
To measure the bidirectional reflectance of the mineral samples across the hemisphere, we utilized the robotic spectro-gonio-radiometer system described in Harms et al. (2017). A photograph of the experimental setup is included in Fig. 6. This system utilizes an ASD FieldSpec FR4 spectro-radiometer at wavelengths ranging from 350 nm to 2500 nm. The robotic arm positions the downward-looking spectroradiometer's foreoptic at emergence angles ranging from nadir to 70 • , and azimuth angles from 0 • -360 • . A laser range finder ensures that the foreoptic is pointed at the same position on the sample throughout each scan. The illumination zenith angle can be varied from nadir to 70 • . The construction of the robotic arm allows radiance measurements to be performed at phase angles of 5 • and greater (smaller phase angles will result in the arm shadowing the sample). We used a 3 • field-of-view foreoptic with an optical scrambler for all measurements. The illumination source is a 750 W tungsten halogen source with an adjustable iris and focusing optics to collimate the light and provide a D.J. Shiltz and C.M. Bachmann Table 1 Geotechnical properties of minerals used in this experiment. Bulk densities and fill factors are measured for minerals prepared according to Fig. 6 uniform illumination field. With this light source, adequate signal-tonoise ratio (SNR) was obtained from 500 nm to 2250 nm for the two minerals.
One of the challenges with measuring the rough mineral samples is the fact that the surface is not uniform across the sample, as seen in Fig. 6. The surface height and slope are pseudorandom variables whose statistics are controlled by each mold's 3D printing process. Therefore, the radiance scattered by a given region of the surface is also a pseudorandom variable. To ensure a representative sample of radiance was obtained at each measurement geometry, we placed the mineral sample on a rotating stage whose center was horizontally offset by three inches from the foreoptic's field of view. The sample was then rotated continuously beneath the foreoptic during the measurement process. At each measurement position on the hemisphere, we recorded 20,000 independent radiance measurements over a duration of 10 min of rotation. This ensured a representative measurement of the radiance scattered by the pseudorandom surface.
To calibrate the measurements, after each scan we inserted a Spectralon™ plaque with NIST-traceable directional-hemispherical reflectance in place of the sample at the same distance from the illumination source and spectro-radiometer foreoptic, and we measured its reflected radiance under identical lighting conditions. Rather than assuming that the Spectralon™ plaque is Lambertian, we utilized an empirical model for the non-Lambertian bidirectional reflectance (Lévesque and Dissanska, 2016) which is accurate to ≈1% out to 70 • zenith in illumination and emergence. This uncertainty was combined with the uncertainty provided in the Spectralon™ calibration certificate, and was then propagated through the reflectance processing chain to provide an overall estimate of the uncertainty of each reflectance measurement. All of the error bars presented in this paper represent one standard deviation of uncertainty.

Smooth surface characterization
For the purpose of this experiment, it is of paramount importance that the smooth surface bidirectional reflectance ( , , ) be known with high accuracy across the entire hemisphere. Otherwise, it would be difficult to discern whether errors in the roughness correction models are due to the models themselves or to errors in the underlying smooth surface reflectance measurement. To this end, smooth samples of each mineral were measured from nadir to 70 • in illumination and emergence and from 0 • -360 • in azimuth, each in 10 • increments. Due to the constraints of the robotic arm and illumination setup, the remaining 20 • in illumination and emergence zenith angles between 70 • − 90 • could not be measured.
One option for filling in these remaining regions would be to fit an analytical bidirectional reflectance model to the measured data. Such models include Hapke's isotropic multiple scattering approximation (IMSA) (Hapke, 2012b) or its derivatives such as SOILSPECT (Jacquemoud et al., 1992). However, we found that neither of these models could fit the measured bidirectional reflectance over the entire hemisphere to the desired accuracy. Since an accurate characterization of the smooth surface bidirectional reflectance is critical to this study, another approach was needed. Therefore, instead of utilizing an analytical model, we estimated the unmeasured portions of the hemisphere by extrapolating the measurements with low order polynomial fits. We then used this extrapolated data set to construct a radial basis function (RBF) interpolation so that ( , , ) could be estimated at any angle. In particular, for every ψ ∈ [0 measurements at = 50 • , 60 • , and 70 • . Finally, the extrapolated estimates were smoothed along each dimension ( , , and ψ) by fitting a cubic polynomial and replacing the estimates with the polynomial interpolation. It should be noted that this smoothing operation was applied only to the extrapolated points, not to the measurements. The original measurements were not adjusted or smoothed in any way.
The results of the extrapolation/interpolation process for the two minerals with illumination zenith angle = 30 • are shown in Fig. 7. As described above, the reflectance at ≤ 70 • is measured directly, while the reflectance at = 80 • and = 90 • is extrapolated from the measurements at each value of ψ. We see in the figure that this procedure results in a physically reasonable bidirectional reflectance profile. No obvious visible artifacts have been introduced by the extrapolation process. As seen in the figure, the two minerals exhibit very different scattering behavior: the olivine is heavily biased toward forward scattering, whereas the quartz has a much more pronounced opposition effect.

Diffusive reflectance estimate
To test Hapke's modification for multi-facet scattering, the diffusive reflectance of each mineral must be estimated at each wavelength. For isotropic scatterers, the diffusive reflectance is defined as where γ = √ 1 − is the albedo factor and is the single scattering albedo of the material (Hapke, 2012b). For anisotropic scatterers, the diffusive reflectance can be written in terms of an effective single scattering albedo * , given by where β is the hemispherical asymmetry factor (Hapke, 2012b). In particular, the effective albedo factor is defined as and the diffusive reflectance for the general case of anisotropic scatterers can be written as The hemispherical asymmetry factor is a property of the effective single particle phase function of the medium. It is defined as the average cosine of the scattering angle: π − . For the particular case of a two parameter Henyey-Greenstein single particle phase function, β = − (Hapke, 2012a). Therefore, to determine the wavelength dependent diffusive reflectance in the general case of anisotropic scatterers, the values of , , and must be estimated at each wavelength.
where is the Ambartsumian-Chandrasekhar function and ( ) is the single particle phase function (Hapke, 2012b). For the two parameter Henyey-Greenstein function, the single particle phase function is written as To perform the inversion, we performed a least squares fit of Eq. (53) to the measured bidirectional reflectance using the differential evolution algorithm (Storn and Price, 1997). The results of the inversion are shown in Fig. 8. It should be noted that in this study, IMSA was used only to estimate the diffusive reflectance by inverting , , and . The IMSA inversion was not accurate enough to model the smooth surface bidirectional reflectance ( , , ) across the entire hemisphere to the accuracy needed for this study, even when the opposition effect terms were included. Instead, the interpolation/extrapolation process detailed in the last section was used to characterize ( , , ).

Validation of single-facet scattering model via Monte Carlo simulation
Next, we used the Monte Carlo simulation presented in Section 4.4 to validate the single-facet scattering portion of our proposed model. To perform the validation, for each rough surface reflectance measurement geometry, we compared the model̃P roposed, single with a corresponding Monte Carlo simulatioñM onte Carlo, single . As seen in Fig. 9, there is very good agreement between the single-facet scattering portion of our proposed model and the Monte Carlo simulation. The only approximation made in the single-facet scattering portion of our model is the treatment of projected shadowing via the term̂(¬ | ¬ ) given in Heitz et al. (2013). Our results shown in Fig. 9 appear to corroborate their earlier findings that̂(¬ | ¬ ) is indeed a good model for the projected shadowing, at least for the roughness scales and view geometries studied in this paper.

Multi-facet scattering
Given the very good agreement between the single-facet scattering portion of our model and the Monte Carlo simulation, it is reasonable to assume that the residual between the measured rough surface D.J. Shiltz and C.M. Bachmann Fig. 8. Inversion of single scattering albedo , phase function parameters and , and the resulting hemispherical asymmetry factor β, effective albedo * , and diffusive reflectance 0 .  Fig. 10. Measurement residual̃for < 90 • along with our Lambertian approximation as a function of * , cos , and while holding the other variables constant. Olivine measurements are represented by green dots and quartz measurements by blue dots. In the left figure, the residual is plotted for = 30 • and = 3 . In the center figure, the residual is plotted for = 3 . In the right figure, the residual is plotted for = 30 • . The best fit value for the fitting parameter is L = 0.19. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 11. Measurement residual̃for all along with our non-Lambertian approximation. Olivine measurements are represented by green dots and quartz measurements by blue dots. The best fit value for the non-Lambertian fitting constant is NL = 6.5. The multi-facet scattering is approximately Lambertian for phase angles < 90 • but gradually increases at higher phase angles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) reflectance and the single-facet scattering model is due primarily to multi-facet scattering. In our measurements, we observed that in almost every case, the measured reflectance was greater than what was predicted by the single-facet scattering model. This observation is consistent with the hypothesis that the missing radiance due to multi-facet scattering is the primary factor affecting the residual.
We found that this residual was approximately Lambertian for phase angles < 90 • . Therefore, we chose to approximate the residual multifacet scattering by a Lambertian surface with an effective Lambert albedo that depends on the material's single scattering albedo as well as the roughness of the surface. Our goal then was to find an approximation for this effective Lambert albedo that satisfies physical intuition and also agrees with measurements.
We define the residual as We observed that the residual is well correlated with the diffusive reflectance 0 as defined in the previous section. Furthermore, we observed that the residual is weakly correlated with the RMS slope . Therefore, our proposed Lambertian multi-facet scattering approximation has the form Proposed, multi, Lambertian = L 0 cos π where L is a constant fitting parameter for the Lambertian multifacet scattering approximation. In this expression, the term cos ∕ π represents the bidirectional reflectance of a perfectly white diffuse surface, and the term L 0 represents the effective Lambert albedo of this surface. Fig. 10 displays the residual as a function of * , cos , and separately for phase angles < 90 • . The solid black line in the figure corresponds with our proposed Lambertian multifacet scattering approximation for the best fit value L = 0.19. As expected, the multi-facet scattering contribution is greater at higher albedo and at higher roughness, which was anticipated in Buratti and Veverka (1985), Shepard and Campbell (1998) and Hapke (2012b). In these references, it was assumed that the primary effect of multi-facet scattering is the filling in of shadowed regions, creating the appearance of a smoother surface. However, our measurements suggest that multifacet scattering is most significant at nadir illumination, i.e. when there are no shadows at all. The most likely explanation is that multi-facet scattering increases the scattered radiance from all facets, whether they are in shadow or not, due to scattering from adjacent facets within their fields of view.
For phase angles greater than 90 • , the Lambertian multi-facet scattering approximation becomes less accurate. In particular, our measurements indicate that the multi-facet scattering becomes biased toward the forward direction for phase angles greater than 90 • . To account for this, we modify the Lambertian approximation with a Gaussian function centered at = π with width π∕2 and a non-Lambertian fitting constant NL : Fig. 11 shows the residual normalized by the Lambertian approximation, where the black line indicates the approximate non-Lambertian correction term for the best fit value of NL = 6.5. While the multi-facet scattering contribution can be roughly approximated by the non-Lambertian model given in Eq. (57), there is still significant variance in the data, especially with regard to the roughness parameter . As shown in Fig. 12, the correlation between the empirical multi-facet scattering model and the measured residual is relatively low ( 2 = 0.38) when compared with the correlation between the single-facet scattering model and the Monte Carlo simulation ( 2 = 0.9998). However, the overall magnitude of the multi-facet scattering is much lower (<10% on average) than the single-facet scattering part, D.J. Shiltz and C.M. Bachmann Fig. 14. Bidirectional reflectance measurements of olivine at 1100 nm ( 0 = 0.29) and = 30 • at various RMS slope values , along with comparisons with our proposed roughness correction and Hapke's single-facet and modified roughness corrections. even for the relatively high albedo minerals used in this study. At lower albedos, the contribution from multi-facet scattering is expected to be even lower.

Final results
Now that both the single-facet and multi-facet portions of our proposed roughness correction model have been specified and validated, the complete model may be written as In this section, we compare the performance of our proposed roughness correction against Hapke's. We begin by examining high-resolution reflectance measurements of the rough mineral surfaces at 1100 nm while illuminated at = 30 • as shown in Fig.'s 13 and 14. For purposes of illustration, the wavelength of 1100 nm was chosen because it is the wavelength at which the difference between the diffusive reflectance of the two minerals is greatest: for quartz, 0 = 0.91 and for olivine, We begin with the quartz (Fig. 13). It is seen that the increase in macroscopic roughness causes a shift in the overall scattering behavior. For the smooth sample, a prominent forward peak and SHOE are both clearly visible (top row of Fig. 13). However, as the roughness increases, the sample becomes progressively more backscattering. For the roughest sample (bottom row of Fig. 13), the forward peak has disappeared entirely. The disappearance of the forward peak can be explained by the effect of shadowing: when viewing the sample in the forward direction, a significant portion of the visible area is in shadow. The broad backward peak is caused by two related phenomena. The first effect may be called macroscopic shadow hiding. When viewing the sample in the backward direction, a significant portion of the shadowed facets are hidden from view. In the case that ψ = 0 • and > , no shadows will be visible to the detector at all. The second effect is due to the preferential orientation of the visible facets toward the detector at these angles. In the backward direction, the visible facets are illuminated closer to their local nadir directions, causing them to appear brighter than they would be if the surface were smooth.
We refer to the combined effect of both of these phenomena as the macroscopic roughness backscattering bias (MRBB). In remote sensing applications, the effect of the MRBB is that a macroscopically rough surface appears to be more strongly backscattering than it would if it were smooth. It should be noted that the photometric effect of the MRBB is qualitatively different than that caused by the SHOE. The SHOE is D.J. Shiltz and C.M. Bachmann  an opposition effect appearing only at small phase angles due to selfshadowing of individual particles and/or clumps of particles (Hapke, 2012b). The MRBB is not an opposition effect, but rather an overall bias toward backward scattering that affects the entire hemisphere. In the case of a very rough surface, a wide backward lobe is observed that is not localized to the opposition direction, as seen in the bottom row of Fig. 13. The effect is particularly pronounced when ψ = 0 and > , since no macroscopic shadows are visible at these angles.
We see that our proposed model (second column) correctly predicts the MRBB effect and the overall shape of the reflectance distribution for each roughness level. We see that Hapke's original single-facet scattering model (third column) does correctly predict a shift toward backward scattering, but it underpredicts the overall reflectance in every direction. This is not surprising, since his original model does not include multi-facet scattering. We also see that the forward lobe is incorrectly preserved, even for the roughest sample. Finally, we see that Hapke's modified roughness model (fourth column) predicts that the reflectance will remain nearly unchanged by the roughness. This is because the diffusive reflectance is high ( 0 = 0.91) at this wavelength, and his modification reduces the roughness parameter ϑ by multiplying it by 1 − 0 , as described by Eq. (12). The net result is that the photometric effect of the roughness has been almost completely neglected by the model.
Next, we examine the olivine in Fig. 14. We again see the MRBB effect, and this time it is even more dramatic. As seen in the first column, the reflectance distribution goes from heavily forward scattering to heavily backward scattering as a result of the roughness. Our proposed model (second column) correctly predicts this behavior. Hapke's single-facet model (third column) also correctly predicts a shift toward backward scattering, but the reflectance is once again underpredicted in every direction. Hapke's modified model (fourth column) does improve this underprediction somewhat. The diffusive D.J. Shiltz and C.M. Bachmann  reflectance of olivine at 1100 nm is 0 = 0.29, and so the effect of the modification in Eq. (12) is to multiply ϑ by 0.71, which in turn increases the reflectance across the hemisphere. However, the MRBB effect predicted in Hapke's modified model is not as pronounced as the measurements indicate it should be.
To more closely examine the performance of our proposed model at a variety of incidence angles and view geometries, we plotted the bidirectional reflectance along transects of the hemisphere, i.e. for fixed ψ, for the two minerals at λ = 1100 nm and using the roughest mold ( = 0.354).
For the scans at = 10 • (Fig's. 15 and 16), our proposed correction and Hapke's correction correctly predict the overall distribution of reflectance. However, we see that Hapke's single-facet roughness correction underestimates the reflectance in every direction. Once again, this is not surprising, since his model does not include multi-facet scattering. As expected, the absolute error is greater for the quartz than for the olivine, since quartz has a higher albedo at this wavelength. Hapke's modified model accounts for this underestimation to some degree. However, it appears to overestimate for the quartz and underestimate for the olivine. Our proposed model correctly accounts for the multi-facet scattering contribution at both albedos. It should also be noted that the difference between the Lambertian and non-Lambertian multi-facet scattering approximations is minimal in these geometries because the measured phase angle does not exceed 80 • .
For the scans at = 30 • (Fig's. 17 and 18), we begin to see qualitative differences in the predicted scattering behavior between the three models. For instance, for the quartz, Hapke's single-facet model predicts a dip in reflectance in the backward direction for > , whereas the measurements and our proposed model do not. While his model underpredicts the reflectance over much of the hemisphere, this underprediction is less pronounced in the forward direction at high phase angles. Hapke's modified model does increase the reflectance D.J. Shiltz and C.M. Bachmann  with respect to his single-facet model, as expected. However, the shape of the distribution is not correct. It still underestimates the reflectance over much of the hemisphere, but overpredicts the reflectance in the forward direction for > 40 • . This may be due to the fact that his modification reduces the effective ϑ value, and so the MRBB effect is not properly modeled. For these geometries, we also begin to see the difference between the Lambertian and non-Lambertian multi-facet scattering approximations at higher phase angles.
Finally, for the scans at = 60 • (Fig's. 19 and 20), we see that while Hapke's single-facet model performs the worst of the three models over much of the hemisphere, it performs surprisingly well in the forward direction. In fact, in the forward direction it out-performs the Monte Carlo single-facet scattering simulation, which is unexpected. For a sufficiently large number of simulations, the Monte Carlo simulation approaches the exact solution for single-facet scattering. Considering the fact that Hapke's single-facet roughness correction is an attempt to approximate the single-facet scattering solution, one would expect that the performance would be worse than that of the Monte Carlo simulation. It appears that Hapke's correction overestimates the singlefacet scattering reflectance in the forward direction for these angles, but that this overestimation corresponds approximately with the effect of multi-facet scattering. Therefore, his modified model greatly overestimates the reflectance in the forward direction, especially for the quartz. Meanwhile, we see that at these higher phase angles, our non-Lambertian multi-facet scattering approximation performs better than the Lambertian one ( NL = 0), such that our proposed model agrees with our measurements across the hemisphere.
To further quantify the performance of the three roughness corrections, we computed the percent error in the reflectance predicted by each model at incidence zenith angles ranging from = 10 • to = 60 • , emergence zenith angles ranging from = 0 • to = 70 • , and azimuth angles ψ = 0 • , 60 • , 120 • , and 180 • at 1751 wavelength bands ranging D.J. Shiltz and C.M. Bachmann  from 500 nm to 2250 nm, for both the olivine and quartz, and for each of the three roughness molds. To see how the model performs across the spectral range of our instrument, we plotted the percent error as a function of wavelength, averaged over all illumination and view geometries, for each mineral as shown in Fig. 21. We see that our proposed model has the lowest reflectance error across this spectral range. The standard deviation of the reflectance error across the different view geometries and roughness levels (indicated by the error bars) is also lower for our proposed model. Since Hapke's single-facet scattering model does not include multi-facet scattering, the reflectance is underestimated across the entire spectrum. It appears that his modified model does correct for this underestimation to an extent, however the variance is greater than that of our proposed model.
We then plotted the root-mean-square error (RMSE) as a function of phase angle , as shown in Fig. 22. As seen in the figure, the RMSE of our proposed roughness correction is an improvement over Hapke's original single-facet scattering model for phase angles less than 90 • , and the performance of the two models is comparable above this threshold. Due to the constraints of our laboratory setup, the models could not be evaluated at phase angles greater than 130 • , and so the performance of the two models in this regime is uncertain. We also see that the performance of Hapke's modified model is superior to the single-facet model for smaller phase angles, but the error becomes much larger at higher phase angles. This can be explained by the fact that the model incorrectly neglects the effects of roughness at high 0 (as was seen in Fig. 13), and this effect is most pronounced at high phase angles. One could theoretically choose to utilize a hybrid model: i.e. use Hapke's modified model for ≤ 80 • and use the singlefacet model for > 80 • . However, if one chooses this strategy, one should ensure that the parameter ϑ is interpreted correctly. As discussed previously, ϑ as defined in Hapke's model is not equal to the mean value of ϑ. Rather, ''ϑ'' as used in Hapke's model should be interpreted in terms of the RMS slope as described by Eq. (29). We note that our proposed model has lower RMSE across the entire range of measured phase angles and also correctly predicts the MRBB effect shown in Fig's. 13 and 14. Based on our measurements, it is expected that our proposed model will perform as well or better than either of Hapke's models in virtually every situation.

Conclusion
In this paper, we derived an alternative to Hapke's macroscopic roughness correction and evaluated the models using rough mineral samples prepared in the laboratory. Our model subdivides rough surface scattering into single-facet scattering and multi-facet scattering. Our single-facet scattering model was derived directly from the definition of bidirectional reflectance of a rough surface. We incorporated a result published by Heitz et al. (2013) to approximate the effect of projected shadowing. Through a Monte Carlo simulation for single-facet scattering, we evaluated our proposed single-facet scattering model, finding good agreement. The simulations demonstrated that the singlefacet scattering component of our proposed model is nearly an exact solution for single-facet scattering. By fitting the residual between the bidirectional reflectance measured in the laboratory and the prediction of the single-facet scattering model, we developed an empirical approximation for multi-facet scattering. We found that the contribution of multi-facet scattering is approximately Lambertian for phase angles less than 90 • , but that for larger phase angles, a non-Lambertian, forward scattering model is more appropriate. We evaluated our proposed model using rough mineral surfaces for phase angles up to 130 • , surface RMS slopes up to 0.354, and single scattering albedos between 0.86 D.J. Shiltz and C.M. Bachmann Fig. 20. Bidirectional reflectance measurements of olivine at = 60 • , λ = 1100 nm, = 0.354, and comparisons with our proposed roughness correction, Hapke's roughness correction, and the single-facet scattering Monte Carlo simulation. and 0.99. Under these conditions, our proposed model demonstrated improved accuracy over Hapke's single-facet roughness correction, as well as his modified roughness correction.
In the course of developing and evaluating our proposed roughness correction model, we identified several important features of rough surface scattering. It has previously been suggested in the literature that the primary effect of multi-facet scattering is the filling in of shadowed regions of the surface, creating the appearance of a smoother surface. However, we demonstrated that, in fact, multi-facet scattering increases the scattered radiance in all directions, even when no shadows are visible from a given view direction. Indeed, the contribution from multi-facet scattering is greatest at nadir illumination, when there are no shadows at all. We also observed that increasing the roughness of a surface produces a more strongly backscattering reflectance distribution due to macroscopic shadow hiding and the preferential orientation of the visible facets when viewed from the backward direction. We refer to this effect as the macroscopic roughness backscattering bias (MRBB) which we showed is distinct from the shadow hiding opposition effect (SHOE).
Although our proposed roughness correction offers an improvement over Hapke's correction over the conditions we studied, there are several areas that warrant future research. As discussed previously, the single-facet scattering component of our proposed model is nearly an exact solution for single-facet scattering. Therefore, the error in the model is due almost entirely to multi-facet scattering, which is  more difficult to model. In this paper, we presented a relatively simple empirical approximation for multi-facet scattering. However, the underlying physical mechanisms behind multi-facet scattering are still not well understood. In particular, the forward scattering behavior that is seen at higher phase angles warrants further investigation. Another possible avenue for future research is to evaluate our proposed model on surfaces with anisotropic roughness statistics, or with statistics that are otherwise not well modeled by the form of ⃗ ( , ) given by D.J. Shiltz and C.M. Bachmann Eq. (21). Our proposed model may be easily adapted to such surfaces by replacing ⃗ ( , ) in Eq. (58) with any desired distribution, rather than the isotropic one given by Eq. (21).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The laboratory data and Python code used to implement these roughness correction models can be found at https://github.com/ djshiltz/macroscopic-roughness-correction. DOI: 10.5281/zenodo.7020667.