Uncertainty Introduced by Darkening Agents in the Lunar Regolith: An Unmixing Perspective

: On the Moon, in the near infrared wavelength range, spectral diagnostic features such as the 1- µ m and 2- µ m absorption bands can be used to estimate abundances of the constituent minerals. However, there are several factors that can darken the overall spectrum and dampen the absorption bands. Namely, (1) space weathering, (2) grain size, (3) porosity, and (4) mineral darkening agents such as ilmenite have similar effects on the measured spectrum. This makes spectral unmixing on the Moon a particularly challenging task. Here, we try to model the inﬂuence of space weathering and mineral darkening agents and infer the uncertainties introduced by these factors using a Markov Chain Monte Carlo method. Laboratory and synthetic mixtures can successfully be characterized by this approach. We ﬁnd that the abundance of ilmenite, plagioclase, clino-pyroxenes and olivine cannot be inferred accurately without additional knowledge for very mature spectra. The Bayesian approach to spectral unmixing enables us to include prior knowledge in the problem without imposing hard constraints. Other data sources, such as gamma-ray spectroscopy, can contribute valuable information about the elemental abundances. We here ﬁnd that setting a prior on TiO 2 and Al 2 O 3 can mitigate many of the uncertainties, but large uncertainties still remain for dark mature lunar spectra. This illustrates that spectral unmixing on the Moon is an ill posed problem and that probabilistic methods are important tools that provide information about the uncertainties, that, in turn, help to interpret the results and their reliability.


Introduction
Reflectance spectroscopy can provide valuable insights into the composition of planetary surfaces. On the lunar surface, the mineral composition is dominated by plagioclase and pyroxenes [1], which show diagnostic absorption bands in the near-infrared wavelength range at 1-µm and 2-µm (e.g., [2]). These spectral features enable us to estimate the abundance of certain minerals or elements. Correlations between spectral parameters and spectra with known composition can be used to create large-scale maps of, for example, compounds of chemical elements, such as FeO or TiO 2 (e.g., [3][4][5]). To determine the abundances of known materials, spectral unmixing can be performed (e.g., [6,7]). With higher spatial and spectral resolution data available in the last decade, such as data from the Moon Mineralogy Mapper (M 3 ) instrument on-board the Chandrayaan-1 spacecraft [8], it is possible to investigate the composition of the Moon in more detail and to gain insight into the evolution and geology of the Moon.
Unmixing describes the process of determining the abundances of a set of constituent materials (endmembers) for a given spectrum employing radiative transfer modeling to describe the interaction of light with the different endmembers. Here, it must be differentiated between a surface where the endmembers build a macroscopic mixture, such that they are spatially separated as in a checkerboard pattern. This type of mixture can be described by a linear superposition of endmember reflectance spectra. The Moon, however, is covered by a layer of grains from crushed rocks, the so-called regolith. This constitutes a typical intimate mixture, where light interacts with the different particles by multiple scattering inside the medium. Intimate mixtures cannot be described by a linear mixture of reflectance spectra in the near infrared. It has, however, been shown (e.g., [6]) that the single scattering albedo (SSA) spectra can be unmixed linearly. Models such as the Hapke model [9] can be used to calculate the SSA spectra, thus transforming a non-linear problem into a linear unmixing problem.
The vast majority of unmixing approaches rely on classical optimization techniques (e.g., gradient descent or Levenberg Marquardt algorithm) to fit the parameters of the model, such as abundances or grain size, to match the measured spectrum. This has been done for several planetary bodies, such as on Mars (e.g., [10,11]) or on the Moon (e.g., [2,12,13]). The endmembers can for example be extracted from the image itself such that they span a simplex. Each spectrum in the image can then be calculated as a linear superposition of the edges of the simplex (endmembers) [6]. While this has the advantage of being applied easily, it has the disadvantage that the composition of these so called image endmembers might not be known. A different approach is to select endmember spectra from a library of laboratory samples of known composition [6]. This way the composition can be directly determined, but the endmembers must be selected carefully and some spectra might still lie outside of the resulting simplex.
Spectra obtained under lunar conditions, however, are systematically different compared to laboratory spectra of earth-analogs measured in the laboratory. Due to the lack of an atmosphere, the lunar regolith is subject to the continuous influence of the space environment. Micrometeorite bombardment and solar wind change the physical and chemical characteristics of the upper layers of the regolith [14,15]. These changes include (1) creation of submicroscopic iron particles (smFe 0 ) [14,15], (2) selective comminution of larger grains into smaller grains, (3) melting of grains into agglutinates [15] and (4) an increase in porosity creating a 'fairy-castle' structure [14]. The smFe 0 particles can be grouped into two categories. Firstly, the smaller nanophase iron (npFe 0 ) particles, which mostly redden the overall spectrum and dampen the absorption bands, but also darken. These particles are about 2-30 nm in diameter and accumulate in the rims of the mineral grains [15]. Secondly, the microphase iron (mpFe 0 ) particles which darken, but do not lead to an increase in spectral slope [16]. All these changes significantly influence the optical properties of the lunar surface and make it substantially more difficult to distinguish between minerals. For example, Sunshine and Pieters [17] and Mustard and Pieters [2] found that, due to low spectral resolution and signal-to-noise ratio multiple solutions might result in similar least squares fits, but are physically implausible.
Therefore, unmixing lunar spectra either rely on image endmembers or mature samples returned by the Apollo missions. A catalog that covers a variety of compositions of the maria is the Lunar Soil Characterization Consortium (LSCC) catalog [18,19]. However, because the simplex is limited to only the main mare compositions and a few highland samples, this catalog cannot cover the entire mineralogy of the lunar surface. For example, the highlands are dominated by plagioclase (90%) [1] and the LSCC catalog cannot accurately describe the typical highland spectra [20]. Further, uncommon compositions such as the spinel-rich regions detected by Pieters et al. [21] cannot be explained by unmixing purely based on the LSCC catalog.
Usually, when unmixing a spectrum, only a single best-fit solution is obtained by using one of the optimization techniques mentioned above. Rommel et al. [13] used a similarity measure between the measured and the modeled spectrum to select the best combination between all possible combinations of endmembers. The differences between best fit solution for a particular combination of endmembers and the measured spectra for many combinations are, however, quite similar. Tiny changes in the spectrum could result in an entirely different combination to be selected as the best solution. Lapotre et al. [22] showed with a variety of experiments that the uncertainties of the endmember abundances and grain size are relatively large and many solutions are equally likely. On the Moon, these uncertainties are likely to be even higher, with the influence of space-weathering. In Rommel et al. [13] the best-fit solution often contained the ilmenite endmember, in contrast to the absence of ilmenite in the actual mixture. Mineral darkening agents such as ilmenite, which is abundant in some lunar maria, add another layer of complexity to the unmixing as they are almost featureless and, therefore, mostly darkened.
In this work, we develop a Bayesian approach to unmixing. In Bayesian inference (e.g., [23]), all parameters of the model are considered probabilistic variables. In this framework, in addition to the point estimates of the best fit, we also estimate the uncertainties of the parameters. Due to the darkening agents, which darken the spectra and dampen the absorption bands, the uncertainties of the mineral abundances are likely to be higher on the Moon than on other planetary bodies with an atmosphere (e.g., [22]). Additional prior information can be conveniently incorporated in the Bayesian framework. To reduce uncertainties of the mineral abundances, estimates of elemental abundances can be included such that the algorithm will prefer solutions that use endmember combinations with similar elemental compositions.

Methods
Generally, the methods employed in this work can be subdivided into two main frameworks. Firstly, the radiative transfer model describes the scattering and absorption behavior of the materials. This framework is built upon the Hapke model [9] and the space weathering model of Wohlfarth et al. [24]. Secondly, the Bayesian framework formalizes the probabilistic perspective of the problem. In the following section, we will describe the choices of the models and their parameters.

Hapke Model
The semi-empirical Hapke model has become the quasi-standard for the analysis of particulate planetary surfaces. Over the years the original model [25] has been refined to also include opposition effects and better approximations for the phase function and Chandrasekhar's H-function to model multiple scattering inside of the medium. The Modified Isotropic Multiple-Scattering Approximation (MIMSA) formulation of the model [26] is as follows: The single scattering albedo ω (SSA) can be calculated by inverting the equation for ω. The cosines of the incidence (Θ i ) and emittance (Θ e ) angles are denoted as µ 0 and µ, respectively. The phase angle g is the angle between the incident and emitted light rays. The phase function is denoted as p(g), the opposition effect is subsumed in the terms for shadow hiding B SH and coherent backscattering B CB . We also include corrections for macroscopic surface roughness as described in Hapke [27], but for the laboratory mixtures used in this work, this effect is not relevant. The multiple scattering term M(µ 0 ,µ) models the scattering between the particles inside the medium and is described in detail in Hapke [26] and not repeated here.
In general, the phase function p(g) can be represented by a sum of Legendre polynomials P n (cos(g)) weighted with the material-specific Legendre coefficients b n : p(g) = 1 + ∞ ∑ n=1 b n P n (cos(g)). ( Specifically, the phase function p(g) of the soil can be approximated with the double Heney-Greenstein function [28]: with the parameter c DHG being a measure for the amplitude of the backscattering lobe to the forward scattering lobe. c DHG > 0 denotes particles that are mainly backscattering. The parameter b DHG is a measure of the width of the lobe. For the lunar regolith, Warell [29] found the standard values for the parameter b DHG and c DHG to be 0.21 and 0.7, respectively. According to Yang et al. [30], ilmenite is more back-scattering while the other lunar minerals are mainly forward-scattering. We adopt the parameters for ilmenite from Yang et al. [30] and set b DHG =0. 19 and c DHG =0. 33 for the ilmenite endmembers. For the other endmembers we use the parameters from Warell [29] to be consistent with the previous work of, for example, Rommel et al. [13]. By expanding the Double Heney Greenstein function the Legendre coefficients can be derived from the parameters b DHG and c DHG as follows [26]: The multiple scattering term M(µ, µ 0 ) in Hapke [26] is derived from the Legendre coefficients of the single scattering phase function p(g). Using only the first 15 coefficients is sufficient, as the b n decreases with increasing n.
Under small phase angles, the surface appears brighter than would be expected. This is accounted for in the model by the correction terms B SH and B CB , denoting the shadow hiding and coherent backscatter effects, respectively [26]. The samples used in this work were acquired under a phase angle of 30°, therefore, we will assume that the opposition effect is negligible, if not mentioned otherwise.
This model is inverted to find the best fit albedo for the measured reflectance values. This procedure is applied to the endmember spectra that are then used according to the forward mixing model.

Forward Mixing Model
The single scattering albedos of the endmembers obtained by inverting the Hapke MIMSA model are subsequently used to calculate the SSA spectra of the mixture. According to Hapke [9] the single scattering albedo is defined as the scattering coefficients (S) divided by the extinction coefficient (E) of the mixture, with N i as the number of particles of particle type i per unit volume and σ i = πa 2 i as the cross-sectional area, where a i denotes the radius of the corresponding particle type [9]. The scattering efficiency is Q S,i = Q E,i ω i , such that the mixing formula becomes: The bulk density of particle type i is denoted M i and the average particle diameter is D i . Assuming that the endmembers have the same average grain size D i and the same extinction efficiencies Q E,i , the fraction can be reduced. The albedo of the mixture can then be expressed as a linear superposition of the endmember SSA spectra ω i weighted with the Hapke coefficients θ em,i , For a mixture of typical lunar minerals, the extinction efficiencies, grain size, and density are constant and approximately equal among the endmembers; therefore, these simplifications will be used to calculate the SSA spectrum of the mixed mineral endmembers. This spectrum is then subsequently used as the input to the space weathering framework to calculate the reflectance spectra of the mixture.
When the mixture directly results from the particular known endmembers used in the unmixing, the sum of the weights should be equal to one (∑ i θ em,i = 1). In this work we are not enforcing this constraint, but prefer solutions where the sum is close to one. By not enforcing the sum-to-one constraint for the fractional abundances θ em,i the sum can be seen as a "catch-all" parameter for all linear, wavelength independent, differences among the endmembers and compared to the examined sample. It is, however, physically implausible to have endmembers with a negative contribution. Therefore, all endmember weights must be larger than zero, that is, θ em,i > 0.
Apart from the mineral composition and maturity of the sample, other physical properties and the laboratory setups themselves also influence the measured spectra. In this work, we use samples and endmembers of the same grain size fraction to remove the influence of the grain size that has a strong effect on the overall spectra (e.g., [31]). The grain size also introduces large uncertainties for immature spectra [22] and would, therefore, influence the dependencies of the parameters we want to examine in this work.
The porosity or compactness [32] of a sample changes the brightness of a sample almost without introducing a change to the continuum slope. Therefore, this parameter introduces a scaling factor on the weights of the endmembers depending on the compactness of the sample.
Depending on the laboratory setup, other factors might influence the measured spectra. For example, the gain of the sensor might be different such that a linear offset is introduced, when comparing the measurements of the same mineral for different laboratories.

Space Weathering
Space weathering poses a major challenge for unmixing mineral abundances. In addition to the darkening and reddening, the diagnostic absorption features are obscured by the optical influence of the space weathering particles. This makes it harder to differentiate between the different minerals. While the band position remains constant the depth of the absorption that could normally be used to estimate the fraction of the constituent minerals is less clear.
Wohlfarth et al. [24] defined a framework based on ab-initio Mie modeling to characterize the effects of the space weathering particles. We include these characteristics in the unmixing framework. It is important to note that extinction efficiencies are not constant over the wavelength and are very different for the mineral endmembers and for each of the iron particles (see Figure 1a). As a result, the fraction in Equation (7) cannot be reduced. The SSA spectrum of the mixed mineral endmembers (ω med ) according to Equation (7) is then used as the starting value for the mixture with the iron particles. It is here assumed that all endmembers have a similar grain size distribution, density and extinction efficiencies. The values are adapted from Wohlfarth et al. [24] and listed in Table 1. The mixture of the iron particles and the medium is then defined as: The phase function of a mixture is also defined by the weighted sum of the phase functions of the components. The phase functions of the soil and the small npFe 0 particles are independent of wavelength, but the larger mpFe 0 particles have a wavelength dependent phase function (see also Figure 2), which is, therefore, given as a vector.    The phase functions are defined by the Legendre coefficients as in Equation (2). These b n are calculated based on the refractive indices, grain size, and the Legendre expansion of the Mie phase function from Fowler [33]. By equating the coefficients it can be shown that the b n of the mixture are also the weighted b n of the soil, npFe 0 , and mpFe 0 particles. The Legendre coefficients of the mixture are then used to calculate the wavelength dependent phase function of the mixture and are additionally used in the multiple scattering term M(µ, µ 0 ). Finally, we obtain reflectance spectra of the mixture, according to Equation (1).

Endmember Catalogs
In order to generate spectra similar to lunar conditions, we need to select endmembers representative of the general mineralogy of the Moon. We are using only mineral endmembers and are not trying to model the influence of agglutinate glasses. This glassy fraction usually has a similar composition to the mineral fraction [34], as it is created endogenously.
The absorption band depth is also not affected by glasses, but by a large amount of smFe 0 in the agglutinates [35]. Agglutinates here are only relevant for the unmixing of the LSCC catalog in Section 3.3, but this should be kept in mind when evaluating the results.
The highlands are dominated by plagioclase (around 90 wt.%) with the remaining fraction mainly consisting of pyroxenes and minor phases of olivine and spinels [1]. The composition of the maria is more diverse. About 30-40 wt.% are pyroxenes and plagioclase still forms a major part of the composition with up to 60 wt.% of the mineral composition. Olivine contributes up to 10 wt.% and the opaque mineral ilmenite is abundant in some mare areas with up to 20 wt.% at maximum [1].
Among these groups of minerals are also variations depending on the elemental composition. For example, Sun and Lucey [36] showed that by generating albedo spectra based on the magnesium number and grain size the abundance estimation of olivine and pyroxenes can be improved. We will use fixed albedo spectra as endmembers that are chosen to be a good representation of a wide variety of possible true endmembers while retaining the most distinct endmembers.
In general, we want to have one representative of each of the major minerals on the Moon. Pyroxenes can be subdivided into low calcium pyroxenes (LCP), mainly orthopyroxenes (OPX), and high calcium pyroxenes (HCP), usually clinopyroxenes (CPX). The LCP pyroxenes exhibit an absorption 0.9-µm and 1.9-µm [37,38]. The HCP pyroxenes can be further subdivided into two spectrally distinct groups. Type A, with an absorption at 0.9-µm and at 1.15-µm, but no absorption around 2-µm [37] similar to olivine. Type B, exhibits absorption bands at 1.05-µm and 2.35-µm [37]. The absorption band positions, however, are not fixed, but vary with Mg 2+ , Fe 2+ , and Ca 2+ content [39,40]. For the unmixing, the distinction between these three types of pyroxenes (HCP, LCP-A, and LCP-B) is, however, beneficial, because we want to minimize the linear dependencies between the endmember spectra.

Bayesian Inference
From experience of using unmixing algorithms based on an exhaustive search of all possible endmember combinations (e.g., [13]), we know that several possible combinations can produce similar results. The solution is, therefore, not unique. There might be multiple local minima and some endmembers that are similar, when they become mature, are even gradually interchangeable, without significantly reducing the quality of the fit. Classical unmixing techniques, that are based on finding the best fit in the least squares sense, ignore these uncertainties and only take a single solution, which by a somewhat arbitrary similarity measure best describes the measured spectrum. By exploring the parameter space more thoroughly we can obtain uncertainties and can also learn more about the dependencies among the model parameters.
The process of estimating the probability distributions of parameters of a model based on a set of data and additional assumptions or knowledge about the data is called Bayesian inference (see e.g., [23]). All parameters of the model are described as probability distributions. The goal, therefore, is to estimate the probability densities of the parameters Θ given some data D. This is called the posterior distribution p(Θ | D). According to Bayes' rule the posterior density can be described as: The data are not changing for one observation, therefore, we can assume that p(D) is constant. Consequently, we can define the unnormalized posterior density [23] as: The probability of the parameters p(Θ) is the prior distribution and can introduce already known information about the distribution of the parameters into the procedure or it can be uninformative. The probability of the data given a set of parameters Θ is denoted as the likelihood (p(D | Θ)). The likelihood is evaluated as a function of Θ and gives us an estimate on how likely the data have originated from the current parameters.
This unnormalized posterior distribution can be calculated analytically if the probability distributions are chosen in a specific way, which might be limiting in practice. With the advancement of computer technology and the introduction of Markov Chain Monte Carlo (MCMC) methods [41], such as the Metropolis-algorithm [42], it is possible to estimate complex probability distributions numerically. The algorithm used in this work to sample from the posterior density distribution is described in the following section.

Sampling Algorithm
The goal of the Monte Carlo simulation is to approximate the posterior distribution. Most of the sampling algorithms as well as our approach are based on the Metropolis algorithm [42], which was subsequently generalized to non-symmetrical proposal distributions and is then called Metropolis-Hastings algorithm [43]. This algorithm describes a random walk through the posterior space such that for a large number of samples, the samples converge to the posterior distribution. Generally, at each step of the Markov Chain a proposal (Θ t ∼ π(Θ t−1 )) is randomly generated based on the current location in parameter space. This proposal is then either accepted or rejected and the current parameters are consequently used for the next proposal. If the proposal is more probable than the current sample, then the proposal is always accepted; if the proposal is less likely than the current sample, it is accepted with a certain probability. Therefore, less probable solutions are also sometimes accepted.
The random walk behavior of the Metropolis algorithm can, however, lead to slow convergence because many proposals are rejected (e.g., [23]) and correlated parameters can make the algorithm much less efficient [44]. Modified versions of the Metropolis algorithm improve the proposals by, for example, estimating the covariance matrix on the previous examples and scaling the proposal distribution accordingly (e.g., [45]) or by additionally proposing an alternative step that is chosen based on the rejection of the first proposal [46]. However, it still remains that these methods tend to be less efficient for non-linearly correlated parameters, due to their inherent random walk behavior [44].
Hamiltonian Monte Carlo (HMC) borrows ideas from Hamiltonian dynamics [47] and introduces a 'momentum' variable. As in Hamiltonian dynamics, the momentum variable and the position are building a dynamic system. The entire sampling process is subdivided into simulations of L leapfrog steps each of length . To calculate the momentum at a position in posterior space we need to calculate the derivative of the posterior density for the local environment. With the advances in automatic differentiation in tools like, for example, Theano [48], which is used in the pymc3 framework [49], it is possible to quickly calculate the derivatives of more complex functions.
The two parameters of HMC (L and ), however, must be tuned by hand and strongly influence the performance of the algorithm. The parameter L determines the length of one simulation of Hamiltonian dynamics and, if chosen too large, the leapfrog steps tend to do a U-turn and return to the starting location [23]. The No U-Turn Sampler (NUTS) [44] is an extension to HMC that automatically stops one simulation of Hamiltonian dynamics when the leapfrog step would not increase the distance to the starting location and, therefore, removes the need to tune the L parameter by hand [44]. This way, the NUTS sampler can effectively explore the posterior space.
This algorithm is likely the best choice for our problem, because we can effectively sample from a posterior distribution where we expect some parameters to be correlated and the target distribution might be multimodal.

Likelihood
The unmixing problem according to Section 2.1 is that the abundances ( θ em ) of the endmembers and the space weathering components npFe 0 and mpFe 0 ( θ smFe ) with known albedo spectra, extinction efficiencies and phase functions must be determined. They are the parameters we estimate with the help of MCMC. The known data that is used to calculate the reflectance spectra are denoted as X. The likelihood/goodness of fit of the current parameter estimates is evaluated with both the reflectance spectrum (r mix ) and the continuum removed spectrum (r cr ) to emphasize the absorption band depths over the absolute fit. The continuum removed spectrum is calculated in a straightforward way such that certain fixed wavelengths are chosen as the points of the continuum (see, Figure 3). Usually, the continuum is defined by the convex hull of the spectrum (e.g., [50]). Calculating the continuum in this way is very time consuming considering that the function and its derivative must be evaluated thousands of times. The convex hull is not a linear function such that it is inefficient to implement it in Theano because the derivatives cannot be determined easily. The indices of the local maxima would be dependent on the current spectrum. Because for each of the current parameters of the Markov Chain the continuum is different and depends on the slope, it cannot be computed prior to the sampling. Both the method using the convex hull and the method using fixed wavelengths used in this work are able to characterize the absorption bands usually found on the Moon.  Therefore, the concatenation of the spectrum and the continuum removed spectrum represents the data S. The output of the model is, therefore, a function of the parameters θ em and θ smFe such that: The likelihood is then evaluated by building the sum over all log likelihoods or the product of all probability density function values at the measured data ( S), given a normal distribution around the modeled data and some variance σ 2 , which is also considered a parameter of the model and, therefore, is also sampled.
K denotes the number of channels. If we know that the model's accuracy is wavelength dependent, it is reasonable to include this information in the likelihood function. Therefore, we introduce a wavelength dependent scaling factor ( ζ) that can be used to emphasize, or de-emphasize certain wavelength channels.

Prior
Bayesian inference offers a convenient way to include a-priori knowledge in the sampling process. As described in Equation (11), the posterior is proportional to the likelihood and the prior p(Θ). For multiple priors, the joint probability is given by the product of the individual probability. Because our parameters are usually limited to the interval of [0, 1], a good choice for the priors on our parameters is the Beta distribution, which is defined by two shape parameters α and β as follows [51]: with Γ denoting the gamma-function. Choosing α and β controls the shape of the Beta distribution. Some examples are shown in Figure 4. For example, for a distribution with α = 2 and β = 5, the mode of the distribution is 0.2, but for a distribution with α = 1 and β = 3, the mode is zero and the probability decreases for increasing values. For values outside the interval [0,1], the probability is zero. The prior distribution on the endmember abundances can be set as follows, with the shape parameters being chosen depending on the application: The space weathering components, namely npFe 0 and mpFe 0 , form only a small fraction of the composition. Therefore, we can choose the prior similar to α = 1 and β = 3 in Figure 4 for each abundance of npFe 0 and mpFe 0 ( θ smFe ) The shape parameters will be chosen such that small values are preferred by the sampler.
We do not enforce the sum to one constraint to account for differences in the laboratory setups of the endmembers as well as porosity and other influences. Nonetheless, it is sensible to set a prior on the sum of the endmember weights θ em : This prior is centered around one such that solutions where the sum of the endmembers is close to one are favored. The variance σ 2 sum is chosen depending on the differences among the endmembers and the investigated samples.
The elemental abundances can be relatively well characterized remotely with, for example, gamma ray spectroscopy (e.g., [52]). Spectral parameters in the UV (e.g., [3,4]) or in the NIR (e.g., [5]) can also be used to estimate the elemental abundances, but should be used with care. For example, according to the TiO 2 maps of Sato et al. [4] the swirl Reiner Gamma is deprived of TiO 2 compared to its surroundings, which is physically implausible, as the optical differences are likely due to a difference in maturity and/or compaction [53]. Setting a prior, however, does not rule out solutions that do not confirm with prior knowledge, as long as the probability of the prior does not become zero. If, for example, the spectrum fits better for a high ilmenite concentration, this solution will be found to be more likely.
The RELAB database provides elemental abundances for many of the samples (e.g., [37]) and the elemental composition of the endmembers from Rommel et al. [13] were determined by electron microprobe analysis. We here use endmembers with well-defined elemental abundances (E em , M × N matrix, for M elemental abundances and N endmembers) and use the estimated elemental abundances of the samples ( e sample , M × 1) as prior information. If the samples are mixed directly from the endmembers with weights c (N × 1) the actual elemental abundances can be calculated as For realistic mixtures, this relationship produces highly non-unique solutions for the mineral weights because, for example, the MgO abundance could be attributed to a pyroxene or an olivine. For each sample the elemental abundances are calculated based on the proposed parameters of the mineral abundances θ em according to: As the beta distribution should be centered around the estimated elemental abundance we want to enforce the mode of that distribution to be equal to the estimated elemental abundance. The mode of the Beta distribution for α > 1 and β > 1 is given by: By increasing α (α ≥ 1), the variance of the distribution decreases. Therefore, we choose a α elem approximately according to the desired variance of that particular element and calculate β elem according to Equation (20) such that: Consequently, we define the prior for the elemental abundances to be: This distribution then effectively acts as a prior on the weights of the mineral endmembers, but there are multiple combinations that fit equally well.

Results
The measured spectrum of a mixture is influenced by several factors, namely (1) composition (2), space weathering, (3) grain size, and (4) porosity. To disentangle the effects, we designed several experiments. The grain size is constant for all experiments and is expected to be known. Firstly, we calculate mature synthetic spectra to remove effects not described by the mixing model itself, except for noise. This can also be interpreted as a general validation of the sampling approach. Secondly, to remove the influence of space weathering and to test the general viability of the mixing model we used laboratory mixtures of known endmember spectra. Thirdly, we apply the full framework to unmix the LSCC spectra with endmember spectra from the RELAB catalog and include space weathering particles. These LSCC spectra can be seen as a good representation of different maturity levels and the general mineralogy of the lunar maria.

Synthetic Mixtures
This experiment is designed to estimate the uncertainties of the mixtures when all influences not covered by the model are removed. Spectra are generated with the mixing model including additional Gaussian noise with a standard deviation of σ n = 0.0008. This corresponds to a signal-to-noise ratio (SNR) between approximately 170 and 400 (with SNR= µ r,ch σ n ) depending on the channel. We selected the endmembers according to the categorizations above (see Section 2.1.4). A representative of each of the classes LCP, HCP Type A, HCP Type B, Olivine, Plagioclase, and Ilmenite is selected from the RELAB catalog (https://pds-geosciences.wustl.edu/ spectrallibrary/default.htm, accessed on 1 November 2020) and the SSA values are obtained by inverting the Hapke model and displayed in Figure 5. All endmembers are samples of the 0-45 µm size fraction. The elemental abundances of the endmembers are listed in Table A1. The generated mixtures are based on a mare composition of the endmembers of 10 wt.% hypersthene, 40 wt.% augite, 5 wt.% olivine, 30 wt.% plagioclase, and 15 wt.% diopside. The ilmenite abundance is varied between 0 wt.% and 16 wt.% and the remaining fraction is normalized, such that the endmember weights again sum up to unity. The abundances of npFe 0 and mpFe 0 is varied between 0 wt.% and 2.2 wt.%, and 0 wt.% and 1.2 wt.%, respectively. Ultimately, 48 synthetic mixtures (mix 0 -mix 47 ) are created. The detailed results for all mixtures are listed in Tables A3 and A4.
Generally, all endmembers are measured under similar conditions. According to these conditions the priors have to be chosen. In our case, we do not want to directly employ knowledge about mineral abundances ( θ em ). Therefore, we choose an uninformative prior, that is: Both shape parameters of the Beta distribution are set to one for the prior distribution, such that the prior is a uniform distribution between zero and one. Because negative weights are physically implausible, these samples are always rejected.
The mixtures are synthetically generated, such that we set σ sum = 0.02 to prefer solutions where the sum of the fractional abundances is close to one. The experiment is carried out once with no prior on the elemental abundances and once with an elemental prior on the abundances of TiO 2 and Al 2 O 3 (see also Table 2). For theoretical abundances below 1 wt.% the variance becomes very low when using the method mentioned above. Therefore, we set α elem = 1 and β elem = 10 for the shape of the prior, making the mode zero and the PDF is given in Figure 6. Table 2. Prior setups for the elemental abundances used during this work. The shape parameter β is then calculated according to Equation (21). If α elem > 1 and the theoretical elemental abundance is below 1 wt.%, the shape parameters will be set to α elem = 1 and β elem = 10.  Figure 6. The probability distribution function of the prior on the elemental abundances for a theoretical abundance below 1 wt.%.

SiO
In general, we set the shape parameters for npFe 0 and mpFe 0 to α smFe = 1 and β smFe = 30. These parameters were tuned by hand such that 80 % of the probability density is distributed in the interval [0 wt.%, 5 wt.%]. Higher values than 5 % are physically implausible. Trang and Lucey [54] found about 2 wt.% of smFe at the maximum for the Moon, excluding the polar regions. The model of Trang and Lucey [54], first versions of which were introduced in Trang et al. [55] and Lucey and Riner [16], differs from the model of Wohlfarth et al. [24], for example, Trang et al. [55] do not consider the influence of the phase function. Therefore, the predicted concentration of npFe 0 and mpFe 0 is not necessarily the same, leading us to set a relatively uninformative prior on each of the smFe abundances. The prior on the abundance of smFe particles is consistent throughout this work and is not changed between experiments.
Generally, the mean predicted solution does not necessarily correspond to the actual mixture coefficients used, but the theoretical values are always part of the solution. Classical unmixing approaches are limited to just one solution and small changes may lead to different optimal solutions, which are also not necessarily equal to the ground truth values.
Ilmenite shows, compared to the mafic minerals, an almost featureless spectrum in the NIR wavelength range. Because of its low albedo, the abundance of ilmenite also dampens the absorption bands and darkens the overall spectrum. This is, to some extent, similar to the spectral influence of mpFe 0 particles. If the absolute reflectance is poorly defined, so is the ilmenite abundance. Consequently, depending on the prior setup and the sum constraint the uncertainties of ilmenite are different. If no prior on the elemental or mineral abundances is set, the abundance of ilmenite is relatively poorly defined. For the mixtures that, in theory, do not contain ilmenite, a majority of the accepted mixtures contain significant amounts of ilmenite (see Figure 7). A combination of ilmenite and other endmembers can, therefore, equally explain the measured spectrum. This effect can be mitigated when a prior is used (see Figure 8). For the fresh mixture (mix 0 ), without any smFe particles, the uncertainty of the spectrum is small (see Figure 9). The mineral uncertainties are small, but for the other mixtures, the uncertainties significantly increase. The majority of posterior density of plagioclase (PLG) is generally lower than the actual abundance, and for olivine (OLV) and/or diopside (CPXA) the majority of posterior density is below the actual value. This suggests that a mixture of olivine/diopside and ilmenite produces very similar results to plagioclase. Olivine and diopside have similar absorptions at about 1-µm, making them interchangeable to some extent when the maturity increases. For example, for mix 15 there are two distinct modes with either high concentration of diopside or olivine. The mineral uncertainties of mix 5 are large, but the spectrum is clearly defined. This highlights that there are several solutions that produce nearly identical spectra. Sometimes, however, the modal abundance values do not sum up to one, because, for example, just one of the distributions is skewed. This can lead to mode predicted spectra outside of the confidence band. The modal abundance values were likely never sampled in this particular combination.
Ilmenite is the only endmember that contains significant amounts of TiO 2 . Including prior knowledge about the abundance of TiO 2 can, therefore, favor solutions that are close to the true ilmenite value. In Figure 8 the violin plots of the acceptable solutions considering the TiO 2 and Al 2 O 3 prior are displayed. The abundance of ilmenite is then much more clearly defined. The uncertainties of the other minerals also mostly decrease, while the uncertainties of the spectra remain similar.
The results for the most mature spectrum with the highest abundance of ilmenite are shown in Figure 10. With the prior, the uncertainties of all endmembers can be reduced. The confidence plots remain very similar. Figure 11 shows the distribution of the differences between theoretical values and mean predicted abundance. By including the prior the differences can be reduced. The mean of olivine is usually above the true value, while the mean of plagioclase and clinopyroxene is usually below the theoretical value. The theoretical abundance of olivine is relatively small compared to the other minerals, consequently, the mean is usually higher. The uncertainties of olivine and clinopyroxenes are high for both with and without the prior. This can be seen in Figure 12. Including prior knowledge about the elements Al 2 O 3 and TiO 2 the uncertainties of all mineral abundances are mitigated. This is despite the fact that these elements are almost depleted in the pyroxenes and olivine.
The abundance of npFe 0 and mpFe 0 is usually well defined, as changes to the abundance of npFe 0 and mpFe 0 can have a strong influence on the modeled spectrum. In Figures 7, 8 and 10 it can be seen that the predicted npFe 0 and mpFe 0 do not change significantly with the introduction of a prior. In Table 3 the correlation coefficients between true smFe abundance and the standard deviation of the minerals is listed. For all minerals it can be seen that the correlation coefficients are positive. Therefore, the higher the abundance of smFe is, the larger the uncertainties become. The uncertainty of ilmenite is, however, uncorrelated to the smFe abundance, even if no prior knowledge is used. Because ilmenite does not show a prominent absorption band such a feature cannot be obscured by the smFe particles. This is visible in Figure 13. While the overall uncertainties decrease with a prior on the elemental abundances, the trend that the uncertainties increase with increasing maturity is visible for both with and without a prior.     Figure 10. Results for the most mature spectrum, which also has the highest ilmenite concentration. The spectrum becomes very dark. The abundance of npFe 0 and mpFe 0 is well defined either with or without the prior information. The uncertainties of the mineral abundances decrease when including the prior. The yellow diamonds in the violin plots mark the true abundances that were used for mixing the endmembers.  Generally, it can be seen that the ilmenite abundance is difficult to determine reliably without prior information, for both immature and mature spectra. For mature spectra, the uncertainties of all mineral abundances become very large, even for an almost perfect case, where influences which are not modeled, are removed. This illustrates that unmixing on the Moon is an ill-posed problem, and additional knowledge about the composition is necessary. The addition of a prior to indirectly include information about the plagioclase and ilmenite abundance is effective in reducing some of the uncertainties.

Fresh Laboratory Mixtures
This experiment is set up to unmix real laboratory mixtures (see also Rommel et al. [13]) from laboratory endmembers. Thus, we select the same endmembers as in Rommel et al. [13] but omit one of the two labradorites (plagioclases) as the two are very similar, and we want to select endmembers that are spectrally distinct. We also omit the diallagite endmember, to reduce the catalog to the same size as for the other experiments and this endmember is only present in one mixture, which is, consequently, also omitted here. The amphibole endmember pargasite (PRG) was selected in Rommel et al. [13] to test whether uncommon minerals can also be detected by the algorithm. We keep this endmember for our analysis even though it is not present on the Moon.
For the inversion of the Hapke model, the parameters are chosen in the same way as in Rommel et al. [13] for comparability. Therefore, the term for the shadow hiding opposition effect is not set to zero, but B SH 0 = 3.1 and h SH = 0.11 are adopted from Rommel et al. [13]. The resulting SSA spectra are shown in Figure 14.
The elemental abundances of the endmembers were determined in Rommel et al. [13] based on electron microprobe analysis in the laboratory (see Table A5) and the modes and shape parameters are determined according to Section 2.2.3. The results were obtained by using no elemental priors on the one hand and using a prior for the TiO 2 abundance on the other hand (see first two rows of Table 2).
The true values of the abundance of submicroscopic iron particles in this experiment is zero. In order to test whether fresh and mature spectra can be reliably distinguished from each other, we still employ the same prior for smFe as for the other experiments.
All endmembers and mixtures were measured in the same laboratory and the mixtures are directly mixed from the endmembers; therefore, the weights of the endmembers should sum up to a value close to one. We set a prior on the sum according to Equation (17) with a standard deviation of σ sum = 0.02.
The detailed results of the unmixing without a prior on the elemental abundances and with a prior on the TiO 2 abundance are listed in Tables A7 and A8, respectively. Generally, the mean predicted results are very similar to those of Rommel et al. [13] considering that the endmembers labradorite A and diallagite were omitted. Even though we did not use a strong prior on the abundance of iron particles, the predicted abundances are almost exclusively below 0.02 wt.% or a fraction smaller than 2 × 10 −4 .
For the version without a prior on the TiO 2 abundance, ilmenite is erroneously predicted in the mixtures of pargasite and plagioclase. For example, Figure 15 shows the violin plots of all similarly well-fitting solutions. It can be seen that ilmenite is predicted to be present, and a solution without ilmenite is not part of the uncertainty. Therefore, a solution with ilmenite fits better than one without ilmenite. If we, however, prefer solutions where the contribution of ilmenite is close to zero, we obtain a better result and the errors in pargasite and plagioclase abundance are reduced. In Figure 16 the confidence plots are plotted. Here it can be seen that the mean and mode predicted spectra fit quite well for both prior setups. By introducing the prior for the TiO 2 abundance, the uncertainties of the fit increase (σ = 0.006 compared to σ = 0.003 for no prior). Therefore, slightly less well fitting solutions are accepted if the ilmenite abundance confirms better with the prior knowledge. In the mixtures containing plagioclase and ferrosilite, plagioclase is confused with pargasite if no prior knowledge is included. The two endmembers are spectrally very similar and are both bright (see Figure 14). In Figure 17, an example violin plot is shown, where the confusion of plagioclase and pargasite can be seen in Figure 17a. Small contributions of ilmenite can already dampen the absorption band of the pargasite endmember, such that the similarity to plagioclase increases. The mean predicted abundance of ilmenite is quite low with 0.48 wt.% for the version without a prior, but still leads to a significantly different solution. The confidence plots in Figure 18 look nearly identical, except that for the TiO 2 prior to the uncertainty of the 0.9 µm absorption band slightly increases. This suggests that both solutions, with and without ilmenite, represent regions of high posterior density. It is, however, very difficult for the sampler to traverse from one region to the other, as the two are far apart and due to the prior on the sum it is unlikely for a limited number of samples that both regions are sampled. However, with a relatively uninformative prior as in Figure 6 this ambiguity can be removed, such that the solutions with no ilmenite are preferred. This trend illustrated by the two examples above, can also be seen in Figure 19. The differences between mean predicted solution and true abundance for the endmembers augite, olivine, and ferrosilite are not improved by including prior knowledge about the TiO 2 abundance. This is likely because of inaccuracies of the data and the simple mixing model and also a problem in Rommel et al. [13]. The inherent simplifications of the Hapke model makes it such that the fractional abundances are not necessarily equal to their actual weight fraction (e.g., [56]). Similar to a regularization term in classical optimizations, we can use the prior to favor certain solutions, without imposing hard constraints. Generally, without prior knowledge the abundance of plagioclase is usually underestimated, while the abundance of ilmenite and pargasite is overestimated. The uncertainties are mostly small, both with and without TiO 2 prior (see Figure 20).

LSCC Samples
In the final experiment, we unmix spectra from the LSCC catalog with laboratory mineral endmembers and the npFe 0 and mpFe 0 particles. The LSCC catalog contains mature samples that were returned from the lunar surface by the Apollo missions and is, therefore, a good representation of lunar mineralogy. We use the <45 µm size fraction spectra for the unmixing, and because there is no data on the mineral abundances of this size fraction we use the mineral abundances from the 20-45 µm size fraction, which correspond to a large portion of the weight, as a reference for evaluating the results. The endmembers are the same samples from the RELAB library listed in Section 3.1. These endmembers do not correspond to the actual constituent minerals and there are additional minerals that are not covered by the endmember catalog, for example, spinel. Additionally, the endmember catalog does not contain agglutinates in comparison to the LSCC samples with large fractions of agglutinates. The ground truth abundances are obtained from the tables of Taylor et al. [18,19]. The abundances of the minerals were grouped into clinopyroxene, orthopyroxene, olivine, plagioclase, and ilmenite. The abundances were then normalized to the mineral fraction because we are not considering glasses. We have seen that, without a prior on the elemental abundances, the uncertainties become very large for mature spectra. Therefore, we will use the same prior once on TiO 2 and Al 2 O 3 abundances, which in the previous experiments reduced the uncertainties effectively for plagioclase and ilmenite, but also for many of the other minerals. Additionally, we will employ a prior on all elements with a relevant abundance on the Moon, which might contribute to distinguishing between minerals (see also Table 2). The elemental abundances of the LSCC samples are listed in Table A9. The spectra were measured under laboratory conditions, but the mixtures are much more complex compared to the laboratory mixtures from Section 3.2. We are, therefore, not enforcing the sum-to-one constraint as strongly, as we set σ sum = 0.1. Because plagioclase and ilmenite have both almost featureless flat spectra, the uncertainties between these two increase significantly the more relaxed the prior on the sum is. Therefore, the prior knowledge about TiO 2 and Al 2 O 3 abundance becomes even more essential.
Detailed results for the unmixing of all LSCC spectra are listed in Tables A10 and A11 for the two experiments with different priors. Generally, the uncertainties are similar between both versions, but for some samples the priors on the other elements help in reducing the uncertainties. This can be seen in Figure 21. Therefore, we will focus on the Ti/Al prior version. In some cases, the plagioclase abundance is underestimated; such an example can be seen in Figure 22. For the violin plots, the two clinopyroxenes are grouped together. The spectrum fits quite well and the true plagioclase abundance is not part of the uncertainty. Therefore, a prior cannot mitigate the differences between mean predicted and theoretical abundance. Instead, olivine and orthopyroxene abundances are overestimated. A similar problem can be seen for 12001 and 15041. These samples are the ones with the highest fraction of agglutinic glasses (12001: 56 Figure 22. Results for LSCC sample 10084 with Ti/Al prior. Yellow diamonds mark the estimated theoretical mineral fractional abundances. Plagioclase is underestimated, while olivine is generally overestimated. This is a mature spectrum with an I s /FeO value of 88, therefore, also the estimated npFe 0 and mpFe 0 abundance is high. The spectrum is generally very dark. (a) Endmember abundances (10084); (b) smFe abundances (10084); (c) Confidence plot (10084). Figure 23 shows the typical results for a mature mare spectrum. Most of the abundances are relatively well defined with the prior on the TiO 2 and Al 2 O 3 abundance. Compared to sample 10084 in Figure 22, there are significantly fewer npFe 0 and mpFe 0 particles predicted to be present. As sample 14260 contains less FeO (9.65 wt.% compared to 14.8 wt.% in 10084) it is likely that the smFe abundance saturates more quickly, such that additional smFe cannot easily be created as quickly anymore. The ilmenite abundance is comparatively small, which might also be a contributing factor to the uncertainties being relatively small. Similarly, for the mature highland spectrum in Figure 24, the iron abundance in the highlands is generally much smaller so that less smFe can be created. According to the I s /FeO parameter, this is also a very mature spectrum and the absorption bands are weakly pronounced. However, the smFe abundance is small compared to the mare samples. In Figure 25, the results for sample 79221, with both the Ti/Al prior and also with the prior on all major elemental abundances, is shown. This sample is rich in ilmenite and also relatively mature, therefore the uncertainties are quite high. In such cases, the results can be improved when additional elemental abundance information is available. For most other LSCC samples the Ti/Al prior is, however, sufficient to constrain the procedure. The npFe 0 and mpFe 0 abundance is not much effected by the change in prior knowledge. In Figure 26 the correlations between I s /FeO or I s and mean predicted npFe 0 or mpFe 0 abundance are displayed. The abundances are only weakly correlated with the I s /FeO maturity index. In contrast the correlation between the ferromagnetic resonance I s and the abundance is very strong. Even for the small dataset this correlation is significant. Because the I s is a measure on the amount of iron particles in the size fraction between 4 nm and 33 nm [57], this suggests that the abundance of npFe 0 can be obtained relatively reliably even with a limited endmember catalog. This difference between I s /FeO and I s also supports the findings of Trang and Lucey [54] that the saturation limit of npFe 0 particles is dependent on the FeO abundance. The saturation limit increases with increasing FeO abundance. Figure 27 shows the relationship between npFe 0 and mpFe 0 abundance in the LSCC samples. Usually, the factor between the abundance of npFe 0 and mpFe 0 abundance is approximately constant, where the amount of npFe 0 is typically twice as high as the amount of mpFe.

Discussion
Generally, the three experiments show that a probabilistic perspective provides important insights into the unmixing process overall. When trying to unmix lunar spectra many factors have to be considered to obtain reliable results. A non-probabilistic approach needs to make many assumptions about the composition and the maturity of the surface to obtain plausible results. Even if these assumptions are accurate only a single best fit solution is obtained. All experiments shown in this work illustrate that there is no unique solution and that a variety of combinations can produce similar results and, thus, a number of different interpretations are possible.
The synthetic mixtures in Section 3.1 were used to estimate the influence of the darkening agents on the overall uncertainties of the mineral abundances. The results showed that with increasing smFe abundance the uncertainties increased. Ilmenite or TiO 2 also contributed to uncertainties in the data, but there appears to be no direct relationship between true ilmenite abundance and uncertainties. However, introducing prior knowledge about the TiO 2 and Al 2 O 3 abundance can significantly reduce the uncertainties of all mineral abundances as well as the errors between mean predicted and true mineral abundances and is, therefore, suitable for the application of the approach to real lunar spectra. For very mature spectra the uncertainties become very large. This illustrates the need for a probabilistic approach to the problem of unmixing spectra for mineral and smFe particle abundances. The abundance of the smFe particles could, however, still be reliably detected for any maturity level. The uncertainties of the minerals all increased with increasing smFe concentration, except for ilmenite. This might be explained by the fact that ilmenite with its featureless spectrum is already difficult to determine accurately if the absolute brightness of the sample is not well defined. Ilmenite, therefore, is not easily distinguished by any absorption bands, which become more difficult to quantify with increasing maturity.
Ambiguities between spectrally similar endmembers such as plagioclase, pargasite, and ilmenite in Section 3.2 can lead to larger errors (up to 57 wt.%) in the best fit solution of the Hapke model and the measured spectra. When these ambiguities are mitigated with using prior information, the errors are usually well below 15 wt.%. All models have to include some simplifications such that they cannot perfectly represent reality. These problems inherent to the mixing model itself usually cannot be predicted by the MCMC sampler. Ambiguities and uncertainties due to the choice of model parameters on the other hand can be quantified and may indicate problems with the choice of endmembers and/or assumptions made about the data. In some cases where plagioclase and pargasite were confused, two local minima far apart from each other in parameter space were not found by the sampler, even though the mean σ of the confidence plots is very similar. Generally, the fresh laboratory spectra could easily be distinguished as such. The abundance of smFe is predicted to be effectively absent.
The final experiment of unmixing the LSCC spectra with RELAB mineral endmembers demonstrated that this approach can be applied to the lunar surface. When the abundance of agglutinitic glasses is high, olivine and plagioclase abundances tend to be poorly defined. The reconstruction is very accurate, but the uncertainties of the mineral abundances are high for mature spectra. This shows that classical optimization techniques are prone to fail on the task of unmixing the spectra with accurate abundances. The abundances of the smFe particles are, however, well defined. The results align with the findings of Trang and Lucey [54] that the npFe and mpFe abundance for mature spectra depends on the FeO abundance. Furthermore, the predicted abundance of npFe and mpFe particles correlates well with the I s parameter of the LSCC samples. For npFe, this can be interpreted to be a sign that the algorithm predicts plausible results. The npFe and mpFe abundances are also strongly correlated, such that the correlation of mpFe might be attributed to this nearly constant factor. According to Lucey and Riner [16] mpFe sized particles are not present in the samples investigated, but the larger particle size is necessary to model the effects of space weathering. Wohlfarth et al. [24] further suggested that clusters or layers of npFe particles that effectively act as a larger mpFe particle are responsible for the darkening. This effect can also be simulated with light scattering theory as in Arnaut et al. [58]. This would suggest that npFe and mpFe (or clusters of npFe) are both created by the same processes. The cluster/layer theory is consistent with the results of the LSCC unmixing.
All endmembers and investigated mixtures have the same average grain size, which usually is not the case for remotely measured datasets. For the Martian surface, Lapotre et al. [22] and Lapotre et al. [59] investigated the influence of the grain size on the uncertainties and showed that the solution to this inverse problem is highly non-unique. Integrating the grain size into the existing framework would add another layer of complexity to the model, making it even more challenging to obtain meaningful results. Therefore, prior knowledge about the grain size would likely be necessary to constrain the procedure. If the grain size varied significantly between endmembers, systematical errors with respect to their relative fraction would be introduced. Some simplifications to the model, however, have to be made for this highly complex inverse problem to be solvable. The effects of the grain size in combination with space weathering components and the other agents that darken the overall spectrum, discussed in this work, can be investigated in possible future research.
The laboratory samples examined in this work all have a relatively high SNR in comparison to, for example, the Moon Mineralogy Mapper data. When the SNR decreases, the standard deviation of the likelihood function is likely to increase, such that more solutions might be considered acceptable. This, in turn, most likely also increases the uncertainties of the fractional abundances. Nonetheless, we can indirectly measure the influence of the SNR based on the average differences between model prediction and observed spectra. In our work, which follows a Bayesian approach, the variance of the differences is a part of the model and is thus also estimated (cf. Equation (13)). However, averaging over neighboring pixels, where the composition is similar, improves the SNR. Due to the high computational effort to sample from these complex distributions, the number of pixels that can be investigated is limited. Therefore, clustering the data and applying the unmixing to the centroids also mitigates the effects of a low SNR and makes computing the abundances over large areas feasible (e.g., [60,61]).

Conclusions
Unmixing can be used to obtain quantitative information about the composition of remote surfaces. Building on the Hapke model (e.g., [9]) and the space weathering model of Wohlfarth et al. [24], we have developed a combined unmixing model. The parameters of the model and their corresponding uncertainties were obtained by Bayesian inference with an MCMC sampler. The results showed that, especially on the Moon, there is no unique solution that would perfectly describe the measured spectra. The absorption bands at around 1-µm and 2-µm are the main features for estimating the abundances. These features are obscured by submicroscopic iron particles, making it harder to quantify the corresponding mineral abundances. Some mineral combinations even become gradually interchangeable, only slightly changing the modeled spectrum. Therefore, additional knowledge is necessary to constrain the procedure. We find that information about the elemental abundances, which can be obtained from different independent data sources, is well suited to helping to distinguish between different minerals and we prefer solutions that correspond better with this prior knowledge. Inherent simplifications of the mixing model can lead to systematic errors. In such cases, the uncertainties might still be small, therefore, it is difficult to predict such problems with the model itself if no ground truth is available. Overall, the uncertainties can be accurately estimated from the MCMC sampler. For very mature spectra, these uncertainties become large, suggesting that a classical optimization approach is destined to fail at consistently predicting accurate results. The abundance of submicroscopic iron particles, responsible for the optical effects of maturity, could reliably be determined with the proposed approach. The results of the unmixing on LSCC spectra showed that the probabilistic unmixing approach is applicable to realistic lunar spectra. The knowledge about the uncertainties inherent to the problem are, however, essential for the interpretation of the results obtained for lunar spectra. Data Availability Statement: RELAB database (https://pds-geosciences.wustl.edu/spectrallibrary/ default.htm, accessed on 1 November 2021). LSCC catalog of lunar returned samples (https://pgi.utk. edu/lunar-soil-characterization-consortium-lscc-data/, accessed on 12 October 2021). Laboratory mineral and mixture spectra from Rommel et al. [13].

Acknowledgments:
We want to thank the RELAB and LSCC teams for making the data to the catalog of spectra available publicly. We also want to thank Daniela Rommel for all the hard work, measuring the spectra and preparing the mixtures in Section 3.2. And finally the authors wish to thank the Department of Geochemistry of the University of Göttingen for providing the mineral samples examined in Section 3.2.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Appendix A.1. Synthetic Mixtures Table A1. Elemental abundances of endmember catalog for experiments in Sections 3.1 and 3.3. The sample spectra and elemental abundances are taken from the RELAB library (https://pds-geosciences.wustl.edu/spectrallibrary/default.htm, accessed on 1 November 2020).   Table A9. Elemental abundances used for the LSCC samples taken from Taylor et al. [18] and Taylor et al. [19]. The elemental and mineral abundances of the 20-45 µm size fraction were used for the priors, but the 0-45 µm size fraction spectra were used for the unmixing. The 0-45 µm size fraction mineral abundances were not available.