<figure><img src="image_1.png" /></figure>

HAL
open science

# Mapping the surface properties of the Asal-Ghoubbet rift by massive inversion of the Hapke model on Pleiades multiangular images

D T Nguyen, S Jacquemoud, A Lucas, S. Douté, C Ferrari, Sophie Coustance, Sebastien Marcq, A Meygret

► To cite this version:

D T Nguyen, S Jacquemoud, A Lucas, S. Douté, C Ferrari, et al.. Mapping the surface properties of the Asal-Ghoubbet rift by massive inversion of the Hapke model on Pleiades multiangular images. Remote Sensing of Environment, 2025, 322, 10.1016/j.rse.2025.114691. hal-04989770

HAL Id: hal-04989770

https://hal.science/hal-04989770v1

Submitted on 13 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Mapping the surface properties of the Asal-Ghoubbet rift by massive inversion of the Hapke model on Pleiades multiangular images

D.T. Nguyen<sup>a</sup>, S. Jacquemoud<sup>a,*</sup>, A. Lucas<sup>a</sup>, S. Douté<sup>b</sup>, C. Ferrari<sup>a</sup>, S. Coustance<sup>c</sup>, S. Marcq<sup>c</sup>, A. Meygret<sup>c</sup>

## ARTICLE INFO

Editor: Jing M. Chen

Keywords:

Hapke photometric model
Surface roughness
Multiangular optical remote sensing
Pleiades satellite
Massive Bayesian inversion
Asal-Ghoubbet rift
Wavelet analysis

## ABSTRACT

A massive inversion of the Hapke model is carried out over the Asal-Ghoubbet rift (Republic of Djibouti) using high-resolution multiangular Pleiades images. This is the first time that such an inversion is performed on Earth over an entire image, previous studies having focused on planetary surfaces. This work addresses challenges such as atmospheric and geometrical corrections of these images to produce parameter maps. The use of fast Bayesian inversion significantly reduces computation times thanks to efficient exploration of the parameter space and leads to improved prediction. The parameters of the Hapke model are also interpreted in terms of surface physical properties, thanks to field measurements. Single scattering albedo is the parameter extracted with the greatest reliability, although its validation is still difficult due to the absence of a simple formula linking it to surface reflectance. Our study reveals a close relationship between photometric roughness and single scattering albedo, indicating that accurate extraction of the former is highly dependent on values of the latter, which must be below 0.8 for reliable estimation. Finally, the correlation between phase function parameters and grain properties depends on surface type and material properties.

## 1. Introduction

Understanding the origin and evolution of land surfaces depends on our ability to accurately determine their physical properties and mineralogical composition. Although in situ measurements provide highly accurate information, their availability is limited, so we must rely on satellite imagery. The mineralogical composition of surfaces is typically extracted from hyperspectral data, using the spectral signatures of the minerals that make up the surface. Physical properties, on the other hand, are often derived from texture analysis of satellite imagery. They can reveal patterns, orientations, and distributions of materials that indicate various surface-forming processes such as volcanic activity, sediment deposition or erosion. However, the effectiveness of texture analysis is inherently limited by the resolution of satellite sensors, which cannot capture details as small as grain size or shape, critical parameters for understanding surface properties.

Photometry is a complementary method to the above techniques that provides valuable information on the smallest scales of surface texture

and properties. The bidirectional reflectance distribution function (BRDF) describes the scattering of a parallel light beam incident in one direction towards another direction in the hemisphere (Nicodemus, 1965; Schaepman-Strub et al., 2006). It captures the intrinsic scattering properties of a surface, enabling other relevant quantities to be derived. However, because it is a ratio of infinitesimal quantities, the BRDF cannot be measured directly and must therefore be determined using the bidirectional reflectance factor (BRF). The inversion of a radiative transfer model on the BRDF/BRF allows to derive the physical properties (surface roughness, grain size and shape, porosity) and, to a lesser extent, the chemical properties (mineralogical composition, moisture content, organic matter content) of the surface.

In recent decades, such an approach has therefore been used in Earth observation and planetary science to analyze surface materials. On Earth, it is mainly applied to data acquired in the laboratory (e.g., Jacquemoud et al., 1992; Lucey, 1998; Helfenstein and Shepard, 2011; Souchon et al., 2011; Pilonget et al., 2016) or in the field (e.g., Pinty et al., 1989; Yang et al., 2011; Bachmann et al., 2017; Eon et al., 2018),

<footnote><sup>a</sup> Institut de physique du globe de Paris, Université Paris Cité, UMR CNRS 7154, 1 rue Jussieu, 75005 Paris, France</footnote>

<footnote><sup>b</sup> Institut de planétologie et d'astrophysique de Grenoble, OSUG-A, 414 rue de la Piscine, Domaine Universitaire, 38400 Saint-Martin-d'Hères, France</footnote>

<footnote><sup>c</sup> Centre National d'Étude Spatiales, 18 avenue Edouard Belin, 31401 Toulouse, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: tridnguyen@ipgp.fr (D.T. Nguyen), jacquemoud@ipgp.fr (S. Jacquemoud), lucas@ipgp.fr (A. Lucas), sylvain.doute@univ-grenoble-alpes.fr (S. Douté), ferrari@ipgp.fr (C. Ferrari), Sophie.Coustance@cnes.fr (S. Coustance), Sebastien.Marcq@cnes.fr (S. Marcq), Aime.Meygret@cnes.fr (A. Meygret).</footnote>
much more rarely to satellite data (Wu et al., 2009; Labarre et al., 2019; Ding et al., 2022). Inversion of BRDF models to retrieve surface physical properties on satellite images presents both challenges and opportunities. In the simplest case, researchers manually select pixels on an image, extract their BRDF, then run inversion schemes to infer model parameters. This approach has been implemented by Labarre et al. (2019) on multiangular Pleiades images (GSD 2 m) acquired over the Asal-Ghoubbet rift and by Ding et al. (2022) on multiangular POLDER images (GSD 6 km × 7 km). Inversion of the Hapke model over extended areas has only been tested on the Moon (Sato et al., 2014) and Mars (Fernando et al., 2016; Kugler et al., 2022). Several reasons have limited the inversion of radiative transfer models on multi-angular images on Earth: (1) Insufficient angular coverage or inadequate geometric configurations, impacting model selection and parameter constraints; (2) The interdependence of model parameters, leading to local minima during inversion; (3) A lack of field measurements for validation; (4) Atmospheric effects affecting the photometric response of surfaces; (5) Important computational resources required to process millions, if not billions, of pixels in several spectral bands. A massive image-scale inversion would provide much wider coverage and a detailed representation of land units.

The aim of this article is to explore the relationship between parameters estimated by massive inversion of the Hapke model and various physical (surface roughness, grain size and shape) and chemical (mineralogical composition) properties of natural surfaces. To this end, we use a unique dataset acquired in the Asal-Ghoubbet rift (Republic of Djibouti) as a case study. While our dataset is specific to this region, the workflow we develop is not constrained to this context and can be applied to any multi-angular remote sensing data, making it accessible for a wide range of users. This includes UAV-acquired imagery (Li et al., 2021; Cao et al., 2023), small satellite formations (Nag et al., 2017) or future missions such as the CO3D satellite (Lebègue et al., 2020), planed for launch in 2025, which are likely to facilitate the acquisition of BRDFs from space. We will first describe the study site (Section 2) and the experimental datasets (Section 3). Next, we introduce the image processing and inversion methods (Section 4). Finally, we compare the estimated parameters with field measurements and discuss their significance (Section 5).

## 2. Region of interest: Asal-Ghoubbet rift

The Asal-Ghoubbet rift (Republic of Djibouti) is a unique natural laboratory for geological and environmental research. It extends for 30 km between Lake Asal, located to the northwest of the rift, and Ghoubbet al-Kharab (or Lake Ghoubbet), a bay located to the southeast of the rift and separated from the Gulf of Tadjoura by a narrow channel (Fig. 1). Lake Asal that lies 155 m below sea level is characterized by thick salt deposits. The rift itself is a dynamic zone, with Quaternary sediments and lava flows associated with significant volcanic and tectonic activity. In particular, the central Fieale volcano, which erupted about 100,000 years ago, and the Ardoukoba volcano, which erupted in 1978, have played a key role in shaping the current topography (Pinzuti et al., 2013; Smittarello et al., 2016).

The Asal-Ghoubbet rift has been the subject of intense research interest, supported by extensive geodetic, topographic and geological surveys, as well as satellite imagery (Stieltjes, 1980; Pinzuti et al., 2013; Smittarello et al., 2016). It has a unique radiometric diversity, ranging from the bright white of saline deposits to the deep black of recent basaltic lava flows, and a wide variety of topographic features, from coarse lava flows to flat sedimentary basins separated by near-vertical faults. These features, along with the hot semi-arid climate of this part of the world, make it an ideal location for remote sensing studies of bare soils by. Due to the near-zero soil moisture and extremely sparse vegetation, photometric signals originate primarily from the soil, making it possible to study its physical properties.

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Asal-Ghoubbet rift (Republic of Djibouti) (11.43°N-11.74°N, 42.34°-42.65°E): a) Pleiades composite image acquired with a solar zenith angle of 37.8°, a viewing zenith angle of 24°, and a relative azimuth angle of 152.8°. The ground sampling distance (GSD) is 2 m. Red dots indicate the location of the study sites, which are shown in yellow letters (see Section 3.2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3. Data collection

We have divided our datasets into two categories: 1) Multiangular Pleiades images covering the entire Asal-Ghoubbet rift, which allow mapping of soil physical properties by massive inversion of the Hapke model; 2) Field data corresponding to twenty study sites in the rift; they include images acquired by drone and by hand, in order to build milli-metric digital elevation models (DEM), as well as soil and rock samples collected in the field.

### 3.1. Pleiades images

Pleiades is a high-resolution optical imaging constellation launched by the Centre national d'études spatiales (CNES) in late 2011 and 2012. It consists of two satellites and can provide daily coverage of any part of the Earth. The two satellites have a panchromatic channel (GSD = 50 cm), and four spectral bands from the visible to the near-infrared (GSD = 2 m). The radiometric accuracy is estimated to be around 4 %. The great agility of the platform allows different acquisition modes (stereo pairs, triplets) to facilitate the construction of metric DEMs using photogrammetry (Gleyzes et al., 2012). An acquisition mode called VIDEO was programmed by CNES over the Asal-Ghoubbet rift on January 26, 2013, allowing 21 images to be taken in a single pass at angles ranging from 55.73° forward to 56.69° backward, with a measurement at 9.73° near nadir (Fig. 2). During the four-minute satellite flyby, the solar zenith angle changed from 37.97° to 37.37°, a variation of only 0.6°. The corresponding phase angles range from 30° to 90°. These values will constrain the version of the Hapke used and the quality of the parameters estimated by the inversion.

The preprocessing of the Pleiades images includes atmospheric corrections, which were performed by the Physics for Optical Remote Sensing Measurement Service at CNES. They were performed using MAJA (MACCS ATCOR Joint Algorithm) (Hagolle et al., 2015; Rouquié et al., 2017), an advanced multi-mission atmospheric correction processor designed to process time series of high-resolution images, such as Sentinel-2 images. The fifth generation ECMWF atmospheric reanalysis
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Polar plot showing the geometric configuration of the Pleiades acquisitions: sun directions (red dots) and viewing directions (blue dots). The radial distance of the dots indicates the zenith angles while the azimuth is measured from the reference direction of the Geographic North Pole. Roman numbers link Pleiades images to their acquisition configuration. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

of the global climate (ERA5) and the Copernicus Atmosphere Monitoring Service (CAMS) reanalysis were used to provide water vapor content and aerosol optical thickness. The effects of clouds, mainly low clouds, and of cast shadows was minimized by filtering out pixels with significant residuals after the inversion of the Hapke model on Pleiades images. These pixels were easily identified and excluded from the estimated parameter maps (Section 5.1).

Accurate image alignment is essential in this work because the BRDF must be extracted on the same pixels. Since the ortho-images provided by CNES could not be fully superimposed, we performed an additional coregistration step. We first applied the Scale-Invariant Feature Transform (SIFT) algorithm (Lowe, 2004) in Orfeo ToolBox (Grizonnet et al., 2017) to identify the tie points between the master image, which has the smallest phase angle and is therefore assumed to be the least distorted, and the other images. We then calculated the average distance between the tie points, and applied a uniform adjustment to align all images to the reference position of the master image. The average distance measured between the master image and the other images is less than or equal to five pixels for 16 images out of 21, which corresponds to a ground distance of 10 m. It should be noted, however, that this approach mainly corrects translational misalignments but no other distortions. In addition, since the area common to all the images is 26.7 km × 20.3 km, over one hundred million pixels would have to be processed at 2 m spatial resolution. The computation time for such a large amount of data would be prohibitive, on the order of several months on a PC. Therefore, we have reduced the resolution of the images by a factor of five (GSD = 10 m), except at experimental sites corresponding to regions of interest of only a few hundred pixels.

The same Pleiades images were used to generate a DEM covering the entire region, with a resolution of 4.5 m (Fig. 3) (Labarre, 2017; Rupnik et al., 2018). This DEM was also co-registered with the images: it is important for the terrain correction process, which will be described in detail in Section 3.2.2. This DEM and all the others mentioned in this article were generated by photogrammetry using MicMac, a free, open-source, multi-view stereo software package developed by the Institut national de l'information géographique et forestière (IGN) and the École nationale des sciences géographique (ENSG) (Rupnik et al., 2017). MicMac was chosen for its great versatility: knowing the camera parameters, such as focal length and pixel size, it is able to generate 3D models from images taken from arbitrary positions (Bretar et al., 2013).

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** DEM of the Asal-Ghoubbet rift generated from Pleiades images at 4.5 m resolution (Labarre, 2017; Rupnik et al., 2018).</figcaption></figure>

## 3.2. Experimental sites

Field data, especially those describing the physical properties of the surface, play a crucial role in validating the results of the inversion. They were collected in the Asal-Ghoubbet rift in February 2016 and partially presented in Labarre et al. (2019). Of the 21 sites studied in this field experiment, only 18 have both ground images and samples in good enough condition to be presented here (Fig. 4). These sites are located in four different areas: i) near Lake Asal (site H), ii) along the western flank of the rift, which shows great variability in terms of geological units, soil types and erosion patterns (sites G, I, K, S, T and U), iii) on the eastern flank north of Ghoubbet al-Kharab (site N), and iv) around the central Fieale volcano, where volcanic deposits and structures predominate (sites A, B, E, F, M, O, P, Q, Y, Z). The vast area of basaltic flow between
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** The eighteen sites of the Asal-Ghoubbet rift, ranked in ascending order of surface roughness. From top left to bottom right: salt bank (site H); clay deposits (site F); sandy fluvial deposits (site K); remodeled hyaloclastic projections (site Z); fluvial deposits (site O); fractured ancient (aphyric) lava (site N); lapilli deposits (site A); altered aphyric lava (site S); powdery limestone rich in mollusk shells (site T); centimetric lava blocks (hyaloclastite) on lacustrine deposits (site E); decimetric lava blocs on lacustrine deposits (site I); small alkaline basalt blocks (site U); remodeled hyaloclastic projections (site Y); decimetric porous lava blocks (site G); fissured lava crust (site M); slab-like basalt porphyry lava (site Q); a'a lava flow (site P); pahoehoe lava flow (site B).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Digital elevation models of site T, reconstructed from a) drone photography; b) photographs taken by hand within a circle; c) photographs taken by hand along a transect. Distances on the x- and y-axes are expressed in meters.</figcaption></figure>
Ardoukóba volcano and Lake Asal is identified as having the darkest surfaces on satellite images and lowest elevations on DEM. However, no field data could be collected in this area due to the difficulty of access.

Each site covers an area between 400 m² and 600 m². Their physical and chemical properties are detailed in Table A1. Surface roughness was determined as the mean slope angle from millimetric DEMs generated by photogrammetry. Reflectance spectra of soil and rock samples were measured in the laboratory. Other surface properties such as compactness, grain size, and shape were assessed by visual inspection of the samples.

### 3.2.1. Ground imagery and high-resolution DEMs

Photographs were used to generate high-resolution DEMs. They were taken by drone and by hand along a transect or within a circle. This allowed DEMs to be obtained at different resolutions, opening up the possibility of studying the multi-scale roughness of the sites. Drone photography is particularly effective for surveying large areas on a scale ranging from a few millimeters to a few centimeters. Hand-captured photos produce sub-millimeter DEMs that allow the study of surface texture and microforms. Transects are obtained by walking along parallel lines in a study site, while circles involve taking photos from multiple positions around a central point. Ground control points (GCPs), easily identifiable on the images, were placed on the ground during acquisition to correct distortions in the DEMs and improve their accuracy. As mentioned above, all DEMs were generated using MicMac.

We supplemented the Labarre et al. (2019) data set, which contains 12 mm-scale DEMs, by generating an additional 6 mm-scale DEMs and 18 sub-millimeter-scale DEMs. Fig. 5 shows the three high-resolution DEMs corresponding to site T.

Represented as a regular grid, these DEMs can be used to calculate the mean slope angle, either by generating randomly oriented facets (Labarre et al., 2017) or by plotting parallel, circular or random height profiles. In this study, we plotted horizontal and vertical linear profiles and determined the mean slope angle for each from statistical distributions. This method requires fewer computational resources than the random facet method, while providing a good estimate of roughness when the number of profiles is sufficient, with the results of the two methods being similar.

### 3.2.2. Soil samples

Soil samples in the form of rock, gravel and powder were collected from the experimental sites (Fig. 6). Their reflectance spectra were measured in the laboratory between 400 nm and 2500 nm using the ASD FieldSpec© 3 spectroradiometer. The light source was an ASD lamp inclined at 15° to the vertical. The pistol grip, equipped with an 8° field-of-view, was positioned 24 cm from the samples to take measurements at nadir. The soils samples were rotated 90° between each series of 10 acquisitions. The 40 spectra obtained were then averaged to reduce the influence of surface heterogeneities, and the final reflectance spectrum resampled to 1 nm spectral sampling interval. These spectra may differ slightly from those obtained in the field, as the surface characteristics of the samples are not perfectly preserved at the time of collection. However, this does not significantly change their representativeness. X-ray diffraction analysis of soil and rock samples was carried out at the Institut de minéralogie, de physique des matériaux et de cosmochimie (IMPMC) to determine their chemical composition. This information is essential for the interpretation of certain parameters of the Hapke model.

## 4. The Hapke model and inversion method

### 4.1. The Hapke model

The photometric model developed by Bruce W. Hapke calculates the reflectance of particulate surfaces such as lunar regolith (Hapke, 1981). It has been modified over the years, with each new version introducing additional terms to improve its accuracy and applicability. The most complete version published by Hapke (2012), which combines all the terms, is based on a number of parameters and functions of varying complexity: single scattering albedo, phase function, porosity, single and multiple scattering components, shadow hiding and coherent backscatter opposition effects. The version to choose and the number of parameters to estimate depend on the radiometric data available. Only measurements at low phase angles (g < 20°), which are not available in our dataset (30° < g < 90°), can constrain the opposition effect and porosity, so the corresponding terms are set to unity. The bidirectional reflectance *r* in Eq. (1) is then calculated from only four parameters: the single scattering albedo (*ω* or SSA), the photometric roughness (Θ)

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Examples of field samples of a) rock (site G), b) lapilli (site A), and c) powder (site T) and corresponding mean reflectance spectrum calculated from 40 spectra acquired using the ASD FieldSpec© 3 spectroradiometer (e, f, g).</figcaption></figure>
occurring in the shadowing function, and the asymmetry (b) and the backscattering (c) parameters of the Heney-Greenstein phase function:

$$r(i, e, g) = \frac{\omega}{4\pi} \frac{\mu_{0e}}{\mu_{0e} + \mu_e} (p(b, c, g) + H(\mu_{0e})H(\mu_e) - 1) S(i, e, g, \bar{\theta}) \quad (1)$$

where *i* and *e* are the solar (incidence) and viewing (emergence) zenith angles, $\mu_{0e}$ and $\mu_e$ are their modified cosines (also called cosines of the effective angle, see Appendix B), and *g* is the phase angle. $H(x)H(y) - 1$ represents the Isotropic Multiple Scattering Approximation (IMSA), where *H* is the Chandrasekhar function. This version of the model is the most commonly used in planetary photometry (e.g., Ceamanos et al., 2013; Fernando et al., 2016; Kugler et al., 2022). The bidirectional reflectance *r* and the bidirectional reflectance factor *BRF* are related by the following equation:

$$\mathrm{BRF}(i, e, g) = \frac{r(i, e, g)\pi}{\cos i} \qquad (2)$$

A detailed description of each parameter is given in Table 1. The global sensitivity analysis performed by Labarre et al. (2017) showed that SSA and the photometric roughness parameter were the two most influential parameters affecting the BRDF. While SSA cannot be measured experimentally, surface roughness can be determined directly or indirectly by various remote sensing techniques from platforms ranging from handheld instruments to satellites to drones. It is important to note that the Hapke model defines roughness as the mean slope angle integrated over all scales (Hapke, 2012). However, research has shown that the photometric response of a surface is primarily influenced by its properties at submillimeter or millimeter scales (Shepard and Helfenstein, 2007; Souchon et al., 2011; Labarre et al., 2017). Phase function parameters provide information about the physical and chemical properties of particulate materials (e.g., McGuire and Hapke, 1995; Souchon

<table><caption>**Table 1**

Physical description of Hapke model parameters.</caption><thead><tr><th>Symbol</th><th>Range</th><th>Parameter</th><th>Description</th><th>Physical significance</th></tr></thead><tbody><tr><td>ω</td><td>[0,1]</td><td>Single scattering albedo</td><td>The ratio of the amount of light at a given wavelength scattered by a representative incremental volume of a medium to the combined amount of light scattered from and absorbed by it</td><td>Particle mineralogical composition, size and microstructure</td></tr><tr><td>θ̃</td><td>[0°,60°]</td><td>Macroscopic roughness</td><td>Mean slope angle, averaged over a size range bounded on the upper end by the angular resolution of the detector and on the lower end by radiative scattering</td><td>Typically roughness at microscale (100–1000 μm)</td></tr><tr><td>b</td><td>[0,1]</td><td>Asymmetry parameter of the Heney-Greenstein phase function</td><td>Anisotropy of the scattering (b &lt; 0.5 broad scattering lobe / b ≥ 0.5 narrow scattering lobe)</td><td>Particle shape, mineralogical composition and internal structure</td></tr><tr><td>c</td><td>[0,1]</td><td>Backscattering parameter of the Heney-Greenstein phase function</td><td>Main scattering direction (c &lt; 0.5 forward scattering / c ≥ 0.5 backward scattering)</td><td>Particle shape, mineralogical composition and internal structure</td></tr></tbody></table>

et al., 2011; Pilorget et al., 2016).

## 4.2. Terrain correction and photometric normalization

The angles of incidence and emergence vary from pixel to pixel in the image depending on the local topography; therefore, terrain correction is a critical step before performing a massive inversion. We determined the slope and orientation for each pixel using the Pleiades DEM (Fig. 3). Given a constant phase angle *g*, we computed a new local illumination and viewing geometric configuration using the law of spherical cosines. Let *i*' be the new local illumination zenith angle, *e'* the new local viewing zenith angle. The geometric correction formulas are:

$$\cos i' = \cos i \times \cos \theta + \sin i \times \sin \theta \times \cos(\phi_i - \varphi), \quad (3)$$

and

$$\cos e' = \cos e \times \cos \theta + \sin e \times \sin \theta \times \cos(\phi_e - \varphi), \quad (4)$$

where *i* and $\phi_i$ are the initial solar zenith and azimuth angles, *e* and $\phi_e$ are the initial viewing zenith and azimuth angles, $\theta$ and $\varphi$ are the slope and orientation of the terrain, and *g* is the phase angle. The local geometric configuration of certain steeply inclined pixels may indicate no direct illumination (shadow) or occultation of the satellite, because the solar or viewing zenith angles may exceed 90°. These represent only a few pixels that have been filtered out.

While the terrain correction assigns a unique geometric configuration to each pixel of the image, the fast Bayesian inversion method described below assumes that all pixels have the same geometric configuration. It is therefore necessary to perform a photometric normalization, which involves calculating the BRDF of each pixel in a user-defined reference configuration ($i''$, $e''$, $g$). To do this, we inverted the Ross-Thick-Li-Sparse (RTLS) model (Wanner et al., 1995) and then ran it with the retrieved parameters to calculate the BRDF in the reference configuration. The RTLS model expresses the bidirectional reflectance of a surface as the sum of three kernels derived from physical considerations that characterize different scattering modes (Eq. 5): i) a Lambertian or isotropic contribution $f_{iso}$, ii) a volume scattering contribution simulated by a set of dispersed facets $f_{vol}K_{vol}$, and iii) a geometric component modeling diffuse scattering and accounting for shading effects $f_{geo}K_{geo}$.

$$r(i, e, g) = r(i', e', g) = f_{iso} + f_{vol}K_{vol}(i', e', g) + f_{geo}K_{geo}(i', e', g) \quad (5)$$

where $f_{iso}$, $f_{vol}$ and $f_{geo}$ are the three parameters to be estimated. This parametric model, originally developed to mimic the bidirectional reflectance of plant canopies, has also been applied to bare soils (Ceamanos et al., 2013). However, some authors have shown that it misrepresents reflectance near the hotspot direction (Lucht et al., 2000; Jiao et al., 2016). The Pleiades BRDFs are not in such a configuration, as the phase angles are generally >30°. Inversion of the RTLS model is performed using the Nelder-Mead simplex method (Nelder and Mead, 1965). Next, we recompute the normalized reflectance factors for each pixel in the user-defined geometric configuration. In our case, it is the same as the one declared in the auxiliary data for the 21 Pleiades images. The new local reflectance factor to be subjected to the massive inversion is then calculated as:

$$r(i'', e'', g) = f_{iso} + f_{vol}K_{vol}(i'', e'', g) + f_{geo}K_{geo}(i'', e'', g) \quad (6)$$

In short, we assume the same geometric configuration ($i,e,g$) for all pixels (Fig. 7A). The terrain correction consists of removing topographic effects by making all pixels perfectly horizontal, without changing their radiometry. This changes the geometric configuration of each pixel (Fig. 7B). Only at the photometric normalization stage, new BRDFs are calculated using the RTLS model, corresponding to a unique geometric configuration for all pixels (Fig. 7C). Finally, we do not invert the Hapke model directly on the Pleiades BRDFs, but on those simulated by the
<figure><img src="image_10.png" /><figcaption>Fig. 7. Simplified graphical representation of terrain correction and photometric normalization for two pixels with different tilt angles ($s_1, s_2$) (and possibly different orientations). Three stages are shown: (A) initial condition (original BRDF, multi-configuration); (B) after terrain correction (original BRDF, multi-configuration); (C) after photometric normalization (simulated BRDF, single configuration). $\vec{n}_1$ and $\vec{n}_2$ are the normal vectors to the facets. $\vec{i}$ represents the illumination vector, $\vec{e}$ the emergence vector, and $\vec{g}$ the phase angle. The marks (') and (°) indicate a change in vector direction at each step. The subscripts (1,2) distinguish the two pixels.</figcaption></figure>

RTLS model after photometric normalization.

### 4.3. Fast Bayesian inversion of the Hapke model

Successful retrieval of Hapke model parameters is closely related to the choice of model inversion technique, due to the nonlinear nature of the model, which can lead to multiple solutions. We used Bayesian inversion, which is based on Bayes' conditional probability theorem and the principle of updated knowledge in the form of posterior likelihood $\rho = p(x, y) \propto p(y|X = x)p(x)$, where $x$ represents the parameters and $y$ the observables. It offers two main advantages over deterministic methods: i) a greater chance of bypassing local minima due to its probabilistic rather than deterministic nature, ii) better control of possible multiple solutions and their uncertainty. The fast Bayesian inversion method developed by Kugler et al. (2022) uses a statistical model known as Gaussian Locally Linear Mapping (GLLiM) (Deleforge et al., 2015). In practice, it approximates the joint distribution $p(X, Y) = (X, Y = F(X) + \epsilon)$ with a parametric statistical model $p\theta$ trained on a learning database ($X_N, F(X_N)$), $n = 1...N$, where $X$ implies uniform sampling over the multidimensional parameter space using the Sobol sequence (Sobol, 1967), $\epsilon$ is Gaussian noise, and $F(X)$ is the expression of the Hapke model. The approach closely related to Gaussian mixtures assumes that the direct and inverse relations between $X$ and $Y$ can be locally approximated by a combination of $K$ local affine transformations. Then, for each vector $Y_{obs}$, which in our case are the bidirectional reflectance factors evaluated for a pixel in the reference geometric configuration, the GLLiM method calculates the approximate posterior $p\theta(X|Y)$ and uses it to predict $X$. Since $p\theta(X|Y)$ is based on a unified set of global parameters, the computational cost is significantly reduced when analyzing a large number of $Y_{obs}$ vectors.

The first step in the inversion is to create a training data set. This is done by generating parameter sets of the Hapke model using Sobol sequences, then running the model to compute surface reflectance and adding centered Gaussian noise. Next, building the GLLiM parametric model $p\theta$ involves two crucial steps: an initialization step, which proposes initial values for the global parameters, and a learning step, which uses the expectation-maximization algorithm (Dempster et al., 1977). The latter iteratively adjusts the parameters to maximize the likelihood of the probabilistic model, thus refining the global parameters. Once the learning phase is complete, the GLLiM algorithm produces a posterior prediction distribution $p\theta(X|Y_{obs})$, which forms the basis for selecting the optimal solution. Standard predictions can be made by averaging the posterior distribution, an efficient method when it is unimodal and the solution is unique, or by identifying the modes of the distribution by merging the components of the Gaussian mixture model into a small number of centroids (typically two or three) if it is multimodal. In addition, Incremental Mixture Importance Sampling (IMIS), known for its simplicity and computational efficiency, can be used to refine the mean and mode predicted solutions. IMIS iteratively samples the true

<figure><img src="image_11.png" /><figcaption>Fig. 8. Fast Bayesian inversion framework.</figcaption></figure>

target probability distribution function $p(X|Y_{obs})$ in a series of updated small sample groups with the proposed distribution set to $p_{GLLiM}(X|Y_{obs})$ around the mean or around the centroids/modes, combines the results, and proposes a new mean and modes. Kugler et al. (2022) have shown that IMIS generally gives better results. In practice, we use different prediction methods and choose the one with the lowest residuals. The entire fast Bayesian framework is illustrated in Fig. 8.

Fast Bayesian inversion is highly efficient for large-scale inversions, while retaining the benefits of Bayesian analysis. It improves prediction accuracy, reduces computation time and provides information about parameter estimates thanks to the posterior prediction distribution. It is also very useful in poorly constrained inverse problems where multiple solutions may coexist (e.g., Fernando et al., 2013; Schmidt and Fernando, 2015; Kugler et al., 2022).

### 4.4. Evaluation of the accuracy of estimates

As mentioned above, estimating a global solution by inversion is challenging due to the nonlinearity of the model. The lack of prior knowledge about parameter values often turns the inverse problem into an attempt to identify an optimal combination of variables that reduces the gap between observations and model predictions through iterative numerical analysis. Therefore, we calculated the root mean square error (RMSE) of the reconstructed BRF as follows:

$$RMSE = \frac{1}{N} \sqrt{\sum_{k=1}^{n} (r_k - r(i_k, e_k, g_k))^2} \quad (7)$$

where $n$ is the number of data points (number of images), $r_k$ is the measured reflectance, and $r$ is the reflectance predicted by the model. Therefore, we select the solution with the lowest RMSE. Fast Bayesian inversion also provides the standard deviation (SD) of each parameter
for each pixel, which is useful for understanding the constraints on each parameter during inversion.

# 5. Results and discussion

In this section, we first assess the quality of the inversion results in terms of model fit to the data, before interpreting the parameter maps and associated error maps.

## 5.1. Overall quality of the inversion

To assess the accuracy of the model's predictions, we consider RMSE values of 0.02 or less to be acceptable (Fernando et al., 2013). This excludes clouds, their shadows and steep slopes where the fit is poor. As shown in Fig. 9, the model reproduces the amplitude and general shape of the BRDF. A comparison between the RMSE map (Fig. 10a) and the slope map generated from the Pleiades DEM (Fig. 3) clearly shows that the RMSE is higher in steep rift areas. RMSE values are generally below 0.01 for most sites (Fig. 10b), indicating good agreement between predicted and observed values. For sites K and S, they vary between 0.1 and 0.12 due to cloud contamination of some images, which produced 2–3 outliers. However, the effect on the inversion is moderate. The fit is also worse at bright sites F, H and T, with RMSE values between 0.15 and 0.2. This suggests that an increase in the BRF corresponds to an increase in the residuals, as shown in Figs. 9 and 10.

Fernando et al. (2016) and Schmidt and Fernando (2015) have demonstrated the importance of phase angles in BRDF analysis. Specifically, they pointed out that phase angles >90° and ranges of variation >50° are essential to effectively constrain the Hapke model parameters. With phase angles ranging from 30° to 90°, our dataset partially meets these requirements, with 80 % of the region meeting the stringent requirement of an RMSE ≤0.02 (Fig. 10a).

## 5.2. Single scattering albedo (ω)

### 5.2.1. Results

Fig. 11a shows the single scattering albedo (SSA) map in the green band obtained by massive Bayesian inversion of the Hapke model on the Asal-Ghoubbet rift. The SSA is well constrained with $SD < 0.015$ almost everywhere in the image (Fig. 11b). The lowest values are observed on basaltic units, but they vary significantly between the inner floor of the rift and its two flanks. The large porphyritic basalt region on the southern shore of Lake Asal, which is predominantly composed of plagioclases, has an average SSA value of ~0.38. The Ardoukôba caldera, characterized by phyric basalt (Fig. 12), has systematically lower SSA values of ~0.26. Similar values are found in the center of the rift, where there is another distinct zone of phyric basalt. The basalt of the Fieale volcano, whose caldera is identified as a porphyric unit on the geologic map, is classified as phyric in our study, with SSA values decreasing progressively from 0.34 to 0.26 between the caldera and Lake Ghoubbet. The SSA values of both rift flanks composed of porphyritic and phytic basalt range from 0.45 to 0.55. This is much higher than in the inner rift.

In contrast to these low values associated with basaltic materials, several bright patches of lacustrine limestone, clearly visible on satellite imagery, have very high SSA values of ~0.9. Nearby basaltic pyroclastic units also stand out with SSA values of ~0.65. They are often associated with lacustrine formations and contribute to the complex geological mosaic of this region. Finally, to the southwest of the rift are several zones characterized by fractured alkaline basaltic flows, particularly hawaiite, with SSA values of ~0.7.

Considering the study sites, Fig. 11c shows consistent results with the SSA map in Fig. 11a: SSA is <0.4 for pure basalt ($P < Q < B < N$), > 0.7 for lacustrine limestone, salt, sand field ($T < H < F$), and in between for mixed sand and lava rock or reworked hyaloclastite ($U < G < Z < M < I < Y < A < S < E < O < K$).

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison between the BRF measured by Pleiades at sites G (RMSE = 0.003), A (RMSE = 0.005) and T (RMSE = 0.014) (circles) and simulated by the Hapke model (triangles) in the four Pleiades spectral bands.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Root mean square error (RMSE) between measured and modeled BRDF a) over the entire image in the Pleiades spectral band centered on 555 nm; b) for the 18 experimental sites in the four Pleiades spectral bands, in ascending order of value.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Estimation of the single scattering albedo a-b) over the Asal-Ghoubbet rift in the spectral band of Pleiades centered on 555 nm and corresponding SD map; c) Over the 18 experimental sites in the four Pleiades spectral bands, in ascending order of value.</figcaption></figure>

### 5.2.2. Discussion

Several surface properties have been reported to have a significant effect on SSA. These include soil chemical composition, grain size and shape, and particle microstructure (e.g., crystals, fractures, pores) (Hapke, 2012). Shepard and Helfenstein (2007) noted that it may be

related to soil compaction. SSA also decreases with increasing absorption coefficient and/or particle size (Fernando et al., 2016; Gao et al., 2021) and it tends to be higher for fine-grained soils than for coarse-grained soils (Johnson et al., 2013). This interpretation must be made with caution, as few surfaces are granular. In addition, it would be
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Relationship between the single scattering albedo (ω) obtained by inversion of the Hapke model on the Pleiades BRDF and the bidirectional reflectance factor (BRF) of soil and rock samples from a) site A, b) site B, c) site F and d) site H measured in the laboratory in the four Pleiades spectral bands.</figcaption></figure>

difficult to assess grain size from SSA without information on the refractive index of these soils. Because basaltic surfaces are generally not granular, low SSA values are mainly due to their chemical composition, in particular the presence of magnesium- and iron-rich minerals, which are ubiquitous in lava flows. The effect of iron concentration on SSA has been extensively documented by Lucey and Riner (2011) and Wang et al. (2017). The porosity of mafic minerals in basalt flows may also play a role. As noted by Hapke (2012), porosity tends to decrease SSA. However, it is difficult to measure in the field. The SSA of sand and salt fields is much higher. For example, the clay deposits at site F, composed mainly of calcite and quartz, and site T, composed of powdery limestone with calcite as the main component, both have high SSA values. The same is true for site H, which is characterized by a salt bank containing halite and gypsum. The common feature of these sites is the presence of crystallized and glassy materials, which generally increase the SSA. In addition, the texture of these materials also influences SSA. Site T is very fine-grained and compact, while site H is more coarse-grained and site F has a mixture of granular and blocky materials. This nuanced analysis of chemical composition and texture allows us to understand how different geological factors affect the photometric properties of these surfaces.

Fig. 11c shows that sites with intermediate SSA deserve special attention. They present a diverse range of mixtures (Fig. 4), but this diversity does not seem to be directly correlated with significant variations in SSA. A more relevant factor seems to be the fraction of bright materials present in the mixtures. Pilorget et al. (2015) have emphasized that the BRDF pattern is mainly influenced by the properties of the most abundant, brightest or most anisotropic grains in the mixture. We note that sites E, K and O are predominantly fluvial or lacustrine deposits, characterized by a significant presence of sand and bright sediments. These geological properties generally result in slightly higher SSA values, reflecting the brighter nature of these materials. Conversely, sites A, G, I, M, S, Y, and Z are predominantly composed of lava rocks, which by their nature have darker materials and therefore lower SSA values.

By extrapolating the results obtained from the 18 study sites, we can interpret the photometric properties of other regions of the rift. For example, we find that the SSA of the rift flanks is influenced by either degradation or mixing with bright materials. This corroborates the geologic map, which indicates that the rift flanks, which date from 87,000 to 220,000 years ago, are older than the inner zone, which dates from 66,000 years ago to the present day (Appendix C). Weathering and erosion processes have likely contributed to the mixing and alteration of the original materials, potentially increasing the SSA through the incorporation of brighter sedimentary components. Similarly, the behavior of pyroclastic units can be compared to that of compact sand or lacustrine deposits, assuming that their composition consists mainly of fine volcanic ash of medium brightness.

### 5.2.3. Comparison between SSA and measured reflectance spectra

A linear relationship is observed between measured reflectance values and estimated SSA values at different wavelengths across sites, but it is not unique as the slopes and intercepts vary from site to site (Fig. 12). Figs. 12c and d show that this relationship becomes non-linear as SSA reaches its upper limit of 1, which could reflect an SSA saturation effect. This suggests that this relationship could be influenced by factors, such as mineral composition, grain size or other environmental conditions and the model limitations. The lack of such information highlights the need for more laboratory or field measurements to interpret remote sensing products.

## 5.3. Photometric roughness ($\bar{\theta}$)

### 5.3.1. Results

We refer to photometric roughness as either $\bar{\theta}$ or simply “roughness”. Most terrains have roughness values between 35° and 40° with a standard deviation of 2 to 3°. Areas with low roughness values are the basalt units and the salt bank west of Lake Asal (Fig. 13). The large porphyritic basalt unit has a wide range of roughness values between 10° and 45°.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Estimation of photometric roughness a-b) over the Asal-Ghoubbet rift in the Pleiades band centered on 555 nm and corresponding SD map; c) Over the 18 experimental sites in the four Pleiades spectral bands, in ascending order of value.</figcaption></figure>

Surface roughness is ~20° on the shores of Lake Asal, ~25° in the region of the Ardoukôba volcano, and between 25° and 30° on the phyric basalt flow in the center of the rift. The basaltic flows of the Fieale volcano have a surface roughness that progressively decreases from 30° to 20° from the highest to the lowest altitudes. On the salt bank of Lake Asal, the roughness decreases from 20° to 5° as one moves away from the shore. These are the lowest values recorded in the Asal-Ghoubbet rift. The lacustrine and pyroclastic basalt units have a more homogeneous and higher roughness (35°-40°), as do the rift flanks (~37°). Note that surface roughness cannot be accurately estimated on steep slopes.

Fig. 13c shows that sites P, H and F have significantly lower photometric roughness (10–20°) than the others (30–40°). Although the roughness should not vary too much within a site, sites P, B, U and N show greater variability than the others. These are the sites where the SSA was low. We also observe that the roughness decreases with wavelength for most of the sites, with the exception of sites P, H and F, previously mentioned as having the lowest roughness.

### 5.3.2. Discussion

As mentioned above, photometric roughness encompasses spatial scales ranging from several times the average interparticle distance (100 to 1000 µm) to a maximum determined by the instantaneous field of view (IFOV) of the detector. In Earth observation, the IFOV defines the ground sampling distance, which can vary from a few decimeters to a

few kilometers. The interpretation of $\bar{\theta}$ is a matter of debate within the scientific community, due to the ambiguous definition of these scales (Shiltz and Bachmann, 2023). However, there is a general consensus that photometric roughness is most sensitive to the millimeter scale and below (Shepard and Helfenstein, 2007; Pilorget et al., 2016; Labarre et al., 2019). We applied the isotropic undecimated wavelet transform (IUWT) (Starck and Murtagh, 2006) to decompose the high-resolution DEMs, so as to disaggregate the roughness components at different scales, and then calculated their contribution. Our results confirm that smaller scales generally have higher roughness (Fig. 14). This is because smaller-scale topographic features contribute more significantly to the overall surface roughness due to their greater frequency and variation in elevation within a given area.

Jacquemoud et al. (1992) and Protopapa et al. (2020) hypothesized that, unlike other photometric parameters, $\bar{\theta}$ does not vary with wavelength. Conversely, Cord et al. (2003) found that when surface reflectance varies with wavelength, photometric roughness is also affected by wavelength. The optical properties of minerals or grains, such as absorption and scattering, can vary as a function of wavelength. Meanwhile, as mentioned earlier, the photometric behavior of natural surfaces is strongly influenced by roughness on the microscopic scale, which is related to absorptivity. When these microscopic structures are opaque, there is a marked photometric contrast between shadowed and directly illuminated areas, resulting in a situation where photometric
<figure><img src="image_17.png" /><figcaption>Fig. 14. Mean slope angle (°) (y-axis) across various scales (meter) (x-axis), calculated from different levels of IUWT decomposition of drone DEMs.</figcaption></figure>

roughness is closely related to geometric roughness. Conversely, when the grains are partially translucent, this contrast is greatly reduced, resulting in a reduced photometric roughness that is less correlated with geometric roughness. This relationship suggests that there may be a general anti-correlation between photometric roughness and SSA for different wavelengths, at least in the solar domain and at most of the sites studied.

The relationship between SSA and photometric roughness in photometric models is key to understanding surface reflectance properties. These two parameters are independently influenced by the physical properties of a surface, but are also interdependent in the way

they affect the BRDF. Shepard and Helfenstein (2007) have shown, through experiments on various soil samples, that SSA and photometric roughness are intrinsically linked, providing empirical evidence of their combined effect on photometric observations. A sensitivity analysis of the Hapke model conducted by Labarre et al. (2017) also showed this coupling, which can explain up to 30 % of the reflectance variability. Figs. 11 and 13 show that SSA values below 0.8 are required to estimate photometric roughness. Sites with these values can still be considered “opaque”, and their photometric roughness corresponds more to the geometric roughness than to the surface material.

<figure><img src="image_18.png" /><figcaption>Fig. 15. Relationship between a) photometric roughness ($\bar{\theta}$) and mean slope angle ($\bar{\theta}_M$). b) ($\bar{\theta} - \bar{\theta}_M$) and single scattering albedo ($\omega$). Only outliers are labeled.</figcaption></figure>
### 5.3.3. Comparison between photometric roughness and mean slope angle

Fig. 15a shows the relationship between photometric roughness estimated by inversion of the Hapke model ($\bar{\theta}$) and the mean slope angle calculated from DEMs at millimeter scales ($\theta_M$). Our results corroborate those of Labarre et al. (2019) and also show a weak correlation. Photometric roughness is larger than the calculated mean slope angle at most sites, with the exception of sites H and P. We have calculated the difference $\bar{\theta} - \bar{\theta}_M$ and plotted it as a function of SSA (Fig. 15b).

According to Hapke (2012), SSA depends on material refractive index but also on grain size, which represents the smallest scale of roughness. The linear relationship observed here could therefore be explained by this submillimeter scale, which has not been measured in the field, however, it needs to be validated by further experiments. On the other hand, sites F, H and T sites with high SSA values show significant deviations. This raises several questions, for which we have no complete answer: What determines such a linear relationship? Why do sites F, H and T deviate from the correlation line? Do sites with very high SSA

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Estimation of the asymmetry parameter *b* and backscattering parameter *c* over the Asal-Ghoubbet rift in the Pleiade band centered on 555 nm (a,c) and corresponding SD maps (b,d); Over the 18 experimental sites in the four Pleiades spectral bands, in ascending order of value (e,f).</figcaption></figure>
follow a different trend? Site M also shows an outlier, with a mean slope
angle estimated to be too low, which could be a local characteristic not
reflecting the geological and topographic conditions of the region.

5.4. Phase function parameters (b and c)

5.4.1. Results

The double Henyey-Greenstein phase function over the Asal-Ghoubbet rift is shown in Fig. 16. The latter generally exhibits broad scattering lobes in both the backward and forward directions. Although both maps are noisy, the maximum SD is only 0.016, which is smaller than for the SSA. Parameters *b* and *c* are mostly in the range of 0.2 to 0.4 and 0.4 to 0.6, respectively. This scattering behavior is consistent across most geological units, making it challenging to identify distinct features. The most notable exception is the salt bank, which features a very high asymmetry parameter (*b* > 0.9), indicating narrow, directional scattering lobes. This geological unit also exhibits values of *c* < 0.1, indicating a predominance of forward scattering. The sandy zone in the middle of the rift also stands out, with values of *b* that can be as low as 0.2, but also close to 1 on some pixels. We also observed a high value for the parameter *c*, which can reach a maximum of 0.7, but it can average 0.5. In pyroclastic areas, *c* values can be very low, reaching 0.15, but averaging 0.3.

The results for the asymmetry parameter at the local scale corroborate the observations at the rift-scale. Most sites have values of *b* < 0.4, with the exception of sites H and F. The latter are very bright and therefore have a high SSA. Sites H and F have the lowest values of *c* and the highest values of *b*, highlighting their distinct photometric responses. Fig. 16c also shows that *c* varies widely between the different sites within its range of variation, with no other significant features emerging. There is no clear variation in these two parameters with wavelength. Volcanic terrains exhibit broad scattering lobes and a balance between forward and backward scattering. In contrast, sites F and H characterized by high albedo and glassy grains, have narrow forward scattering lobes.

5.4.2. Discussion

McGuire and Hapke (1995) studied particle scattering as a function of physical properties: smooth, bright and spherical particles tend to have narrower scattering lobes (high values of *b*) and greater forward scattering (low values of *c*). Conversely, irregular and rough particles or those with internal scatterers tend to have wider scattering lobes (low values of *b*) and greater backward scattering (high values of *c*). Plotting *b* versus *c* produces an “L-shape” or “hockey-stick” curve (Fig. 16) that represents the transition from rough, opaque and irregular grains to smooth, clear and spherical grains. This plot highlights the influence of particle structure on light scattering properties. Hartman and Domingue (1998) found a behavior independent of particle size, suggesting that the observed effects are due to the overall structure of the particle rather than its size. Cord et al. (2003) and Shepard and Helfenstein (2007) have suggested that the scattering patterns are more sensitive to the organization of the grains on the surface, such as density or micro-roughness, than to the intrinsic optical properties of the particles. The work of Lucas et al. (2023) on the effect of soil compactness on reflectance has confirmed this observation. To better understand how light is scattered by a surface, Souchon et al. (2011) performed BRF measurements on natural granular surfaces composed of volcanic materials that may vary in grain size (microns to millimeters), shape, surface texture, and mineralogy. They showed that these surfaces scatter light mainly forward, even though they contain only a small fraction of isolated translucent monocrystals or fresh glass. This revealed new modes of light scattering and allowed the *b* and *c* parameter space to be extended. The

<figure><img src="image_20.png" /><figcaption>Fig. 17. Relationship between phase function parameters estimated by inversion of the Hapke model at the 18 sites studied. Texts in italic correspond to groups of sites and designate their main constituents, while other texts are the physical properties at the relative position on the graph (after Souchon et al., 2013; Fernando et al., 2016).</figcaption></figure>

Asal-Ghoubbet rift contains granular surfaces (reworked volcanic ma-
terials, sand, and salt) as well as rocky surfaces, lava flows and mixtures
of all these terrain types.

Consistent with the experimental results of Souchon et al. (2011), most sites show irregular, rough, opaque, and hollow grains on granular surfaces (0.2 < *b* < 0.4 and 0.4 < *c* < 0.6), which is typical of pyroclastic material. Fig. 17 shows three groups based on the mineral composition and location of the sites. Group 1 comprises lava flows and hyaloclastite. Sites in this group are characterized by hard, blocky lava flows with a wide range of backscatter parameters (from 0.2 to 0.9) and low asymmetry parameters (*b* < 0.3). This suggests materials with moderate to high internal density, consistent with surfaces composed of rough, opaque particles and hollow grains. Group 2 includes calcareous and lacustrine deposits, consisting mainly of sand and gravel. The position of this group on the graph indicates a structure consisting of agglomerated fine particles, typical of sites I and T. However, some sites show a loose combination of rocks and gravels (e.g., sites A, K, M). There is significant overlap between groups 1 and 2, corresponding to sites containing a mixture of the two materials, characterized by low to intermediate values of *b* and *c*. Group 3 is composed of bright salt and clay, notably occupying a unique region on the graph with very high *b* values ( > 0.8) and low *c* values ( < 0.3). This group, which includes sites F and H, has surface properties similar to fresh opaque glass or isolated single crystals (Souchon et al., 2013). The mineral composition of these sites (F with calcite, quartz, albite, zeolite, palygorskite, smectite and talc, and H with halite and gypsum) supports their inclusion in this group. The diagram in Fig. 17 is useful for a general understanding of the physical properties of surfaces. For a higher level of accuracy, interpretation becomes complicated when the sites represent a mixture of materials where the percentage contribution of each material is unknown.

6. Conclusion

We have performed a comprehensive inversion of the Hapke model
over the Asal-Ghoubbet rift (Republic of Djibouti) using high-resolution,
multi-angular Pleiades imagery. To our knowledge, this is the first large-scale inversion of its kind on Earth, with most previous efforts limited to localized areas of interest. Our extensive coverage provides valuable insights into the landscape, including the spatial arrangement of geological units. This broad-scale inversion technique facilitates the detection of new surface features that might be missed by more localized analyses. Beyond simple data processing, this approach offers exciting potential for developing innovative methods for analyzing multi-angular photometric data. The capability is integrated into a framework that applies the Hapke model.

Our study also advances the interpretation of Hapke model parameters by linking them to field measurements. While single scattering albedo (SSA) remains the most reliable parameter to extract, as demonstrated in other studies, its validation presents a significant challenge. There is no simple formula relating SSA to surface reflectance from a single measurement configuration, indicating the need for further methodological development or alternative approaches to bridge the gap between theoretical models and observables. In addition, our results confirm the close relationship between photometric roughness and SSA, emphasizing that accurate extraction of photometric roughness may depend on SSA values. As for the phase angle parameters, their correlation with grain properties varies depending on the surface type and the physical and optical properties of the materials, but the general relationship holds for most sites.

This work could be improved in several ways. First, data acquisition should prioritize measurements from the principal plane, where incidence, emergence, and normal vectors align, to maximize the retrieval of photometric information (Schmidt and Bourguignon, 2019). Unfortunately, this is challenging for satellite data. Selecting a plane that is approximately 30° offset from the principal plane, or using other geometries that better sample the diffusion lobes of the phase function (Kugler, 2021), may be beneficial. Second, inverting the Hapke model on BRDFs with a GSD of 2 m would allow maximum detail to be preserved, although computational time constraints must be considered. Finally, a more in-depth analytical approach and field measurements are needed to understand the correlation between photometric roughness and SSA derived from the inversion. This will involve not only validating model results with field data, but also developing new methods to better

integrate these data sets. These steps will enable us to refine our methods, leading to more accurate predictions and a better understanding of surface properties.

### CRediT authorship contribution statement

D.T. Nguyen: Methodology, Formal analysis, Writing – original draft. S. Jacquemoud: Supervision, Methodology, Conceptualization, Writing – review & editing. A. Lucas: Supervision, Methodology, Writing – review & editing. S. Douté: Software, Methodology, Writing – review & editing. C. Ferrari: Methodology, Writing – review & editing. S. Coustance: Resources, Methodology, Writing – review & editing. S. Marcq: Resources, Methodology, Writing – review & editing. A. Meygret: Resources, Methodology, Writing – review & editing.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Stephane JACQUEMOUD reports financial support was provided by Université Paris Cité. Stephane JACQUEMOUD reports a relationship with Université Paris Cité that includes: employment. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

Nguyen Tri Dung's PhD thesis was funded by the STEP'UP Doctoral School of Université Paris Cité. This work was supported by the CNES in the framework of the CAROLInA (Characterization of Multi-Scale Roughness using OpticaL ImAgery) and SURFACEs (Etude de la rugosité de surface par radiométrie multi-angulaire) projects funded by the CNES and by the Labex UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001) in the frame of the SOILSPECT project. We would like to thank Erica Ucellatori (University of Pavia) and Rodolphe Marion (CEA) for their help in measuring the reflectance spectra of the soil samples.

## Appendix A

<caption>Table A1</caption>

Geographical coordinates, geological description and mineral composition of 18 sites.

<table><thead><tr><th>Site</th><th>Latitude N</th><th>Longitude E</th><th>Description</th><th>Mineral composition</th></tr></thead><tbody><tr><td>A</td><td>11.581021</td><td>42.493072</td><td>Lapilli deposits</td><td>Anorthite, Diopside, Hematite, Smectite</td></tr><tr><td>B</td><td>11.594598</td><td>42.472655</td><td>Ropy pahoehoe lava</td><td></td></tr><tr><td>E</td><td>11.593649</td><td>42.491789</td><td>Centimetric lava blocks on lacustrine deposits</td><td>Anorthite, Quartz, Calcite, Augite, Albite, Amphibole, Palygorskite, Talc, Smectite</td></tr><tr><td>F</td><td>11.592512</td><td>42.495872</td><td>Clay deposits</td><td>Calcite, Quartz, Albite, Zeolithe, Palygorskite, Smectite, Talc</td></tr><tr><td>G</td><td>11.562228</td><td>42.440316</td><td>Decimetric porous lava blocks</td><td></td></tr><tr><td>H</td><td>11.619922</td><td>42.396506</td><td>Salt bank</td><td>Halite, Gypsum</td></tr><tr><td>I</td><td>11.563368</td><td>42.448571</td><td>Decimetric lava blocs on lacustrine deposits</td><td></td></tr><tr><td>K</td><td>11.538943</td><td>42.489818</td><td>Sand, fluvial deposit</td><td>Diopside, Anorthite, Smectite</td></tr><tr><td>M</td><td>11.589644</td><td>42.503011</td><td>Fissured lava crust</td><td>Diopside, Anorthite, Chlorite, Smectite, Palygorskite</td></tr><tr><td>N</td><td>11.604637</td><td>42.522814</td><td>Fractured ancient (aphyric) lava</td><td>Diopside, Anorthite, Amphibole, Chlorite, Smectite</td></tr><tr><td>O</td><td>11.584279</td><td>42.513356</td><td>Fluvial deposit</td><td>Anorthite, Diopside, Smectite</td></tr><tr><td>P</td><td>11.575606</td><td>42.511158</td><td>A'a lava</td><td></td></tr><tr><td>Q</td><td>11.561884</td><td>42.513632</td><td>Slab-like basalt porphyry lava</td><td>Diopside, Smectite</td></tr><tr><td>S</td><td>11.541699</td><td>42.481100</td><td>Altered aphyric lava</td><td></td></tr><tr><td>T</td><td>11.558583</td><td>42.455869</td><td>Pulverulent limestone, rich in mollusk shells</td><td>Calcite, Muscovite, Mg-rich Calcite, Quartz, Palygorskite, Smectite</td></tr><tr><td>U</td><td>11.545372</td><td>42.453302</td><td>Small alkaline basalt blocks</td><td></td></tr><tr><td>Y</td><td>11.586405</td><td>42.488147</td><td>Remodeled hyaloclastic deposits</td><td>Anorthite, Diopside, Smectite, Mg-rich Calcite, Halite, Palygorskite</td></tr><tr><td>Z</td><td>11.583280</td><td>42.489402</td><td>Reshaped hyaloclastic deposits</td><td></td></tr></tbody></table>
## Appendix B. Extra components of the Hapke model

* Two-parameter Henyey-Greenstein phase function

$$p(b, c, g) = (1-c) \frac{1-b^2}{(1+2b \cos g + b^2)^{3/2}} + c \frac{1-b^2}{(1-2b \cos g + b^2)^{3/2}} \quad (B.1)$$

We used a slightly different backscattering parameter $c = (1 + c')/2$, where $c$ is the new backscattering parameter and $c'$ is the old one from Hapke (2012). This was done to reset the parameter range (Table 1) from [-1,1] to [0,1].

Chandrasekhar $H$ function

$$r_0 = \frac{1 - \sqrt{1 - \omega}}{1 + \sqrt{1 - \omega}} \quad (B.2)$$

$$H(x) = \left( 1 - \omega x \left( r_0 + \frac{1 - 2r_0 x}{2} \ln \frac{1+x}{x} \right) \right)^{-1} \quad (B.3)$$

* Modified cosine

$$\mu_{0e} = \begin{cases} \chi(\bar{\theta}) \left( \cos i + \sin i \tan \bar{\theta} \frac{\cos \phi E_2(\bar{\theta}, e) + \sin^2(\phi/2) E_2(\bar{\theta}, i)}{2 - E_1(\bar{\theta}, e) - (\phi/\pi) E_1(\bar{\theta}, i)} \right) & \text{when } i < e \\ \chi(\bar{\theta}) \left( \cos i + \sin i \tan \bar{\theta} \frac{E_2(\bar{\theta}, i) - \sin^2(\phi/2) E_2(\bar{\theta}, e)}{2 - E_1(\bar{\theta}, i) - (\phi/\pi) E_1(\bar{\theta}, e)} \right) & \text{when } i \ge e \end{cases} \quad (B.4)$$

with $\phi = [\phi_1 - \phi_e]$

$$\mu_e = \begin{cases} \chi(\bar{\theta}) \left( \cos e + \sin i \tan \bar{\theta} \frac{E_2(\bar{\theta}, e) - \sin^2(\phi/2) E_2(\bar{\theta}, i)}{2 - E_1(\bar{\theta}, e) - (\phi/\pi) E_1(\bar{\theta}, i)} \right) & \text{when } i < e \\ \chi(\bar{\theta}) \left( \cos e + \sin i \tan \bar{\theta} \frac{\cos \phi E_2(\bar{\theta}, i) + \sin^2(\phi/2) E_2(\bar{\theta}, e)}{2 - E_1(\bar{\theta}, i) - (\phi/\pi) E_1(\bar{\theta}, e)} \right) & \text{when } i \ge e \end{cases} \quad (B.5)$$

$$\chi(\bar{\theta}) = \frac{1}{\sqrt{1 + \pi \tan^2 \bar{\theta}}} \quad (B.6)$$

$$E_1(\bar{\theta}, x) = \exp\left(-\frac{2}{\pi}\cot\bar{\theta}\cot x\right) \quad (B.7)$$

$$E_2(\bar{\theta}, x) = \exp\left(-\frac{1}{\pi}\cot^2\bar{\theta}\cot^2 x\right) \quad (B.8)$$

* Shadowing function

$$\eta(x) = \chi(\bar{\theta}) \left( \cos x + \sin x \tan \bar{\theta} \frac{E_2(\bar{\theta}, x)}{2 - E_1(\bar{\theta}, x)} \right) \quad (B.9)$$

$$S(i, e, \phi, \bar{\theta}) = \begin{cases} \frac{\mu_e}{\eta(e) \eta(i)} \frac{\chi(\bar{\theta})}{1 - f(\phi) + f(\phi)\chi(\bar{\theta}) (\mu_0 / \eta(i))} & \text{when } i < e \\ \frac{\mu_e}{\eta(e) \eta(i)} \frac{\chi(\bar{\theta})}{1 - f(\phi) + f(\phi)\chi(\bar{\theta}) (\mu / \eta(e))} & \text{when } i \ge e \end{cases} \quad (B.10)$$

with $\mu_0 = \cos i$

$$f(\phi) = \exp\left(-2\tan\frac{\phi}{2}\right) \quad (B.11)$$
<figure><img src="image_21.png" /><figcaption>Appendix C. Geological map of the Asal-Ghoubbet rift (after Stieltjes, 1980; Pinzuti et al., 2013)</figcaption></figure>

## Data availability

The authors do not have permission to share data.

## References

Bachmann, C.M., Eon, R.S., Ambeau, B., Harms, J., Badura, G., Griffo, C., 2017. Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones dunes. J. Appl. Remote. Sens. 12 (1), 012005. https://doi.org/10.1117/1.JRS.12.012005.

Bretar, F., Arab-Sedze, M., Champion, J., Pierrot-Deseilligny, M., Heggy, E., Jacquemoud, S., 2013. An advanced photogrammetric method to measure surface roughness: application to volcanic terrains in the piton de la Fournaise, Reunion Island. Remote Sens. Environ. 135, 1–11. https://doi.org/10.1016/j.rse.2013.03.026.

Cao, H., You, D., Ji, D., Gu, X., Wen, J., Wu, J., Li, Y., Cao, Y., Cui, T., Zhang, H., 2023. The method of multi-angle remote sensing observation based on unmanned aerial vehicles and the validation of BRDF. Remote Sens. 15 (20), 5000. https://doi.org/10.3390/rs15205000.

Ceamanos, X., Douté, S., Fernando, J., Schmidt, F., Pinet, P., Lyapustin, A., 2013. Surface reflectance of Mars observed by CRISM/MRO: 1. Multi-angle approach for retrieval of surface reflectance from CRISM observations (MARS-ReCO). J. Geophys. Res.: Planets 118 (3), 514–533. https://doi.org/10.1029/2012JE004195.

Cord, A., Pinet, P., Daydou, Y., Chevrel, S., 2003. Planetary regolith surface analogs: optimized determination of Hapke parameters using multi-angular spectro-imaging laboratory data. Icarus 165, 414–427. https://doi.org/10.1016/S0019-1035(03)00204-5.

Deleforge, A., Forbes, F., Horaud, R., 2015. High-dimensional regression with Gaussian mixtures and partially-latent response variables. Stat. Comput. 25 (5), 893–911. https://doi.org/10.1007/s11222-014-9461-5.

Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. B. Methodol. 39 (1), 1–22. https://doi.org/10.1111/j.2517-6161.1977.tb01600.x.

Ding, A., Ma, H., Liang, S., He, T., 2022. Extension of the Hapke model to the spectral domain to characterize soil physical properties. Remote Sens. Environ. 269, 112843. https://doi.org/10.1016/j.rse.2021.112843.

Eon, R.S., Bachmann, C.M., Gerace, A.D., 2018. Retrieval of sediment fill factor by inversion of a modified Hapke model applied to sampled HCRF from airborne and satellite imagery. Remote Sens. 10 (11), 1758. https://doi.org/10.3390/rs10111758.

Fernando, J., Schmidt, F., Ceamanos, X., Pinet, P., Douté, S., Daydou, Y., 2013. Surface reflectance of Mars observed by CRISM/MRO: 2. Estimation of surface photometric properties in Gusev Crater and Meridiani Planum. J. Geophys. Res.: Planets 118 (3), 534–559. https://doi.org/10.1029/2012JE004194.

Fernando, J., Schmidt, F., Douté, S., 2016. Martian surface microtexture from orbital CRISM multi-angular observations: a new perspective for the characterization of the geological processes. Planet Space Sci. 128, 30–51. https://doi.org/10.1016/j.pss.2016.05.005.

Gao, C., Xu, M., Xu, H., Zhou, W., 2021. Retrieving photometric properties and soil moisture content of tidal flats using bidirectional spectral reflectance. Remote Sens. 13 (7), 1402. https://doi.org/10.3390/rs13071402.

Gleyzes, M., Perret, L., Kubik, P., 2012. PLEIADES system architecture and main performances. Int. Arch. Photogramm. Remote. Sens. Spat. Inf. Sci. XXXIX-B1, 537–542. https://doi.org/10.5194/isprarchives-XXXIX-B1-537-2012.

Grizonnet, M., Michel, J., Poughon, V., Inglada, J., Savinaud, M., Cresson, R., 2017. Orfeo Toolbox: open source processing of remote sensing images. Open Geospat. Data Softw. Stand. 2 (1), 15. https://doi.org/10.1186/s40965-017-0031-6.

Hagolle, O., Huc, M., Villa, Pascual D., Dedieu, G., 2015. A multi-temporal and multi-spectral method to estimate aerosol optical thickness over land, for the atmospheric correction of FormoSat-2, LandSat, VENUS and Sentinel-2 images. Remote Sens. 7 (3). https://doi.org/10.3390/rs70302668, 2668–269.

Hapke, B., 1981. Bidirectional reflectance spectroscopy: 1. Theory. J. Geophys. Res. 86 (B4), 3039–3054. https://doi.org/10.1029/JB086iB04p03039.

Hapke, B., 2012. Theory of Reflectance and Emittance Spectroscopy. Cambridge University Press. https://doi.org/10.1017/CBO9781139025683, 529 pp.

Hartman, B., Domingue, D., 1998. Scattering of light by individual particles and the implications for models of planetary surfaces. Icarus 131, 421–448. https://doi.org/10.3847/1538-3881/ab5e82.

Helfenstein, P., Shepard, M.K., 2011. Testing the Hapke photometric model: improved inversion and the porosity correction. Icarus 215 (1), 83–100. https://doi.org/10.1016/j.icarus.2011.07.002.

Jacquemoud, S., Baret, F., Hanocq, J.F., 1992. Modeling spectral and bidirectional soil reflectance. Remote Sens. Environ. 41 (2), 123–132. https://doi.org/10.1016/0034-4257(92)90072-R.

Jiao, Z., Schaaf, C.B., Dong, Y., Román, M., Hill, M.J., Chen, J.M., Wang, Z., Zhang, H., Saenz, E., Poudyal, R., Gatebe, C., Bréon, F.M., Li, X., Strahler, A., 2016. A method for improving hotspot directional signatures in BRDF models used for MODIS. Remote Sens. Environ. 186, 135–151. https://doi.org/10.1016/j.rse.2016.08.007.

Johnson, J.R., Shepard, M.K., Grundy, W.M., Paige, D.A., Foote, E.J., 2013. Spectrogoniometry and modeling of Martian and lunar analog samples and Apollo soils. Icarus 223 (1), 383–406. https://doi.org/10.1016/j.icarus.2012.12.004.
Kugler, B., 2021. Analyse massive d'images multi-angulaires hyperspectrales de la planète Mars par régression inverse de modèles physiques. PhD thesis. Université Grenoble Alpes, 149 pp. https://www.theses.fr/2021GRALM021.

Kugler, B., Forbes, F., Douté, S., 2022. Fast Bayesian inversion for high dimensional inverse problems. Stat. Comput. 32, 31. https://doi.org/10.1007/s11222-021-10019-5.

Labarre, S., 2017. Caractérisation et modélisation de la rugosité multi-échelle des surfaces naturelles par télédétection dans le domaine solaire. PhD thesis. Université Sorbonne Paris Cité, 228 pp. http://www.theses.fr/2017USPCC231.

Labarre, S., Ferrari, C., Jacquemoud, S., 2017. Surface roughness retrieval by inversion of the Hapke model: a multiscale approach. Icarus 290, 63–80. https://doi.org/10.1016/j.icarus.2017.02.030.

Labarre, S., Jacquemoud, S., Ferrari, C., Delorme, A., Derrien, A., Grandin, R., Galludin, M., Lemaître, F., Métois, M., Pirot-Deseilligny, M., Rupnik, E., Tanguy, B., 2019. Retrieving soil surface roughness with the Hapke photometric model: confrontation with the ground truth. Remote Sens. Environ. 225, 1–15. https://doi.org/10.1016/j.rse.2019.02.014.

Lebègue, L., Cazala-Hourcade, E., Languille, F., Artigues, S., Melet, O., 2020. CO3D, a worldwide one one-meter accuracy DEM for 2025. In: International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XLIII-B1-2020, pp. 299–304. https://doi.org/10.5194/isprs-archives-XLIII-B1-2020-299-2020.

Li, L., Mu, X., Qi, J., Pisek, J., Roosjen, P., Yan, G., Huang, H., Liu, S., Baret, F., 2021. Characterizing reflectance anisotropy of background soil in open-canopy plantations using UAV-based multiangular images. ISPRS J. Photogramm. Remote Sens. 177, 263–278. https://doi.org/10.1016/j.isprsjprs.2021.05.007.

Lowe, D.G., 2004. Distinctive image features from scale-invariant keypoints. Int. J. Comput. Vis. 60 (2), 91–110. https://doi.org/10.1023/B:VISI000029664.99615.94.

Lucas, S., Ribardiere, M., Pacanowski, R., Barla, P., 2023. A micrograin bsdf model for the rendering of porous layers. In: SIGGRAPH Asia 2023 Conference Papers, 40, pp. 1–10. https://doi.org/10.1145/3610548.3618241.

Lucey, P.G., 1998. Model near-infrared optical constants of olivine and pyroxene as a function of iron content. J. Geophys. Res.: Planets 103 (E1), 1703–1713. https://doi.org/10.1029/97JE03145.

Lucey, P.G., Riner, M.A., 2011. The optical effects of small iron particles that darken but do not redden: evidence of intense space weathering on mercury. Icarus 212 (2), 451–462. https://doi.org/10.1016/j.icarus.2011.01.022.

Lucht, W., Schaa, C.B., Strahler, A.H., 2000. An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Trans. Geosci. Remote Sens. 38 (2), 977–998. https://doi.org/10.1109/36.841980.

McGuire, A.F., Hapke, B.W., 1995. An experimental study of light scattering by large, irregular particles. Icarus 113 (1), 154–155. https://doi.org/10.1006/icar.1995.1012.

Nag, S., Hewagama, T., Georgiev, G.T., Pasigot, B., Aslam, S., Gatebe, C.K., 2017. Multispectral snapshot imagers onboard small satellite formations for multi-angular remote sensing. IEEE Sensors J. 17 (16), 5252–5268. https://doi.org/10.1109/JSEN.2017.2717384.

Nelder, J.A., Mead, R., 1965. A simplex method for function minimization. Comput. J. 7 (4), 308–313. https://doi.org/10.1093/comjnl/7.4.308.

Nicodemus, F.E., 1965. Directional reflectance and emissivity of an opaque surface. Appl. Opt. 4 (7), 767–775. https://doi.org/10.1364/AO.4.000767.

Pilorget, C., Fernando, J., Ehlmann, B.L., Douté, S., 2015. Photometry of particulate mixtures: what controls the phase curve? Icarus 250, 188–203. https://doi.org/10.1016/j.icarus.2014.11.036.

Pilorget, C., Fernando, J., Ehlmann, B.L., Schmidt, F., Hiroi, T., 2016. Wavelength dependence of scattering properties in the VIS–NIR and links with grain-scale physical and compositional properties. Icarus 267, 296–314. https://doi.org/10.1016/j.icarus.2015.12.029.

Pinty, B., Verstraete, M.M., Dickinson, R.E., 1989. A physical model for predicting bidirectional reflectances over bare soil. Remote Sens. Environ. 27 (3), 273–288. https://doi.org/10.1016/0034-4257(89)90088-6.

Pinzuti, P., Humler, E., Manighetti, I., Gaudemer, Y., 2013. Petrological constraints on melt generation beneath the Asal rift (Djibouti) using quaternary basalts. Geochem. Geophys. Geosyst. 14 (8), 2932–2953. https://doi.org/10.1002/ggge.20187.

Protopapa, S., Olkin, C.B., Grundy, W.M., Li, J.Y., Verbiscer, A., Cruikshank, D.P., Gautier, T., Quirico, E., Cook, J.C., Reuter, D., Howett, C.J.A., Stern, A., Beyer, R.A., Porter, S., Young, L.A., Weaver, H.A., Enrico, K., Dalle Ore, C.M., Scipioni, F., Singer, K., 2020. Disk-resolved photometric properties of Pluto and the coloring materials across its surface. Astron. J. 159 (2), 74. https://doi.org/10.3847/1538-3881ab5e82.

Rouquié, B., Hagolle, O., Bréon, F.M., Boucher, O., Desjardins, C., Rémy, S., 2017. Using Copernicus atmosphere monitoring service products to constrain the aerosol type in the atmospheric correction processor MAJA. Remote Sens. 9 (12), 1230. https://doi.org/10.3390/rs9121230.

Rupnik, E., Daakir, M., Pirotot, Deseilligny, M., 2017. MicMac – a free, open-source solution for photogrammetry. Open Geospat. Data Softw. Stand. 2 (1), 14. https://doi.org/10.1186/s40965-017-0027-2.

Rupnik, E., Pirotot-Deseilligny, M., Delorme, A., 2018. 3D reconstruction from multi-view VHR-satellite images in MicMac. ISPRS J. Photogramm. Remote Sens. 139, 201–211. https://doi.org/10.1016/j.isprsjprs.2018.03.016.

Sato, H., Robinson, M.S., Hapke, B., Denevi, B.W., Boyd, A.K., 2014. Resolved Hapke parameter maps of the moon. J. Geophys. Res.: Planets 119 (8), 1775–1805. https://doi.org/10.1002/2013JE004580.

Schaepman-Strub, G., Schaepman, M.E., Painter, T., Dangel, S., Martonchik, J., 2006. Reflectance quantities in optical remote sensing—definitions and case studies. Remote Sens. Environ. 103 (1), 27–42. https://doi.org/10.1016/j.rse.2006.03.002.

Schmidt, F., Bourguignon, S., 2019. Efficiency of BRDF sampling and bias on the average photometric behavior. Icarus 317, 10–26. https://doi.org/10.1016/j.icarus.2018.06.025.

Schmidt, F., Fernando, J., 2015. Realistic uncertainties on Hapke model parameters from photometric measurement. Icarus 260, 73–93. https://doi.org/10.1016/j. icarus.2015.07.002.

Shepard, M.K., Helfenstein, P., 2007. A test of the Hapke photometric model. J. Geophys. Res.: Planets 112 (E3). https://doi.org/10.1029/2005JE002625. E03001.

Shiltz, D.J., Bachmann, C.M., 2023. An alternative to Hapke's macroscopic roughness correction. Icarus 390, 115240. https://doi.org/10.1016/j.icarus.2022.115240.

Smittarello, D., Grandin, R., De Chabalier, J.B., Doubre, C., Deprez, A., Masson, F., Socquet, A., Saad, I.A., 2016. Transient deformation in the Asal-Ghoubbet rift (Djibouti) since the 1978 diking event: is deformation controlled by magma supply rates? J. Geophys. Res. Solid Earth 121 (8), 6030–6052. https://doi.org/10.1002/2016JB013069.

Sobol, I.M., 1967. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comput. Math. Math. Phys. 7 (4), 86–112. https://doi.org/10.1016/0041-5553(67)90144-9.

Souchon, A.L., Pinet, P.C., Chevrel, S.D., Daydou, Y.H., Baratoux, D., Kurita, K., Shepard, M.K., Helfenstein, P., 2011. An experimental study of Hapke's modeling of natural granular surface samples. Icarus 215 (1), 313–331. https://doi.org/10.1016/j.socquet.2011.06.023.

Souchon, A.L., Besse, S., Pinet, P.C., Chevrel, S.D., Daydou, Y.H., Josset, J.L., D'Uston, L., Haruyama, J., 2013. Local spectrophotometric properties of pyroclastic deposits at the Lavoisier lunar crater. Icarus 225 (1), 1–14. https://doi.org/10.1016/j. icarus.2013.03.004.

Starck, J., Murtagh, F., 2006. Astronomical Image and Data Analysis. Springer-Verlag. https://doi.org/10.1007/978-3-540-33025-7, 338 pp.

Stieltjes, L., 1980. Carte géologique du rift d'Asal (dépression Afar, est-africain), BRGM (Orléans), CNRS (Paris), 1/50 000, 88 × 114 cm.

Wang, Z., Wu, Y., Blewett, D.T., Cloutis, E.A., Zheng, Y., Chen, J., 2017. Submicroscopic metallic iron in lunar soils estimated from the in situ spectra of the Chang'E-3 mission. Geophys. Res. Lett. 44 (8), 3485–3492. https://doi.org/10.1002/2017GL072652.

Wanner, W., Li, X., Strahler, A.H., 1995. On the derivation of kernels for kernel-driven models of bidirectional reflectance. J. Geophys. Res. Atmos. 100 (D10), 21077–21089. https://doi.org/10.1029/95JD0237.

Wu, Y., Gong, P., Liu, Q., Chappell, A., 2009. Retrieving photometric properties of desert surfaces in China using the Hapke model and MISR data. Remote Sens. Environ. 113 (1), 213–223. https://doi.org/10.1016/j.rse.2008.09.006.

Yang, G.J., Zhao, C.J., Huang, W.J., Wang, J.H., 2011. Extension of the Hapke bidirectional reflectance model to retrieve soil water content. Hydrol. Earth Syst. Sci. 15 (7), 2317–2326. https://doi.org/10.5194/hess-15-2317-2011.