Retrieving the characteristics of slab ice covering snow by remote sensing

We present an effort to validate a radiative transfer model previously developed, and an innovative bayesian inversion method designed to retrieve the properties of slab ice covered surfaces. This retrieval method is adapted to satellite data, and is able to provide uncertainties on the results of the inversions. We focused in this study on surfaces composed of a pure slab of water ice covering an optically thick layer of snow. We see sought to retrieve the roughness of the ice/air interface, the thickness of the 5 slab layer and the mean grain diameter of the underlying snow. Numerical validations have been conducted on the method, and showed that if the thickness of the slab layer is above 5 mm and the noise on the signal is above 3%, then it is not possible to invert the grain diameter of the snow. On the contrary, the roughness and the thickness of the slab can be determined even with high levels of noise up to 20%. Experimental validations have been conducted on spectra collected from laboratory samples of water ice on snow using a spectro-gonio-radiometer. The results are in agreement with the numerical validations, and show 10 that a grain diameter can be correctly retrieved for low slab thicknesses, but not for bigger ones, and that the roughness and thickness are correctly inverted in every case.


Introduction
Various species of ice are present throughout the solar system, from water ice and snow on Earth to nitrogen ice on Triton (Zent et al., 1989), not to forget carbon dioxide ice on Mars (Leighton and Murray, 1966).The physical properties of the cover also have an impact on the energy balance; for example, the albedo depends on the grain size of the snow (Dozier et al., 2009;Negi and Kokhanovsky, 2011;Picard et al., 2009;Mary et al., 2013), on the roughness of the interface (Lhermitte et al., 2014), on the presence or lack of impurities and the physical properties of impurities (Dumont et al., 2014).The study and monitoring of theses parameters is key to constraining the energy balance of a planet.
Radiative transfer models have proven essential for retrieving such properties (Zege et al., 2008;Negi and Kokhanovsky, 2011) and their evolution at a large scale, and different families exist.Ray-tracing algorithms, such as those described in Picard et al. (2009) for snow, Pilorget et al. (2013) for compact polycrystalline ice or Muinonen et al. (2009) for particulate media such as rough ice grains in an atmosphere, simulate the complex path of millions of rays into the surface.Such modelings are booming due to the positive comparison between models and exact calculations (e.g., Muinonen et al., 2012;Mishchenko et al., 2015).Analytical solutions of the radiative transfer in homogeneous granular media have been developed, for example, by Shkuratov et al. (1999) and Hapke (1981).They are fast, but when the surface cannot be described statistically as a mono-layer, they must be combined with another family of techniques such as discrete ordinate methods like DISORT (Stamnes et al., 1988).These methods have been widely studied on Earth snow (Carmagnola et al., 2013;Dozier et al., 2009;Dumont et al., 2010;Painter and Dozier, 2004) and other planetary cryospheres (Appéré et al., 2011;Eluszkiewicz and Moncet, 2003), modeling a granular surface.Compact polycrystalline ice has, however, been recognized to exist on several ob-jects: CO 2 on Mars (Kieffer and Titus, 2001;Eluszkiewicz et al., 2005), N 2 on Triton and Pluto (Zent et al., 1989;Eluszkiewicz and Moncet, 2003) and probably SO 2 on Io (Eluszkiewicz and Moncet, 2003), as suggested by the very long light path lengths measured, over several centimeters to decimeters (Eluszkiewicz, 1993;Quirico et al., 1999;Douté et al, 1999Douté et al, , 2001)).In particular, the Martian climate is mostly controlled by a seasonal CO 2 cycle that results in the condensation and sublimation of deposits constituted of a layer of nearly pure CO 2 ice up to 1 m thick, possibly contaminated with H 2 O ice and regolith dust (Leighton and Murray, 1966;Eluszkiewicz et al., 2005).The monitoring of the evolution of these deposits' composition would bring key pieces of information on the water ice cycle on Mars and on the dust transport to the surface.
Compact slabs have very different radiative properties from closely packed granular media, and radiative transfer models have been developed to study their characteristics (e.g., Mullen and Warren, 1988;Jin et al., 1994;Perovich, 1996;Jin et al., 2006) in the case of sea or lake ice.We developed an approximated model (Andrieu et al., 2015d) that has the capability of being able to model a layer of ice covering a surface with radically different optical properties, for instance a different refractive index, unlike its predecessors.It was originally designed for the study of Martian seasonal ice deposits, using massive spectro-imaging datasets, such as the OMEGA (Bibring et al., 2004) or CRISM (Murchie et al., 2007) datasets.For this purpose, it is semi-analytic and implemented to optimize the computation time.However, the algorithm was built to be adaptable to any other spectroscopic data, from terrestrial water ice laboratory measurements, as is the case of this work, to the study of SO 2 ice on Io or N 2 ice on Pluto using the NIMS (Carlson et al., 1992) and RALPH (Reuter et al., 2009) datasets respectively.
In the present article, we test the accuracy of this approximated model on laboratory spectroscopic measurements of the bidirectional reflectance distribution function (BRDF) of pure water ice on top of snow.At the same time, we present an innovative Bayesian inversion method that was developed to retrieve the properties of solar system compact ice using satellite spectro-imaging data.In this paper, the term "inversion" is used and means "solving the inverse problem".The slabs that are studied contain no impurity, and the surface properties we seek to retrieve are the thickness of the ice, the roughness of the ice-air interface and the grain diameter of the underlying snow.The main goals of this work are (i) to test the ability of the model to reproduce reality and (ii) to propose an inversion framework able to retrieve surface ice properties, including uncertainties, in order to demonstrate the applicability of the approach to satellite data.
We present a set of spectro-goniometric measurements of different water ice samples put on top of snow using the spectro-radiogoniometer described in Brissaud et al. (2004).Three kind of experiments were conducted.First, the BRDF of only one snow layer was measured, and then it was mea-sured again after adding a slab ice layer at the top.The objective was to test the effect of an ice layer at the top of the snow on the directivity of the surface.Second, the specular lobe was closely investigated, at high angular resolution, at the wavelength of 1.5 µm, where ice behaves as a very absorbing media.Finally, the bidirectional reflectance was sampled at various geometries on 61 wavelengths ranging from 0.8 to 2.0 µm.In order to validate the model, we did qualitative tests to demonstrate the relative isotropization of the flux.We also conducted quantitative assessments by using a Bayesian inversion method in order to estimate the sample thickness, the surface roughness and the snow grain diameter from the radiative measurements only.A comparison between the retrieved parameters and the direct independent measurements allowed us to validate the model.
The inversion algorithm that is tested is based on lookup tables (LUTs) that minimize the computation time of the direct model.The solution is formulated as a probability density function, using Bayesian formalism.This strategy is very useful for analyzing hyperspectral images.The thickness of ice estimated from the inversion is compared to real direct measurements.In addition, the power distribution in the specular lobe, which is determined by the roughness of the surface, is adjusted to demonstrate that the model is able to reasonably fit the data with a consistent roughness value.

Description of the model
The model by Andrieu et al. (2015d) is inspired by an existing one described in Hapke (1981) and Douté and Schmitt (1998), which simulates the bidirectional reflectance of stratified granular media.It has been adapted to compact slabs, contaminated with pseudo-spherical inclusions, and a rough top interface.In the context of this work, we suppose a layer of pure slab ice, overlaying an optically thick layer of granular ice, as described in Fig. 1.The roughness of the first interface is described using the probability density function of orientations of slopes defined in Hapke (1984).This distribution of orientations is fully described by a parameter θ , which can be interpreted as a mean slope parameter at the surface, in the case of small angles.The ice matrix is described using its optical constants and its thickness.
Figure 2 illustrates the general principle of the model.The simulated bidirectional reflectance results from two separate contributions: specular and diffuse.The specular contribution in the model is estimated from the roughness parameter, the optical constants of the matrix and the apertures of the light source and the detector.In practical applications, the optical constants of the ice matrix and the optical apertures are known parameters.The specular contribution for a given geometry only depends on the orientation of slopes at the surface, which is fully determined by θ .The total reflection coefficient at the first rough interface is obtained by integrating specular contributions in every emergent direction, The Cryosphere, 10, 2113-2128, 2016 Illustration of the radiative transfer in the surface medium.Anisotropic transits are represented in red.F is the incident radiation flux, R spec and R Diff are the specular and diffuse contributions to the reflectance of the surface respectively, r s is the Lambertian reflectance of the granular substrate and R 0 and T 0 are the total reflection and transmission factors of the slab layer respectively.A prime indicates an anisotropic transit.The reflection and transmission factors are different in the cases of isotropic or anisotropic conditions.The granular and slab layers are artificially separated in this figure to help the understanding of the coupling between the two layers.Top: illustration of the reflections and transmission at the first interface, used in the calculations of R spec and the determination of the amount of energy injected into the surface.z is the normal to the surface, W f the local normal to a facet, i and e are the incidence and emergence angles respectively and e f is the local emergence angle for a facet.Each different orientation of a facet will lead to a different transit length in the slab.A more detailed description can be found in Andrieu et al. (2015d).
at a given incidence.This gives the total amount of energy transmitted into the system constituted of the contaminated slab and the substrate.The diffuse contribution is then estimated through solving the radiative transfer equation inside this system under various hypotheses.The following considerations are made.(i) The first transit through the slab is anisotropic due to the collimated radiation from the source, and due to there being an isotropization at the second rough interface (i.e., when the radiation reaches the semi-infinite substrate).For the refraction and the internal reflection, every following transit is considered isotropic.(ii) The geometrical optics is valid.The reflection and transmission factors of the layers are obtained using an analytical estimation of the Fresnel coefficients described in Chandrasekhar (1960) and Douté and Schmitt (1998), as well as a simple statistical approach, detailed in Andrieu et al. (2015d).The contribution of the semi-infinite substrate is described by its singlescattering albedo.Finally, as the slab layer is under a collimated radiation from the light source, and under a diffuse radiation from the granular substrate, the resulting total bidirectional reflectance is computed using adding-doubling formulas (Stamnes et al., 1988;Douté and Schmitt, 1998;Van de Hulst, 2012).
In this work, the radiative transfer model described in Andrieu et al. (2015d) is used to simulate the reflectance factor spectrum of a pure slab of water ice covering a layer of snow, as represented in Fig. 1.This spectrum may vary in terms of three parameters: (i) the roughness θ of the slab ice surface, which characterizes the power distribution in the specular lobe; (ii) the thickness h of the slab ice layer, which determines the absorption in the ice layer and (iii) the grain diameter s of the snow, which determines the absorption in the snow.

Spectro-radiogoniometer
The bidirectional reflectance spectra were measured using the spectro-radiogoniometer from Institut de Planétologie et d'Astrophysique de Grenoble, fully described in Brissaud et al. (2004).We collected spectra in the near-infrared spectrum at incidences ranging from 40 to 60 • , emergence angles from 0 to 50 • and azimuth angles from 0 to 180 • .The sample is illuminated with a 200 mm large monochromatic beam (divergence < 1 • ), and the near-infrared spectrum covering the range from 0.800 to 4.800 µm is measured by an InSb photovoltaic detector.This detector has a nominal aperture of 4.2 • , which results in a field of view on the sample that has a diameter of approximately 20 mm.The minimum angular sampling of illumination and observation directions is 0.1 • , with a reproducibility of 0.002 • .In order to avoid azimuthal anisotropy, the sample is rotated during the acquisition.The sample rotation axis may be very slightly misadjusted, result-F.Andrieu et al.: Retrieving the characteristics of slab ice covering snow by remote sensing ing in a notable angular drift on the emergence measured up to 1 • .

Ice BRDF measurements
The ice samples were obtained by sawing artificial columnar water ice into circular sections 20 cm in diameter.These sections were put on top of an optically thick layer of compacted snow that was collected in Arselle, in the French Alps.The spectral measurements were conducted in a cold chamber at 263 K.However, the ice and the snow were unstable in the measurement's environment, due to the dryness of the chamber's atmosphere.The grain size of the snow showed evolution, and the thickness of a given slab showed a decrease of 0.343 mm day −1 .Each sample needed an acquisition time of 10 h.For each measurement, the ice slab was sliced, and its thickness was measured in five different locations.It was then set on top of the snow sample, and this system was put into rotation in the spectro-radiogoniometer for the measurement.The sample completes a full rotation (10 s) during the measurement of the reflectance at one wavelength and one geometry.As the surface is not perfectly planar, the measured thickness is not constant.This results in an 2σ standard deviation in the measurement of the thickness that ranges from 0.54 to 2.7 mm in our study, depending on the sample.

Specular contribution
The specular reflectance was measured on a slab sample on top of Arselle snow, 12.51 mm thick.This sample is described as sample 3 in the next paragraph.The illumination was at an incidence angle of 50 • , and 63 different emergent geometries were sampled, ranging from 45 to 55 • in emergence and from 170 to 180 • in azimuth.The specular lobe is sampled every 1 • in emergence and azimuth, and within 47 and 53 • in emergence and 175 and 180 • in azimuth.

Ice on snow diffuse reflectance spectra
The diffuse contribution was measured on three samples of different slab thickness.The three thicknesses were measured on different locations of the samples with a caliper before doing the spectro-goniometric measurement, resulting in h 1 = 1.42±0.47mm, h 2 = 7.45±0.84mm and h 3 = 12.51± 2.7 mm, respectively, for samples 1, 2 and 3, with errors at 2σ .Sixty-one wavelengths were sampled ranging from 0.8 to 2.0 µm.Spectra were collected on 39 different points of the BRDF for the incidence, emergence and azimuth angles: 40, 50, 60 • , [0, 10, 20] and 0, 45, 90, 140, 160, 180 • , respectively.This set of angles results in only 39 different geometries because the azimuthal angle is not defined for a nadir emergence.
Diffuse reflectance spectra of natural snow only were also measured before putting a slab on top of the snow.The objective was to estimate the effect of a slab layer on the BRDF.

Method
We designed an inversion method aimed at massive data analysis.This method consists of two steps: first, the generation of a synthetic database that is representative of the variability in the model, and then the comparison with actual data.To generate the synthetic database, we used optical constants for water ice at 270 K.The 7 K difference between the actual temperature of the room and the temperature assumed for the optical constants has a negligible effect.We combined the datasets of Warren and Brandt (2008) and Schmitt et al. (1998), making the junction at 1 µm, the former set for the shorter wavelengths and the latter for the wavelengths larger than 1 µm.
In order to validate the model on the specular reflection from the slab, we chose to use the reflectance at 1.5 µm, where the ice is very absorptive.Figures 7 and 8 clearly demonstrate that there is a negligible diffuse contribution in geometry outside the specular lobe from the sample with a 12.51 mm thick pure slab.Thus, the roughness parameter θ is the only one impacting the reflectance in the model.We chose to invert this parameter first and validate the specular contribution.
We then focused on the validation in the spectral domain, for the diffuse contribution.We used the estimation of the roughness parameter θ obtained earlier and the spectral data in order to estimate the slab thickness and the grain size of the snow substrate.To do this, we assumed that the roughness was not changing significantly enough to have a notable impact on diffuse reflectance from one sample to another.This assumption is justified by the fact that the different columnar ice samples were made the same way, as flat as possible, and the low value of θ retrieved as discussed in the next section.It is confirmed by the results of Sect.4.2, which suggest a very low roughness, as expected.Such low roughness parameters have a negligible influence on the amount of energy injected into the surface.

Inversion strategy
The inversion consists in estimating the model parameters m (i.e., the slab thickness, the roughness parameter, the snow grain diameter) from the models F (m)) (the reflectance simulations) that best reproduce the data d (the reflectance observations).Looking for these parameters sets is called the inverse problem.Tarantola and Valette (1982) showed that this inverse problem can be mathematically solved by considering each quantity as a probability density function (PDF).In nonlinear direct problems, the solution may not be analytically approached.Nevertheless, it is possible to sample the solutions' PDF with a Monte Carlo approach as shown in Mosegaard and Tarantola (1995), but this solution is very time-consuming.
The actual observation is considered as a priori information on the data ρ D (d) in the observation space D. It is as-sumed to be an N-dimension multivariate Gaussian distribution G(d mes , C), centered on d mes with a covariance matrix C. The measurements at any given wavelength/geometry are supposed to be independent from each other, as each measurement of one wavelength, at one geometry is done individually.The matrix C is thus assumed to be diagonal and its diagonal elements C ii are σ 2 1 , . .., σ 2 N , with σ i being the standard deviations corresponding to the uncertainties of each measurement.The a priori information on model parameters ρ M (m) in the parameters space M is independent of the data and corresponds to the state of null information µ D (d) if no information is available on the parameters.We consider a uniform PDF in their definition space M. The a posteriori PDF in the model space σ M (m) as defined by Bayes' theorem (Tarantola and Valette, 1982) is where k is a constant and L(m) is the likelihood function, where θ (d | m) is the theoretical relationship of the PDF for d given m.We do not consider errors of the model itself, so Therefore, the likelihood is simplified into and in the case of uniform a priori information ρ M (m), the a posteriori PDF is This expression is explicitly where the superscript t is the transpose operator that applies to (F (m) − d mes ).The factor k is adjusted to normalize the PDF.The mean value of the estimated parameter can be computed by and the standard deviation, In order to speed up the inversion strategy but keep the advantage of the Bayesian approach, we choose to sample the parameter space M with regular and reasonably fine steps, noted i.The likelihood for each element is (8) The derivation of the a posteriori PDF with such formalism for specular lobe inversion and for spectral inversion is explained in the next sections.

Specular lobe
To study the specular lobe, we have to consider the whole angular sampling of the spot as a single data measurement.Similar to the "pixel" (contraction of picture element), we choose to define the "angel" (contraction of angular element) as a single element in a gridded angular domain.Interestingly, angel also refers to a supernatural being represented in various forms of glowing light.A single angel measurement could not well constrain the model, even at different wavelengths.Instead, a full sampling around the specular lobe should be enough, even at one single wavelength.We chose a wavelength where the diffuse contribution was negligible in order to simplify the inversion strategy.We chose to focus on the 1.5 µm wavelength, as we measured a penetration depth lower than 1 mm and thus much lower than the thickness of the used sample for this channel.We first generated a synthetic database (lookup table), using the direct radiative transfer model.We simulated spectra in the same geometrical conditions, for a 12.5 mm thick ice layer over a granular ice substrate constituted of 1000 µm wide grains.These two last parameters are not important since the absorption is so high in ice, such that the main contribution is from the specular reflection, and the diffuse contribution is negligible (the penetration depth inside a water ice slab at the 1.5 µm wavelength is lower than 1 mm).
The sampling of the parameter space, i.e., the lookup table, must correctly represent the variability of the model according to its parameters.For this study, we sampled the roughness parameter from 0.1 to 5 • with a constant step d θ = 0.01 • .We used a likelihood function L defined in Eq. ( 8), where d sim and d mes are n geom -element vectors, with n geom the number of angels (63 in this study).They respectively represent the simulated and measured reflectance at a given wavelength in every geometry.C is a n geom × n geom matrix.It represents the uncertainties of the data.In this case, we considered each wavelength independently, thus generating a diagonal matrix, containing the level of errors given by the technical data of the instrument described by Brissaud et al. (2004).It corresponds at this wavelength to 2 % of the signal.The roughness parameter θ returned by the inversion is described by its normalized PDF: The value θ (i) in the database with the highest probability (maximum likelihood) corresponds to the sampling step that is closest to the maximum of the a posteriori PDF.If the PDF is close to a Gaussian, then the solution to the inverse problem can be estimated by its mean, and associated standard deviations, We give error bars on the results that correspond to two standard deviations, and thus a returned value for θ that is

Diffuse spectra
When out of the specular lobe, the radiation is controlled by the complex transfer through the media (slab ice and bottom snow).The experimental samples were made of pure water slab ice, without impurity.We generated the lookup table for every measurement geometry at very high spectral resolution (4.10 −2 nm), and then downsampled it at the resolution of the instrument (2 nm).We sampled the 17 085 combinations of two parameters for the 39 different geometries: p 1 , the thickness of the slab from 0 to 20 mm (noted i = [1, 201]) every 0.1 mm (noted dp 1 ), and p 2 , the grain size of the granular substrate from 2 to 25 µm every 1 µm and from 25 to 1500 µm every 25 µm (noted j = [1, 85] and the corresponding dp 2 (j )).The parameters' space is thus irregularly paved with dp(i, j ) = dp 1 .dp 2 (j ).
For the inversion, we used the same method as previously described, with a likelihood function L that is written as in Eq. ( 8).Two different strategies were adopted.First, we inverted each spectrum independently.Thirty-nine geometries were sampled (described in Sect.3.2), and thus we conducted 39 inversions for each sample.This time d sim and d mes are thus respectively the simulated and measured spectra.Therefore d sim and d mes are n b -element vectors, where n b is the number of bands (61 in this study) and C is a n b × n b matrix.As previously done (see Sect. 4.2), we considered each wavelength independently, thus generating a diagonal matrix, containing the level of errors given by the technical data of the instrument given by Brissaud et al. (2004).The error is a percentage of the measurement, and thus C is different for every inversion.
Secondly, we inverted the BRDF as a whole, for each sample.For this method, d sim and d mes are the simulated and measured BRDF respectively and are thus n b × n geomelement vectors (2379 in this study), where n b is the number of bands (61 in this study) and n geom is the number of geometries (39 in this study) sampled.C is a n b × n geom × n b × n geom diagonal matrix, containing the errors of the data.We represent the results the same way as previously, but there are two parameters to inverse.For the sake of readability, we plot the normalized marginal probability density function for each parameter.We present here the general method for the inversion of n p = 2 parameters: the slab thickness and the grain size of the substrate.The PDF for the two parameters p is described by For a given parameter p 1 , the marginal PDF of the solution is with L (i) = j L(i, j )dp 2 (j ).The value p 1 (i) in the database with the highest probability (maximum likelihood) corresponds to the sampling step that is closest to the maximum of the a posteriori PDF.The marginal PDF can be described by the mean, and the associated standard deviation, As for the roughness parameter, we give error bars on the results that correspond to two standard deviations, and thus a returned value for p 1 that is

Numerical validations of the inversion method
In order to numerically validate the inversion method described above, two kinds of tests were conducted.First, we applied Gaussian noise and inverted every spectrum in the synthetic spectral database.We show that with a negligible noise, the parameters are always correctly retrieved with negligible uncertainties, and as the level of noise on the data increases, so do the uncertainties on the results.Secondly, we generated spectra for parameters that were not sampled in the database and tried to recover their characteristics successfully.On Fig. 4 each curve corresponds to a the normalized sum of 1000 a posteriori PDFs for the grain diameter of the underlying snow resulting from 1000 random noise draws of the same 2 % level.We name this normalized sum of PDF a "stack" in the remainder of this article.Figure 4a is obtained for a low slab thickness of 1 mm.In this case, the grain diameter of the snow can be correctly estimated: the PDF are The Cryosphere, 10, 2113-2128, 2016 centered on the correct value and the dispersion suggests an a posteriori uncertainty lower than the retrieved value.When the thickness of the slab layer increases, so does the a posteriori uncertainty on the estimation of the grain diameter.For a slab thickness of 5 mm (Fig. 4b), the a posteriori uncertainty is of the same order as the estimated value, meaning that the grain diameter cannot be retrieved.The grain diameter of the snow thus cannot be retrieved for slab thicknesses greater than 5 mm.
Figure 5a represents the stack of 1000 a posteriori PDFs for the thickness of the ice layer.These PDFs do not depend on the grain diameter of the snow, but only on the thickness itself and the level of noise.It shows that the thickness can be estimated, in the experimental conditions (2 % noise level), with an uncertainty of 2 % for lowest thicknesses to 5 % for highest ones.All obtained a posteriori PDFs for the thickness were very close to a Gaussian function.We were thus able to sum them up by their means and standard deviations, allowing us to plot, for example, the uncertainty of the thickness estimation as a function of the thickness that we want to estimate (Fig. 5b) and of the level of noise on the data (Fig. 6). Figure 5b shows the uncertainty (at 2σ ) on the estimation of the thickness of the slab layer as a function of the thickness itself, in the experimental conditions described by Brissaud et al. (2004), which means a 2 % noise level on the signal.This relative uncertainty does not depend on the thickness in the range of values tested.The low values for thicknesses below 1 mm are an effect of the discretization in the LUT: the thickness has been sampled every 0.1 mm.Below 1 mm, this sampling step is large relatively to the values itself and ranges from 10 to 100 %.The relative uncertainty that we expect to be about 5 % is then no longer measurable, and the value drops to 0.
Figure 6 shows the evolution of the a posteriori uncertainties for the estimations of thicknesses and grain diameters as a function of the noise level.For the grain diameters, a slab thickness of 2 mm has been used.The results show that with very low noise, i.e., lower than 0.5 %, the a posteriori uncertainties of the results are of the same order of magnitude, even for the grain diameters.When the level of noise increases, the uncertainties of the thickness estimations increase in the same proportions (Fig. 6b), unlike the uncertainties of grain diameters (Fig. 6a) that increase drastically with the noise level.The uncertainties of the grain diameters seem to be saturated for high noises.This effect is only an edge effect due to the size of the LUT: the dispersion of the a posteriori PDFs cannot get bigger than the range of values tested.
With the level of noise at 2 % as expected for the measured spectra (Brissaud et al., 2004), a posteriori uncertainties are expected to be about 5 % on the thickness, and should be lower than 50 % for the grain diameter for low thicknesses.This means that the method should be able to retrieve thicknesses with an uncertainty that corresponds to the level of noise, but cannot retrieve the grain diameters of the snow when the ice layer above is thicker than 5 mm.

Impact of a slab on the BRDF
Figure 3 shows the reflectance factor (the ratio between the bidirectional reflectance I /F of the surface and the reflectance of a perfectly Lambertian surface) vs. phase angle (angle between incident and emergent directions) of the snow and the snow covered with an ice slab 1.42 mm thick (sample 1).It illustrates the two most notable effects of a thin layer of slab ice on top of an optically thick layer of snow.The most intuitive effect is to lower the level of reflectance: it is due to absorption during the long optical path lengths in the compact ice matrix as the dependance of the reflectance on the www.the-cryosphere.net/10/2113/2016/The Cryosphere, 10, 2113-2128, 2016  phase angle is strongly attenuated by the addition of the ice layer.The second effect is that the radiation is more Lambertian than that of snow only.These data give credit to the first hypothesis of isotropization of the radiation formulated in the model (see Sect. 2).The description of the bottom granular layer as Lambertian, defined only by its single-scattering albedo, may be considered simplistic, but this dataset shows that a thin coverage of slab ice, even on a very directive material such as snow, is enough to strongly flatten the BRDF.

Specular lobe: roughness retrieval
We performed the inversion taking into account 63 angel measurements, but for the sake of readability, Fig. 7 represents only the reflectance in the principle plane.The shapes and the intensities in Fig. 7a are compatible, but the measurement and simulation are not centered at the same point.
The simulation is centered at the geometrical optics specular point (emergence 50 • and azimuth 180 • ), whereas the measurement seems to be centered around an emergence of 50.5 • .This could be due to slight misadjustment of the rotation axis of the sample in the instrument.This kind of misadjustment is common, and can easily result in a notable shift up to 1 • of the recorded measurement geometries.We simulated different possible shifts in this range, and found a maximum likelihood for the simulation represented in Fig. 7b for a shift of 0.5 • in emergence, as suggested by the first plot in Fig. 7a, and 0.2 • in azimuth.The measurements and the simulations corresponding to the maximum likelihood are represented in Fig. 8.The shape and the magnitude of the specular lobe are very well reproduced.Both lobes show a small amount of forward asymmetry.This asymmetry is not due to the sampling as it is also present when the simulation is not shifted (see the red curve in Fig. 7).It is due to an  increase in the Fresnel reflection coefficient when the phase angle increases for this range of geometries.Figure 9 shows the a posteriori PDF for the parameter θ .The simulation corresponding to the maximum likelihood was obtained with θ = 0.43 • .The inversion method gives a result with a close to Gaussian shape at θ = 0.424 • ± 0.046 • .Unfortunately, we have no direct measurement of θ because the experimental determination of the roughness of a sample requires a digital terrain model (DTM) of its surface.These DTMs are measured with a laser beam, and are thus very difficult to obtain for ice samples, as multiple reflections create false measures.Still, we find a low value, which is consistent with the production of slabs of columnar ice in the laboratory that are very flat, but still imperfect as described in the dataset.The average slope is compatible with a long-wavelength slope at the scale of the sample, demonstrating that the microscale was not important in our case.Indeed, for a sample that has a length L, a 1σ standard deviation on the thickness h can be attributed to a general slope ϑ = arctan h L due to a small error in the parallelism of the two surfaces of the slab.In the case of sample 3, L = 20 cm and h = 1.35 mm result in ϑ = 0.39 • , which is compatible with the roughness given by the inversion.We thus think that what we see is an apparent roughness due to a small general slope on the samples, and that the roughness at the surface is much lower than this value.Moreover, the value retrieved by the inversion is very well constrained as the probability density function is very sharp.This means that we have an a posteriori uncertainty on the result that is very low.The quality of the reproduction of the specular lobe by the model suggests that the surface slope description is a robust description despite its apparent simplicity.In particular, one single slope parameter is enough to describe this surface.

Example of individual geometries
To reproduce diffuse reflectance we used the results obtained with the specular measurements and assumed that the roughness of the samples did not change much between the experiments.The range of variations in roughness should be negligible in the spectral analysis.We simulated slabs over snow, with the grain diameter of the substrate and the thickness of the slab as free parameters.Figure 10 represents three examples of measured and simulated reflectance spectra at the maximum likelihood for three different geometries.We also represented the mismatch between these simulations and the observations.We find an agreement between the data and the model that is acceptable.Nevertheless, there seems to be a decrease in quality in the fits as the thickness increases.Figure 11 shows an example of the marginal PDF for the three samples that are associated with the previous fits.The thickness is well constrained as the marginal a posteriori probability density functions are sharp and very close to a Gaussian function.However, the grain diameter of the substrate seems to have a limited impact on the result since it is only slightly constrained.The marginal PDFs for the grain diameter of the substrate are broad, and thus the a posteriori relative uncertainties in the result are very high.Unfortunately, we have no reliable measurement of the grain diameter of the substrate, as it is evolving during the time of the measurements.Numerical tests show that the snow grain diameter is not accessible for slab thicknesses above 5 mm.The a posteriori PDFs for samples 2 and 3 are then not to be interpreted.

Results for 39 geometries
Figure 12 shows the measurements and the final result of the inversion of the thickness for the three samples, and for 39 measurement geometries independently.The data and the model are compatible.Still, the thickness of sample 1 is slightly overestimated.This may reveal a sensitivity limit of the model.The thickness of sample 3 seems underestimated.This could be partly due to the duration of the measurement: the slab sublimates as the measure is being taken.Moreover, the specular measurements were performed on that sample, increasing the duration of the experiment even more.The inversion points in Fig. 12 are sorted by increasing incidence and, for each incidence, by increasing azimuth.There seems to be an influence of the geometry on the returned result: it is particularly clear for sample 2. The estimated thickness  tends to increase with incidence and decrease with azimuth.This effect disappears for large thicknesses (sample 3).

Full BRDF inversion
Figure 13 shows the measure and the simulation corresponding to the maximum likelihood at the λ = 1.0 µm wavelength when conducting the inversion on the whole BRDF dataset for each sample.The relatively flat behavior of the radiation with the phase angle is reasonably well reproduced.The quality of the geometrical match increases with the thickness of the sample.This is consistent with the fact that a thicker slab will permit a stronger isotropization of the radiation.It is also consistent with the disappearance of the geometrical dependence on the estimation of large thicknesses noted in Fig. 12.The values of thicknesses returned by the inversion are displayed in Fig. 14a.They are also compatible with the data, and the results are close to the one given by independent inversions of each geometry (see Figs. 11 and 12).The grain diameter returned (see Fig. 14) for sample 1 is lower, but is compatible with the one given by independent measurements.For samples 2 and 3, the PDFs are not interpreted, as the grain diameter cannot be constrained by the method.

Discussion
The two main goals of this work were (i) to develop and validate an inversion method that is adapted to the treatment of massive and complex datasets such as satellite hyperspectral datasets, and (ii) to partially validate a previously developed radiative transfer model.The first criterion is the speed of the whole method, including the direct computation of the LUT and the inversion.The lookup tables used for this project were computed in 150 s for the roughness study (1763 wavelengths sampled, 30 933 spectra) and ∼ 2.5 h for the thickness and grain diameter study (33 186 wavelengths sampled, 666 315 spectra).The inversions themselves were performed in less than onetenth of a second for specular lobe and independent specwww.the-cryosphere.net/10/2113/2016/The Cryosphere, 10, 2113-2128, 2016 tral inversions, and 2 s for BRDF-as-a-whole inversions.Every calculation was computed on one Intel CPU with 4 GB RAM.It has to be noted that once the lookup table has been created, an unlimited number of inversions can be conducted.This means that this method satisfies the speed criterion for the study of massive and complex datasets.For inversions over very large databases, the code has been adapted to GPU parallelization.It is also possible to increase the speed of the calculation of the lookup tables by means of multi-CPU computing.This Bayesian method has been designed to deal with the particular case of a direct model for which the computation time per simulation decreases with the size of the database.In our case, the direct radiative transfer algorithm is slow to simulate only one spectrum (∼ 1 s), but becomes fast for the calculation of large spectral databases (75 spectra per second in the case of the grain diameter and thickness study in this work).This particularity makes the usual Bayesian approaches such as Markov chain Monte Carlo methods inefficient, and in contrast, makes the method presented in this paper efficient.
A second aspect is the reliability of the inversion method, regardless of the direct model.For a given level of measurement errors, the user shall know the quality of the retrieval of any parameter.The Bayesian statistics in our method allowed us to determine that the thicknesses estimated in this work were reliable, with a 5 % uncertainty.Moreover, for the radiative transfer model used in this work (see Sect. 2), and in the experimental conditions described in Sect.3, we could determine on synthetic cases (see Sect. 4.4) that a 5 % uncertainty of ice thickness estimation should be expected, and that the grain diameter of the underlying snow could not be determined for ice thicknesses higher than 5 mm.The experimental results on the thickness were in agreement with these estimations.
A third point to be discussed is the capability of the model to reproduce reality.Section 5 showed that every thickness estimation was in agreement with independent measurements.This means that the modeling of radiative transfer in the ice layer is satisfactory, and that the thickness can be determined only using spectral measurements.However, this is not the case for the estimations of the grain diameter of the snow.Indeed, when the ice layer is thicker than 5 mm, our synthetic study predicts that it cannot be retrieved.Still, the results obtained on experimental data for slab thicknesses greater than 5 mm (blue and green curves in Figs.11 and  14) showed a posteriori PDFs for the grain diameters with surprisingly low standard deviations compared to what was obtained for synthetic data.The experimental results favor situations in which the geometrical optics hypothesis that is fundamental in the radiative transfer model is no longer valid.This shall not be interpreted as a result on the grain diameter, as the synthetic test showed that it was unaccessible.These low a posteriori uncertainties shall rather be interpreted as a compensation effect: a behavior that cannot be reproduced by the model may be approached by the most extreme values tested.In our case, small grain diameter results, even if they are either not realistic or not in agreement with the model's hypothesis, will produce an effect in the simulation that reproduces the data better than the other cases.
This work show that the radiative transfer model and the inversion method tested are adapted to retrieve the characteristics of an ice slab overlaying a granular layer.In particular, they are adapted to the study of the Martian CO 2 ice deposits, but also to the study of other planetary compact ice such as nitrogen ice on Pluto or Triton, or SO 2 ice on Io.However, our results show that this method is adapted to study the granular material underneath only in the most favorable cases, when the uncertainties of the data are lower than 1 %, or when the absorption in the slab layer is weak.In the general case of a slab ice layer covering a granular material, the retrieval method used in this work is not adapted to the study of the bottom layer.

Conclusions
The aim of this present work is to validate an approximate radiative transfer model developed in Andrieu et al. (2015d) using several assumptions.The most debated one is that the radiation becomes Lambertian when it reaches the substrate.We first qualitatively validated this assumption with snow and ice data.We then quantitatively tested and validated our method using a pure slab ice with various thicknesses and snow as a bottom condition.The thicknesses retrieved by the inversion are compatible with the measurements for every geometry, demonstrating the robustness of this method to re-  trieve the slab thickness from spectroscopy only.The result given by the inversion of the whole dataset is also compatible with the measurements.We also validate the angular response of such slabs in the specular lobe.Unfortunately, it was not possible to measure the microtopography in detail to compare it with the retrieved data.Nevertheless, we found very good agreement between the simulation and the data.In future work, an experimental validation of the specular lobe and roughness should be addressed.
The large uncertainties in the grain diameter inversion demonstrate that the bottom condition is less important than the slab for the radiation field at first order, as predicted by the synthetic tests conducted.The inconsistency between the a posteriori PDFs of the grain diameters for experimental data and numerical tests stresses that synthetic tests must be F. Andrieu et al.: Retrieving the characteristics of slab ice covering snow by remote sensing performed in order to determine which quantities can be retrieved or not in the context of the study, and to precisely calculate the expected uncertainties.
The comparison of the a posteriori uncertainties on the thickness of the slab and of the grain diameter of the snow substrate illustrates the fact that those uncertainties depend both on the constraint brought by the model itself and on the uncertainty introduced into the measurement, which only the Bayesian approach can handle.The use of Bayesian formalism is thus very powerful in comparison with traditional minimization techniques.We propose here a fast and innovative method focusing on massive inversions, and we demonstrated that it is adapted to remote sensing spectro-imaging data analysis.The radiative transfer model used in this study was proven appropriate to study the superior slab layer, but not the bottom one, unless the top layer is thin (thinner than 5 mm in our case).The whole method is thus adapted to study the top slab layer of a planetary surface using satellite hyperspectral data, for instance Martian seasonal deposits, that are constituted of a slab CO 2 ice layer resting directly on the regolith.

Data availability
The experimental spectra used in this work are publicly available on the SSHADE (previoulsy GhoSST) database.The following links are permanent.

Figure 1 .
Figure 1.Scheme of the surface representation in the radiative transfer model applied to the laboratory measurements.h represents the slab thickness and θ represents the mean slope parameter used to describe the surface roughness.

Figure 3 .
Figure 3. (a) Reflectance factor at a wavelength of λ = 1.4 µm vs. phase angle for snow only (black crosses) and the same snow but covered with a 1.42 ± 0.27 mm water ice slab (red squares).(b) Same data but normalized by the value at a phase angle α = 20 • .

Figure 4 .
Figure 4. Normalized stacks of 1000 a posteriori PDFs for the grain diameter of the snow, when conducting the inversion on synthetic data, with added random noise.The legends indicate the values for the grain diameter used to create the synthetic data.(a) The ice layer is 1 mm thick.(b) The ice layer is 5 mm thick.

Figure 5 .
Figure 5. (a) Normalized stacks of 1000 a posteriori PDFs for the thickness of the slab ice layer, when conducting the inversion on synthetic data, with added random noise.The legends indicate the values for the thickness used to create the synthetic data.(b) The a posteriori uncertainty (at 2σ ) on the thickness estimation as a function of the slab thickness.

Figure 6 .
Figure 6.(a) A posteriori uncertainties at 2σ on the grain diameter as a function of the noise standard deviation, for a 2 mm thick ice layer.(b) A posteriori uncertainties at 2σ on the thickness as a function of the noise standard deviation.

Figure 7 .
Figure 7. (a) Measured (black) and simulated (red) reflectance factor at 1.5 µm in the principal plane for an incidence angle of 50 • .The specular lobe measured is not centered at 50 • .(b) Measured (black) and simulated (red) reflectance factor at 1.5 µm in the principal plane for an incidence angle of 50 • .We simulated a small misadjustment of the sample, resulting in a shift of the observation of 0.5 • in emergence and 0.2 • in azimuth.

F
Figure 8. Measured and simulated reflectance factor around the specular geometry at 1.5 µm for an incidence angle of 50 • .The simulation was computed assuming the determined shift of 0.5 • in emergence and 0.2 • in azimuth.

Figure 9 .
Figure 9.The a posteriori probability density function for the roughness parameter θ , noted P θ .The inverted value at 2σ is θ = 0.424 ± 0.046 • .The simulation corresponding with the highest likelihood is obtained for θ = 0.43 • .

Figure 10 .
Figure 10.Reflectance factor spectra for the measure and the simulation at the maximum likelihood, and for the geometry for which maximum likelihood was highest, for each sample: at incidence 40 • , emergence 10 • and azimuth 140 • for sample 1 (thickness: 1.42 mm) (a); at incidence 40 • , emergence 20 • and azimuth 45 • for sample 2 (thickness: 7.45 mm) (b); and at incidence 60 • and emergence 0 • for sample 3 (thickness: 12.51 mm) (c).The thicknesses indicated were measured before putting the sample into the spectro-radiogoniometer, and the errors are given at 2σ .The absolute differences are shown in blue on each graph.

Figure 11 .
Figure 11.Marginal a posteriori probability density functions for (a) the thickness of the slab P {p 1 (i)} and (b) the grain diameter of the snow substrate P {p 2 (j )} for the three samples, and for the geometries described in Fig. 10.

Figure 12 .
Figure 12. Results of the inversions and measurements with error bars at 2σ for samples 1, 2 and 3, and for the 39 different geometries of measurement.The inversion points (in red) are sorted by incidence (three values), and each incidence is then sorted by azimuth (13 values: one for emergence 0 • and six each for the 10 and 20 • emergences).

Figure 14 .
Figure 14.Marginal a posteriori probability density functions for (a) the thickness of the slab P {p 1 (i)} and (b) the grain diameter of the snow substrate P {p 2 (j )} for the three samples.