Bedrock sculpting under an active alpine glacier revealed from cosmic-ray muon radiography

Mountain glaciers form landscapes with U-shaped valleys, roche moutonées and overdeepenings through bedrock erosion. However, little evidence for active glacial carving has been provided particularly for areas above the Equilibrium Line Altitude (ELA) where glaciers originate. This is mainly due to our lack of information about the shape of the bedrock underneath active glaciers in highly elevated areas. In the past years, the bedrock morphology underneath active glaciers has been studied by geophysical methods in order to infer the subglacial mechanisms of bedrock erosion. However, these comprise surveys on the glaciers’ surface, from where it has been difficult to investigate the lateral boundary between the ice and the bedrock with sufficient resolution. Here we perform a muon-radiographic inspection of the Eiger glacier (Switzerland, European Alps) with the aid of cosmic-ray muon attenuation. We find a reach (600 × 300 m) within the accumulation area where strong lateral glacial erosion has cut nearly vertically into the underlying bedrock. This suggests that the Eiger glacier has profoundly sculpted its bedrock in its accumulation area. This also reveals that the cosmic-ray muon radiography is an ideal technology to reconstruct the shape of the bedrock underneath an active glacier.

surveys and gravity measurements, ground penetrating radar surveys, and topographic modelling to reconstruct the bedrock topography of formerly glaciated areas or underneath active glaciers [20][21][22][23][24][25][26][27] . Despite the progress of such technologies, the bedrock morphology along the sides of a glacier, particularly in remote high-elevation alpine areas, has been hardly constrained, because these methods are mainly performed from above the glaciers' surfaces. In addition, most of the alpine cirques are hardly accessible particularly in their accumulation areas where they originate. Where surveys were possible, the resolution of the data decreases rapidly towards the glaciers' bases and lateral sides. This is related to (1) the presence of non-consolidated sediment deposits on the glacier margin and (2) the steep dip of the bedrock surface underlying most of the surveyed glaciers. In addition, in the case of temperate glaciers, high fluid pressure ratios lower the resolution of seismic and radar surveys 28,29 , thereby thwarting a high-resolution reconstruction of the bedrock shape underneath modern glaciers.
In the present work we apply the cosmic-ray muon radiography technology to map the lateral margin of the Eiger glacier ( Fig. 1) situated in Switzerland (Central European Alps at 46°34′05″N latitude and 7°59′56″E longitude). This method is based on the high penetration power of cosmic-ray muons that hit the Earth's surface continuously 30 , where the attenuation rate of the muon flux mainly depends on the density of the traversed material 31 . The idea dates back to 1955 with the measurements of the thickness of rock above an underground tunnel 32 , followed by the pyramid's inspection by L. Alvarez in 1970's 33 , and it has spread to multiple fields such as volcanology 34,35 , geology 36 , and the non-destructive inspection of reactors 37 in the past decade. In 2015, Nishiyama et al. 31 launched a pilot survey for the application of such a technology in an Alpine environment using emulsion films as muon detectors. These authors 31 successfully mapped a small portion of the bedrock underneath the uppermost part of the Aletsch glacier (Central European Alps) over an area of c. 50 × 100 m 2 . Based on the success of this pilot work, the present study aims at imaging of a much larger region where a glacier originates and where active erosion and bedrock sculpting are likely to occur.
Muon radiography requires muon detectors to be placed at an elevation lower than the surveyed target because of the downward going path of cosmic rays. We benefited from the railway tunnel of the Jungfraubahn, which runs through the bedrock surrounding and beneath the headwall reaches of the Eiger glacier (Fig. 1). In 2017, we installed muon detectors made of emulsion films at three sites along this tunnel (sites ES, TA and TB, Fig. 1). Emulsion films, which have been applied in various experiments on fundamental physics 38,39 , record the trajectories of incoming muon particles, which can be observed with optical microscopes as sequences of tiny silver grains after chemical development 38 . The microscopic analysis of the films allowed the attenuation pattern of muon intensity to be measured 38 , which in turn gives the information on the density of overburden material 31 . The observed density information is then quantitatively converted to the location of the boundary separating the low-density ice from the high-density bedrock 31 , providing a unique and high-resolution reconstruction of the bedrock topography at the bottom and lateral margins of the Eiger glacier.

Results
Morphology of the eiger glacier. We explored the morphology of the Eiger glacier, which is the target of our survey (Fig. 1). The nearly 2.1 km 2 large glacier originates on the western flank of Mt. Eiger at approximately 3600 m a.s.l., from where it flows and terminates at an elevation of approximately 2400 m a.s.l. over a length of c. 2 km (Fig. 1). The bedrock geology of the study area comprises a suite of NW-dipping recrystallized www.nature.com/scientificreports www.nature.com/scientificreports/ micritic limestones that display a penetrative horizontal foliation and that overly polymetamorphic gneisses (inset Fig. 1b) 40 . The bedrock has thus a constant fabric and comprises the same lithology within the entire glacial catchment 40 . The mean value of the bulk density of the bedrock in this region is 2.68 g cm −3 , determined from 16 bedrock samples 31 , and its standard deviation is 0.02 g cm −3 . The Eiger glacier consists of four morphologic domains I-IV ( Fig. 1): (I) a headwall reach forming a concave amphitheatre-shaped cirque with c. 50° steep flanks (including the headwalls), (II) a relatively flat domain, and (III) a middle segment where a prominent bedrock ridge forces the ice flow to diverge over c. 250 m distance before converging again farther downslope. There, the bedrock ridge, which most likely corresponds to a "roche moutonée" (i.e. a distinctive bump of 10-100-m scale typical of glacially-scoured landscapes with smooth abraded slopes facing upglacier, and abrupt slopes on its downglacier face), is exposed in the middle of the glacier, and the ice surface in segment (iii) is rugged and dissected by several transverse crevasses. Segment (III) also evidences a drastic change in ice-flow direction from South-West to West. The exposed bedrock ridge separates the Eiger glacier in two parts: a slip-off, steep margin situated on the inner northern-side of the glacier (referred hereafter to the "inner glacier"; Fig. 1), and a cut-bank, flatter glacier part on its southern outer margin (referred hereafter to the "outer glacier"; Fig. 1). The down-slope end of segment III also lies within the altitudinal range of the modern ELA, which is situated at c. 3000 m a.s.l. 41 . Finally, segment IV is characterized by a terminal lobe with longitudinal crevasses and a till.
A closer inspection of the glacier morphology discloses further details, particularly along segment II where the Eiger glacier starts to separate into two branches. Downslope of segment II where the ice-flow bifurcation begins, the ice-surface strongly differs between the inner and outer glaciers. The inner glacier has a straight ice flow, with a relatively small ice surface, and the surface elevation rapidly drops from 3400 to 3100 m a.s.l. over a short distance of c. 650 m, yielding in an average surface slope of c. 25°. Contrariwise, the outer glacier follows a large bend where the ice-flow orientation changes by 90° within a few hundred meters. Along this bend, the ice-surface elevation drops of about 300 m over a distance c. 1000 m. Interestingly, this glacier segment is also the location where cuts into its bedrock wall are visible, forming two secondary concave niches α and β with a shape that is similar to an incipient cirque ( Fig. 1). While the upper niche α is situated slightly upstream the bend and only weakly developed, the lower niche β is located within the bend and characterized by a distinct concavity. The ice surface maintains a constant dip of c. 10° above the 90°-direction bend, and then steepens to c. 22° farther downslope after the flow direction has changed to the West. At this point, transverse crevasses suggest that the ice is under extension, while the absence of any crevasses along the flatter segments implies that the glacier is under compression. Visual inspection shows that the outer glacier has a surface area of c. 300,000 m 2 and is thus >200% larger than the inner glacier ( Fig. 1). Figure 2 shows the direction of arriving muons at the three detector sites (ES = Eismeer Station, TA = Tunnel site A and TB = Tunnel site B; Fig. 1). In these diagrams, the direction is represented as azimuth and elevation angles (θ x , θ y ). The number of observed muons are 2.3, 5.3 and 7.9 × 10 3 for ES (exposure time: 106.8 days, effective area: 1512 cm 2 ), TA (164.9 days, 1296 cm 2 ) and TB (106.8 days, 1080 cm 2 ) sites, respectively. The population of the recorded muons is clearly anti-correlated with the thickness of the bedrock and the ice along the straight muon trajectories (Fig. 2). Specifically, the shadows on each plot where the data density is low (sparsely distributed data points) coincide with the regions where the obstacle thickness is greater than 1 km. This indicates that most of the muons from these directions were absorbed in the thick mountain edifice. We binned these registered muons in a rectangular histogram (tanθ x , tanθ y ) and converted the data into a normalized flux (cm −2 sr −1 sec −1 ). Figure 3 illustrates how the muon flux data is used to infer the density of the traversed material. First, we focus on the flux data of muons, which passed purely inside the bedrock without crossing the glacier. Such data are available from the TA detector ( Fig. 3: open circles, and Fig. 1 for location of detector). As seen in Fig. 3, the flux decreases nearly exponentially with the increase of the bedrock thickness, and its attenuation depends also on the zenith angle of the muons (represented by the colours of the data points). These features of the muon flux attenuation in the bedrock are well reproduced by the simulation for different zenith angles (45°, 60° and 75°) under the assumption of a uniform bulk rock density of 2.68 g cm −3 ( Fig. 3: solid curves) 31 . This agreement assures that the variety of the bulk density inside the bedrock is small, as is observed with the standard deviation of the density sampling (0.02 g cm −3 ). On the other hand, trajectories where muons cross partially the glacial ice before arriving at the detectors yield flux values that are by a factor of up to 2 higher than these theoretical calculations ( Fig. 3: solid circles). This excess of the muon flux is due to the lower bulk density of ice compared to that of bedrock, thus allowing a higher muon transmission. We benefitted from these differences in muon attenuation between ice and rock: the magnitude of the muon flux attenuation recorded at the different tunnel sites can be quantitatively interpreted in terms of the ice-bedrock proportion along the muons' trajectories 31 , which is then converted into the position of the boundary between the basal ice and its underlying bedrock. The three detectors surrounding the glacier provided information about the bedrock shape from different perspectives. The bedrock positions were first estimated for individual detectors, and the results are combined and re-sampled into a three-dimensional representation of the interface geometry (see Method section).

Muon attenuation pattern.
Reconstruction of the bedrock shape. We determined the bedrock shape underneath the Eiger glacier over a 600 m-long (NE-SW) and 300 m-wide reach (NW-SE) (Fig. 4). This corresponds to the segment where the glacier surface is relatively flat and where it dips at c. 10° (segment II on Fig. 1) just below the >30°-50° steep headward reach (segment I on Fig. 1). The elevation resolution of the inferred bedrock surface ranges from 10 m (1σ) in vertical dimension along the relatively flat segment (which is very close to the TB detector), to about 30 m (1σ) in the headward reach where statistics are deteriorated due to the thick mountain edifices. The inferred glacier thickness is typically around 50 m (down-slope segment II on Fig. 1) to 100 m (up-slope segment I on www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 1). Along the thalweg axis, the bedrock surface continuously dips in the downslope direction, and neither an overdeepening nor transverse bedrock knobs could be observed (Figs 4 and 5).
The cross-sectional views (Fig. 5) present quantitative information on the bedrock shape beneath the Eiger glacier, particularly on the southern outer side of the bedrock ridge. In the up-slope reach (segment II), represented by cross-section (a) on Fig. 5, the bedrock steeply dips at a constant angle >50° from the lateral valley wall above the glacier down to 50 m beneath the ice surface. There, the transition to the flat glacier base is gradual without any distinct break-in-slope. The situation markedly changes farther down-slope where the boundary between the aerial and subglacial lateral margin is characterized by a sharp break-in-slope (sections (b) and (c) on Fig. 5). There, the dip angle of the bedrock surface changes from c. 50° above the glacier to sub-vertical at the ice contact. In addition, cross-section (c) suggests that some undercutting may occur along this segment (Fig. 5c). At c. 50-60 m depth, the cross-sectional base of the subglacial bedrock is flat and nearly horizontal, resulting in a L-shaped cross-sectional geometry of the bedrock surface underneath the ice. This is also the area where the ice surface starts to bend from a South-West to a West-oriented flow direction together with the development of a niche (α) similar to an incipient cirque (Figs 1 and 5). Farther downslope of cross-section (c), the cross-sectional geometry changes again. Above the glacier, the bedrock surface dips at 40°−50° and steepens by a few degrees at the ice surface, thus forming a break-in-slope where the bedrock plunges beneath the glacier (sections (d) and (e) on Fig. 5). Along these sections, the base of the glacier is flat to slightly tilted towards the glacier's outer side, and the transition to the lateral margin is much more gradual than for the upslope reaches (sections (b) and (c) on Fig. 5).
While the outer glacier is covered by many muon trajectories (Figs 1 and 5), the inner glacier on the northern side of the bedrock ridge is out of the field views of the detectors. Therefore, the elevation of the bedrock underneath the inner glacier has to be constrained through indirect evidence. Here, we use information from the elevation of the bedrock ridge itself (Fig. 5) to infer the elevation of the bedrock underneath the inner glacier. In particular, the bedrock depth along the outer glacier is situated at the same elevation, or even higher, than the ice surface of the straight inner glacier at least along cross-sections 'd' and 'e' in Fig. 5. This suggests that the underlying bedrock underneath the steep and straight inner glacier is likely to be at a lower elevation than the reconstructed bedrock below the outer glacier, at least along a section perpendicular to the ice flow.

Discussion
We verify our subglacial bedrock reconstruction within the Eiger headward region (segments I and II) following two pieces of evidence. First, the reconstructed bedrock positions agree with the mapped hillslope walls above the lateral glacier margin (Fig. 5). Second, the reconstructed bedrock positions also match with the domain where the bedrock is partially exposed. This is mainly the case for the bedrock ridge that separates the inner from the outer glacier (cross-sections in Fig. 5). These observations are supporting evidence for a reliable performance of the muon-flux attenuation analysis and the related calibration for our study site (see Methods). The muon www.nature.com/scientificreports www.nature.com/scientificreports/ radiography method thus yields quantitative information on the ice thicknesses (c. 10-30 m precision) of remote Alpine glaciers, if suitable detector sites such as natural bedrock galleries or tunnels are available underneath the target glacier. We thus propose that the technology of cosmic-ray muon radiography with emulsion detectors offers a quantitative tool that complements established geophysical methods 42 . In the case of the Eiger glacier, the steep bedrock slopes below the glacier even sharpen the density contrasts in the retrieved muon-derived images, while other methods (such as seismic or radio-echo soundings) often fail to yield information with acceptable resolution for similar bedrock geometries. Moreover, our present survey suggests that this novel technology works well in remote and harsh environments with the advantage of the emulsion films as muon detectors. The emulsion films showed an excellent performance in the remote and cold environment during the entire measurement period (spanning several months) without electricity or regular maintenance. This passive and robust feature of the measurement device cannot be attained by electronic detectors such as digital detectors including plastic scintillation and gaseous detectors, which have so far been applied to muon radiography in many prior works [43][44][45][46][47] . www.nature.com/scientificreports www.nature.com/scientificreports/ With respect to the Eiger glacier and its erosional mechanisms, the most prominent result of our work is the reconstructed shape of the lateral margin along the outer glacier. In the upper reach of our surveyed area (section (a) on Fig. 5), the shape of the bedrock-ice interface is characterized by a generally smooth transition from the lateral margin to the base of the glacier. In the reach where the glacier begins to curve and where a bedrock niche points to the occurrence and possibly ongoing development of an incipient cirque, the bedrock underneath the lateral glacier margin appears much steeper than the exposed hillslope (sections (d) and (e) on Fig. 5) and almost vertical-to-overhanging in some reaches (sections (b) and (c) in Fig. 5). Such bedrock-slope variations over short distances (few tens to hundred meters) require a strong lateral and headward erosional component to shape the nearly vertically-oriented margin of the subglacial bedrock. Interestingly, this is also the glacier segment where no crevasses are observed on the ice surface. This implies that the glacier is likely under compression along this reach, which provides a precondition for erosion on the lateral margin of a glacier to occur. Along the inner glacier reach, vertical glacial carving has been more efficient than along the outer curved glacier. We based this inference on the observation that the subglacial bedrock under the inner glacier is situated at lower elevations than along the outer curved reach.
The erosional mechanisms of an alpine glacier, related efficiencies and controls thereof have been explored through both empirical observations and numerical investigations over the past decades. This resulted in the overall notion that the local erosional capacity of an Alpine glacier depends on its flow velocity raised to some power, modulated by different controlling factors such as subglacial hydrology, fracture spacing of the underlying bedrock and debris concentration in the basal ice 4,6,7,[48][49][50] . It has also been proposed that the ice-flow velocity and the basal shear stresses are closely related, where shear increases with flow velocity raised to some power 6 . This likewise suggests that subglacial erosion and thus the glacial impact on mountainous landscapes exponentially increase with the basal shear stress at the ice-bedrock interface. Under warm-based glaciers, shear stresses and erosion rates can also be modulated by the subglacial hydrology 7 . In the case of the Eiger glacier, the curved geometry of the glacier reaches along the outer part of segments II and III may additionally result in a highly asymmetric pattern of shear stresses and flow velocities due to the curvature. This is explained by the contrasts in the slopes between the inward and outward segments of a curved glacier. As a consequence, the flow velocity and stress centre line, commonly situated in the middle of a straight glacier, may be shifted toward the inner side of a bend, where the amount of this shift increases with the curvature of the ice flow 51 . We use these mechanisms to explain the lower depth at which the bedrock has been carved along the straight and steep inner reach of the Eiger glacier (Fig. 5).
While the larger bottom shear stresses along the straight inner glacier could be invoked to interpret the lower elevations of the bedrock topography (see cross-sections of Fig. 5), these mechanisms alone are not capable of explaining the nearly vertical lateral dip of the ice-bedrock interface along the southern margin of the outer glacier. Such an erosional component, operating on the side of a glacier, requires substantial shear stresses on the lateral ice margin. This suggests that glacial erosion depends not only on the ice thickness and the slope, but also on the momentum of the flowing ice particularly where glaciers appear under compression (segment II, Fig. 1). Ice flow thus shows similarities with a fluid when considering the erosional mechanisms: both apparently have a vertical and a lateral component where concave banks can experience stronger erosion. This L-shape sculpting acts to widen the glacial thalweg in the lateral direction with a strong over-steepening of the lateral margin. After the glacier retreats, such oversteepened cliffs would be exposed above the glacier and collapse due to the loss of mechanical support or activated hillslope/fluvial processes 52 . This could be a reason why we do not observe this morphology in deglaciated areas. Alternatively, the oversteepened bedrock could represent the headwall of an incipient cirque, which additionally feeds, or starts to feed, the Eiger glacier from the southwest. We base this inference on the plan-view shape of the ice margin, which is nearly straight along segment I (Fig. 1) and then curves towards the headwall in the middle of segment II. Accordingly, the vertical to nearly over steepened glacier margin, which was imaged through the cosmic-ray muon radiography (sections (b) and (c) of Fig. 5) better reflects the erosional work along the cirque headwall through backward erosion.
The search for the controls on cirque wall retreat has received much attention in the past years 12,13,16 , yet with diverging conclusions, which range from bedrock shattering through freeze-thaw cycles 53 , to bedrock carving in response to rotational flow of ice, and to quarrying 16 . In addition, using numerical simulations where glacial erosion has been treated as a function of basal sliding, McGregor et al. 16 proposed that cirques are preferentially formed hundreds of meters below the ELA, and several studies suggested that the position of the cirque floor may reflect the averaged ELA position over multiple glacial-interglacial cycles 16,54,55 . Because our results highlight the occurrence of erosional carving at the inferred incipient cirque, and since the modern ELA (and thus also the long-term averaged ELA) is lower than the elevation of this incipient cirque, we suggest that glaciers are also capable of efficient headwall scouring above the ELA. In addition, the bedrock geometry underneath the cirque does correspond to the shape that is expected if a rotational flow, paired with shear along the lateral margin, is responsible for the occurrence of cirque wall retreat, supporting current views on the erosional mechanisms of glaciers and cirque evolution.
As a summary, we have illustrated a successful application of the cosmic-ray muon radiography, where we reconstructed the bedrock topography under an active glacier over an area of several hundreds of meters. Our results suggest that this technology is capable for mapping the interface between ice and bedrock at a resolution of 10-30 m, suitable for this study. It also allows the imaging of steep and nearly overhanging boundaries between ice and bedrock, which would not be depicted with other geophysical methodologies. From a geomorphological perspective, the results of our survey imply that active glacial backward erosion does occur in the accumulation area above the modern ELA (situated at c. 3000 m a.s.l. 41 ) and that this mechanism can be invoked to explain the formation and the long-term evolution of glacial cirques. Further high-resolution constrains on the bedrock topography for high-elevation glaciers will allow us to improve our understanding on the erosional mechanisms at work in cirque environments where glaciers originate. www.nature.com/scientificreports www.nature.com/scientificreports/ Methods emulsion detectors and analysis. Double-side coated emulsion films were used as muon detectors 38 .
The thicknesses of the plastic base and emulsion layers were 180 μm and 60 μm, respectively. The emulsion gels, which were produced by the Nagoya University 56 , were poured onto the plastic base at the underground facility at the University of Bern. The produced films were transported to the experiment sites and installed in the detector frames mounted onto the wall of the railway tunnel. In the detector frames, five or six films were layered, and 2-mm-thick stainless steel plates were inserted between each adjacent film 57 . The location and size of the detectors, as well as the exposure time are reported in Table 1.
The films were chemically developed after extraction and scanned with the automated readout microscopes at the Laboratory for High Energy Physics, at the University of Bern. The microscope consists of a CMOS camera, a motorized microscope stage and illumination 38 . This allows to take tomographic profiles of silver grains in the emulsion layers and to measure the position and direction of incident charged particles 38 . The procedure is as follows: (i) a sequence of silver grains in each emulsion layer is selected as muon track candidates (micro-tracks); (ii) micro-tracks from two layers of a film, which penetrate the plastic base in a straight way, are selected (base-tracks); (iii) base-tracks aligned in consecutive films are identified as muon trajectories, here defined by at least three base-tracks out of the total of five (TB) or six films (ES and TA). The processes (i) and (ii) were performed by means of an in-house developed software 58 , while (iii) was performed with the FEDRA software 59 . Figure 2 shows the incident angle distribution for the reconstructed tracks, where the angular resolution of the tracks is 3-10 milliradian. The number of detected muons was then determined in rectangular bins and converted into a muon flux (cm −2 s −1 sr −1 ) by normalizing with respect to the scanned area, the exposure time and the solid angle. The inefficiency of the films was estimated to ~10% at base-track level, which leads to the muon detection efficiency of 95-99%. We performed a correction of the measurement results taking the detection efficiency into consideration.
Flux simulation and calibration. Before reconstructing the bedrock shape, a calibration of the flux attenuation analysis was performed by comparing the observed flux of muons through pure rock with the corresponding simulated flux. TA site provides such data (18 bins from forward direction and 60 bins from backward). The simulation is performed assuming a rock density ρ_rock of 2.68 g cm −3 , which is the average bulk density of the local limestone, determined from 16 samples (2.68 ± 0.02 g cm −3 ) collected inside the railway tunnel and at the surface 31 . The calculation of the muon flux is based on the energy spectrum model of cosmic-ray muons and the interaction of muons with matter. We employed the spectrum model proposed by Tang 60 and the muon range for standard rock tabulated by Groom 61 . The topography of the mountain and the glacier surface is taken from a 2 m mesh digital elevation model (Swisstopo© with elevation uncertainty of <3 m) to derive the length of the muon trajectories from the topographic surface to the tunnel sites (L). Since this obstacle length varies within the bin due to the steep topography of the mountain, the bin is further divided into small hundred bins so that the roughness within the subdivided bin can be negligible. For each subdivided bin, the density-length traversed by muons is calculated by multiplying the rock density (ρ_rock) and the length of the muon trajectory (L). The minimum energy of muons (E_min) needed to penetrate this density-length can be looked up from the range table 61 . The muon flux is obtained by integrating the energy spectrum from E_min to infinity. The muon flux values calculated for hundred subdivided bins are then averaged out and set to the representative value of the original rectangular bin. Figure 3 shows the observed flux as a function of rock thickness and zenith angle. The observed data at the TA site generally agrees with the theoretical curves. The ratio between the observed and simulated flux is found to be 0.94 ± 0.04 (stat), independent of the zenith angle and the rock thickness. A small deviation of this ratio (0.94) from one can be regarded as a bias of this observation due to the systematic uncertainties in the detector efficiency and the muon energy spectrum model. The statistical fluctuation (4%) is comparable to the effect of the intrinsic variability of the bedrock density (1% from rock sample measurements) on the resultant muon flux (3%). In a further analysis, therefore, this difference is calibrated when comparing the observed and simulated flux by multiplying the simulated flux by this factor (0.94).
Bedrock shape reconstruction. We followed the approach by Nishiyama et al. 31 upon reconstructing the shape of the bedrock underneath the Eiger glacier. The observed muon flux for each bin yields the average bulk density <ρ> of the material confined in the viewing range of the bin. Since the confined region is a mixture of bedrock and ice, the fraction of the bedrock (x) is related to the bulk density of a bedrock (ρ_rock) and that of an ice component ρ_ice = 0.85 g.cm −3 by Huss 62 through the relationship: www.nature.com/scientificreports www.nature.com/scientificreports/ _rock x _ice (1 x) <ρ > = ρ ⋅ + ρ ⋅ − Once x is obtained for each bin, the boundary position can be plotted at a distance Lx from the detector position. Here, a constant density is inferred for the ice because snow and firn exhibit a lower density in only the uppermost 10 m of the glacier 62 , with significant density changes only occurring in the varying snow cover layer or in the presence of crevasses 63 , which are almost absent along the reconstructed glacier reach (Fig. 1). The resulting uncertainty is in the order of the DEM error 63 and thus negligible compared to the scale of our observation. The three-dimensional representation and error estimation are performed as follows. First, the observed muon flux dataset is multiplied to hundred synthetic datasets by adding a statistical fluctuation for each bin. Specifically, a random variable following a Gaussian distribution with a standard deviation of N is added when the number of muons in the bin is N. Subsequently, the ice-rock boundary position is calculated for each flux data in every synthetic dataset and plotted in a 3D space. In the end the plotted point clouds are re-sampled in a cylindrical coordinate system with the axis parallel to the glacial flow direction (south-west) with a division of Δz = 40 m along the axis and Δφ = 12° in azimuth. For each divided cylindrical bin (z, φ), the best position of the ice-bedrock boundary is given by taking the average of the radial coordinates of the points. The error of the boundary position is given by their standard deviation (yellow bands in Fig. 5). A synthetic data reproduction was introduced so that the magnitude of the statistical fluctuations is properly taken into account after re-sampling. Cylindrical bins, which contain less than 10 points, were neglected as lack of statistics and resolving power.