Introduction

Avalanches, crevasses or bad weather conditions have been the cause of many mountaineer fatalities on glaciers in the European Alps over the past two centuries. Since 1992, eight people have disappeared on Swiss glaciers (personal communication from U. Mosimann, Swiss Alpine Club, 2013). Owing to the ice dynamics, their corpses will reappear on the glacier surface at some time in the future. In June 2012, human bones were found at the surface of Aletschgletscher, Switzerland, and were identified as being those of three brothers who went missing in March 1926 (Police Report, 2012: Polizeirapport, Kanton Wallis, unpublished and confidential). Since they were found in the ablation area at ˜2400 m a.s.l., they are likely to have disappeared in the accumulation area, several kilometres upstream. Indeed, any object buried in the accumulation area of a glacier is transported via ice flow along a trajectory and reappears at the surface of the ablation area some time later. In the case of fatally injured people this has been observed several times in the past (Reference Krämer, Lehner and ZinkKrämer and others, 1988; Reference Ambach, Ambach, Tributsch, Henn and UnterdorferAmbach and others, 1992; Reference LoreilleLoreille and others, 2010).

The trajectories of ice particles in a glacier have already been computed by glacier modellers, though not for the same purpose as that of this paper. Most of the time, trajectories are computed in order to date the ice. For instance, a stream-function formulation is added to the ice-flow model by Reference Stolle and KilleavyStolle and Killeavy, (1986), while a dating advection equation is used by Reference Mügge, Savvin, Calov and GreveMügge and others, (1999), Reference DurandDurand and others, (2007), Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, (2007) and Martin and others, (2009). In a further application, the accumulation of meltwater is computed by integrating strain heating along trajectories starting at the surface in the accumulation area (Reference Aschwanden and BlatterAschwanden and Blatter, 2005). The forensic aspects of some cases of corpses released from glaciers were studied by Reference Ambach, Ambach, Tributsch, Henn and UnterdorferAmbach and others, (1992). In that study, the authors give a coarse estimate of the length of the trajectory based on a single assessment of the averaged velocity. To our knowledge, no study has ever applied an entire model to the determination of the trajectory of a set of objects, or corpses, in ice.

The goal of this paper is to reconstruct the trajectory along which the corpses of the brothers who disappeared on Aletschgletscher in 1926 were transported, using the time-dependent reconstructed velocity field of the full-Stokes glacier model (Reference Jouvet, Huss, Funk and BlatterJouvet and others, 2011). By integrating the transient ice velocity field backward-in-time, starting from the location of their reappearance at the glacier surface in 2012, we estimate the location where the three brothers are most likely to have disappeared on the glacier in 1926. This paper is organized as follows. First, we give a short summary of what is known about the disappearance of these three mountaineers. Then we present our method for reconstructing the trajectory from the model data of Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011). Following that, we present our results, discuss their sensitivity with respect to the sources of error and provide some indications as to what could have happened to the mountaineers in March 1926. Finally, we offer a perspective on potential future applications of our methodology.

Historical Facts and Data

Grosser Aletschgletscher is the largest glacier in the European Alps. It consists of three sizeable accumulation basins: Grosser Aletschfirn (West), Jungfraufirn (Nord-West) and Ewigschneefeld (Nord), which merge at Konkordiaplatz (Fig. 1). Downstream of the glacier is a 10 km long curved tongue. Ever since 1912, when Jungfraujoch became accessible by railway, Aletschgletscher has been an attractive destination for mountaineers (Reference RossbergRossberg, 1983).

Fig. 1. (a) Map of Grosser Aletschgletscher and the surrounding area. (b) Overview of the uppermost part of Grosser Aletschgletscher. The stars indicate Hollandiahütte on the left, and the location where the skeletons were found on the right.

In the early morning of 4 March 1926, three brothers and a fourth person left the village of Kippel, in the Lötschental, to climb to the Hollandiahütte refuge (Fig. 1), and arrived there around noon. At the refuge they talked to a group of tourists and said they planned to go for a hike up to Konkordiaplatz in the afternoon and return to the refuge in the evening. During the afternoon, the weather deteriorated considerably and a heavy snowstorm broke out, which lasted for 3 days. As the group did not return to the valley the next day as expected, a search party set out in the days that followed. Despite many search operations during the next few weeks, no trace of the group was found. It was assumed that they had died somewhere in the vicinity. Seasonal mass-balance measurements at Jungfraujoch (located ˜5 km from Hollandiahütte, at 3450 m a.s.l.) confirm heavy snowfalls at the beginning of March (Reference Huss and BauderHuss and Bauder, 2009). More importantly, those measurements show that melting during the summer did not remove the snow that fell after 4 March. This means that the bodies were very likely buried immediately and did not reappear after the snowstorm, so it is not surprising they were not found by the search party.

On 27 June 2012, 86 years later, human skeletons (later identified as those of the three brothers) and the remains of old mountaineering equipment were found on Aletschgletscher at the coordinate

in the Swiss referential (Swiss geodetic datum CH1903), which is ˜5 km downstream of Konkordiaplatz on the glacier tongue (Fig. 1). The fourth group member is still missing.

Method

In this paper, (x, y) denote the horizontal coordinates in the Swiss referential while z marks the vertical coordinate; b (x, y) is the bedrock elevation at point (x, y), s (x, y,) is the top surface elevation at point (x, y) and time t, and u(x, t) is the velocity field at point x = (x, y, z) and time t. The trajectory of the particle passing at time on the time interval , containing , is the setx

Where , solves the ordinary differential equation:

with the symbol ‘ indicating the time derivative. Given a trajectory , we call the haed of the trajectory the pair (x(t^h), t^h, such that

By changing ‘inf’ into ‘sup’ in the above equation, we define the foot of the trajectory, which is denoted (x(t ^f), t ^f). Intuitively, the head, (x(t ^h), t ^h), and the foot, (x(t ^f), t ^f), of the trajectory correspond to the position and date of the entrance and the exit of the particle in the glacier, respectively. With this notation, our problem takes the following mathematical form: If we know the foot of the trajectory,

(x(t ^f), t ^f) (corresponding to where and when the skeletons were found), then we want to compute the trajectory and deduce its head, (x(t ^h), t ^h) (corresponding to where and when the bodies were buried).

To compute the trajectory solving Eqn (2), we need the velocity fields over the time period [t ^h, t ^f]. For that, we use the results of the model of Aletschgletscher from 1926 to 2012 given by Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011). For convenience the main features of this model are recalled in the next paragraph.

The motion of ice was described by the nonlinear Stokes equations based on Glen’s law. Ice was assumed to be isothermal, as records of englacial temperature show that Grosser Aletschgletscher is mainly temperate. Consequently, the rate factor, A, in Glen’s law is constant. Sliding of ice over the bedrock below 2400 m a.s.l. was described by a nonlinear Weertman’s law with c as a sliding coefficient, with ice assumed to be fixed to the bedrock elsewhere. This is due to the fact that more sliding is expected in the lowest regions of the glacier, where meltwater production increases. Parameters A = 100 MPa^-3 a^-1 and c = 23 300 MPa^-3 m a^-1 were found to minimize the discrepancy between modelled and measured surface velocities (Table 2). We calculated the surface mass balance in daily time-steps, using a distributed accumulation and temperature index melt model and climate data. The mass balance was included in the model through a transport equation, which describes the mass conservation. Finally, a decoupling algorithm allowed the Stokes problem and the transport equation to be solved using different numerical methods. First, the nonlinear Stokes equations were solved on a fixed, unstructured mesh consisting of tetrahedrons, using a fixed-point algorithm and the finite-element method. Second, the transport equation was solved on a fixed, structured grid of smaller cells, using the volume-of-fluid method. Further details of the numerical methods are given by Jouvet and others (2008, 2009).

The model proved to be capable of reconstructing Aletschgletscher in the 20th century. Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011) initialized the model once in 1880 and ran it until 1999. Here we use a model for the period 1880–2012, re-initialized in 1926, 1957, 1980 and 1999, years for which digital elevation models (DEMs) are available (Reference Bauder, Funk and HussBauder and others, 2007). Such re-initializations prevent cumulative errors from increasing with the length of the modelling period.

Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011) supply the following approximations: the velocity field, u _a , and the upper elevation, s_a , each half-year, or equivalently for the years, a, in {1926, 1926.5, 1927, 1927.5,…,2011.5, 2012}. Note that the u _a are only available at some nodes of an unstructured mesh of volume {x = (x, y, z) ϵ R³, b(x, y) ≤ z ≤ s_a (x, y, a)}, while the s_a are available at the centres of the cells of a regular and rectangular horizontal grid. Using u _a instead of u, x ₂₀₁₂ from Eqn (1) instead of , 2012 instead of , [and 1926, 2012] instead of , we solve the system of Eqn (2), using a fourth-order Runge–Kutta method to obtain an approximation of the trajectory, x(t, x ₂₀₁₂, 2012), along which the corpses were transported. Finally, x(1926, x ₂₀₁₂, 2012) provides an estimate of the location where the corpses were buried in 1926.

Results

Once the trajectory, x(t, x ₂₀₁₂, 2012), was computed, we obtained the following point of immersion:

The distance from x ₁₉₂₆ to the surface in 1926 was

Figure 2 shows an aerial view of the trajectory, x(t, x ₂₀₁₂, 2012), on a topographic map from 1926 to 2012. Figure 3a shows the elevation of the trajectory, the bedrock below the trajectory and the surface above the trajectory with respect to time. Figure 3b and c show the velocity, u(x(t), t), and the position along the flowline, , of the particle following the trajectory with respect to time, respectively.

Fig. 2. Aerial view of the modelled trajectory on a topographic map. The black curve represents the trajectory, x(t, x ₂₀₁₂, 2012), while the black dots indicate the positions of the trajectory in 1926, 2012 and every decade between 1930 and 2010. The rectangle indicates the confidence of the reconstructed location, x(1926, x ₂₀₁₂, 2012), and the star indicates the Hollandiahütte refuge.

Fig. 3. (a) Elevation of the trajectory x(t) = x(t, x ₂₀₁₂, 2012) (i.e. z (t)), of the bedrock b (x (t), y (t)) and of the upper surface s (x)t, y (t), t) with respect to time, t. (b) Velocity across the trajectory, u(x(t), t). (c) Position along the flowline,

It is remarkable that the elevation of the final upstream position of the trajectory (1926 in Fig. 3a) is only 26 m below the top surface elevation. If the model were free of any numerical and physical errors, then the trajectory would exactly reach the surface in 1926. As a consequence, the observed mismatch of 26m is the result of all the approximations in the model. Thus, the depth of x ₁₉₂₆ below the top surface is a valuable indicator in assessing the quality of the model. Moreover, this indicator is global, since it accounts for the entire simulation period and a large part of the ice domain. Regarding the accuracy of the numerical model used to solve the Stokes equations by Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011) (the elements of the mesh are ˜30m high) the accuracy (3 to 5 m) of the DEM of the surface elevation in 1926 (Reference Bauder, Funk and HussBauder and others, 2007), the mismatch can be considered small. That there is only a slight discrepancy further underlines the validity of the model.

According to Figure 3, the corpses were transported ˜10.5 km and lost ˜800 m in altitude from 1926 to 2012 (mean velocity 122 m a^-1). As expected, the corpses were regularly advected deeper into the ice until reaching a maximum depth of ˜250 m in 1980, for a total ice thickness of ˜600 m at that location. With such a column of ice cover, the bodies experienced a pressure of ˜2 MPa, which is corroborated by the deformation observed in the bones that resurfaced in 2012. In 1980, the elevation of the top surface above the bodies was ˜2750 m a.s.l., which is consistent with the equilibrium-line altitudes occurring on Aletschgletscher in the 1980s and 1990s (Reference Huss, Bauder, Funk and HockHuss and others, 2008). After 1980, the path of the corpses was directed mostly towards the glacier surface. According to Figure 2, the trajectory turned to the right at the confluence of Konkordiaplatz from 1980 to 1990, and accelerated at the same time from ˜100 to 200 m a^-1 (Fig. 3).

Our result relies on a number of model approximations, which might contribute to an error in the evaluation of Eqn (3). The source of error are:

1. 1. error regarding the location of x ₂₀₁₂;

2. 2. errors of the surface elevations in 1926 and 2012;

3. 3. mass-balance model errors;

4. 4. ice-flow model errors;

5. 5. discretization errors in solving the ice flow, the mass-conservation equations and Eqn (2);

6. 6. error in the date of 2012 being the year that these corpses emerged at the surface.

For item 2, an error of the surface in 2012 has significantly more impact on the reconstruction of the trajectory than an error of the surface in 1926, since the starting time of the integration is 2012. Additionally, the DEM for the surface in 1926 is known and was used to initialize our model, such that the accuracy of the model was rendered optimal (±3 to 5 m; Reference Bauder, Funk and HussBauder and others, 2007). As a consequence, only the error for the surface in 2012 was investigated. Further, errors of the mass-balance model (item 3) affect only the top surface and have a limited effect on the trajectories of deep ice particles. Thus, the effects of such errors are more pronounced when the particle gets close to the surface. For instance, if the melting is overestimated in our model, then the modelled ice particle would reach the glacier surface too soon. However, the particle was close to the s urface only during the extreme years of 1926 and 2012. For these years, the errors of surfaces have already been taken into account in item 2, so we dropped item 3 from the sensitivity analysis. Finally, we also dropped the effects of the discretization errors (item 5), since they are likely to be negligible compared with the modelling errors (item 4).

For the other sources, we analyzed the sensitivity of our result by modifying a parameter in the model or the location of x ₂₀₁₂ and determining the resulting shift in the upper emergence point of the trajectory. To perform such an analysis, it was convenient to define local frames oriented to the lateral and longitudinal direction of the trajectory. We let (or be the local frame at the head (or foot) of the trajectory, such that x ^h (or x ^f) was the location of the head (or foot), was the projection onto the horizontal plane of the tangent at x ⁱ to the trajectory oriented with the time variable eⁱ _lat was horizontal and orthogonal to . If we perturbed the model parameter by changing (A, C) or the foot of the trajectory, x^f _, by

With , then the haed of the trajectory, x^h, become

With (Fig.4).We defined the distance to the surface (dts) functions as:

Fig. 4. Definition of the local frames at the upstream and downstream ends of the trajectory.

Which provided a more natural indicator than δ^h _z for the vertical deviation. In the following, we estimated the effects of any perturbation on

We first quantify the error in the location of x ₂₀₁₂. The record was made by the police using a conventional GPS device with an accuracy of ˜10 m, and the remains were found roughly within a 10 m wide square (Police Report, 2012: Polizeirapport, Kanton Wallis, unpublished and confidential). For this reason, we recomputed the trajectories, taking as starting locations four points in the vicinity of x ₂₀₁₂: two 15 m forward and backward, and two 15 m to each side (Table 1). As a result, the lateral deviation of the final location of the trajectory never exceeded 57 m and the longitudinal deviation never exceeded 92 m.

The uncertainties of the surface in 2012 and their impact on the reconstructed trajectory are now evaluated. This surface relies on the DEM of the year 1999 (which is ± 1 m accurate; Reference Bauder, Funk and HussBauder and others, 2007), and 13 years of modelling. As a consequence, the uncertainty can be directly evaluated by computing the difference between the elevation of the modelled surface at (x ₂₀₁₂, y ₂₀₁₂) and the exact altitude, z ₂₀₁₂, which is known. It was found that the surface is 8 m above z ₂₀₁₂:

To evaluate the consequences of this error on Eqn (3), we recomputed the trajectory by taking x ₂₀₁₂ + δ_z e _z with δ_z = –4, 4, 8 m as starting locations (Table 1). As a result, the end of the trajectory obtained by starting exactly from the modelled surface (δ_z = 8) deviated significantly, by 33 m laterally and 260 m longitudinally, while its distance to surface, dts _δ , was quite high, 47 m. Although it is expected that two particles that are vertically close to each other may have followed very different trajectories and have distant origins, this deviation was raised by a numerical artefact. Indeed, the top surface of the glacier in the volume-of-fluid method is represented by a thin diffusive layer (Reference Jouvet, Picasso, Rappaz and BlatterJouvet and others, 2008, 2009), which might artificially slow down the particles at the surface. To remove such an artefact, it is necessary to start the integration slightly below the surface. It was found that the trajectories obtained by doing so (with δ _z = 4, 0, – 4 m) are less sensitive to the aforementioned artefact, such that δ _z = ±4 m can be used as bounds in the sensitivity analysis. As a result, the final locations of the trajectories never deviated more than 7 m laterally or 106 m longitudinally, while their distance to surface, dts_δ, remained small.

Since we used the Stokes model, which is the most advanced ice-flow model (Reference Greve and BlatterGreve and Blatter, 2009), no uncertainty due to mechanical simplifications was introduced. However, the model was isothermal. As a consequence, the main source of error was related to the rate factor, A, and the basal sliding, c. Extending the work of Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011), for each A = A_i in

we searched for the optimal parameter c = c_i that minimizes the difference between modelled and averaged measured surface velocities, and found

Among the pairs {(A_i, c_i ), i = 1, 2, 3, 4, 5}, the optimal parameters were (A, c) = (A ₃, c ₃) with a root-mean-square error (RMSE) between measured and simulated velocities of ˜19 m a^-1. However, other pairs also proved to be realistic with ˜25, ˜19, ˜23 and ˜31 m a^-1 RMSE for (A ₁, c ₁), (A ₂, c ₂), (A ₄, c ₄) and (A ₅, c ₅). These pairs represent two kinds of velocity profile: (A ₁, c ₁) and (A ₂, c ₂) account for more sliding and less ice shearing, while (A ₄, c ₄) and (A ₅, c ₅) account for more ice shearing and less sliding (no sliding for (A ₅, c ₅)). Thus, the measured surface velocities do not determine (A, c) uniquely, as they can be reproduced with various ratios of sliding and shearing. Since we know that the corpses of the mountaineers that entered the glacier in 1926 went through the location x ₂₀₁₂ in 2012, we can corroborate the choice of (A, c) with these additional data. For all pairs (A_i, c_j ), the model was run from 1926 to 2012 and the trajectory x(t, x ₂₀₁₂, 2012) was integrated backward in time to obtain the distance to the surface, dtsδ (Table 2), and the final location, x(1926, x ₂₀₁₂, 2012) (Fig. 5). Since x ₂₀₁₂ might have been above or below the modelled surface in 2012, the integration of the trajectory was not initialized at x ₂₀₁₂, but slightly below at its vertical projection on the modelled surface. Table 2 and Figure 5 clearly indicate two categories of parameters. First, the diagonal pairs {(A_i, c_i ), i = 1,…,5} led to trajectories ending in the same basin (Fig. 5). Importantly, the dtsδ and the RMSE were simultaneously the smallest for these pairs (Table 2). This shows that the diagonal pairs {(A_i, c_i ), i = 1,…,5} that were calibrated using velocity measurements, were confirmed by the new data. However, since all dtsδ were small, this indicator cannot be used to choose between the diagonal pairs. In contrast, all other pairs {(A_i, c_j ), i ≠ j} (except (A ₃, c ₁)) led to trajectories ending either in, or at the entrance to, another basin (Fig. 5). Additionally, most of the final positions of trajectories were still deep below the glacier surface, as indicated by the high dtsδ (Table 2). These results confirm that non-diagonal pairs must be excluded. Interestingly, Figure 5 shows that an incorrect parameterization of the ice-flow model changes the result dramatically. Since rate factors A between 80 and 100 MPa^-3 a^-1 are usually recommended in the literature (Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998; Reference GudmundssonGudmundsson, 1999) for comparable mountain glaciers, we exclude the pair (A ₅, c ₅) and restrict our consideration to

Fig. 5. Aerial view, similar to Figure 2. The black dots labelled (i, j) indicate the horizontal position of the final point of the trajectory obtained with model parameters (A_i, C_j ) of Table 2.

Table 1 shows the deviations of the final location of the trajectory obtained with model parameters (A ₁, c ₁), (A ₂, c ₂) and (A ₄, c ₄). If (A ₁, c ₁) (or (A ₂, c ₂)) were chosen instead of (A ₃, c ₃), then the location of the corpses in 1926 would be ˜680 m (or 230 m) downstream. Indeed, since no sliding was assumed above 2400 m a.s.l., reducing A made the ice flow globally slower. In contrast, the parameters (A ₄, c ₄) led to faster flow and to the location of the corpses in 1926 ˜530 m upstream. Note that the deviation is mostly longitudinal, since the lateral deviation never exceeds 175 m.

The skeletons were discovered in June immediately after winter snowmelt and there is no evidence that they had reached the surface earlier. Human bodies are usually preserved (by mummification) when confined in ice (Reference Ambach, Ambach, Tributsch, Henn and UnterdorferAmbach and others, 1992). Since the mountaineers died in March, it is very likely that the bodies were immediately covered by snow and thus remained frozen and intact until they reached the surface. Regarding the time required for the bodies to decompose and the bones to erode, it was established that the corpses were exposed at the surface of the glacier for several seasons. Additionally, the human remains had little chance of being found immediately they came to the surface, since their location, x ₂₀₁₂, lay on the west side of the tongue and far from the east-side summer route usually taken by mountaineers. Nevertheless, Aletschgletscher is popular during the summer season and the west side is also visited from time to time. Thus, it seems unlikely that the corpses emerged more than 10 years prior to 2012 (and even more unlikely it could have been 20 years earlier). To explore the effects of this unknown period at the head of the trajectory, we performed several experiments assuming time intervals of 5, 10, 15 and 20 years between the emergence of the corpses at the surface and their discovery. For instance, if the corpses had emerged in year Y instead of 2012, the bodies would very probably be close to

where x _Y = x(Y, x ₂₀₁₂, 2012) is the point of the trajectory in Y (Fig. 2). Taking Y equal to 2007, 2002, 1997 or 1992, we recomputed the trajectory by taking (), (), (), and () as the foot of the trajectory, integrating backward-in-time over the remaining time period, and checked the deviation of the final location (Table 1). We found that the lateral deviation never exceeded 27 m, while the longitudinal deviation never exceeded 303 m. Note that the distance to the surface, dts_δ, is slight in all cases, so this cannot be used as an indicator to distinguish between assumptions.

Assuming that the bodies of the mountaineers appeared on the surface after 2002, the interval of confidence for x ₁₉₂₆ given in Eqn (3) (obtained by totalling all sources of uncertainty) is [-648, 974Š m longitudinally and [-60, 247] m laterally. Thus, the area where the corpses were probably buried is a 1622 m long and 307 m wide rectangle (Fig. 2), which represents ˜2.5% of the Grosser Aletschfirn basin area and ˜0.6% of the entire Aletschgletscher area.

Probable Causes of Death

Based on the uncertainty range of the immersion location, three causes of death are conceivable: (1) a fall into a crevasse, since a highly crevassed zone can be seen in the upstream area of the uncertainty rectangle of Figure 2; (2) avalanche burial, since the surface topography slope immediately above the likely trajectory head is > 30˚; and (3) exhaustion/freezing to death. Cause (1) is improbable, because three people falling into a crevasse at the same time is virtually impossible. Cause (2) is also extremely unlikely, because in such a case the bodies would have been spread out over a large area, which is not compatible with the small area in which the remains were found. Therefore we conclude that (3) is the most probable cause of death in the case of the three brothers. Additionally, we note that their last position was not on the path joining the refuge to Konkordiaplatz (where they presumably were coming from), but ˜1 km to the north of the refuge (Fig. 2). This supports the theory that the victims became disoriented due to a lack of reference points in conditions of poor visibility. What happened to the fourth group member remains a mystery. Since no trace of him was found in the remains, he probably separated from the group and died elsewhere.

Conclusions and Perspectives

Using the model results published by Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011) and the location where the human remains were found in 2012, we reconstructed the trajectory of the corpses of three brothers who disappeared in 1926 on Aletschgletscher. In particular, the position of their immersion into the glacier – where they probably died – could be determined with great confidence: the area is a ˜1600 m long and ˜300 m wide rectangle, which represents only ˜0.6% of the entire Aletschgletscher area. In addition, the resulting trajectory offers a new and reliable validation of the model presented by Reference Jouvet, Huss, Funk and BlatterJouvet and others, (2011). Indeed, the elevation of the obtained submergence point of the computed trajectory in 1926 is only 26 m beneath the surface. This small discrepancy attests to the ability of the model to reproduce reality. In contrast to local measurements normally used for model verification purposes, the distance to the surface is global in time and space as an indicator of the model quality, since it accounts for a substantial part of the glacier and a long modelling period. Finally, we have shown that it was crucial to accurately calibrate the ice-flow parameters. As illustrated in Figure 5, an incorrect parameterization changes the immersion location from one basin to another, since a bifurcation of ice flows occurs downstream. Our findings proved to be useful in a ˜100 year old investigation. Corroborated with historical facts, our estimate of the location where the brothers died contributes to a better understanding of what might have caused their deaths. Considering all possible causes, we conclude that the loss of all points of reference produced by severe weather conditions very likely led to the death of the mountaineers by exposure.

In this paper, we have reconstructed the trajectory by integrating backward in time a time-dependent velocity field, starting from the known emergence point of the human remains. It would be straightforward to begin with the head and to predict the trajectory by integrating forward-in-time. This would open the method to another application: to predict when and where the bodies of other mountaineers or certain objects will reappear at the surface, when the location and the date of disappearance are known. For instance, the method could also be used to determine the trajectory of the Dakota aircraft of the US Army that crash-landed on the Gauligletscher, Bernese Alps, Switzerland, in November 1946 and is expected to emerge in the coming decades on the top surface of the glacier tongue. In view of the current global warming situation, bodies and objects submerged within Alpine glaciers can be expected to emerge from the ice faster than originally anticipated. A method for establishing the position of various bodies in a glacier, with an estimated place and time of reappearance, would be helpful for identification purposes, and is of immediate interest to the police in charge of such investigations.