SMAUG v1.0 – a user-friendly muon simulator for transmission tomography of geological objects in 3D

Knowledge about muon tomography has spread in recent years in the geoscientific community and several collaborations between geologists and physicists have been founded. As the data analysis is still mostly done by particle physicists, we 15 address the need of the geoscientific community to participate in the data analysis, while not having to worry too much about the particle physics equations in the background. The result hereof is SMAUG, a toolbox consisting of several modules that cover the various aspects of data analysis in a muon tomographic experiment. In this study we show how a comprehensive geophysical model can be built from basic physics equations. The emerging uncertainties are dealt with by a probabilistic formulation of the inverse problem, which is finally solved by a Monte Carlo Markov Chain algorithm. Finally, we 20 benchmark the SMAUG results against those of a recent study, which however, have been established with an approach that is not easily accessible to the geoscientific community. We show that they reach identical results with the same level of accuracy and precision.


Introduction
Among the manifold geophysical imaging techniques, muon tomography has increasingly gained the interest of geoscientists 25 during the course of the past years. Before its application in Earth sciences, it was initially used for archaeological purposes. Alvarez et al. (1970) used this method to search for hidden chambers in the pyramids of Giza, in Egypt; an experiment which was recently repeated by Morishima et al. (2017), as better technologies have continuously been developed. Other civil engineering applications include the monitoring of nuclear power plant operations (Takamatsu et al., 2015) and the search for nuclear waste repositories (Jonkmans et al., 2013) as well as the investigation of underground tunnels (e.g. Thompson et 30 al., 2020;Guardincerri et al., 2017). A serious deployment of muon tomography in Earth sciences has only begun in the past https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License.
Other experiments have been performed in order to explore the geometry of karst cavities in Hungary (Barnaföldi et al., 35 2012) and Italy . Further studies (see also review article by Lechmann et al., 2021a) have been conducted by our group to recover the ice-bedrock interface of Alpine glaciers in central Switzerland (Nishiyama et al., 2017;.
The core component of every geophysical exploration experiment is formed by the inversion, which might be better known to other communities as fitting or modelling. This is where the model parameters are found, which best fit the observed data. 40 Up until now, this central part has mostly been built specifically to meet the needs of the experimental campaign at hand. On the one hand this approach has the advantage of allowing the consideration of the peculiarities of particle detectors, their data processing chain, and other models involved (e.g. the cosmic ray flux model). On the other hand, when every group develops a separate inversion algorithm, the reconstruction of the precise calculations performed in the data analysis procedure becomes a challenge. For a researcher who is not familiar with the intricacies of inversion, this might even be tougher. We 45 thus see the need for a lightweight programme that incorporates a structured and modular approach to inversion, that also allows users with little inversion experience to familiarise themselves with this rather involved topic. This programme can be used to directly analyse experimental data in a stand-alone working environment, and the modules and theoretical foundations can be adapted, customised, and integrated into new programmes. For this reason, the code is built in the programming language Python as to facilitate the exchange between researchers and to enhance modifiability. Moreover, the 50 source code is freely available online (see code availability section below).
To facilitate the further reading of our code, we introduce the reader at this point to our benchmark experiment, to which we will refer on multiple occasions throughout this work. The experimental campaign is explained in detail in Nishiyama et al. (2017) and thus we will resort to a description of the experimental design at this point. In the Nishiyama et al. (2017) study we aimed at we aimed at recovering the ice-bedrock interface of an Alpine glacier in central Switzerland. Figure 1 shows 55 that we had access to the railway tunnel of the Jungfraubahnen, where we installed three detectors. In our measurement we recorded muons from directions that consisted purely of rock and others where we knew that the muons must have crossed rock and ice (see the two cones in Fig. 1). From the former it was possible, together with laboratory measurements, to determine the physical parameters of the rock more precisely. Conversely, we utilised the measurement of the latter to infer the interface between rock and ice underneath the glacier. Finally, we will also use the results of that experiment (Nishiyama 60 et al., 2017) to verify our new algorithm in the present study.

Inversiona modular view 70
The goal of every muon tomography study is essentially to extract information on the physical parameters of the radiographed object through a measurement of the cosmic ray muon flux and an assessment of its absorption as the muons cross that object. In geological applications these objects are almost always lithological underground structures such as magma chambers, cavities, or other interfaces with a high-density contrast. The reconstruction of the geometry of such structures can only be achieved if the measured muon data is compared to the results of a muon flux simulation. As stated 75 earlier, this is the basic principle of the inversion procedure. However, the aforementioned "muon flux simulation" is not just a simple programme, but it consists of several physically independent models that act together. Taking a modular view, we will call these models "modules" from here on, as they will inevitably be part of a larger inversion code. We have visualised the components that are necessary to build an inversion in Fig. 2.
The first of the modules is the input module for the experiment results, which also considers the detectors that were used in 80 the experiment. Typical detector setups include nuclear emulsion films (e.g. Ariga et al., (2018), cathode chambers (e.g. Oláh et al., 2013), scintillators (e.g. Anghel et al., 2015) or other hardware solutions. Although the detailed data processing chain may be comprehensive, the related output almost always comes in the form of a measured directional (i.e. from various incident angles) muon flux, which will be the input to the inversion scheme. Here, we primarily work with the premise that the muon flux data and the associated errors are given. The corresponding errors can then be furnished to the code by means 85 of an interface. The simulation module on the other hand, consists of four autonomous components (see Fig. 2). First, a cosmic ray muon flux model is necessary, which describes the muon abundance in the atmosphere and which is generally dependent on the muon energy, its incidence angle, and the altitude of the detector location. Lesparre et al. (2010)  compare various muon flux models that may be incorporated into an extensive simulation. Second, it is necessary to model the spatial distribution of the detectors as well as the initial distribution of the lithologies. Related pre-existing software 90 solutions mainly comprise GIS-and geological 3D-modelling applications, that excel at capturing and compiling geological information from various sources (e.g. field, maps, etc.) into a spatially organised database. Third, the lithologies consisting of different minerals have to be translated to a set of parameters, which are a necessary input for the subsequent physical simulation. This can be done by a rock model (e.g. Lechmann et al., 2018), which considers the effects of the mass density as well as the average atomic mass and charge of the rock as a function of its mineralogical composition. Lastly, the muon 95 fluxes at the detector sites have to be simulated by means of a muon transportation model, which calculates all physical processes by which a muon loses kinetic energy while travelling through matter. The particle physics community has a great variety of particle simulators, the most prominent being GEANT4 (Agostinelli et al., 2003), a Monte Carlo based simulator.
These have the advantage that stochastic processes resulting in energy loss are simulated according to their probabilistic occurrence -an upside that has to be traded off for longer computation times. In contrast to obtaining the full energy loss 100 distribution, lightweight alternatives often resort to the calculation of only the mean energy loss. The solution of the resulting differential equation can even be tabulated, as has been done by Groom et al. (2001). The interplay of these four submodules allows for the simulation of muon fluxes at the detector sites that are mostly located in an underground environment.  consists of a model for rocks, detectors, the cosmic ray flux and a particle physical model on how muons lose energy upon travelling through rocks. These models allow for a synthetic data set to be computed, which will be compared to the actual measured data from the experiment. An optimisation problem then solves for the best set of parameters.

The need for a consistent inversion environment
The sole combination of the aforementioned four submodules does not fully justify the need for a new software, as cosmic 110 ray flux models as well as rock models can also be programmed within existing Monte Carlo simulators such as GEANT4 (Agostinelli et al., 2003) or MUSIC (Kudryavtsev, 2009). Unfortunately, the application of such a Monte Carlo approach requires a rather good understanding of programming and nuclear physics processes. Thus, it might prove time-consuming to programme a specific code. Moreover, these codes are often written in a specialised programming language such as C++.
Third, the compatibility between different modules (e.g. cosmic ray flux and energy loss) may be severely hampered, if the 115 programme interfaces are not taken into consideration. It might be even worse if the two modules are written in two different programming languages. In addition, one has to carefully evaluate the benefit of such an undertaking, especially if the resulting code will most likely be tailored only to a specific problem. We thus see the need for a versatile, user-friendly simulator, which allows users not only to quickly perform the necessary calculations without the need of additional coding, but also tailor the individual models to custom needs. A new simulator can be more useful if an inversion functionality is 120 already included. As can be seen in Fig. 2, the inversion compares the simulated flux data with the measured ones. It also attempts to reduce the discrepancy between measurements and simulations by optimising the parameters in the simulation, namely material density and the thickness distribution of the overlying materials. This results in a density-or structural rock model, which best reproduces the measured data. As the mathematical optimization in muon tomography generally is nonlinear, one has to employ nonlinear solvers or even Monte Carlo techniques. This circumstance encourages us to work 125 with a lightweight version of a muon transport simulator, because a nonlinear inversion of Monte Carlo simulations, although mathematically preferable, is computationally prohibitive. This allows us to make use of methods from the Bayesian realm, that thrive when measurements from different sources have to be combined into a single comprehensive model. With the code presented in this paper, we aspire to make muon tomography accessible to a broader geoscientific community, as the know-how in this field is mainly concentrated in particle physics laboratories. We want to provide the 130 tools for Earth scientists, or users that are mainly focused on the application of the method, so that they can perform their own analyses.
In this contribution we present our new code, SMAUG (Simulation for Muons and their Applications UnderGround), that allows a broader scientific community to plan and analyse muon tomographic experiments more easily, by providing them with data analysis and inversion tools. Specifically, we describe the governing equations of the physical models, and the 135 mathematical techniques that were used. Chapter 2 depicts how the muon flux simulation is conducted by its submodules and how a muon flux simulation is performed. Chapter 3 then dives into the inversion module and explains how the parameters of the inferred density/rock-model can be estimated based on measured data. This chapter includes a description of the model and data errors and an explanation on how a subsurface material boundary can be constructed. Chapter 4 provides a short overview of the program, explaining which functionality can be found in which source code. Chapter 5 140 discusses the model's performance based on the data that we collected in the framework of an earlier experimental campaign https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License.
(see supplement of Nishiyama et al., 2017). Chapter 6 then concludes this study by outlining a way of how this code can be developed further to fit the needs of the muon tomography and geology community.

The forward model: Muon flux simulation
In geophysical communities this part is generally known as the forward model, i.e. a mathematical model which calculates 145 synthetic data for given "model" parameters. In muon tomography experiments this forward model consists of different physical models which are serially connected.

Cosmic ray flux model
The nature of the data used in muon tomography generally consists of several counts within a directional bin, defined by two polar and two azimuthal angles. Additionally, the measurement is taken over a defined period of time, as well as over a given 150 extent within the detector area. The simulated number of muons, in the i-th bin, can be calculated by this integral, (1) Here, T denotes the exposure time interval, A the detector area, Ω the solid angle of the bin and E the energy range of the muons that were able to be registered by the detector. There are various differential muon flux models, also referred to as the integrand in Eq.
(1), that can be employed at this stage. Lesparre et al. (2010) provide a good overview on the different flux 155 models, which can broadly be divided in two classes. On the one hand there are theoretical models, which capture the manifold production paths of muons and condense them in an analytical equation, e.g. the Tang et al., (2006) model. They contrast, on the other hand, with empirical models that were generated by fitting formulae to the results of muon flux measurements. The model of Bugaev et al. (1998) falls into this category, with later adjustments for different zenith angles (Reyna, 2006) and altitude (Nishiyama et al., 2017), which are also utilised in this study. The details of the formula are 160 explained in Appendix A. The evaluation of Eq. (1) is rather cumbersome as strictly speaking several of the integration variables depend on each other. We may facilitate the calculation by considering that the differential muon flux model is only dependent on energy, E, and zenith angle, whereas the effective area, Δ , , is solely dependent on the orientation of the bin. This is the case because muons do not necessarily hit the detector perpendicularly, such that the effective target area is usually smaller. By averaging over the zenith angle and keeping the bin size reasonably small, we may approximate 165 Eq. (1) by where Δ is the exposure time and Δ , is the effective detector area. ΔΩ is the solid angle, ̂ and ̂ are the mean azimuth and zenith angle of the i-th bin, respectively. , describes the energy needed for a muon to enter the detector. Δ , has to be scaled by the cosine of the angle between the bin direction and the detector facing direction, 170 which can be calculated using the formula for a scalar product, where ⃗ is the normal vector to the detector surface and (̂,̂) is the mean vector of muon incidence within the i-th bin, both of which can be chosen to feature unit length. Evaluating the scalar product in spherical coordinates, Eq. (4) 175 Here, and are the zenith and azimuth angles of the detector facing direction. It is important to note that except for , all variables in Eq.
(2) are predetermined by the experimental setup (Δ , ΔA) as well as by the data processing (̂, ̂) , such that the number of muons , can be interpreted as a function of one variable, , only.

Muon transportation model
Since muons permanently lose energy when travelling through matter, they also need a certain amount of energy to enter the 180 detector. This energy, , , was introduced in Eq.
(2) and is called the cutoff energy. If the detector is now positioned underground, the muons have to traverse more matter to reach the detector and consequently need a higher initial energy to reach the target. For this purpose, we introduce the new variable 0 , which refers to the energy needed to penetrate the detector, and we reinterpret , as the minimum energy that is required to traverse the matter and to be registered at the detector. For the goal of studying the interactions between particles and matter, physicists regularly use energy loss models. 185 We base our calculations in large parts on the equations of Groom et al. (2001), where the energy loss of a muon along its path is described by an ordinary differential equation of 1 st order, In Eq. (5), denotes the density of the traversed material, and and are the ionisation loss and radiation loss parameters respectively. The radiation loss parameter groups the effects related to the bremsstrahlung, , the pair-production, 190 , and the photonuclear interactions, ℎ , where ( , ) = ( , ) + ( , ) + ℎ ( , ).
Each of the radiative process is, in turn, calculated through where ∈̃= {bremsstrahlung, pair-production, photonuclear} is the set of radiative processes, is Avogadro's number, 195 is the atomic weight of the traversed material, is the fractional energy transfer, and ⁄ the differential cross-section  Groom et al. (2001). The only exception in Eq. (6) is , which is calculated after GEANT4 (Agostinelli et al., 2003). We selected the solution of these latter authors because it is computationally less time consuming. As the two results agree within 1 %, we deem it acceptable to exchange the two differential cross-sections. 200 Because Eq. (5) describes the energy loss in response to the interaction with a single-element material, certain modifications have to be made to make it also valid for rocks, which in this context represent a mixture of minerals and elements. In this case, the modified equation takes an equivalent form to Eq. (5) when replacing , , with their mixture counterparts We show in Appendix B3 how the rock model (explained in Ch. 2.3) can be used to determine these quantities.
By applying a change of variables to Eq. (8), i.e. ′ = − , the energy loss equation can be transformed to an energy gain equation. This has the advantage of being much easier to solve than the "final value problem" in Eq. (8). We can reorganise Eq. (8) into an initial value problem by setting the initial energy to 0 , In this context 0 is the minimal energy needed for a muon to penetrate the detector, which can be influenced by the detector design. Equation (9) is a well-investigated problem that can be solved by numerous methods. In our work we employ a standard Runge-Kutta integration scheme (see for example Stoer and Bulirsch, 2013), with a step size of 10 cm. As a result, it is now possible to write the cut-off energy in a functional form, where 215 , = ( , , ) . (10) Here (⋅) is the function that returns the Runge-Kutta solution of Eq. (9) for defined thicknesses of materials, , with densities and compositional parameters . Thickness and density are allowed to be vectors, as there may be more than just one material. In this case, the final energy, after the muon has passed through the first segment of materials, is the initial energy for the second segment, etc. In order to speed up the computationsespecially the calculation of the pair production 220 cross-section, which includes two nested integrations. In particular, a log-log table of muon energy vs. radiation loss parameters is produced, from which the b-values, see Eq. (7), can be interpolated. We justify this approach because the radiative losses are almost linear in a log-log plot, as can be seen in Fig. 33.1 of Tanabashi et al. (2018, p.447)

Rock model
Equation (10) shows that for the calculation of the cut-off energy two types of material parameters are required, which are the material density and its average composition . The pre-tabulated values from Groom et al. (2001), however, include only pure elements as well as certain compounds. To extract the relevant parameters in a geological setting, a realistic rock model is needed. In an earlier work  we have shown how an integrated rock model can be 230 constructed and how the physical parameters for a realistic rock can be retrieved. In the present work we generally use the same approach, apart from a few aspects. First, we measured the average material density directly in the laboratory, using various techniques which are explained in detail in Appendix B1. Second, in order to be able to compare the results of this study with the ones in the Nishiyama et al. (2017) publication, we consider a rock composition that corresponds to a densitymodified standard rock. This is applicable, as the rock in the study of Nishiyama et al. (2017) is mostly of granitic/gneissic 235 origin, with thicknesses rarely larger than 200 m, with the consequence that the differences are negligible. However, as the inclusion of compositional data is a planned feature for a future version of our code, we decided to include the theoretical treatment in this work. Hence, all equations are tailored to include the statistical description of such data. Compositional data for whole rock samples, which can be scaled to outcrop scale, are usually presented in one of two forms, the first being the measurement results of X-Ray Diffraction (XRD). This kind of data yields the mineral phases within a rock. Unfortunately, 240 XRD is a rather time-consuming method. This is the reason why in muon tomographic experiments researchers often resort to a bulk chemical analysis of the rock, which is the second form of compositional data. This type of data is usually the output of dedicated X-Ray Fluorescence (XRF) measurements, describing the bulk rock composition by major oxide fractions. We note here that by the absence of information on the spatial distribution of mineral phases within a rock, we implicitly infer a homogeneous mixture of elements within the rock itself, which is thus different from our previous work 245 . From a particle physics perspective this does not pose a real problem as the difference to a mixture of minerals is rather small. Nevertheless, we lose the power to obtain meaningful inferences that could be drawn if compositional information is being considered. As the present work aims to infer positions and uses material parameters as constraints, we can accept this drawback. Details on how compositional parameters are derived from XRF measurements, including an example, can be found in Appendix B2, and an explanation of the related influence on the energy loss equation 250 in Appendix B3.

Spatial models of detectors and materials
In addition to the above explained physical models, we may also utilise available spatial data for our purposes. In this context, the use of a digital elevation model (DEM) of the surface allows the visualisation of the position of the detectors relative to the surface, as well as the spatial extent of the bins. Additionally, it allows us to determine the location where 255 these bins intersect with the topographic surface. As a first deliverable, we can draw conclusions on which bins consist of how many parameters. For example, if we know that the detector is located underground and that there is ice at the surface, https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License.
we can already infer the existence of at least 2 materials (rock and ice). For this purpose, we wrote the script "modelbuilder.py", which allows the user to attach geographic and physical information to the selected bins. This process of building a coherent geophysical model is needed for the subsequent employment of the inversion algorithm to process all the 260 data.

The inverse model: A Bayesian perspective
As stated in the Introduction, we solve the inversion by using Bayesian methods. This needs an explanation as to why we chose this way and not another. First, the equations in Ch. 2 enable us to calculate a synthetic dataset for fixed parameter values. There, one can see that the governing equations constitute a nonlinear relationship between parameter values and 265 measured data. Despite this being of no particular interest in the forward model, the estimation of the parameters from measured data is rendered more complicated. Among muon tomographers, linearised versions have been extensively used with deterministic approaches (e.g. Nishiyama et al., 2014;Rosas-Carbajal et al., 2017), which are successfully applicable when the density or the intersection boundaries are the only variables. When deterministic approaches are viable, they efficiently produce good results. Descent algorithms or, generally speaking, locally optimising algorithms, offer a valid 270 alternative, as they could cope with the nonlinearity of the forward model, while including all desired parameters. Even though these algorithms suffer from possible non-uniqueness solutions (i.e. the solution depends heavily on the starting model, possibly yielding multiple solutions), the main problem is the calculation of the derivatives of the forward model with respect to the parameter values. The analytical calculation of the derivatives is enormously tedious because the cut-off energy results from a numerical solver of a differential equation, as can be seen in Eq. (9) & (10). Unfortunately, numerical 275 derivatives do not produce better results, because they might easily produce artefacts, which are hard to track down. This is especially true if the derivative has to be taken from a numerical result, which is always slightly noisy. In that case the differentiation amplifies the "noise", resulting in unreliable gradient estimates. A good overview over deterministic inversion methods can be found in Tarantola (2005).
The reasons stated above and our goal to include as much information on the parameters as possible nudges us towards 280 employing probabilistic methods. Those approaches are also known as Bayesian methods. The main feature that distinguishes them from the deterministic methods described above is the consistent formulation of the equations and additional information in a probabilistic way, i.e. as probability density functions (pdf). This allows us to (i) incorporate, for example, density values that were measured in the lab (including its error), (ii) set bounds on the location of the material interface, or (iii) define a plausible range for the composition of the rock. All these changes act on the pdf of the respective 285 parameter and naturally integrate into the Bayesian inversion. Readers may find the book of Tarantola  The flexibility of being able to include as much information on the parameters as we consider useful comes at the price of having to solve the inversion in a probabilistic way. This can either be done using Bayes' Theorem and solving for the pdfs 290 of the parameters of interest, or if the analytical way is not possible by employing Monte Carlo techniques. As the presence of a numerical solver renders the analytical solution impossible, we resort to the Monte Carlo approaches. In the following sections we guide the reader through the various stages of how such a probabilistic model can be set up, how probabilities may be assigned, and how the inversion can finally be solved.

Probabilistic formulation of the forward model 295
The starting point for a probabilistic formulation is denoted by the equations that were elaborated in Ch. 2. These deterministic equations need to be upgraded into a probabilistic framework, where their attributed model and/or parameter uncertainties are inherently described. In the following paragraphs we describe how each model component can be expressed by a pdf before the entire model is composed at the end of this subchapter. The model is best visualised by a directed acyclic graph (DAG), i.e. see Kjaerulff and Madsen (2008), that depicts which variables enter the calculation at what point. For our 300 muon tomography experiment this is visualised in Fig. 3. In the following the pdfs are denoted with the bold Greek letter to differentiate them from normal parameters.

Muon data
The data in muon tomography experiments are usually count data, i.e. a certain number of measured tracks within a directional bin, which has been collected over a certain exposure time and detector area. As the measured number of muons is always an integer, we may model such data by a Poisson distribution, where denotes the measured number of muons in the i-th bin and is the poisson parameter in the same bin, which can be interpreted as mean and variance of this distribution. Equation (11) may be rewritten in terms of a flux, by where is the muon flux in the i-th bin and Δ = Δ , * Δ * ΔΩ i (13) 320 is the exposure, in which Δ , is the effective total detector area from Eq. (4), Δ is the exposure time and ΔΩ is the solid angle.

Flux model
The next step is to set up a probabilistic model for the muon flux. First, we observe that "flux" is a purely positive parameter, i.e. ∈ [0, ∞). Thus, it is natural to model it by a lognormal probability distribution if estimates of mean and variance are 325 readily available. The uncertainty on the muon flux is generally taken around 15% of the mean value. As it is possible, by Eq. (2), to calculate a flux for a given cut-off energy, which we interpret as the mean of the non-logarithmic values, the parameters of the lognormal distribution (i.e. , 2 ) may be expressed by ) 2 ) = ln(1.0225) and 330

Energy loss model
The energy loss model has multiple sources of errors that have to be taken into account. Most notably, the relative errors on 335 the different physical cross-sections are given by Groom et al. (2001) as = 6 %, = 1 %, = 5 %, ℎ = 30 % . As it is not clearly stated as to what this error relates to, i.e. one or more standard deviations, we interpret an error like = 6 % as: "within a factor of 1.06", which can be written as where ∈ = {ionisation, bremsstrahlung, pair-production, photonuclear}. Dividing this inequality by and taking the 340 logarithm yields Thus, we may attribute a Gaussian pdf in the log-space for a "log-correction factor, " by setting its mean to zero and its standard deviation to ln(1 + ), i.e.
With a change of variables, using the Jacobian rule as explained in Tarantola (2005), we get the lognormal pdf for the correction factor. The pdf for the cross-section uncertainty ( ) can now be written as a product of the four different pdfs described by Eq. (20) as the errors of the physical cross-sections are stochastically independent from each other.
The calculated energy loss depends also on the material parameters and subsequently on their uncertainties. However, these will be explained in detail in Ch. 3.1.4. A last error enters by the numerical solution of the ordinary differential equation, Eq.
(9). We decided not to model this error, as its magnitude is directly controlled by the user (by setting a small enough step length in the Runge-Kutta algorithm) and thus can be made arbitrarily small. Lastly, we assume that all the errors in the 355 energy loss model are explained by uncertainties in the cross sections as well as in the material parameters. Although this assumption is rather strong, since it excludes the possibility of a wrong model, we argue that this approach works as long as https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. the variation in these parameters can explain the variation in the calculated cut-off energy. If this requirement is met, we may model the pdf for the energy loss model as a delta function, Due to the presence of the delta function in Eq. (22), this integral is solved analytically resulting in where the parameters are given by (27) Please note that M( , , , ) = ( ( , , , )) describes the combined parts of the forward model that include the energy loss and the integrated flux calculation, which is basically a composition of functions.

Rock model 375
The density model can take different forms of probability densities (see Appendix B1), such as normal, lognormal, uniform, etc. For either form, it is possible to describe it by a generic function ( ), which is a short version for a multidimensional pdf, i.e. ( , ) if the i-th cone is known to consist of two segments with two specific densities. Equivalently, the pdf for the composition (see Appendix B2) is either fixed or a multidimensional Gaussian distribution in the space of log-ratios.
Thus ( ) can be split up to ( , ), like in the example above. Generally, we may assume that in our problem j 380 different materials exist. We note here that the description of the composition is already probabilistic. However, the inversion in that case is not functional and works only with the mean values of the multidimensional Gaussian. The support for composition inversion is planned for a future version of the code. The situation for the thicknesses of the segments, ( ), within the i-th cone presents itself in a similar way as for the compositions (e.g., Fig. 1). As the total material thickness is known (detector position and digital elevation models are 385 given), the sub-space containing the thickness parameter is endowed with the same mathematical structure as the one containing the composition parameter (i.e. one sum constraint), if the cone consists of more than just one segment. One can therefore safely assume that the thickness parameters live in a log-ratio space, within which we a-priori possess no information about the parameters. Thus, we attribute the thickness parameters a multidimensional uniform distribution within the log-ratio space. 390

The Joint probability density function
With the help of the DAG, introduced in Fig. 3, it is now straightforward to factorise the joint probability distribution for the whole problem, as their structure is equal. This results in or equivalently (and this will also be of a much better use later on) the log joint pdf 395 where the prefix " " denotes the logarithm of the pdf. This has the benefit of reducing the size of numbers that the code has to cope with. Moreover, many computational statistics packages already have this feature included, which renders it easy to use.
Equation (28) depicts the full joint pdf. However, the relations between the parameters, as shown by the DAG (see Fig. 3), 400 classify this model as a hierarchical model (Betancourt and Girolami, 2013). The key characteristic of such models is their tree-like parameter structure, i.e. the measured number of muons is related to the thickness or the density of the material by the flux parameter only, which "relays" the information. A central problem of such models is the presence of a hierarchical "funnel" (see Fig. 2 & 3 of Betancourt and Girolami, 2013), which renders it very difficult for standard MCMC methods to adequately sample the model space. In high-dimensional parameter spaces this problem exacerbates even more. 405 Our aim to provide a simple and easy-to-use program somewhat contradicts this necessity of a sophisticated method (which inevitably requires the user to possess a strong statistical background). As the main problem is the rising number of parameters, it should be possible to mend the joint pdf by imposing thought-out simplifications.
We first get rid of the flux parameter, as for our problem it merely is a nuisance parameter. This is an official term for a parameter in the inversion which is of no particular interest but still has to be accounted for. Here specifically, we integrate This effectively reduces to a problem where the new likelihoods have to be calculated (as is given) or fully, where and are given by Eqs. (26) and (27), respectively. This integral is not solvable analytically but can be evaluated 420 by numerical integration schemes. The likelihood has a maximum when the Poisson and the log-normal pdfs fully overlap.
Interestingly, this directly shows the trade-off between the flux model and the data uncertainty. Usually, we want to measure enough muons so that the statistical counting error is smaller than the systematic uncertainty of the flux model (i.e. the width of the Poisson pdf is smaller than the width of the log-normal pdf). This can be controlled directly by the exposure of the experiment, via a larger detector area, a coarser binning, or a longer exposure time. 425 This marginalisation roughly halves the number of parameters, but there is still another simplification, which we may use.
Many muon tomography applications deal with a two-material problem, while there may also be measurement directions where only one material is present. If we conceptually split those two problems and solve them independently, it is possible to further reduce the number of simultaneously modelled parameters. In the study of Nishyiama et al. (2017), the results of which we will use later, these two cases encompass bins where we measured only rock and others where we know there is 430 ice and rock. The joint pdf for rock bins subsequently is which leaves the problem effectively with only a handful of parameters. Solving Eq. (33), for the rock density we retrieve ̃( ), the posterior marginal pdf for the rock density. We refer the reader to Ch. 3.2 for the details of how to solve this inverse problem. Theoretically we could also retrieve ̃( ), but this would require the detector to be positioned within the 435 glacier (see chapter 5.2 and Fig. 1 for a  For the second problem, we can interpret ̃( ) as the new prior pdf for the rock density. At this point we employ one last simplification by assuming that the parameters between different cones are independent form each other. This is a rather strong presumption, which must be justified. The main mathematical problem lies in consideration of the hierarchical nature 440 of the density parameter, which is the same for each cone and therefore not independent in different cones. We, however, https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License.
argue that in cones with two materials, there are more parameters than in bins with only rock, such that we may expect the posterior pdf of the rock density of these second kind of models to be less informative than the posterior rock density pdf of Eq. (33). This, in turn, means that the posterior rock density pdf of the two-material model largely equals the prior one if we select the posterior of the first kind of models as the prior of the second kind of models. The same is valid for the 445 composition and the cross-section error parameters . As long as this assumption is valid, we may decompose the joint pdf into independent joint pdfs for each cone

Solution to the inverse problem
Usually in Bayesian inference, the goal is to calculate the posterior pdf, given the measured data, i.e. the quantity This can be interpreted as the inferences one may draw on the parameters in a model given measured data. The denominator 455 on the right-hand side of Eq. (35), also called the "evidence", can be rewritten as the data marginal of the posterior, i.e.
This is the starting point of Monte Carlo Markov Chain (MCMC) methods.

The Metropolis-Hastings algorithm
The basic MCMC algorithm, which we also use in this study, is the Metropolis-Hastings (MH) algorithm (Hastings, 1970;465 Metropolis et al., 1953), which allows for the sampling of the joint pdf to obtain a quantitative sample. We note, however, that many different MCMC algorithms exist for various purposes and that the MH has no special status except for being comparatively simple to use and implement. An example of another MCMC algorithm in muon tomography can be found in Lesparre et al. (2017). The authors used a simulated annealing technique on the posterior pdf in order to extract the https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. maximum a posteriori (MAP) model. As every simulated annealing algorithm has some type of MH-algorithm at its core, we 470 directly use the MH-algorithm in its original form such that we not only retrieve a point estimate but a pdf for the posterior parameter distribution. The algorithm is explained in detail by Gelman (2014), such that we only provide a short pseudocode description.

Algorithm 1 (Metropolis-Hastings):
(1) Draw a starting model, ⃗⃗ 0 = ( ,0 , 0 , 0 , x ⃗ 0 ), by drawing ,0 , 0 , 0 , x ⃗ 0 from their respective prior pdfs and 475 determine the log-pdf value of this model The advantage of this algorithm, compared to a "normal" sampling, lies in its efficiency. It is often not possible, or even reasonable, to probe the whole model space, as the largest part of the model space is "empty", where the pdf-value of the posterior is uninterestingly small. The fact that regions of high probability are scarce, and this becomes worse in high dimensional model spaces, is known as the "curse of dimensionality" (Bellman, 2016). MCMC algorithms (including the 490 here presented MH-algorithm) allows one to focus on regions of high probability, and therefore we are able to construct a reliable and representative sample of the posterior pdf. We again refer to Gelman (2014) for a discussion of why the MHalgorithm converges to the correct distribution and why we may use samples that were gained this way to estimate the posterior probability density.

Assessing convergence, mixing, and retrieving the samples 495
The above stated advantages, however, come at a price. First and foremost, we must ensure that the algorithm advances fast enough, but not too fast, through the model space. This is mainly controlled by the proposal distribution (0, 2 Σ), which is taken to be a multivariate Gaussian distribution. Ideally, the covariance matrix of the proposal distribution Σ is equal to the covariance structure of the posterior pdf. We acknowledge that at the start of the algorithm one generally has no idea what this looks like, but we assume that a combination of the prior variances is a reasonable starting point. After a certain number 500 https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. of steps, it is possible to approximate the covariance matrix of the proposal distribution with the samples taken up to at this point.
A second crucial point is the presence of a warm-up period. The starting point, which usually lies in a region of high prior probability, does not necessarily lie in a region of high posterior probability. The time it takes to move from the latter to the former is exactly this warm-up. This can usually be visualised by a trace plot, e.g. Fig. 4, in which the value of a parameter is 505 plotted against the number of iterations. After this warm-up phase, the algorithm can be run in operational mode and "true" samples can be collected. As in a Markov Chain the actual sample is dependent on the last one, we need a criterion to argue that the samples created in that way really represent "independent" samples. Qualitatively, we may say that if the Markov Chain forgets the past samples fast enough, then we may sooner treat them as independent from each other. Gelman (2014) suggests that in order to 515 assess this quantitatively, multiple MH-chains could be run in parallel and statistical quantities within and between each chain are analysed. For a detailed discussion thereof, we refer the reader to Appendix C.
Once a satisfying number of samples has been drawn from the posterior pdf, a marginalisation of the nuisance parameters can be done by looking at the parameters of interest only. These samples may then be treated like counts in a histogram, i.e. distributional estimates, or simply the interesting statistical moments, such as mean and variance, can be obtained. 520 https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License.

Construction of the bedrock-ice interface
The main analysis programme allows us to export all parameters either as a full chain dataset, where every single draw is recorded, or as a statistical summary (i.e. mean and variance). Both are then converted to point data, i.e. (x, y, z)data. For the subsequent construction of the interface between rock and ice (see section 5.2 and Fig. 1 for the presentation of the test experiment) we only need the full-chain point-data. In the present study we restrict ourselves to a probabilistic description 525 until the bedrock positions within a cone. It would also be possible to treat the bedrock construction within a Bayesian framework, however this would go beyond the scope of this study and is therefore left for a future adaption of the code.
Nonetheless, in order to construct a surface, we rely on deterministic methods, which are explained in detail in what follows.

Interpolation to a grid
The "modelviewer.py" routine is able to read datasets from different detectors (which are saved as JSON-files) and computes 530 for each cone the statistic, which the user is interested in (see "sigma" entry in program). Thus, it is possible to use the mean or, for example, the +1 position of each cone. From here onwards this point cloud is named and contains one interface position (x, y & z coordinates) per cone. These are shown as triangles (▲) in Fig. 5.
As a second step, the programme interpolates this point cloud in a bilinear way to a rectangular grid with a user specified cell size, Δ . This grid can be described by a matrix ∈ ℝ × , where and are the number of rows and columns (i.e. the 535 number of y-and x-cells, needed to cover the whole grid). The procedure is similar to the bilinear interpolation of Lagrangian markers (that carry a physical property) to a (fixed) Eulerian grid in geodynamical modelling (see Gerya, 2010, p. 116), with the difference that our physical property is the height of the ice-bedrock interface (Fig. 1).
We could also have fitted a surface through the resulting point cloud. However, by formulating this surface as a matrix we gain access to the whole machinery of linear algebra. Moreover, can directly be interpreted as a rasterised DEM, which 540 can be easily loaded and visualised in any GIS software. Thus, from a modular design perspective we think that the matrix formulation has more advantages than drawbacks. The bilinear interpolation is shown in more detail in Fig. 5.
where the weights, , are given by 550

Damping & Smoothing
The concept of damping usually revolves around the idea to force parameters to a certain value (e.g. in deterministic inversion by introducing a penalty term in the misfit function for deviations from that value). From a Bayesian viewpoint this would be accomplished by setting the prior mean to a specific value. In our code we implemented this idea by allowing the 555 user to read a DEM and a "damping weight" to the code (see "fixed length group" in code). The programme effectively computes a weighted average between the bedrock positions within the cones and a user defined DEM. The higher the chosen damping weight, the more the resulting interface will match the DEM, when pixels overlap. The matrix formulation also enables us to use a further data processing technique without much tinkering. As geophysical data are often quite noisy, a standard procedure in nearly every geophysical inversion is a smoothing constraint. This 560 effectively introduces a correlation between the parameters and forces them to be similar to each other. From a Bayesian perspective, we could have achieved this correlation by defining a prior covariance matrix of the thickness parameters, such that neighbouring cones should have similar thicknesses (which makes sense as we do expect the bedrock-ice interface to be relatively continuous; Fig. 1). As we work with independent cones in this study, we leave the exploration of this aspect open for a future study. Nevertheless, we offer the possibility in our code to use a smoothing on the final interpolated grid. This is 565 achieved by a convolution of a smoothing kernel, (see Appendix D for details), with the surface matrix , which results in a smoothed surface matrix = * . (40) Please note that the * operator in Eq. (40) denotes a convolution. In index notation the advantage of the linear algebra formalism becomes clear, as can be expressed by 570 The user is free to choose the number of neighbouring pixels, , across which the programme performs a smoothing. As a related matrix we use an approximation to a Gaussian kernel, which corresponds to a Gaussian blur in image processing.
Whereas "smoothing" is a general term used in the geophysical community where such a process is used with an aim to force a correlation of parameters, in our case where the parameters describe a surface (Fig. 1), the convolution effectively 575 smooths the surface, i.e. removes small-scale variations.
Finally, we added a checkbox to our code to allow it to change the order of the damping and smoothing operations.
Sometimes when a strong damping is necessary, this may result in rather unsmooth features at DEM boundaries, such that it makes sense to perform a smoothing only afterwards.

Main modules of SMAUG 580
Our toolbox, SMAUG, contains several subprograms, which are executed separately. This allows the user to inspect This subroutine allows the user to create their own material that will be used in the subsequent model builder. The user may choose a density (either from data or directly insert mean and standard deviations) and a composition (also either from data or from the list of Groom et al., 2001).

DATA_BINNING.py
As the name suggests, this subroutine is used to spatially bin the recorded track data. The bin data (i.e. the output hereof) is then fed to the model builder.

MODEL_BUILDER.py 595
The model builder takes the bin data and the materials as inputs and allows the user, with help of DEMs, to allocate data and materials to certain cones. This is basically the spatial setup of the model. The resulting model file is then provided to the inversion code.

INVERSION.py 600
This is the main module in SMAUG, providing the functionality to perform a MCMC algorithm on the probabilistic model created with MODEL_BUILDER.py. It also includes several analysis tools to assess MCMC performance.

MODEL_VIEWER.py
The model viewer allows us to visualise the interface results, obtained and exported by INVERSION.py. It also has the 605 functionality to dampen and smooth the resulting surfaces.

Model Verification
In this section we present examples of how the model can be employed, what it calculates and how the output is structured. 610 We proceed by verifying, in a first step, that the physical models employed in this work yield results which are numerically consistent with the results of calculations from other studies. We will compare our results with reported values from literature in Ch. 5.1. We do this because we do not change the parameters of the flux model (except the height scaling, which has been verified by Nishiyama et al. (2017), and since the energy loss calculations is based on equations that stem from different frameworks. 615 In a second step, we will benchmark the results obtained by this code from real data against previously published results. For this purpose, we will re-analyse the raw data from Nishiyama et al. (2017). This is a study (see also Fig. 1) that was conducted in the Central Swiss Alps in a railway tunnel that featured a glacier (part of the Great Aletsch glacier) above. Our goal there was to estimate the thickness of the overlying glacier and thus the subsurface structure of the ice-bedrock interface. The respective calculation and discussion thereof are presented in Ch. 5.2. 620

Verification of energy loss calculations
The energy loss model that we use in our code generally reproduces the literature values well (below 1%) across the different energy loss processes and relevant energies. In Fig. 6 we present the energy loss calculations for each energy loss process (i.e. ionisation, bremsstrahlung, pair-production and photonuclear interactions) across energies from 10 MeV to 100 TeV for silicon. 625  We note that the energy losses by ionisation are reproduced very well over the entire energy range. We also note that the relative error on the radiative energy losses is rather large below 10 GeV. This does, however, not introduce a major bias, because below this energy, radiative energy losses are negligible compared to ionisation losses, as can be seen in Fig. 6. 640 Furthermore, the related errors are in an acceptable range at the energy level at which radiative losses begin to become noticeable (i.e. around 100 GeV). This can be seen in Fig. 7, in the sense that the total relative error remains well bounded within 0.5 %. In the ionisation domain (i.e. below 100 GeV) the total relative error is dominated by the ionisation relative error, whereas above this energy level the relative errors on radiative losses start to prevail. A close-up of this energy range is given in Fig. 8. 645 There are different sources and circumstances that contribute to the error in the different energy losses processes. The scatter of the relative ionisation-loss error around 0 with a rather small deviation can be viewed as simple rounding errors. The errors on the radiative processes, however, seem to be of a more systematic nature. We explain this behaviour through a https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. different numerical integration scheme in Eq. (7), which tends to systematically under-/overestimate the true value, especially when the integrand comprises exponential functions. Whereas we used a Double Exponential Integration scheme 650 (see Takahasi & Mori, 1974), the integration scheme from Groom et al. (2001) is not discernible. However, as the relative errors on the processes of energy loss remain well within the theoretical uncertainties, (see Ch. 3.1.3), we consider, that our calculation accurately reproduces the literature values for elements. The above calculations were performed for pure silicon. The respective figures for other four important elements in the Earth's crust (Al, Fe, Ca & O) can be found in Appendix E. Those elements are, however, not representative for any real material encountered in geological applications. For this reason, we compiled the same computations for four selected, geologically important compounds (SiO2, CaCO3, Standard Rock, ice) that are also shown in Appendix E. We summarise, 660 that with the exception of Standard Rock, all calculations yield results that are similar to the silicon calculation above. The discrepancy for Standard Rock stems from its inconsistent definition, with respect to the different parameters. In particular, the "Standard Rock" according to Lohmann et al. (1985) has an atomic number Z of 11 (i.e. sodium) and an atomic weight A of 22, which yield the characteristic parameter values of 〈 ⁄ 〉 = 0.5 and 〈 2 ⁄ 〉 = 5.5 respectively. Note that Groom et al.
(2001) list sodium as the only constituent of a standard rock. However, this material cannot be modelled by any mixture of 665 pure elements, as common sodium consists of one neutron more and thus has a higher atomic weight (i.e. = 23). Consequently, the use of standard sodium would lead to different characteristic parameter values, i.e. 〈 ⁄ 〉 = 0.478 and 〈 2 ⁄ 〉 = 5.263, thus leading to an inconsistency. This is often conveyed by the phrase that standard rock "is not-quitesodium" (Groom et al. 2001, p.203). In order to circumvent this problem, we advocate the exchange of 11 . With this change, standard rock does not need any more special treatment and can be calculated in a way that is consistent to all other compounds. Furthermore, the relative error between the tabulated values and our modified calculation 675 falls in line with the calculations for the other compounds and elements (Figs. 9 and 10).

Verification of the bedrock-ice interface reconstruction
In this part we test the presented reconstruction algorithm on previously published data. For this purpose, we compare our 685 calculations to the ones already published in the study by Nishiyama et al., (2017), where the goal was to measure the interface between the glacier and the rock, in order to determine the spatial distribution of the rock surface (also below the glacier). We could access the Railway Tunnel to install the muon detectors beneath the ice. A situation sketch is shown in Figs. 1 and 11.
The results shown below (Figs. 12 -14) represent the bedrock-ice interface interpolated to an 8-metre grid, which was first 690 damped (weight 8) and then smoothed (2 grid pixel). We assess the goodness of fit according to the three cross-sections (East, Central, West), that are shown in Fig. 10. The crosscuts are nearly perpendicular to the train tunnel and roughly 40 apart from each other. Figures 12 to 14 show the three cross-sections in detail. In every plot, we also indicate the solution from Nishiyama et al. (2017). Please note that we added a systematic error of 2 to the uncertainty planes, as the DEM we are working with has itself an uncertainty of ± 2 . The dash-dotted lines mark thus the most conservative error estimate. 695 Moreover, we highlighted the parts of the cross-section that had either been damped to the bedrock DEM or that have been solely resolved by the measurement (see "damping marker" in Fig. 11).   80 ) can be explained as a smoothing artefact. Towards the end of the profile, the decreasing data coverage becomes evident as the uncertainties rise. This effect can also be seen in the jagged behaviour of the interface curves around 100 to 120 , hinting at the effect where the interpolation has occurred with few data.  of 65 ). Here we used the same DEM and aerial photograph as in the previous study. This means that newer versions might 720 be available that show more bedrock (due to the glacial retreat as a response to global warming). The eastern profile is shown in Fig. 14. One sees that the results from this study are internally consistent. The surface from 730 the previous study plunges down earlier with respect to the surfaces calculated here. This may in fact be a damping effect, as the bedrock-ice interface from Nishiyama et al. (2017) has not been constrained to the bedrock (via damping) and thus plunges down before the damping mark at ~72 . Still, the two surfaces agree within 5 , which we consider as acceptable.

740
All together the performance of the whole workflow, which is shown in this study, produces results, which are similar to the ones published in the previous study (Nishiyama et al., 2017). We use the results of this comparison to validate the base of our code.

Conclusion 745
In this study we have presented a model that allows us to integrate geological information into a muon tomography framework. The inherent problem of parameter estimation has been formulated in a probabilistic way and solved accordingly. The propagation of uncertainties thus occurs automatically within this formalism. We also considered approaches including DAGs or the simplex subspace of compositions which could be helpful to the muon tomography community while tackling their own research. We condensed these approaches in a modular toolbox. This assortment of 750 python programmes allows the user to address the subproblems during the data analysis of a muon tomography experiment.
The programmes are modular in the sense that the user can always access the intermediate results, as the files are mostly in a portable format (JSON). Thus, it is perfectly possible to only use one submodule of the toolbox while working with an own codebase. As every "tool" is embedded in a GUI, the programme is made accessible without the need to first read and consider several thousand code lines. Furthermore, we have shown that the results we obtain with our code are largely in 755 good agreement with an earlier, already published experiment. The small deviations may be attributed to data analysis subtleties.
We would like to stress that this work is merely a foundation upon which many extensions can be built when it is used in other applications as well. Future content might, for example, include a realistic treatment of multiple scattering and the inclusion of compositional uncertainties in the inversion, for which we laid out the basis in this study.
where is the zenith angle of the incident muon. It is important to note that the parameter values in Eq. (A1) are changed to 0 = 0.2455, 1 = 1.288, 2 = −0.2555, 3 = 0.0209 and = 0.00253. In order to include height above sea level as 770 an additional parameter, Hebbeker and Timmermans (2002) proposed to model the altitude dependency as an exponential decay, which modifies Eq. (A2) into The scaling height, ℎ 0 , is usually to be taken as ℎ 0 = 4900 + 750 −1 * , where , is the momentum of the incident muon in * −1 . However, as this formula is only valid up to an altitude of 1000 m above sea level, Nishiyama 775 et al. (2017) adapted it to ℎ 0 = 3400 + 1100 * * −1 * * cos( ). This was done in order to fit the energy spectrum up to 4000 m above sea level. This formula is now valid for momenta above 3 * −1 , zenith angles between 0° and 70° and an altitude below 4000 above sea level.

B1 -Density model 780
The density distribution of a lithology can be determined through various methods. In our work, we constructed a density model by analysing various rock samples from our study area in the laboratory. Two experimental setups were employed to gain insight into the grain, skeletal as well as the bulk density of the rocks. Grain and skeletal density were measured by means of the AccuPyc 1340 He-pycnometer, which is a standardised method that yields information on the volume. Bulk density values were then determined based on Archimedes' principle, where paraffin coated samples were suspended into 785 water (ASTM C914-09, 2015; Blake and Hartge, 1986).
Every sample = 1, … , (usually the size of a normal hand sample) has been split up into smaller subsamples = 1, … , , that were measured. The bulk density of the i-th subsample can be calculated by , where 2 , , denote the density of water, paraffin and the thread that was used to dip the sample into the liquid, 790 respectively. , , , , , , , describe the mass of the sample, the paraffin coating, the thread and the apparent mass of all three components suspended in water. , , , can then be simply obtained through as , , , denote the mass of the sample including thread and paraffin coating on one hand and , , only the mass of the sample and the thread on the other hand. Further, the mass of the thread is given by 795 The maximal precision of the reading is estimated at ±5 * 10 −5 , and the commonly ignored effects regarding buoyancy in air has been estimated to introduce an error on the order of ±2 * 10 −4 . This error has been attributed to all direct mass measurements. Moreover, because small pieces of material may detach from the sample upon attaching the thread to the sample and during the paraffin coating, we set an error of ±2 * 10 −2 to all measurement results. The variables in Eq. (B1) 800 are strictly positive values. Following Tarantola (2005) we model these "Jeffreys parameters" by lognormal distributions, as they inherently satisfy the positivity constraint. Because Eq. (B1) does not simply allow a standard uncertainty propagation, the script "subsample_analysis.py" performs a Monte Carlo simulation for each subsample and attributes a final lognormal probability density function to the resulting histogram. Figure B1 illustrates such an example, where the calculation has been performed for subsample JT-20-1. 805 We have found 10'000 draws per subsample to be sufficient to retrieve a solid final distribution. However, this parameter can easily be changed in the script, depending on the user's preference of precision/speed. From this point onwards we may work with a Gaussian distribution as Fig. B1 assures us that a normal pdf describes the results of the Monte Carlo simulation rather well.  The determination of the grain and skeletal densities is simpler than the bulk density measurements because the corresponding method consists of a mass and a volume measurement, respectively. The density formula reads then simply 815 / , = . (B4) The question remains as to how it is possible to construct a pdf that represents the knowledge about the whole lithology.
There are two possible methods that can be readily employed at this point. The first, following largely Tarantola (2005), performs a so-called disjunction of the pdfs that corresponds to an averaging of all subsample pdfs. As Vermeesch (2012) points out, even though this might seem a "sensible strategy at first glance", there might be some problems with this method. 820 The main problem lies in the small error on the subsamples, such that the variation between different subsamples may be larger than their attributed errors. This would not be a problem if enough subsamples could be measured, such that the resulting lithology pdf might be sampled correctly. On the other hand when one is faced with a situation where data is rather scarce, then the approach of Tarantola (2005) would result in a rather spikey pdf that would be hard to handle. For this reason we adopted the methodology of Vermeesch (2012) where the lithology pdf is estimated by a kernel density 825 https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. estimation. The main difference lies in the fixed "bandwidth" of the subsample distributions. We refer to Vermeesch (2012) for more details and an in-depth discussion on this problem.
The kernel density estimation has the advantage, that only the mean values of the subsamples have to be processed as the bandwidth is determined from the spread of the subsample means. Following the methodology of Vermeesch (2012) we end up with a pdf like the one visualised in Fig. B2. We could at this point use the kernel density estimated pdf for further 830 calculations. However, for simplicity we approximate the kde with a normal distribution and intend to add support for the kde in a later code version.

B2 -Composition model
We have seen in Eq. (9) that the material density parameter enters the energy loss calculations rather directly. Contrariwise, the compositional model affects the energy loss equations much more subtly through the average { ⁄ } and { 2 ⁄ } values and mean excitation energies that need to be calculated for the entire lithology. Likewise, information on the weight 840 percentages of the main elements within the rock is required for the quantification of the radiation loss term.
Although a modal mineral analysis (e.g. the quantitative determination of mineral volumes) is preferable and can be treated according to Lechmann et al. (2018), its execution is a rather time-consuming effort. This is the reason why compositional data in muon tomography experiments predominantly consist of XRF-data, which show the abundance of major oxides within the rock. We describe here a method to incorporate such type of information in a probabilistic way thereby following 845 Aitchison (1986). Compositional data are usually available in the form of Table B1, which presents an excerpt of four samples for illustration purposes. We refer to the excel sheet in the supplementary material of the present work for the full data. There are several challenges to this kind of data. First, the parameters (i.e. the oxide percentages) can take a value between 0 and 1. This means that normal as well as lognormal distributions are not suitable to describe these parameters. Second, the requirement that the sum of all parameters has to ideally equal 1 poses a constraint on this parameter space, which 855 effectively reduces the number of independent parameters by one. Third, due to measurement uncertainties, this sum is never exactly one. Spaces, which have this unit sum condition can be viewed as a simplex, e.g. if we had three compositional parameters, the simplex would be a 2-dimensional surface (i.e. a subspace) in this 3-dimensional parameter space. The last issue, of not summing up exactly to 1, can be remedied by projecting each sample dataset back to the simplex (Aitchison, 1986, p. 257-860 261). This works only if the measurement imprecisions are not too large, which works well for the examples in Table B1.
With respect to the energy loss calculation, it is preferable to decompose the oxides into elements, which can be done by following formula where and denote the molar mass mass of the i-th element and the j-th oxide, is the j-th datum in the column and 865 is the number of atoms of the i-th element within the j-th oxide. The two transformations are visualised in Table B2.
https://doi.org/10.5194/gmd-2021-342 Preprint. Discussion started: 2 November 2021 c Author(s) 2021. CC BY 4.0 License. In order for the data to be in a statistically convenient form, Aitchison (1986) suggests to further transform the data in Table   B2 by first forming a ratio with an arbitrary element (in the list) and then taking the logarithm. For the exemplary dataset this is shown in Table B3. The rationale behind this transformation is as follows. The division by an arbitrarily present element effectively transforms the space into an N-1-dimensional open space, where the parameters (i.e. ratios) may have values between 0 and ∞. The subsequent application of the logarithm further changes the space, such that the new parameters can have values between −∞ and ∞. This results in so-called log-ratios, which should ideally be following a multivariate normal distribution. As a consequence, we can calculate the mean log-ratio vector across all samples as well as its corresponding covariance matrix, 880 which completely describes the multivariate normal distribution. In addition to these statistical parameters, the script "compo_analysis.py" outputs a graph that plots for all samples an order statistic, , (see Aitchison, 1986). This enables us to visualise how different the data is from a multivariate normal distribution. If equal, they should fall on the red line, shown in  Table B3 is only an excerpt). Each subplot checks for marginal normality. Oxygen is the denominator variable (arbitrarily chosen) and does thus not appear in the plot. for the denominator element (here oxygen). In Eqs. (B6) and (B7) the denote the log-ratios from Table B3 and is the total number of elements (in Table B2).

B3 -Energy loss equation for rocks
As stated in Eq. (7)  The radiation loss term, however, must be calculated as a weighted radiation energy loss over all elements. This means that the average can be written in a rather concise form, . (B10)

Appendix D -Construction of the smoothing kernel
As stated in the main text, the user specifies the number of neighbouring pixels s to smooth over. The main idea is to construct a roughly Gaussian smoothing kernel by approximating it with a binomial distribution. With help of the binomial 940 coefficient we can construct a vector of weights with = (2 * + 1) entries. The weight vector is then given by = 1 2 2 * ( − 1 ) , with ∈ {0, . . , − 1}. It is now possible to create a matrix by forming the dyadic product of ⃗⃗ with itself, i.e.
As an example, we show how a smoothing kernel that smooths over two neighbouring pixels (i.e. = 2) is constructed. This is incidentally also the smoothing kernel we used to construct our ice-bedrock interface. The weight vector in this case is given by