## 1.

## Introduction

In March 2015, a joint field experiment and airborne measurement campaign was conducted in the Algodones Dunes, California. The Algodones Dunes has been selected by NASA as a potential site for satellite intercalibration, and the overall goal of the experiment was to better characterize the spectral variability of the dunes. The experiment was a multiinstitutional effort involving NASA, South Dakota State University, the University of Arizona, the University of Lethbridge, and the Rochester Institute of Technology (RIT). These dunes have significant topographic relief. This fact, combined with natural variations anticipated within a large dune system, such as Algodones, made the direct measurement of the bidirectional reflectance distribution function (BRDF), which characterizes the angular variation in scattered light from the surface for a given illumination geometry, of paramount importance.

Within the context of the larger team effort, our goals were to provide representative hyperspectral BRDF measurements at a number of sites within the dune system and to characterize the geophysical properties of the surface that might influence observed variations in BRDF at these locations and throughout the dune system. With a longer-term view of improving representative BRDF of the Algodones Dunes system, we captured sediment samples from multiple sites within the dune system, both where our BRDF field measurements occurred as well as at other sites that represented surface variations in geophysical properties of the dunes where it had not been possible to bring our field goniometer system within the time constraints of the field experiment. Our plan was to bring these samples back to the laboratory and deliberately manipulate their geophysical properties to mimic variations that might be expected within the Algodones Dunes system. Thus, we report here our laboratory BRDF measurements in which we have manipulated individual geophysical variables such as sediment density and grain size of the sediment samples as well as the illumination geometry. We describe the accuracy of using these BRDF libraries from manipulated specimens in comparison with NASA Goddard’s light detection and ranging (LiDAR), hyperspectral, and thermal (G-LiHT)^{1} imagery. We also report on our efforts to invert radiative transfer models based on the work of Hapke to directly estimate sediment fill factor from the observed BRDF.

Representative BRDF libraries play a critical role in calibration of both airborne and satellite systems. For example, these type of libraries have been used within scene simulation models, such as DIRSIG,^{2} which can predict the expected scene top-of-atmosphere radiance. DIRSIG currently incorporates standard BRDF models, such as the Ross–Li^{3} model. While models such as the Ross–Li model have been widely used in NASA satellite data products, they have a limited number of available parameters to model variations in expected BRDF of the surface. Another goal of this work, therefore, is to undertake analysis of more complex models, such as those developed by Hapke,^{4} for future inclusion in DIRSIG and use in satellite and airborne sensor calibration.

The rest of the paper is organized as follows: Sec. 2 outlines methods used to acquire data in the field and in the laboratory where samples were manipulated. In this section, we also outline approaches to radiative transfer modeling and inversion based on Hapke’s model.^{4} In Sec. 3, we report comparisons between the spectral libraries acquired in the field and developed in the laboratory to airborne G-LiHT imagery acquired over the dune system. In this section, we also examine the results of our approaches to geophysical parameter inversion. These approaches rely on a comparison of BRDF data acquired at two different illumination geometries using the invariance of the single scattering properties of the surface. In this context, we consider the problem of inversion of the Hapke model and a variant of this model to obtain the fill factor, and compare to laboratory measurements of density, which are expected to correlate with the fill factor. Finally, in Sec. 4, we provide a summary and conclusions and describe the next steps needed to eventually transition the approach to retrieval of the fill factor directly from the imagery time series.

## 2.

## Methods

## 2.1.

### Field Measurements of the Angular Dependence of Hyperspectral Reflectance at the Algodones Dunes

During the March 2015 campaign, our team collected BRDF data at multiple sites within the dune system (Fig. 1). Due to time limitations, our measurements were collected from six locations using the goniometer of the Rochester Institute of Technology (GRIT).^{5}^{,}^{6} The GRIT positioned the fiber optic of an ASD spectrometer with 5 deg fore-optic over a set of hemispherical positions between 0 deg and 65 deg in zenith and 0 deg and 360 deg in azimuth. Although the aperture size of the sensor is relatively small, measurements such as these are sometimes referred to as hemispherical conical reflectance factor (HCRF) measurements.^{7}^{,}^{8} The use of the term conical reflects the finite size of the sensing aperture, and “hemispherical” indicates that illumination comes from both direct and diffuse sources in an outdoor setting. Our ASD FR4 spectrometer measurements onboard the GRIT used a 5-deg fore-optic.

Due to the limited number of days available for measurement, we also collected sediment samples at a number of sites within the dune system to capture the typical variation in geophysical parameters. These samples were brought back to the lab where additional testing, described below, was performed to manipulate the geophysical properties of the sample to characterize the potential impact on the HCRF and to compare with radiative transfer models. Relevant parameters include the grain size distribution, sediment density, sediment composition, moisture content, surface roughness, and local topography. Radiative transfer models depend on parameters such as these, and, for example, the fill factor of the sediment (closely related to density) directly parametrizes approximate solutions to the radiative transfer equation (RTE),^{4} which we will review in greater detail below.

Field spectral measurements were accompanied by a set of geotechnical measurements, which included sediment density, grain size distribution, and moisture content, as well as mechanical properties of the sediment, such as shear strength. For example, sediment density was measured with a sand cone density apparatus (Fig. 1). Samples collected on the dune were also analyzed to determine moisture content and grain size distributions using, respectively, a drying oven and a mechanical sieve shaker, both brought to the site in our mobile laboratory. To determine grain size distributions, sediments were dried in the oven prior to sorting into size fractions with the mechanical sieve shaker. At sites where GRIT HCRF measurements were made, all geotechnical measurements were completed. Additional samples were collected from other areas to provide examples of sediment variation found across the dune system. These were returned to our laboratory for use in laboratory-based experiments using GRIT (Fig. 1) and later a second-generation hyperspectral goniometer system known as GRIT-two (GRIT-T), shown in Fig. 5.^{9} In these experiments, individual variables such as grain size distribution and density were manipulated, along with the illumination geometry, and used in a series of laboratory biconical reflectance factor (BCRF)^{7}^{,}^{8}^{,}^{10} measurements.

The primary sites where field HCRF measurements were taken with GRIT in tandem with geotechnical measurements were on the western side of the Algodones Dunes system (Fig. 1). Additional sites where samples were collected for later laboratory use were on the northern edge, central-northern region, and central western edge of the dune system, as shown in Fig. 1. In each case, efforts were made to collect samples from the top, bottom, and middle of individual dunes. Grain size and density variations were observed across sites and even across positions within an individual dune. The next section describes the geotechnical measurements conducted in tandem with the GRIT hyperspectral HCRF measurements.

## 2.2.

### Geotechnical Measurements in the Algodones Dunes

Geophysical properties vary considerably across the dune system. At the sites where we measured HCRF, as well as at those sites where samples were collected for later laboratory analysis, grain size distributions show some significant differences. The average size factor, $X=\frac{\pi D}{\lambda}$, where $D$ represents the average particle diameter and $\lambda $ is the wavelength of light of an image band, varied by as much as a factor of four across all collected samples, while density varied by as much as 43% (Fig. 2). These density measurements were conducted in the field with a sand cone apparatus according to an American Society for Testing and Materials (ASTM) standard.^{11} Grain size distributions were determined onsite by drying the sediments in a drying oven over a 24-h period, and then sieving the sediment in a Humboldt mechanical sieve shaker for a 30-min period. Moisture content of the sediment was determined using the weights of the sample before and after drying.

## 2.3.

### Laboratory Measurements of the Angular Dependence of Hyperspectral Reflectance: Controlling Geophysical Parameters

In a laboratory setting, we conducted a series of BRDF measurements using sediments brought back from the Algodones Dunes site. Since grain size distribution and density both play an important role in the observed BRDF, measurements focused on manipulating these individual parameters. In addition, the illumination geometry was also changed to provide a wide variety of overall measurement scenarios.

To achieve samples of varying density, samples were pluviated using an apparatus in our laboratory following an ASTM standard (Fig. 3).^{11} At the top of the pluviation apparatus, a canister with a trapdoor releases sediment, which falls through a fixed set of coarse sieves placed at a specific drop height above the sample holder. The process is designed to mimic Aeolian deposition. In many of the experiments, we prepared samples over a wide range of densities obtained by varying the drop height from the coarse sieves to the sample holder in the pluviation device.

For the laboratory studies, BRDF measurements were conducted using two different generations of the goniometer system, GRIT^{5}^{,}^{6} and GRIT-T.^{9} In a laboratory setting, BRDF measurements are sometimes referred to as BCRF^{7}^{,}^{8}^{,}^{10} measurements when there is a finite aperture for both the single directional illumination source and the sensor. The term “factor” signifies the referencing of sample radiance data to the radiance of a Lambertian standard. In our experiments, Spectralon^{®} panels, which approximate a Lambertian surface,^{12} were used as the reference standard. These panels were calibrated by LabSphere. Figure 4 shows some examples of BCRF of Algodones Dunes sediment samples for varying density.

## 2.4.

### Radiative Transfer Models for Retrieval of Sediment Geophysical Properties: Fill Factor

Absorption and scattering from sediment surfaces have been modeled with radiative transfer models in both planetary astronomy and Earth remote sensing applications,^{4}^{,}^{13}14.15.^{–}^{16} and Hapke’s approximate solution to the RTEs has been widely used.^{17}18.19.20.21.22.23.24.25.^{–}^{26} The isotropic multiple scattering approximation (IMSA) is an approximate solution to an RTE based on the Chandresekhar–Ambartsumian method of invariance.^{4} In the Hapke treatment of the RTE, five orders of scattering are considered, and the single scattering contribution through the single scattering phase function is treated in an exact form; however, for multiple scattering contributions, the phase function is assumed to be isotropic.^{4} The resulting approximate solution (IMSA) has the form

## (1)

$$r({\theta}_{i},{\theta}_{e},g)=K\frac{w(\lambda )}{4}\frac{1}{{\mu}_{i}+{\mu}_{e}}(p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]+\{H[\frac{{\mu}_{i}}{K},w(\lambda )]H[\frac{{\mu}_{e}}{K},w(\lambda )]-1\left\}\right)[1+{B}_{C0}{B}_{C}(g,K,\lambda )].$$^{4}

^{,}

^{15}

## (2)

$$H[\frac{{\mu}_{i}}{K},w(\lambda )]=\frac{1}{1-w(\lambda )\frac{{\mu}_{i}}{K}[{r}_{0}+\frac{1-2{r}_{0}\frac{{\mu}_{i}}{K}}{2}\mathrm{ln}(\frac{1+\frac{{\mu}_{i}}{K}}{\frac{{\mu}_{i}}{K}}\left)\right]},$$^{4}

^{,}

^{14}In this expression, ${B}_{S0}$ is a constant that must be determined in an inversion procedure. The functional form of the opposition effect ${B}_{S}(g,K,\lambda )$ is discussed further below. The term ${B}_{C0}{B}_{C}(g,K,\lambda )$ describes the coherent backscatter opposition effect (CBOE), which is observed over only a very narrow range of phase angles, typically $\lesssim 2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$, where coherent addition of scattering pathways in equal and opposite directions can occur.

^{4}

^{,}

^{15}${B}_{C0}$ is another free constant that must in principle be determined, while ${B}_{C}(g,K,\lambda )$ describes its functional form, which, similarly to the SHOE, is large at small phase angles and decreases with increasing phase angle. In this paper, we neglect the CBOE since our measurements do not cover phase angles this small. Note that in Eq. (1), factors of ${\mu}_{i}$ and $1/\pi $ are removed since our measurements involve reflectance factors (BCRF and HCRF), which are the ratio of radiance values obtained from the sediment sample and a Spectralon

^{®}white reference, an approximately Lambertian standard.

^{4}

^{,}

^{7}

^{,}

^{8}

^{,}

^{10}

^{,}

^{12}

Importantly, a geophysical parameter, the sediment fill factor $\varphi $ (the fraction of volume occupied by particles), explicitly parameterizes the approximate IMSA solution since the nonlinear porosity factor $K=K(\varphi )$ is a function of the fill factor $\varphi $, and, for equant particles, it has the form^{4}^{,}^{16}

## (3)

$$K\approx -\frac{\mathrm{ln}(1-1.209\text{\hspace{0.17em}}{\varphi}^{2/3})}{1.209\text{\hspace{0.17em}}{\varphi}^{2/3}}.$$Thus, inversion of the IMSA solution provides an explicit estimate of an important geophysical parameter, the sediment fill factor. In this work, our goal was to obtain a relative fill factor through inversion of the radiative transfer model. We wanted to demonstrate that the resulting retrieved fill factor correlates consistently with the directly measured sediment density measured in the laboratory, even if we are not sure of its absolute magnitude.

A second geophysical property that is also implicitly present in the IMSA solution is the sediment grain size distribution. Within the IMSA solution, this plays a role in the term related to the SHOE. Hapke has derived an approximate form for the SHOE, which takes the form^{4}^{,}^{14}

^{4}One such example is a unimodal case, having the form $r{e}^{{-}^{r/\overline{r}}}$, where $r$ is the particle radius and $\overline{r}$ is the average radius. This grain size distribution leads to a SHOE width parameter, ${h}_{S\text{\hspace{0.17em}\hspace{0.17em}}}={(\frac{3}{8})}^{\frac{3}{2}}K(\varphi )\varphi $.

^{4}

^{,}

^{27}Grain size distributions have been the focus of earlier studies, one of which included a retrieval of a weighted average grain size from a desert site.

^{28}

The actual shape of distributions found in nature may be more complicated. Figure 2 shows several examples of the grain size distributions derived from our sediment samples acquired in the Algodones Dunes. Several of these show grain size distributions with more than one mode. In many of the examples shown in the figure, the secondary mode is a smaller peak of fine sediments, such as fine sand and/or silt. However, there is at least one example from the sediment samples that has a complicated set of multiple peaks (sample 1202-M-02-B). For sample 1202-M-02-B, the full distribution extends beyond the range shown in Fig. 2, and this is addressed later in our discussion of laboratory experiments.

While a significant portion of the distributions in the samples was unimodal, many were not, and since multimodal distributions are quite common, we have modeled the width of the SHOE peak in this study as being proportional to $K(\varphi )\varphi $, with the proportionality constant being a free parameter $\epsilon $ that must be optimized in our inversion of the IMSA solution

As each free parameter is introduced into an inversion procedure, the optimization process grows progressively more complicated. Thus, adding the free parameter $\epsilon $, while logical from a modeling perspective, adds to the complexity of the model residual surface being searched, especially for a model as nonlinear as IMSA. In order to avoid this pitfall, we have taken a different approach, building upon an earlier workflow that we previously developed.^{6}We have added some important refinements, which allow the approach to be more successful. Another variation of this approach, which we will introduce later, takes advantage of the explicit form of the SHOE width found in Eq. (5).

In that initial exploration,^{6} we inverted the Hapke IMSA solution solving Eq. (1) for the product of the single particle phase function and the SHOE correction

## (6)

$$p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]=\frac{4}{Kw(\lambda )}({\mu}_{i}+{\mu}_{e})r({\theta}_{i},{\theta}_{e},g)-\left\{H\right[\frac{{\mu}_{i}}{K},w(\lambda )\left]H\right[\frac{{\mu}_{e}}{K},w(\lambda )]-1\}.$$## (7)

$${\mathrm{min}}_{\varphi ,w(\lambda )}({\{{p}_{1}(g,\lambda )[1+{B}_{S0}{B}_{S,1}(g,K,\lambda )]-{p}_{2}(g,\lambda )[1+{B}_{S0}{B}_{S,2}(g,K,\lambda )]\}}^{2}).$$## (8)

$${\mathrm{min}}_{\varphi ,w(\lambda )}[{\{{p}_{1}(g,\lambda )[1+{B}_{S0}{B}_{S,1}(g,K,\lambda )]-{p}_{2}(g,\lambda )[1+{B}_{S0}{B}_{S,2}(g,K,\lambda )]\}}^{2}\phantom{\rule{0ex}{0ex}}+\alpha w{(\lambda )}^{2}+\beta K{(\varphi )}^{2}].$$One subtlety of the minimization approach, described in Eqs. (7) and (8), is that the set of sample phase angles in the two illumination conditions is not guaranteed to be the same, and we elaborate here on some of the technical details required for this to successfully converge to a meaningful solution. An overall flow diagram of the processing steps is shown in Fig. 6.

In the inversion process, the two sets of BCRF scans must be interpolated to a common phase angle grid to facilitate the difference calculation specified in Eq. (8). Although we tried a number of different interpolation schemes, we found that Steffen’s algorithm^{29} provided a stable, well-behaved interpolation of the BCRF data. Additional refinements, which decrease the overall residual, include the use of a low-pass filter (Sovitsky–Golay filter)^{30} during the interpolation onto the common phase angle grid. In the optimization procedure, we typically used a Nelder–Mead simplex method,^{31} although we also evaluated the Levenberg–Marquardt algorithm,^{32} which did not produce appreciably different results.

A new refinement, which we introduce in this paper, involves an additional constraint that improves the convergence of the model to a consistent solution. This constraint is motivated by the goal that the resulting fill factor retrieved should be consistent across wavelengths. This is, of course, a physical requirement of our approach, but without insisting on this constraint, the processing steps described in Eqs. (6) and (8) and illustrated in Fig. 5 are essentially wavelength agnostic. The system of equations can be applied independently to each wavelength, and therefore there is a potential to obtain solutions that do not report a consistent fill factor $\varphi $ across all wavelengths. In contrast, we expect to obtain a unique value of $w(\lambda )$ for each wavelength $\lambda $. To ensure this constraint, the Nelder–Mead optimization step is set up with an additional constraint. We want to constrain the individual estimates obtained for the fill factor at a particular wavelength ${\varphi}_{\lambda}$ so that they do not differ from the average

by more than a specified amount $\delta $ To implement this constraint on the optimization procedure, we must first run the optimization over all wavelengths $\lambda $ [Eqs. (8) and (6)] without the constraint found in Eq. (10). We then repeat the same steps, but this time the Nelder–Mead simplex optimization step enforces the constraint of Eq. (10). We can reiterate this constraint if desired and let the size of $\delta $ be given by $\delta =\delta (t)$, where $\delta (t)$ is a decreasing function of the iteration step $t$.We also explored another variation of this inversion process based on an earlier proposed modification to the Hapke model, in which we reintroduced directional dependence into the multiple scattering term.^{5} The inclusion of this factor was not derived in a rigorous manner; however, by insisting that the single particle phase function and SHOE correction be inserted in front of the multiple scattering term, this approach can be thought of ensuring that the final scattering event in a multiple scattering pathway has the directional dependence of a single particle, while the other scattering events are assumed to average out and not contribute to the directional dependence. The form used in the earlier work was^{5}

## (11)

$$p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]=\frac{4}{Kw(\lambda )}({\mu}_{i}+{\mu}_{e})r({\theta}_{i},{\theta}_{e},g)-{\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}.}\left[H\right(\frac{{\mu}_{i}}{K}\left)H\right(\frac{{\mu}_{e}}{K})-1].$$^{5}then the optimization of Eq. (8) proceeded with ${\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}}$ being used on the right-hand side of Eq. (11). In this paper, however, we use a more refined version of this approach that considers some of the underlying parameters that determine the form of $p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]$ explicitly in separate optimization steps. ${\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}}$ in Eq. (11) is also replaced by $\eta p(g,\lambda )$ in the present version, which requires us to separate out the functional form of $p(g,\lambda )$ as part of our processing chain. The processing steps that we describe are outlined in Fig. 7.

The initial steps of this inversion process follow a similar approach outlined in Eqs. (6) and (8). This provides us with an initial estimate of the fill factor $\varphi (t=0)$ and single scattering albedo $w(\lambda ,t=0)$. These values are then passed through Eq. (6) to obtain a first estimate for both illumination conditions: ${\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}.,\mathrm{1,2}}$. Here the subscripts “1” and “2” refer to the two illumination conditions. These estimates are interpolated onto the common phase angle grid during each step of the initial optimization as before, and then with the initial estimates obtained, $\varphi (t=0)$ and $w(\lambda ,t=0)$, a separate optimization step is used to find an explicit parameterization of ${\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}.,\mathrm{1,2}}$ in terms of the free parameters ${B}_{S0}$, $\epsilon $ [Eq. (5)], and the coefficients of the three-parameter Henyey–Greenstein function ${b}_{1},{b}_{2}$, and $c$^{4}

## (12)

$${p}_{HG}(g,\lambda )=\frac{1+c}{2}\frac{1-{b}_{1}^{2}}{{(1-2{b}_{1}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}g+{b}_{1}^{2})}^{3/2}}+\frac{1-c}{2}\frac{1-{b}_{2}^{2}}{{(1+2{b}_{2}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}g+{b}_{2}^{2})}^{3/2}}.$$## (13)

$${\mathrm{min}}_{{B}_{S0,\epsilon ,{b}_{1},{b}_{2},c}}[{{(\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}.,1}-{\{{p}_{HG}(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{fwd}})}^{2}+{{(\{p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{est}.,2}-{\{{p}_{HG}(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]\}}_{\mathrm{fwd}})}^{2}].$$## (14)

$$p(g,\lambda )[1+{B}_{S0}{B}_{S}(g,K,\lambda )]=\frac{4}{Kw(\lambda )}({\mu}_{i}+{\mu}_{e})+\eta p(g,\lambda )\left[H\right(\frac{{\mu}_{i}}{K}\left)H\right(\frac{{\mu}_{e}}{K})-1].$$## 3.

## Results

## 3.1.

### Field HCRF and Laboratory BCRF Measurements: Impact of Geophysical Parameters

Above, we noted that density and grain size distributions can both impact measured BCRF and HCRF data. In previous studies, we have observed strong density-dependent effects, and the trend for these effects has been observed to depend somewhat on the relative optical contrast of the underlying mixture of materials in the granular sediment in combination with the illumination zenith angle, which determines the degree of multiple scatter within the medium, and thus the potential for darker optically contrasting minerals to diminish multiple scatter.^{5}^{,}^{6}^{,}^{33} Grain size effects also play an important role, since they are not entirely independent from another free variable that affects BCRF and HCRF, namely surface roughness. Although not modeled in the present work, Hapke has also developed a correction to the IMSA model for macroscopic surface roughness.^{4}^{,}^{34} These corrections are applicable to the majority of the grain sizes found in our distributions, which all fall within the geometric optics region, with the exception of very fine clay particles, which typically have size $\lesssim 2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. In the laboratory studies in which grain size was explicitly manipulated as an independent variable, this directly impacted the apparent surface roughness. Algodones sediment sample 1102-MS-02 from the north-central region of the dune system (Fig. 1) was prepared in the laboratory, where different subsets of the grain size distributions were separated through mechanical sieve shaking (Fig. 3). Sieved fractions were included or removed, and these new samples with manipulated grain size distributions (Fig. 8) formed the basis of several experiments in which density was also manipulated through pluviation (Fig. 3). Grain size manipulations resulted in the apparent differences in surface roughness (Fig. 8).

In general, altering the grain size distribution to include significant gaps (case 3 of Fig. 8) or measuring from the original distribution that already had significant gaps and therefore surface roughness tended to produce results showing the greatest change in the angular distribution when different illumination zenith angles were chosen. We have observed trends like this in the past. Samples with the most pronounced roughness (case 3 in Fig. 8) tend to exhibit the greatest backscatter when the light source is more oblique (the 65-deg illumination case, Fig. 9, top row, last column), similar to what we have observed in previous studies.^{27} When illumination angles are close to nadir, the amount of diffuse multiple scatter often significantly increases (second row of Fig. 9), as observed in previous studies.^{5}

Density of granular sediments also plays an important role in observed BCRF and HCRF.^{4}^{,}^{5}^{,}^{33} In general, reflectance tends to increase with density within certain conditions and only up to a maximum density. In our past studies, we have observed that the degree of multiple scatter present and the degree of optical contrast of constituents both play a role in the observed trend.^{33} The trend of increasing reflectance with density tends to be for cases when there is less multiple scatter, as may be expected often at more oblique illumination angles.^{5} At angles close to nadir, where greater multiple scatter may be present, if the sediment mixture contains fractions that are smaller and darker minerals, then increasing the density may lead to a decrease in overall reflectance^{5} since densifying the material will tend to force the smaller, darker particles to more optimally occupy the available pore space, shutting down multiple scattering pathways.^{33} When illumination is closer to nadir and therefore more likely to produce diffuse multiple scatter, this leads to a decrease in reflectance.^{5} At oblique angles with less multiple scatter, if the trend does follow the predictions of the Hapke model, then the reflectance should increase up to a critical maximum at which point particles are so close together that they begin to behave like larger single particles, and reflectance tends to decrease.^{4} We have observed this trend in the past under these types of conditions.^{5} The overall variance between the spectra may vary with illumination angle and thus degree of multiple scatter. In Fig. 9, we see greater diffuse multiple scatter when the illumination zenith angle is closer to nadir at 25 deg compared to the oblique case of 65 deg, and in Fig. 10, we see that the spectral plots of the full set of HCRF measurements for the 20-deg zenith angle illumination also have greater overall spectral variance compared to the more oblique case at 40 deg.

The relative amount of forward to backward scatter may also change with density,^{5} and we see this to some extent in Fig. 10, row 3, where the extent of the backward scatter and relative amount of the backscatter compared to the amount of forward scatter changes with density for the 450-nm wavelength. At the longer SWIR wavelength (fourth row of Fig. 10), the main difference with increasing density is a more gradual change in the extent of the backscatter region (especially along the azimuthal dimension), which grows larger with increasing density. These effects with density may be stronger when the illumination zenith angle is even larger, as we have observed in earlier work.^{5}

## 3.2.

### Matching Field and Laboratory HCRF and BCRF to G-LiHT Hyperspectral Imagery

A spectral library was created from all the field hyperspectral HCRF measurements recorded with GRIT during the March 2015 field experiment at Algodones as well as the hyperspectral BCRF (using GRIT and GRIT-T) of samples brought back from the dunes and prepared as described earlier. The total size of the spectral library generated for this study included 13,519 reflectance spectra. Although to date, the laboratory experiments have not exhausted all the collected sediment samples, the laboratory BCRF makes an important contribution to the representation of the spectral characteristics of the dune system (Figs. 11 and 12). The figures portray a Euclidean spectral match between the spectral library and the G-LiHT imagery after the library has been resampled to the spectral grid of the G-LiHT imagery. The figures show that manipulating geophysical parameters of field samples brought back to the lab can obtain better matches over significant parts of the imagery acquired with G-LiHT during the March 2015 field campaign. In the sequences of laboratory experiments, we systematically explored varying the sediment density, grain size distribution, and illumination zenith angle of the source relative to the sediment surface. All of these libraries match parts of the G-LiHT imagery better than field measurements in at least some part of the image subsets shown in these figures, and especially in Fig. 12, where roughly half of the scene is better matched by the laboratory-manipulated spectral BCRF libraries. Note that in Fig. 12, the manipulation of grain size to prepare three very different types of distributions followed by pluviation to produce different densities had a highly significant role in the spectral matching shown in Fig. 12. Likewise, density effects are highlighted at the top of the dune ridges in this scene subset. This is consistent with what we observed in our field geotechnical data, which showed differences in density between the bottom and top of a dune, for example, between sites BC01-01 and BC01-02 (Fig. 2). Note that in both cases in Figs. 11 and 12, none of the GRIT field measurement locations were within the boundaries of the scene.

These results are significant because it means that laboratory manipulations of field samples can be used to successfully augment the representative spectral library for the field site. Since the Algodones Dunes is envisioned as an ongoing vicarious calibration site, this means that scene simulation packages such as DIRSIG will have more reliable spectra as input.

## 3.3.

### Radiative Transfer Model Inversion Using Laboratory Biconical Reflectance Factor Measurements

We used the optimization procedure for inverting the original Hapke model described in Sec. 2.4. For this particular analysis, a set of five densities was created for the sample from site 0903-T-B (Fig. 1) using the air pluviation technique previously described (Fig. 3). These densities ranged from 1.59 to $1.82\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{g}/{\mathrm{cm}}^{2}$. The retrieved fill factor was then compared against the measured density of the material in the sample holder. The regression between retrieved fill factor and measured density is shown in Fig. 13 for each of the two approaches that we described in Sec. 2.4.

We found that the first iteration of the constraint equation provided the best-correlated estimate with direct laboratory measurements of sediment density. At this stage of processing, the VNIR and SWIR parts of the spectrum give slightly different estimates of the fill factor, and it is clear that we should interpret the results only as a measure of the relative fill factor rather than an absolute measure of this parameter. Correlations of the retrieved fill factor with directly measured density were comparable for the VNIR, with ${R}_{\mathrm{VNIR}}^{2}=0.73$ for the inversion of the original IMSA model and ${R}_{\mathrm{VNIR}}^{2}=0.71$ for the modified model, which includes an angular dependence in multiple scattering. In the SWIR, the estimates were less correlated with directly measured density, with ${R}_{\mathrm{SWIR}}^{2}=0.50$ for inversion of the original IMSA model and ${\mathrm{R}}_{\mathrm{SWIR}}^{2}=0.60$ for the modified model. The weaker correlation in the SWIR is likely due to lower SNR in this spectral region, although that could be improved by recording a larger number of spectral samples at each position. However, it is clear that the modified model was better correlated in this region of the spectrum. In both cases, the estimates of the single scattering albedo (Fig. 14) appear consistent with each other for all wavelengths and across all densities. In only one case, for the density of $1.79\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{g}/{\mathrm{cm}}^{2}$, the solution obtained by the original IMSA model exhibited more noise (a spike in the spectral shape $\sim 700\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$) in the VNIR, which is not observed in the modified model.

## 4.

## Conclusions

We have seen that field-measured HCRF with the GRIT hyperspectral goniometer system is well matched to the hyperspectral imagery collected during the Algodones Dunes campaign. As additional samples brought back to the laboratory have been studied and manipulated for BCRF studies with GRIT and GRIT-T, we have seen that they provide a better match to different parts of the dune system where GRIT measurements were not obtained while on site. It is an important point that laboratory measurements, which manipulate the geophysical parameters, such as density and grain size distribution, can provide important data that augment field BRDF spectral libraries to make the library more representative of the dune system in general. Because it is likely that NASA will use the Algodones Dunes as an ongoing vicarious calibration site for satellite systems, these observations are particularly important. This means that more representative data can be created to augment costly field campaigns and be useful in scene simulation/prediction packages such as DIRSIG, which are important for calibration and rely on representative BRDF as a foundational input to the overall processing chain. In addition, the data collected during this campaign have also served the dual purpose of allowing us to develop and test radiative transfer models that are more sophisticated than those being used, such as Ross–Li, in packages, such as DIRSIG.

Beyond the significance for calibration, radiative transfer models like those developed by Hapke (and variants described here) are important for geophysical parameter inversion from imagery. In this work, we used the samples brought back from Algodones to develop and test the idea of using BCRF measurements from two illumination geometries to invert an important geophysical parameter of the radiative transfer model, namely the sediment fill factor, which is ubiquitous in the approximate solution to the radiative transfer equation developed by Hapke. The laboratory-based results that were used in the inversion allowed us to demonstrate an inversion of the fill factor that strongly correlates with measured sediment density, an observation that suggests the validity of the inversion process. This does not confirm the absolute scale of the fill factor obtained; however, it does confirm the expected trend of relative fill factor with respect to directly measured density. In the future, a source of potential improvement to our approach would be to include a more sophisticated constraint in the spectral domain, as even in imposing the relatively weak constraint across wavelength used in this paper, convergence of the inversion and consistency across wavelengths was greatly improved.

To adapt this approach that we have taken for geophysical parameter inversion to a practical strategy where this can be accomplished for imagery, such as G-LiHT, there are still some technical hurdles to overcome, but the initial results are promising. The additional hurdles include the fact that the number of relative geometries over which the surface is both illuminated and viewed for a given position on the ground will depend on factors, for example, such as satellite revisit schedule or length of campaign and design of imaging flight patterns and times of day for airborne studies. The success of the approach that we described depended on the ability to interpolate an estimate of the single scattering phase function over the same phase angles with SHOE correction and a comparison of a densely sampled set of views over a hemisphere for two different illumination cases. In an imagery time series, the combination of illumination direction and view direction in a time series may be more complex; however, one point is clear: the better represented the range of phase angles is in the imagery time series, the more likely the interpolation and comparison can be successful. In the G-LiHT imagery flown over Algodones, sets of lines were flown along a particular heading, with successive lines having a significant percentage of overlap. Since the lines were short, this meant that at least several looks at the same point on the ground occurred from multiple view directions with roughly similar illumination conditions. On successive days, other orientations were flown, providing a diversity of potential geometries, and AM and PM lines were flown on each flight day so that the illumination conditions also varied. Thus, comparison of pixels from the same point on the ground among the various flight lines could be used to compare and interpolate the single scattering properties needed in our approach for two or more different illumination geometries. Likewise, in a satellite time series with frequent enough revisit time to ensure sites are sufficiently stable during the acquisition of the time series, a similar strategy might be followed. The essential point of interpolating an estimate of the single scattering properties from observations and modeled multiple scattering appears to be feasible.

## References

## Biography

**Charles M. Bachmann** received his AB from Princeton University (1984) and his ScM (1986) and PhD (1990) from Brown University. He was a research physicist (1990–2013) at the Naval Research Laboratory, serving as head of the coastal science and interpretation section in the remote sensing division (1993–2013). He has been Frederick and Anna B. Wiedman Chair at the RIT Chester F. Carlson Center for Imaging Science since 2013. His research emphasizes coastal hyperspectral remote sensing.

**Rehman Eon** received his BSc degree in mathematical physics and chemistry from Viterbo University in 2015. He is currently pursuing a PhD in imaging science at the Rochester Institute of Technology. His research interests include the use of optical remote sensing for the assessment of earth sediments and vegetation, radiative transfer modeling, and sensor calibration of Earth-observing systems.