{'fraction_non_alphanumeric': 0.036529559582065864, 'fraction_numerical': 0.013807697394524534, 'mean_word_length': 4.576106194690266, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The effect of rock composition on muon tomography measurements

In recent years, the use of radiographic inspection with cosmic-ray muons has spread into multiple research and industrial fields. This technique is based on the highpenetration power of cosmogenic muons. Specifically, it allows the resolution of internal density structures of largescale geological objects through precise measurements of the muon absorption rate. So far, in many previous works, this muon absorption rate has been considered to depend solely on the density of traversed material (under the assumption of a standard rock) but the variation in chemical composition has not been taken seriously into account. However, from our experience with muon tomography in Alpine environments, we find that this assumption causes a substantial bias in the muon flux calculation, particularly where the target consists of high {Z2/A} rocks (like basalts and limestones) and where the material thickness exceeds 300 m. In this paper, we derive an energy loss equation for different minerals and we additionally derive a related equation for mineral assemblages that can be used for any rock type on which mineralogical data are available. Thus, for muon tomography experiments in which high {Z2/A} rock thicknesses can be expected, it is advisable to plan an accompanying geological field campaign to determine a realistic rock model.


Introduction
The discovery of the muon (Neddermeyer and Anderson, 1937) entailed experiments to characterise its propagation through different materials. The fact that muons lose energy proportionally to the mass density of the traversed matter (see Olive et al., 2014) inspired the idea of using their attenuation to retrieve information on the traversed material. This was first done by George (1955) for the estimation of the overburden upon building of a tunnel, and then later by Alvarez et al. (1970) to search for hidden chambers in the pyramids in Giza (Egypt). In a related study, Fujii et al. (2013) employed this technology to locate the reactor of a nuclear power plant. Recently, Morishima et al. (2017) successfully accomplished quest of Alvarez's team in the Egyptian pyramids.
Besides these applications, which have mainly been designed for archaeological and civil engineering purposes, scientists have begun to deploy particle detectors to investigate and map geological structures. In recent years, this has been done for various volcanoes in Japan (Nishiyama et al., 2014;Tanaka et al., 2005Tanaka et al., , 2014, including the Shinmoedake volcano (Kusagaya and Tanaka, 2015), the lava dome at Unzen (Tanaka, 2016) and most recently the Sakurajima volcano (Oláh et al., 2018). Further experiments have been conducted in the Caribbean, in France (Ambrosino et al., 2015;Jourde et al., 2013Jourde et al., , 2015Lesparre et al., 2012;Marteau et al., 2015) and in Italy on Etna (Lo Presti et al., 2018) and Stromboli (Tioukov et al., 2017). Recently, Barnaföldi et al. (2012)

used
A. Lechmann et al.: The effect of rock composition on muon tomography measurements this technology to examine karstic caves in the Hungarian mountains. Our group is presently carrying out an experimental campaign in the central Swiss Alps for the purpose of imaging glacier-bedrock interfaces .
Inferences about subsurface structures from observed muon flux (i.e. the number of recorded muons normalised by the exposure time and the detector acceptance) necessitate a comparison of the measurement data with muon flux simulations for structures with various densities. Such a simulation consists of a cosmic-ray muon energy spectrum model and a subsequent transportation of these muons through matter. The former describes the abundance of cosmic-ray muons for different energies and zenith angles at the surface of the Earth. This has been well documented in literature (see, for example, Lesparre et al., 2010). The differences between models and experimental data, and hence the systematic model uncertainty, can be as large as 15 % for vertical muons (Hebbeker and Timmermans, 2002). On the other hand, the attenuation of the muon flux is assumed to depend only on the density of the traversed material. In this context, however, potential effects of its chemical composition have not been taken into account specifically. Instead, previous works employ a certain representative rock, so-called "standard rock", for which the rate of muon energy loss has been tabulated (e.g. Groom et al., 2001).
The origin of this peculiar rock type can be traced back to Hayman et al. (1963), Miyake et al. (1964), Mandò and Ronchi (1952) and George (1952), who gave slightly different definitions of its physical parameters (mass density ρ, atomic weight A and atomic number Z). A comprehensive compilation thereof can be found in Table 1 of Higashi et al. (1966). Various corrections to the energy loss equation were then added in the framework of follow-up studies, which particularly include a density effect correction (see, for example, Sternheimer et al., 1984). Richard-Serre (1971) listed data relevant for muon attenuation for (i) soil from the CERN (European Organization for Nuclear Research) premises near Geneva (Switzerland), (ii) molasse-type material (e.g. Matter et al., 1980) and (iii) a "rock" that is equivalent to the one from Hayman et al. (1963). These latter authors assigned additional energy loss parameters to this particular rock type, which were similar to those of pure quartz. Lohmann et al. (1985) then adjusted these parameters to energy loss variables for calcium carbonate (i.e. calcite) and gave the standard rock its present shape. In summary, this fictitious material consists of a density of crystalline quartz (i.e. ρ qtz = 2.65 g cm −3 ), a Z and A of 11 and 22, respectively (which is almost sodium), and density effect parameters that have been measured on calcium carbonate. However, when the material's Z and A differ greatly from standard rock parameters as for carbonates, basalts or peridotites, a substantial bias would be introduced to the calculation of the muon flux. Such a situation is easily encountered in geological settings such as the European Alps where igneous intrusions, thrusted and folded sedimentary covers and recent Quaternary deposits are found in close vicinity (e.g. Schmid et al., 1996). Currently, our collaboration is performing a muon tomography experiment in the Jungfrau region, in the central Swiss Alps, aiming at imaging the glacier-bedrock interface Nishiyama et al., 2017). There, we face a variety of lithologies ranging from gneissic to carbonatic rocks that have a thickness larger than 500 m . In this context, it turned out that the analyses based on the standard rock assumption might cause an over-or an underestimation of the bedrock position in the related experiment. Such an uncertainty arising from the chemical composition of the actual rock has to be reduced at least to the level of the statistical uncertainty inherent in the measurement as well as in the systematic uncertainty of the muon energy spectrum model.
To achieve this, we investigate how different rock types potentially influence the results of a muon tomographic experiment. We particularly compare the lithologic effect on simulated data with standard rock data to estimate a systematic error that is solely induced by a too-simplistic assumption on the composition of the bedrock.

Rock types
In this study, we chose 10 different rock types that cover the largest range of natural lithologies, spanning the entire range from igneous to sedimentary rocks. The simplest rocks have a massive fabric in the sense that they do not exhibit any planar or porphyritic texture. Typical lithologies with these characteristics are igneous rocks or massive limestones (not sandstones, as they might have a planar fabric such as laminations and ripples). Exemplary thin sections of granite and limestone are shown in Fig. 1. Note that rocks featuring strong heterogenic, metamorphic textures are not treated in the framework of this study for simplicity purposes and will be the subject of future research. Also, for simplicity purposes, we do not consider spatial variations in crystal sizes in our calculations (i.e. a porphyritic texture). We justify this approach because a related inhomogeneity is likely to be averaged out if one considers a several-metre-thick rock column. Additionally, the rock is considered to consist only of crystalline components; i.e. glassy materials such as obsidian have to be treated separately. Porous media can be approximated by assigning one of the constituents as air or (in the case of a pore fluid) water. This is explicitly done for the case of arkoses (10 % air) and arenites (11 % air). We compare the energy loss of muons in these rocks and hence the resultant muon flux attenuation depending on depth with those of the standard rock. The analysed lithologies, together with their relevant physical parameters, are listed in Table 1. Among these parameters, {Z/A} and {Z 2 /A}, i.e. the ratio of the atomic number (and its square) to the mass number averaged over the entire rock, are most relevant to the energy loss of muons (Groom et al., 2001). The former is almost proportional to the ionisation energy loss that occurs predominantly at low energies, whereas the latter is mostly proportional to the radiation energy loss that becomes dominant for muons faster than their critical energy at around 600 GeV. The volumetric mineral fractions of these 10 rocks can be found in Appendix A.

Cosmic-ray flux model
We perform our calculations with the muon energy spectrum model proposed by Reyna (2006), at sea level and for vertical incident muons. This model describes the kinetic energy distribution of the muons before they enter the rock. The calculation of the integrated muon flux after having crossed a certain amount of material is done in two steps. First, the minimum energy required for muons to penetrate a given thickness of rock is calculated considering the chemical composition effects (see Sect. 2.3). Afterwards, the energy spectrum model, dF /dE, is integrated above the obtained minimum energy (which we call from here on "cut-off energy", E cut ) to infinity, i.e.
The integration is necessary, as most detectors, which have been used for muon tomography, record only the integrated muon flux. As already stated in the introduction, we attribute a systematic uncertainty of ±15 % to the integrand dF /dE. All the calculations in this work have been verified with another flux model (Tang et al., 2006) and are presented in the Supplement.

Muon propagation in rocks
As soon as muons penetrate a material, they start to interact with the material's electrons and nuclei and lose part of their kinetic energy. The occurring processes can be categorised into an ionisation process, i.e. a continuous interaction with the material's electrons, and radiative interactions with the material's nuclei (i.e. bremsstrahlung, electron-positron pair production and photonuclear processes), which are of a stochastic nature. All these processes are governed by the material density ρ and the atomic number Z and atomic weight A (see Groom et al., 2001, for details). Our general strategy for the calculation of the energy loss in a rock is to use its decomposition into energy losses for the corresponding minerals. Accordingly, the energy loss of muons travelling a unit length, dE/dx, in a rock can be described by a volumetrically averaged energy loss through its mineral constituents: where ϕ j is the volumetric fraction of the j th mineral within the rock. The derivation of Eq.
(2) can be found in Appendix B.
In order to exploit this abstraction efficiently, we have to assume a homogeneous mineral distribution within the rock. This is a strong simplification, considering, for example, effects related to a local intrusion, tectonic processes like folding and thrusting, or spatial differences in sedimentation patterns. These concerns can be addressed through averaging over a large enough volume. Figure 1 shows two typical thin sections from rock samples of our experimental site that exhibit crystal sizes well below 4-5 mm. As muon tomography for geological purposes generally operates at scales of 10-1000 m, it is safe to assume that small-scale variations are averaged out. Thus, the term on the right-hand side of Eq. (2), i.e. the energy loss across each mineral, can be written as where a and b are the ionisation and radiative energy losses across a given mineral, respectively. These two parameters are in turn calculated by averaging the contribution of each element (i.e. atom) constituting the mineral by their mass (see Eqs. B5 to B15 in Appendix B for details). The density of the minerals, ρ mineral , is estimated from its crystal structures (see Appendix A for more detailed instructions). Once the energy losses are obtained for all minerals, each contribution is summed up according to Eq. (2). The energy loss within the rock can then be expressed in a similar way, as in Eq. (3) (for a detailed discussion, we refer to Appendix B): Again, the values {a} and {b} indicate the averaged ionisation and radiative energy losses across the whole rock, respectively. Equation (4), an ordinary nonlinear differential equation, is usually given as a final value problem; i.e. we know that the muon, after having passed through the rock column, still needs some energy to penetrate the detector, E det . This can be turned into an initial value problem by reversing the sign of Eq. (4) and defining the detector energy threshold as initial condition.
The problem has been transformed into the one of finding the final energy, the cut-off energy, E cut , after a predefined thickness of rock. This is a well-investigated problem, for which a great variety of numerical solvers are available. In this work, we employ a standard Runge-Kutta integration scheme (see, for example, Stoer and Bulirsch, 2002). The energy loss equations are subject to systematic uncertainties, mainly because the experimentally determined interaction cross-sections have an attributed error. According to Groom et al. (2001), the error on ionisation losses is "mostly smaller than 1 % and hardly ever greater than 2 %". These authors also state that, in the case of compounds, the uncertainties might be thrice as large. Therefore, we considered an ionisation loss uncertainty of ±6 % as appropriate for our calculations. The errors on the cross-sections of bremsstrahlung, pair production and photonuclear interactions are ±1 %, ±5 % and ±30 %, respectively. Appendix C shows in detail how we propagated these errors to the cut-off energy, E cut .

Figures 2 and 3
show the muon flux simulations as a function of rock thicknesses up to 2 km for igneous and sedimentary rocks, respectively. The depth-intensity relation is described by a power law, as it is the integration of the differential energy spectrum of muons, which also follows a power law.
To better visualise the difference between the fluxes after having passed these 10 rock types and the standard rock, we report the ratio between fluxes calculated after the different materials and that after the standard rock in Fig. 4: The attenuation of the muon flux expectedly depends predominantly on the rock density, as we can see in Figs. 2 to 4. Rocks exhibiting a high material density result in a larger muon flux attenuation than lithologies with a lower density. This, however, only depicts the overall differences, including density and compositional variations, between real and standard rocks. In this regard, Groom et al. (2001) apply an explicit treatment of density variations of known materials. Thus, the flux data can be simulated for a standard rock with the exact density as its real counterpart. Such a density normalisation enables us to isolate the compositional influence on the computed data. Figures 5 and 6 show the muon flux simulations for each rock compared to a density-    Table 1 as a function of rock thickness.  Table 1 and a density-modified standard rock. The mean flux is indicated by a bold line and 1σ bounds are indicated by the shaded area.
normalised standard rock, and Fig. 7 summarises this information by representing the ratio between muon fluxes after passing through real rocks and the muon flux after passing through a density-normalised standard rock. It is important to note that the standard rock muon flux in each flux ratio has been normalised with respect to the density of the original rock (i.e. the peridotite is compared to a standard rock of density ρ = 3.340 g cm −3 ; the limestone is compared to a standard rock of density ρ = 2.711 g cm −3 ). One notices that the flux ratios are rather close together, mainly within   Table 1 as a function of rock thickness.
2.5 % of the standard rock flux, before they start to diverge towards larger (dolomite, shale and arenite) and smaller (igneous rocks, arkose, limestone and aragonite) flux ratios beyond 300 m thickness of penetrated rock. Even though the errors on the fluxes are relatively large and sometimes even overlap with the standard rock fluxes, the propagated errors on the flux ratios remain well bounded near their means. This effect is due to the correlation of the errors in the numerator and the denominator in Eq. (6). A detailed discussion of how uncertainties have been propagated is presented in Appendix C.

Discussion
The differences in the calculated muon flux illustrated in Figs. 2 and 3 become even more pronounced in Fig. 4, where the fluxes are compared to the case where cosmic fluxes are attenuated by a standard rock. One notices a direct correlation with material density. This is reinforced by the fact that the granite (Fig. 2) has the same density as the standard rock (2.650 g cm −3 ) and shows an overall similar flux magnitude as the standard rock, i.e. a flux ratio of 1. This can be explained by Eq. (4), as the energy loss is almost directly proportional to the density, while the presence of density in the ionisation loss term (i.e. {a(E, ρ, A, Z)}) is negligible com-pared to this factor. Thus, if the rock flux data are compared to a standard rock with equal density, this effect should be removed, and one is left with the composition difference only. A closer look at Fig. 7 reveals that the muon fluxes for every rock below 300 m do not depart more than 2.5 % from their respective density-modified standard rock flux. The chemical composition effect can thus be considered negligible when compared to the systematic uncertainty originating from the muon flux model. We explain this through the dominance of the ionisation energy loss in this thickness region. Muons that penetrate down to 300 m of rock are still slow enough to predominantly lose their kinetic energy for the ionisation of the rock's electrons. As the number of electrons per unit volume is given by the product, ρ rock · {Z/A}, ionisation losses are proportional to this term. When comparing a density-normalised standard rock with a real rock, the only difference can emerge from the second part, i.e. {Z/A}. According to Table 1, these values do not change more than 1 % with respect to each other.
When the rock thicknesses become larger than 300 m, the flux ratios start to exceed ±2.5 % and the ratio patterns diverge. This corresponds to the point where radiative losses start to become the dominant energy loss processes. The latter are interactions of the muon with the nuclei of the atoms within the rock and its cross-section is mainly proportional to the square of the nucleus' charge (i.e. {Z 2 /A}). Hence, rocks that exhibit a lower {Z 2 /A} value than a standard rock (e.g. dolomite, arenite and shale) attenuate the muon flux less (i.e. flux ratio > 1), while all igneous rocks as well as limestone, aragonite and arkose, that have a higher {Z 2 /A} value, attenuate the muon flux more, which results in a lower flux ratio.
The above results reflect only the most striking connections to the chemical composition of a rock. In reality, however, the nature of muonic energy loss processes is much more complex than the shape of the flux ratios in Fig. 4 below 300 m suggests. The actual ionisation energy loss, Eq. (B27), is an interplay of the mean excitation energy {I }, i.e. the mean energy needed to ionise a material's electrons, the material density ρ rock , {Z/A} and various correction terms that depend on these parameters. These additional factors are also responsible for the nonlinear behaviour of the flux ratios between 100 and around 600 m, as effects from radiative losses start to become significant. However, as the resulting differences due to these processes remain smaller than 2.5 %, a detailed discussion of these matters falls beyond the scope of this paper.
As we see above, the muon flux calculation is significantly biased when one employs the standard rock assumption and thus neglects the effect of the chemical composition, especially when the thickness of the rock is beyond 300 m. This systematic error would then later turn into an over-or an underestimation in the assessment of density structures. We can roughly estimate the error on a thickness estimation of a cer-tain structure by employing the following formula: Here, x SR (F ) and x Ro (F ) denote the thickness of standard rock and a real rock, respectively, needed to attenuate the cosmic-ray muon flux to F . This is possible because the flux, as a function of rock thickness, is a strictly decreasing function. The domain of this function ranges from zero to infinite thickness, where its image takes the values from the initial flux, F 0 , to zero. On these two sets, the function is a bijection, and therefore an inverse function, x(F ), exists. Although its functional form might be unknown, it is still possible to interpolate between the simulated points. For our rocks, this is shown in Fig. 8.
As an example, in the case where the target is 600 m thick and made of limestone (ρ = 2.711 g cm −3 ), the standard rock assumption underestimates the flux by 7 %-8 % and thus overestimates the thickness by around 15 m or 2.5 %. The same is valid for basalt and aragonite.
The above discussion concentrates on calculations of the mean values of model parameters. A full description encloses also the propagation of their uncertainties. The rather large error bounds on single flux calculations stem from the uncertainties in the flux model and in the interaction crosssections. However, by taking a ratio, i.e. Eq. (6), of quantities with correlated errors, the resulting uncertainty on the ratio tends to cancel out. If the errors were propagated by linear operations, they would even cancel out perfectly. The small error bars which are still present in Figs. 4, 7 and 8 can be seen as effects of the nonlinearity in the differential equation, Eq. (5).
Because this is a pure sensitivity study, we cannot offer distinct measurements to verify our predictions. The reason for this is mainly because dedicated experimental campaigns have not yet been conducted, and thus such data are not available. We suggest that future studies in this field will address the composition issue and try to experimentally constrain this theoretical model. Nevertheless, our inferences are based on the same conceptual framework that has already been used for other materials, including standard rock. As a result of this, we find significant differences if the rock parameters are changed, especially for rock thicknesses larger than 300 m.

Conclusions
Our results suggest that it is safe to use the standard rock approximation for all rock types up to thicknesses of ∼ 300 m, as the flux ratio will mainly remain within 2.5 % of the standard rock flux, which generally lies within the cosmic-ray flux model error. However, we also find that beyond these thicknesses the use of the standard rock approximation and its density-modified version could lead to a serious bias. This mainly concerns basaltic and carbonate rocks. The flux error for these rock types increases with growing material thickness. It can be extrapolated that the errors grow even further beyond 600 m of material thickness up to a point where any inference based upon this approximation becomes difficult. This is, however, a thickness range where muon tomography becomes increasingly hard to perform, as lower fluxes have to be counterbalanced by larger detectors and longer exposure times.
In order to account for the composition of rock, it is advisable to undertake a geological study of the region alongside the muon tomography measurements, especially when faced with basaltic rocks or carbonates, which includes at the least the analysis of local rock samples. Auxiliary data could comprise rock density measurements (i.e. helium pycnometer or buoyancy experiments), chemical composition and mineralogical information (i.e. X-ray diffractometry/fluorescence measurements) as well as microfabric analyses (i.e. mineral and fabric identification on thin sections). This additional information may help to constrain solutions of a subsequent inversion to a potentially smaller set. The use of additional information, such as spatial information in the form of a geological map or a 3-D model of the geologic architecture, is strongly encouraged, because it might greatly improve the state of knowledge about the physical parameters that are to be unravelled.
Data availability. All data necessary to reproduce our results are included in the paper.
Mineralogical information is available from Tables A1 and A2 in Appendix A.
The equations to calculate the physical parameters for the different rock types and the simulated fluxes are listed in the main body of the text as well as in Appendix B.

Appendix A
To estimate the mineral density, we assume that it can be calculated by dividing the mass of the atoms within the crystal unit cell by the volume of the latter (see, for example, Borchardt-Ott, 2009): In this equation, M is the total molar mass of one mineral "formula unit", Q is the number of formula units per unit cell, and V Unit Cell is the volume of the unit cell. The latter is calculated by the volume formula of a parallelepiped: Equation (A2) can be rewritten as Here, a, b, c denote the unit cell vectors and their lengths; · is measured in Ångströms, i.e. 10 −10 m, whereas α, β, γ are the internal angles between those vectors. These six parameters can be looked up for each mineral in the crystallographic literature (e.g. Strunz and Nickel, 2001). The volumetric percentages of the minerals that constitute the 10 investigated rock types are shown in Tables A1 and A2. They were chosen as a reasonable compromise from literature values (e.g. Best, 2003;Tuttle and Bowen, 1958;Folk, 1980).

Appendix B B1 Energy loss in elements
The average spatial differential energy loss can be written in a rather simple form (Barrett et al., 1952): Here, ρ, A and Z denote the mass density, atomic weight and atomic number of the penetrated material, while E is the kinetic energy of the penetrating, charged particle, and x is the position coordinate. The function a (E, ρ, A, Z) in Eq. (B1) is the differential ionisation energy loss that accounts for the ionisation of electrons of the penetrated material. In the case of incident muons (i.e. electric charge q µ = −1C and mass m µ = 105.7 MeV c −2 ), the relationship expressed in Eq. (B1) takes the form In this equation, βγ are the relativistic factors and are therefore a function of the kinetic energy E. The constant m e denotes the mass of the electron, and c is the speed of light. Q max is the highest possible kinetic recoil energy of scattered electrons in the medium, while K is a constant incorporating information about the electron density. The function δ (ρ, Z, A) is a correction factor, which considers the mechanisms where the material becomes polarised at higher muon energies, with the consequence that the energy loss is weaker (Sternheimer, 1952). The last term in Eq. (B2) is another correction factor, which considers bremsstrahlung from atomic electrons (not the incident muon, which would be the term in Eq. 3) that also appears at higher muon energies. A more detailed explanation of this equation and its parameters can, for example, be found in Olive et al. (2014). In contrast to Eq. (B2), the function b (E, A, Z) describes all the radiative processes that become dominant at higher velocities (above ∼ 600 GeV c −1 for muons). This term includes energy losses due to bremsstrahlung, electron-positron pair production as well as photonuclear interactions. These different contributions can be written independently of each other: b (E, A, Z) = b brems (E, A, Z) + b pair (E, A, Z) Each process in Eq. (B3) is computed by integrating its differential cross-section with respect to every possible amount of transferred energy: Here, N A is the Avogadro number and ν = ε/E the fractional energy loss (whereas ε is the absolute energy loss) for this process. Specific cross-sections for bremsstrahlung (Kelner et al., 1995, 1997), photonuclear (Bezrukov and Bugaev, 1981 and pair production (Nikishov, 1978) energy losses are used by Groom et al. (2001) for the calculations of their tables. As this pair production cross-section involves the calculation of many computationally extensive dilogarithms, an equivalent cross-section (Kelner, 1998;Petrukhin, 1969, 1971), which is used in GEANT4 (Agostinelli et al., 2003) by default, is used in our study.

B2 Energy loss in minerals
Since the above equations are valid for pure elements, adjustments are needed for compounds (e.g. minerals) and mixtures thereof (e.g. rocks). Generally, it is advised to use the physical parameters for materials that have already been measured (see Seltzer and Berger, 1982, for a compilation). However, except for calcium carbonate (i.e. calcite) and silicon dioxide (i.e. quartz), no other minerals have been investigated. This also means that there is no standard approach available for considering natural rocks. Fortunately, for such materials, a theoretical framework has been proposed (see, for example, Appendix A of Groom et al., 2001). The basic idea is to consider the compound as a single "weighted average" material and the energy loss therein as a mass-weighted average of its constituents' energy loss: The weights w i are calculated according to the atomic weights A i of the elements: and can then be used to calculate an average Z/A value: Equivalently, the average Z 2 /A value can be calculated according to One more change must be made to the ionisation loss Eq. (B2) in order to appropriately account for the change in the atomic structure that emerged due to chemical bonding of the elementary constituents. This is reflected in a modified mean excitation energy I , which can be calculated by taking the exponential of the function which is basically a weighted geometric average of the elementary mean excitation energies: One has to pay attention that the mean excitation energies for some elements, I i , can change quite significantly when they are part of a chemical bond. A guideline to address this issue can be found in Seltzer and Berger (1982). Equations (B7) to (B10) are still a consequence of Eq. (B5). However, there is one term in the function δ(ρZ/A) in Eq. (B2) that is calculated differently from Eq. (B5). This concerns the logarithm of the plasma energy of the compound, which for an element is given by (e.g. Olive et al., 2014) ln ω p = ln 28.816 · ρ Z A .
According to Eq. (B5), the plasma energy for a compound should be calculated the same way as the mean excitation energy in Eq. (B9). However, Sternheimer et al. (1982) and Fano (1963) explicitly advise us to use the arithmetic mean within the logarithm when dealing with an atomic mixture (i.e. a molecule), yielding ln ω p = ln 28.816 ρ mineral Z A .
This results in the modified ionisation energy loss: a(E, ρ mineral , A, Z) The radiation loss for the compound, on the other hand, is only a linear combination of the radiation losses of its elementary constituents, Eq. (B3), yielding The resulting Eq. (B15), has now the same form as the energy loss Eq. (B1) for elements and can be solved accordingly.

B3 Energy loss in rocks
To obtain an energy loss equation for rocks, a similar procedure as for forming minerals through the assembly of elements can be applied. Starting from Eq. (B5), we consider the energy loss for a rock as the mass-weighted average of the energy losses of its mineral constituents: where q j are the mass fractions of the j th mineral within the rock, analogous to Eq. (B6), Using dχ = ρ · dx, Eq. (B16) then takes the following form: Multiplying both sides with ρ rock and applying the definition of the density, ρ = m/v, that can also be written as v = m/ρ, Eq. (B19) becomes If one sets ϕ j = v mineral,j /v rock , the volumetric fraction of the j th mineral within the rock, Eq. (B20) transforms into the compound equation for rocks: Analogously to the mineral case, we can now define new average energy loss parameters for the rock, beginning with its overall density: The average {Z/A} is given by and similarly, the average {Z 2 /A} can be calculated according to The rock's mean excitation energy is The only difference between the rock calculation and the mineral calculation enters in the calculation of the plasma energy. While in the mineral case we were advised to use Eq. (B11) instead of what would naturally follow from the weighted average in Eq. (B5), we prefer to use the weighted average, Eq. (B21), for the case of rocks. The reason for this lies in the fact that the density effect operates on a nanometric scale, whereas minerals generally have sizes between several micrometres and a few centimetres. In the case of a mineral compound, the molecular structure comprises also a nanometric scale. These parameters can then be rearranged into an ionisation loss term for a rock: Like Eq. (B14), the radiative losses can be rewritten as a weighted average of the mineral radiative losses: Equations (B27) and (B28) can then be joined together to form again a similar term to Eqs. (B1) and (B15): the energy loss equation for rocks. We want to stress that the starting point of the derivation of the energy loss equation for rocks is a mass averaging of mineral energy losses. Therefore, the mass-averaging approach is inherently included in this approach. In fact, mass-averaging and volumetric averaging are two equivalent descriptions of the same problem. For the mass-averaged formulae, we refer to the Supplement to this paper.

Uncertainty propagation
The first step in our uncertainty treatment includes a propagation of the interaction cross-section errors (σ a = ±6 %, σ b brems = ±1 %, σ b pair = ±5 %, σ b photonucl = ±30 %) to the cut-off energy, i.e. by solving the differential equation Eq. (5). Generally, a higher cross-section yields a higher cutoff energy, as the muon needs more initial kinetic energy, which it then loses on the way, and vice versa. In order to estimate a lower and an upper error bound on the cut-off energy, E cut , we use a conservative approach. This means that the lower cut-off energy error bound is calculated by setting all cross-sections to their lower 1σ bound and running the simulation with these modified values. The upper error bound is calculated accordingly. Of course, this overestimates the effective error; however, if our calculations remain valid within this conservative error, then they can also be trusted with a conventional error.
The second step is the estimation of the error regarding the integrated flux. Here, we need to propagate the errors through Eq. (1) to the simulated flux. There are two different errors present at this stage. The first one includes the error on the lower integration boundary, i.e. E cut , which has just been calculated above. The second error addresses the integrand, i.e. the flux model. Figure C1 visualises the concept behind the propagation of these two errors. The simulated flux error is equivalent to the error which is made by calculating the area under the graphs. We estimate the lower error bound on the simulated flux (i.e. smallest area) by taking the upper error bound on E cut and the lower error bound on dF /dE. Similarly, the upper error bound on the simulated flux (i.e. largest area) is calculated by setting E cut to its lower error bound and dF /dE to its upper error bound. Again, this is a conservative approach, which we justify with the same rationale as above. The last step addresses the propagation of the simulated flux errors to the flux ratio in Eq. (6). Here, we can make use of the fact that the errors in both simulations are perfectly correlated. In other words, if we knew the errors on all affected quantities in one simulation, we would instantaneously know the corresponding values for any other simulation. This allows us, for example, to calculate the upper error bound on the flux ratio by dividing the upper error bound of the simulated flux in the numerator by the upper error bound of the simulated flux in the denominator. The same is valid for any other constellation of errors, including the lower error bound and the mean.
Author contributions. AL, FS and AE designed the study. AL developed the code with contributions by MV. AL performed the numerical experiments with support by RN. DM and AL compiled geological data. AA, TA, PS, RN and CP verified the outcome of the numerical experiments. AL wrote the text with contributions from all co-authors. AL designed the figures with contributions by DM. All co-authors contributed to the discussion and approved the manuscript.