A COMPOSITE MODEL FOR REFLECTANCE AND POLARISATION OF LIGHT FROM GRANULATE MATERIALS

Many natural land surfaces, such as sand or snow, consist of densely packed grains, often covered by dust, water droplets, contaminated with other materials such as possible oil leaks, hoar frost, and can also be internally cracked, porous, and heterogeneous. Most scattering models ignore these complications, but here a more detailed approach is taken to test all these effects. The current model is composed of three techniques: 1) Monte Carlo-based electromagnetic volume integral equation technique for non-spherical wavelength scale dust particles, 2) Monte Carlo ray tracing for stochastic-shaped grains much larger than the wavelength, with optional point scattering from dust cover, internal inclusions, and liquid surface layer, in a layer of an optical depths of few units, and 3) adding-doubling to combine smaller layers into an arbitrary, thick and vertically inhomogeneous medium. The model allows the medium to be built in a modular way, and after initialisation, rather complicated layered structures can be computed quickly and flexibly. The computed results are compared against experimental measurements of snow and sand. The model agrees with measurements usually within the measurement accuracy (∼ 0.05). The scattering is observed to depend significantly on grain size, shape, orientation, composition, fine structures, dust, and some other properties that need to be defined. Both, measurement and modelling, require much deeper attention to these properties.


INTRODUCTION
Many Earth and planetary surfaces are covered by granulate material, such as snow, sand, dust, and regolith. For many purposes, one needs to model the scattering from such surfaces (Muinonen et al., 2019, Martikainen et al., 2018. Several things complicate the modelling: the particles are packed densely compared to the free path, they can be sintered together, there can be multi-scale horizontal and vertical structures, layers can be thick and scattering order high, larger particles may be covered or contaminated by fine dust, and there are also structures in the wavelength scale. Most current models are streamlined for fast and smooth operation and skip many of these complexities. It is not even clear how carefully each factor needs to be modelled.
The aim of this work is to present a physical model that could be used in numerical experiments and analysis of the polarised scattering from simple granulate layered media composed of homogeneous irregular grains with optional dust cover. Specifically we focus here on samples of which we have lot of experimental data and can thus use the model to analyse the data and data to validate models. For this purpose the main model output is the bidirectional reflectance factor [matrix] and albedo as a function of measurable target properties.
Below, we first review existing techniques and experiments. Then we present the model. Next, we compare against experimental data from FGI's Reflectance Library. Finally, we make conclusions and recommendations for future work.
Exact electromagnetic scattering techniques can only be used in limited cases, such as spheres, slabs, and small, simple objects. * Corresponding author In general, frequency-domain methods, at least the interior or the surface of the scatterer, must be discretised in a sufficiently dense grid. This typically leads to very big matrices. A rectangular grid, as used in the discrete dipole approximation (DDA), needs to be of high density to reproduce the shapes well (Yurkin, Hoekstra, 2007, Draine, Flatau, 1994, Lumme, Rahola, 1994. Recently, methods have generalised to much more complex scatterers (Muinonen et al., 2019, Muinonen et al., 2017, Markkanen et al., 2015. The scattport lists over 100 scattering codes already (Hellmers, Wriedt, 2009, SCATTPORT, 2020 Ray-tracing assumes the radiation can be discretised using a finite number of rays crossing through the medium using different paths. The rays propagate directly in space, reflect or refract from flat surfaces, and scatter to selected or random directions from diffuse scatterers. This requires that all effective length scales of the medium be much larger than the wavelength. Smaller-scale effects can be included using diffuse volume or surface scattering. Ray-tracing has been applied to solve scattering from many kinds of particles and media , Peltoniemi, Lumme, 1992. This technique allows 3-dimensional structures to be modelled well. Ray-tracing can be extended towards physical optics and coherence effects (Muinonen, 1989, Muinonen, 2004, Hesse et al., 2009, Muinonen et al., 2012. Unfortunately, ray-tracing is a slow method, especially for thick multiply-scattering layers.
Radiative transfer approximation assumes independent volume scattering. This allows fast computations, especially in horizontal symmetry. However, 3-D structures, dense packing, clustering and interference are usually excluded or oversimplified. For unpolarised plane parallel media, the discrete or-dinate method is most used (Stamnes et al., 1988). Another interesting approach is the approximate asymptotic radiative transfer technique by (Kokhanovsky et al., 2018, Kokhanovsky, Zege, 2004. Polarisation can be conveniently included using the adding-doubling technique (de Haan et al., 1987, Stammes et al., 1989. The limitations of radiative transfer assumptions can be somewhat relaxed by initialising the technique using inputs from less approximate techniques. High-quality experimental data to validate models are still rare. Proper validation needs both maximally simplified experiments to study different parts of the models, and realistic experiments to test their true applicability. In both cases, all model input variables should be known with a high degree of accuracy. Natural particles are very hard to characterise, even statistically. Artificial particles provide better control, e.g., (Ulanowski et al., 2003), though they of course lack many realistic features. Scattering from a stream of particles has been measured, e.g., by (Sasse, Peltoniemi, 1995, Muñoz et al., 2010, and by individually levitating particles by (Maconi et al., 2018, Muinonen et al., 2019, Maconi et al., 2020. Back-scattering from snow and regolith has been measured by (Piironen et al., 2000, Kaasalainen et al., 2002, Näränen et al., 2004, Kaasalainen et al., 2005, Kaasalainen et al., 2006. FGI has a large database of different measurements of bidirectional reflectance factors available that can be used for model validation: , Wilkman et al., 2016, Zubko et al., 2019, Peltoniemi et al., 2020, Xu et al., 2009. Another well-managed library is SPECCHIO by (Hueni et al., 2009), and in planetary field SSHADE (Schmitt et al., 2018).

MODELS
We assume the particles can be divided to two classes: larger particles that are visually detectable and where geometric optics applies -here called grains, and smaller particleshere called dust -that are visible only in quantities and need electromagnetic modelling. We assume the grains are packed densely and the optional dust covers the grain surfaces. The shapes of both particles are modelled as randomly rough ellipsoids obeying log-normal distribution. We assume that the materials are statistically plane-parallel, with minor 3D-effects, and separable to semi-homogeneous layers. There can be small topological roughness on the top layer. We model sintering by allowing grains to overlap. The grains can be either fully randomly oriented, or flat side up. The grains may also have internal scatterers such as airbubbles or contaminants.
The model is composed of three parts: ray-tracer for thin densely packed media (Peltoniemi, Lumme, 1992, Peltoniemi, 2007, adding-doubling radiative transfer programme (Peltoniemi, 1993) for thick layered media, and electromagnetic volume integral equation scattering solver for small dust particles (Peltoniemi, 1996). All components have been upgraded significantly, thus more detailed description follows here.

Definitions
The bidirectional reflectance factor (BRF for short, or R in equations) is defined as the ratio of the reflected light intensity of a given target to an ideal Lambertian reflector with a spherical albedo of 1.0 under the same incident irradiation (Nicodemus Figure 1. Definition of the angles used in surface reflectance work: and ι are the zenith angles of the emergent (Observer) and incident (solar) radiation, respectively (shorthand µ = cos and µ0 = cos ι are also used). φ and φ0 are the corresponding azimuth angles. The phase angle α is the angle between the observer and the Sun. The principal plane is fixed by the solar direction and the surface normal, while the cross plane is a vertical plane perpendicular to the principal plane. et al., 1977et al., , Hapke, 1993. The BRF can be presented as where F0 is the incident collimated flux and I the reflected radiance; ι and φ0 are the zenith and azimuth angles of incidence, and φ are the zenith and azimuth angles of emergence, µ = cos , µ0 = cos ι, α is the phase angle that is defined as the angle between the source and observer. (Fig. 1).
When the polarisation of radiation matters, I must be a Stokes vector [I, Q, U, V ] and R a 4 × 4 Muller matrix.
The degrees of linear polarisation are defined here as

Ray-tracing
The ray-tracing part is based on (Peltoniemi, Lumme, 1992, Peltoniemi, 2007, with many fixes and improvements. This assumes densely packed irregular particles with Fresnel reflections and refractions from the surface. The particles are characterised by their mean semiaxii (a1, a2, a3), radius standard deviation σ, and slope standard deviation ρ. The particles can be randomly oriented, either fully isotropically, or azimuthally isotropically with horizontal/vertical asymmetry. Additionally, particles can contain air bubbles inside, or small dust on the surface. These can be modelled either as a similar Gaussian geometric object as the larger particles, or point scatterers by a given scattering matrix, computed here either using BHMie (Bohren, Huffman, 1983) or below described Monte-Carlo electromagnetic technique, but any other discretised phase matrix can be used.
The top surface of the layer can contain the small roughness of a length scale of a few particles in diameter. When the medium is dark or thin, all scattering can be computed using this ray-tracer, but when the medium is bright and thick, the order of scattering can grow very large, making it particularly slow for ray-tracing. In that case, ray-tracing is only used to compute small layers that are combined using adding-doubling (see Subsection 2.4). For this purpose, the reflected and transmitted radiation is recorded in a polar grid used in the adding-doubling process, i.e., Gaussian quadrature points for the zenith angles, and regular intervals for the azimuths, with a few extra points in the principal plane and around selectable points of interest.
Single and multiple scattering are recorded separately, to allow back scattering corrections in the adding-doubling process and improved accuracy for single scattering. Typically, simulation for one medium configuration takes about one hour or more, depending on the angular resolution and desired accuracy.

Monte Carlo volume integral equation technique
The small wavelength scale dust can be modelled using many electromagnetic techniques. Here, a substantially revised version of (Peltoniemi, 1996) is used, simultaneously testing some other ways to do things.
The electromagnetic wave equation in the presence of a scatterer of volume V and refractive index m, with an incidentharmonic wave E0 of wavelength 2π, including boundary conditions and the scattering condition at infinity, can be presented in the form of an integral equation where χ(r) = m(r) 2 − 1 and G is the dyadic Green function To solve the equation numerically, the electromagnetic field E inside the scatterer is discretised using base functions fi Mathematically most rigorous way to solve this integral equation is least squaring (L2-norm), and in all our tests it yielded best convergence and accuracy against Mie results, but it demands time and memory with larger particles. Thus, here, the integral equation is solved using the method of weighted residuals. Multiply Eq. 6 by a set of test functions hj(r) and integrate over the volume of the scatterer to get a set of linear equations: where Q(ρ) = −1/3+O(ρ 2 ), and O contains some second and higher order corrections for a finite size volume Vρ of radius ρ around the singularity. Here, ρ is set so small that all O terms are insignificant. A and B are the matrix representations, and the separation of A into two parts is probably evident.
After testing with entire domain base functions (best for small and easy particles), pulse functions (as used in DDA), and linear base functions, we selected 64-point Hermitian tricubic functions in a rectangular grid as base functions. They provide a good compromise between matrix sizes, sparsity, and continuity up to the first derivative. The actual particle shape is modelled independently of the field grid. The grid is extended outside the particle as far as the interpolation rule requires, and the external values outside the particle are extrapolated from the nearest internal points linearly. Thus, boundaries behave smoothly, and no restriction for the particle shape is needed.
In the traditional Galerkin method, the complements of base functions are used as test functions, yielding some economy in evaluation, but we try different options again, and instead use plane waves hj(r) = e iK j ·r , because they evaluate quickly, and are continuous, rather complete, and allow semi-analytic treatment of the singularity: where r = r − r = RΩ, R = R(Ω) is the distance to surface (in convex particles, this separates occasionally to parts), κ = (1 + Kj · Ω), R0 is an arbitrary constant, for example, the radius of the largest internal sphere and the r-integral is evaluated using Taylor series at small values and a tabulated numerical integral at larger values. Additionally, the κ → 0 case can be evaluated using Taylor series. The oscillating angular integral is evaluated numerically using a randomly oriented Lebedev grid.
The other volume integral in Eq. 9 -note the change of integration order -is evaluated using a Monte Carlo grid providing full flexibility for arbitrary shapes and structures. The number of integration points NMC must be rather right, typically ∼ 100× the number of expansion (NX ) and test functions (NT ). The number of test functions must be at least the number of expansion/unknowns, but small overkill improves accuracy, if orthogonality or independence is compromised by irregular shapes and mathematically violent procedure. The method works adequately with a base function grid density as low as 4 points per wavelength (2π|m|), which is much less than with DDA, and reduces memory limitation to the point of being mostly insignificant, even on a normal laptop.
Orientation average is computed by vectorising the incident field (=right hand side of equation). When orientation average is not needed, the base function grid can be asymmetrically optimised, providing a speed increase of at least a 50% speed increase.
We have tested this against spheres computed using the BHMie programme by (Bohren, Huffman, 1983). For the size parameter x in the range 0.0001 to 10/|m|, the extinction and scattering cross-sections can be computed with an accuracy of 0.01 to 0.0001 (in most cases, similar or much better than ADDA). With larger refractive indices, the absorption coefficient and internal fields remain more difficult, in 1 to 50% inaccuracy, due to spurious fluctuations in need of a lot of integration points to become smooth. If only scattering is interesting, orientationaveraged scattering matrices can be computed for close to x = 20 particles with moderate refractive index at 24 hours in a single laptop processor.
Non spherical shapes are modelled using a spherical grid of radii and interpolating between. Random shapes are generated using discrete log-normal statistics for the grid points.

Adding-doubling
The adding-doubling part originates from (Peltoniemi, 1993). The method is selected because of its flexibility for variegating layering and full polarisation.
The medium (snow pack) is divided into quasi-homogeneous layers. Each layer is initialised by a small sub-layer from the ray-tracing part or single scattering. These sub-layers are then doubled to full layer thickness, and layers are added to full size medium. Different initialisation layers can be averaged or interpolated together to simulate different materials: R eff = Riwi etc.. Similarly, at least first-guess-accurate grain size distribution can be forked by scaling initialisation layers up and combining with doubled layers at several steps. Thus, a timeconsuming ray-tracing part is needed only to initiate the computations, and arbitrarily, many new combinations can be built during the adding/doubling phase.
To optimise the process, single and multiple scattering are treated separately. Here, single and multiple scattering are considered just methodological, not physical, and single scattering means the part that is expressed in an analytic single-scattering formula, and everything else is included for multiple scattering. (R = R1 +RM , R1 = (1−exp(−h/µ+h/µ0)P(α)Γ(), Phase matrix P expanded as Legendre series). Thus, multiple scattering parts can be computed using a smaller grid, and single scattering terms at higher accuracy. The single scattering part is corrected for dense packing shadow effects Γ. The original work (Peltoniemi, 1993, Jämsä et al., 1993 also included many surface roughness and internal heterogeneity effects, but here they are disabled because they need more testing to fit into the current package. Further, we now use x = µ 2 -discretisation (dx = 2µdµ), making plots look nicer, and fitting into experimental data of limited zenith range that is less divergent.
The adding part takes some seconds or minutes, depending on the level of heterogeneity and packing corrections.

COMPARISONS WITH EXPERIMENTAL DATA
FGI's Reflectance Library contains over 500 measurements of different samples. We have selected a few different cases for comparison. The data are measured using the Finnish Geodetic Institute's Field Goniospectropolarimeter FIGIFIGO. The estimated uncertainty of BRF is about 2-5% . The polarisation has never been calibrated, and the inaccuracy is not known. A guess is at least 0.02 absolutely, and 5% relatively, but there is the possibility of even larger systematic or casual errors.
The comparisons are made in the principal plane, using a visual part of the spectrum, integrated over 550 to 750 nm, which is the least noisy region.

Snow
To get a visual overview of the effects of various model configurations, and a first guess for inversion, a set of different snow configurations is shown in Fig.2. The staring point is "plane snow", which is 20% dense, and grains that are randomly oriented ellipsoids (semiaxii 0.4, 0.5, 0.69 mm) with a moderately rough surface (σ = 0.01, ρ = 0.1). "Needle snow" is more elongated by shape (0.2,0.2, 1.0). "Dense snow" is 40% dense. "Fluffy snow" contains small ripples of hoar on the surface, modelled electromagnetically, a size below a few microns. "Bubbly snow" contains small geometrically modelled air bubbles inside. In "Flat snow", grains are oriented on the flat side. All layers are 64 cm deep. It is obvious that all these variations change the scattering properties. Grain shape affects the forward/backward ratio significantly, dense packing increases scattering. Flatter orientation especially strengthens middle scattering angles. Hoar and bubbles increase polarisation in the middle and backward directions. Next, the model is compared against three measured cases: new snow, old refrozen snow, and volcanic sand-doped snow, Fig.  3. A most passing case was first selected from Fig. 2, based on measured or estimated snow properties, and fine tuned for a nicer fit.
New snow was measured in the Tähtelä research station area, in Sodankylä near the Norsen mast, at a small opening in a pine forest. The snow is loose and starting to metamorphose from flakes to needles. The "Needle snow" from Fig. 2 seems to best fit the measurements and estimated snow properties. For the model, semi axii 0.2 mm, 0.2 mm, 1.0 mm are used, with small roughness (σ = 0.01, ρ = 0.1), and density of 0.1. The model agrees to data well within the spectrum of measurement accuracy, although BRF could be better.
Old, slightly melted and refrozen snow was measured at the Sodankylä airport , Svensson et al., 2016. Snow thickness was about 70 cm and the surface was rather flat. Grain mean diameter was around mm. Grains were rounded and slightly elongated, but no detailed shape data are available. For modelling, the main axes 0.4-0.5-0.625 mm were assumed, oriented flat side up, and a density of 0.20 (∼ Flat snow in Fig. 2). The model slightly overestimates the forward spike, maybe due to some details of grain shape, orientation or surface structure. Polarisation agrees within the measurement uncertainty.
From the same campaign, dirty snow was found with volcanic sand deposited on a snow surface. The dark stuff caused some melting and compressing to the snow. Thus, snow density was set at 0.25, but otherwise, the same parameters were used. Additionally, dark particles of a refractive index of 1.77+0.001i, grain diameter of 0.5 mm, roughness σ = 0.01, ρ = 1.0, and density 0.02 were assumed. Agreement is still strong.

Sand
We compare against a beach sand sample measured in laboratory 7. May 2010, using unpolarised QTH light source (Fig. 4). The surface was made as flat as manually possible. Thickness was several cm, visually sufficient to cover the bottom. Grain size was below 1 mm, and the sand was fully dry. In the model, the real part of the refractive index was set to mR = 1.55, with mean grain radius 0.2 mm, roughness σ = 0.01, ρ = 0.5. The imaginary part was varied in a range from 10 −9 to 0.01 to find the best match, which happened around 3 · 10 −6 . However, the imaginary part scales with the size, and cannot be fully uniquely fixed by the data.
As seen in Fig 4, the reflectance and polarisation are adequately reproduced by the model in the almost full zenith angle range. Some discrepancies are observed near the backward and forward scattering directions where the measurement errors are typically larger due to experimental constrains.

CONCLUSIONS
We have developed a multi-step model for spectro-polarised reflectance from a particulate medium, and compared that against measurement data.
All model components are greatly improved from previous versions and provide some interesting features. The volume integral code proves that one can also compute the difficult integrals using the Monte Carlo technique, mixing different base and test functions, and get excellent accuracy. Even though losing to ADDA by computing time hugely, the memory usage is only a fraction of the usage of ADDA, and the accuracy is consistent and converging with all realistic refractive indices. This may be useful for applications needing orientation and refractive index averaging or inversion. The ray-tracer allows good experimenting with irregular particles covered by stuff. Spectral effects and size distributions can be combined and scaled well. The adding-doubling part integrates the ray-tracing and single scattering results to a fast flexible package with more advanced corrections than plain adding-doubling.
The model can predict the bidirectional reflectance factor close to the the range of measurement uncertainty (∼ ±0.05) in nearly full spectral and angular range. The samples themselves can vary ±20% in many properties in a few minutes, making fitting complicated. There are plenty of parameters to adjust further for fine tuning. Also, the polarisation is mostly well within the measurement uncertainty.
The next steps in modelling are to macroscopically re-implement rough surfaces and internal heterogeneities as in (Jämsä et al., 1993, Peltoniemi, 1993, and play more with grain shapes, orientation, and dust cover. The cyclic coherence effects must be included to effectively predict the near back-scattering phenomena (Muinonen et al., 2012). Further, the forward and specular coherences of the surface dust must be considered, because they really matter inside and between particles. This model frame provides a good base for development and numerical experiments, but it is too slow and complicated for production work.
The measurement accuracy, especially for the polarisation, must be improved to a 0.1% level. Further, the old challenges still remain: how to find targets that are sufficiently simple and controllable that models can be tested, and how to get a realistic sampling of the complex nature. Even the simplest targets are too complex to be characterisable with sufficient rigour. All the measured natural targets are strongly biased towards the easiest ones that one can find, and the real mess nearby is not counted.