Research article 19 Aug 2019
Research article  19 Aug 2019
Glacier thickness estimations of alpine glaciers using data and modeling constraints
 ^{1}Institute of Geophysics, ETH Zurich, Zurich, Switzerland
 ^{2}Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, Zurich, Switzerland
 ^{1}Institute of Geophysics, ETH Zurich, Zurich, Switzerland
 ^{2}Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, Zurich, Switzerland
Correspondence: Hansruedi Maurer (hansruedi.maurer@erdw.ethz.ch)
Hide author detailsCorrespondence: Hansruedi Maurer (hansruedi.maurer@erdw.ethz.ch)
Advanced knowledge of the ice thickness distribution within glaciers is of fundamental importance for several purposes, such as water resource management and the study of the impact of climate change. Ice thicknesses can be modeled using ice surface features, but the resulting models can be prone to considerable uncertainties. Alternatively, it is possible to measure ice thicknesses, for example, with groundpenetrating radar (GPR). Such measurements are typically restricted to a few profiles, with which it is not possible to obtain spatially unaliased subsurface images. We developed the Glacier Thickness Estimation algorithm (GlaTE), which optimally combines modeling results and measured ice thicknesses in an inversion procedure to obtain overall thickness distributions. GlaTE offers the flexibility of being able to add any existing modeling algorithm, and any further constraints can be added in a straightforward manner. Furthermore, it accounts for the uncertainties associated with the individual constraints. Properties and benefits of GlaTE are demonstrated with three case studies performed on different types of alpine glaciers. In all three cases, subsurface models could be found that are consistent with glaciological modeling and GPR data constraints. Since acquiring GPR data on glaciers can be an expensive endeavor, we additionally employed elements of sequential optimized experimental design (SOED) for determining costoptimized GPR survey layouts. The calculated cost–benefit curves indicate that a relatively large amount of data can be acquired before redundant information is collected with any additional profiles, and it becomes increasingly expensive to obtain further information.
Estimating the amount of the glacier ice around the globe is crucial, for example, for sealevel predictions, securing fresh water resources, designing hydropower facilities in highalpine environments and predicting the occurrence of glacierrelated natural hazards. For estimating the overall glacier ice mass and its local distribution, (i) knowledge of the glacier outline, (ii) its surface topography and (iii) the underlying bedrock topography is required. The first two quantities can be observed with aerial and satellite imagery, but the bedrock topography is more difficult to determine.
The conceptually simplest option includes drilling boreholes through the glacier ice (e.g., Iken, 1988). This approach offers groundtruth information, but only a very sparse observation grid can be obtained with realistic efforts. Therefore, geophysical methods have been employed for obtaining more detailed information. Due to the very high electrical resistivity of glacier ice and the relatively high electromagnetic impedance contrast between ice and bedrock material, groundpenetrating radar (GPR) techniques, also referred to as radioecho sounding, have been the primary choice for such investigations (e.g., Evans, 1963). GPR data can either be acquired using groundbased methods (e.g., Watts and England, 1976), or, more efficiently, using fixedwing airplanes (e.g., Steinhage et al., 1999) or helicopters (e.g., Rutishauser et al., 2016).
Despite the powerful capabilities of modern GPR acquisition systems, it is still beyond any practical limits to acquire spatially unaliased 3D data sets. GPR data are therefore collected only along a sparse network of profiles, which leaves considerable uncertainties in the regions between the profiles.
To address this problem, glaciological modeling techniques have been established to relate observable surface parameters to the thickness distribution of ice. One of the earliest concepts was published by Nye (1952). He established a simple relationship between the surface slope and ice thickness. During the past decades, more sophisticated ice thickness modeling techniques have emerged rapidly. Various glaciological constraints, such as mass conservation and/or the relation between basal shear stress and ice thickness, were considered (e.g., Farinotti et al., 2009; Huss and Farinotti, 2012; Clarke et al., 2013; Linsbauer et al., 2012; Morlighem et al., 2011). See Farinotti et al. (2017) for a more complete review of most of the approaches published to date.
Due to inaccuracies of the observed data (GPR measurements, surface topography, etc.) and/or inadequacies in the modeling approaches, modeled ice thicknesses cannot be expected to be perfect. This can be considered by formulating ice thickness estimation as an optimization problem, in which the discrepancies between observed and predicted data are minimized (e.g., Morlighem et al., 2014). In this contribution, we follow an approach similar to Morlighem et al. (2014) but with a different implementation. We introduce the general framework of Glacier Thickness Estimation (GlaTE), with which modeling and data constraints can be combined in an appropriate fashion. After introducing the underlying theory, we demonstrate the performance of the GlaTE inversion procedure with three case studies. In the second part of the paper, we employ elements of GlaTE to address the experimental design problem. Here, we seek a data set to be measured that offers maximum information content at minimal costs. For that purpose, we consider sequentially optimized experimental design (SOED) techniques (e.g., Maurer et al., 2017). The paper concludes with a critical review of potential problems and shortcomings of GlaTE and the associated SOED procedures, and we outline options to address these issues and propose useful extensions of the methodology.
2.1 Theory
The basic idea of GlaTE inversions is to combine observable data with glaciological modeling constraints. A key feature of the algorithm includes appropriate consideration of the uncertainties associated with both constraints. All constraints are formulated, such that they can be integrated into a single system of equations, which can be solved with an appropriate solver.
The first type of constraints includes the GPR data. They can be written in the form of
where h^{est} is a vector including the unknown (estimated) ice thicknesses at M locations (typically defined on a regular grid R on a glacier) and G is a N^{GPR}×M matrix with ones in its main diagonal and zeros everywhere else (N^{GPR}= number of available GPR data points; M= number of elements in h^{est}). The vector h^{GPR} of length N^{GPR} includes the GPRbased thickness estimates. Since the GPR data usually do not coincide with the grid points of R, the values h^{GPR} are obtained by interpolating or extrapolating the GPR data to the nearest grid points of R.
Next, we consider glaciological modeling constraints. In principle, any of the algorithms proposed in the literature can be employed. Here, we closely follow the approach described in Clarke et al. (2013). Input data include a digital terrain model (DTM, defined on R) and the glacier outline.
First, the glacier area is subdivided into socalled flow sheds using the MATLAB TOPOToolbox (Schwanghart and Kuhn, 2010). The subsequent procedure is applied to each flow shed individually (see comments in Clarke et al., 2013, for more information on the flow shed subdivision).
Next, the apparent massbalance is applied, defined as follows:
with $\dot{\mathit{b}}$ being the mass balance rate and $\partial \mathit{h}/\partial t$ the thickness change rate, which is either determined by measuring $\dot{\mathit{b}}$ and $\partial \mathit{h}/\partial t$ or computed via the condition
where Ω_{G} denotes the glacier area (see Farinotti et al., 2009 for more details). In the next step, the flow sheds are partitioned into a prescribed number of elevation zones D_{i} (i=1… number of elevation zones), for which the ice discharge Q_{i} through its lower boundary is computed using
where ${\mathrm{\Omega}}_{{D}_{i}}$ is the area of zone D_{i}. Following Clarke et al. (2013), the basal shear stress τ can then be obtained via the relationship
The parameters n, ρ, g, and A denote the exponent of Glen's flow law, ice density, gravity acceleration, and creep rate factor, respectively (e.g., Cuffey and Patterson, 2010). The factor ξ denotes the creeping contribution (relative to basal sliding) to the ice flux $\left(\mathrm{0}<\mathit{\xi}<\mathrm{1}\right)$, and q is the specific ice discharge ${q}_{\mathrm{i}}={\stackrel{\mathrm{\u203e}}{Q}}_{\mathrm{i}}/{l}_{i}$, where l_{i} is the length of the lower boundary of D_{i} and ${\stackrel{\mathrm{\u203e}}{Q}}_{\mathrm{i}}$ is the average of Q_{i} within D_{i}. Likewise, the angle ϕ represents the surface slope averaged along the lower boundary of D_{i}.
As outlined in Kamb and Echelmeyer (1986), the physics of ice flow can be incorporated into the modeling procedure by applying “longitudinal averaging” of the shear stress (i.e., along the flow direction). We apply this procedure to the results obtained with Eq. (5). Finally, the ice thicknesses ${\widehat{\mathit{h}}}^{\mathrm{glac}}$ (glac stands for glaciological modeling constraints) are obtained using
where τ^{∗} denotes the basal shear stress after longitudinal averaging.
Some of the parameters in Eq. (5) may be subject to considerable uncertainties. For example, the parameter ξ is often poorly known, and it is not guaranteed that the values of the parameters A and n, usually taken from the literature, are accurate. Typically, n is reasonably well constrained but A can vary over orders of magnitudes. Therefore, the overall magnitudes of ${\widehat{\mathit{h}}}^{\mathrm{glac}}$ may be significantly overestimated or underestimated. This can be considered with an additional factor α_{GPR}, yielding
α_{GPR} can be computed with an optimization procedure that minimizes ${\u2225{\mathit{h}}^{\mathrm{GPR}}{\mathit{\alpha}}_{\mathrm{GPR}}{\widehat{\mathit{h}}}^{\mathrm{glac}}\u2225}^{\mathrm{2}}$.
The correction factor α_{GPR} accounts for some inadequacies of Eq. (5), but it is still possible that there are systematic differences between h^{GPR} and h^{glac}. To avoid the resulting inconsistencies, we consider not the absolute value of h^{glac} but instead the spatial gradient of ∇h^{glac} as glaciological constraints, resulting in
where L is a difference operator of dimension M×M .
Further constraints can be imposed via the glacier boundaries that can be determined from aerial or satellite images or ground observations. They are considered in the form of the following equation:
where B is a M×M matrix with ones at appropriate places in its main diagonal.
Depending on the discretization of the glacier models (i.e., the discretization of R), the constraints described above may allow the resulting system of equations to be solved unambiguously. However, in most cases, there will still be a significant underdetermined component, that is, there will be many solutions that explain the data equally well. This requires regularization constraints to be applied (e.g., Menke, 2012). A common strategy for regularizing such problems is to follow Occam's razor, which identifies the “simplest” solution out of the many possible solutions (Constable et al., 1987). Here, we define “simplicity” in terms of structural complexity, that is, we seek a smooth model. This can be achieved via a set of smoothing equations of the form
where S is a M×M smoothing matrix.
All the constraints can now be merged into a single system of equations:
where the parameters λ_{1} to λ_{4} allow a weighting according to the confidence into the individual contributions. The dimension of the system of equations in Eq. (11) can be very large, but the matrices G, L, B and S are all extremely sparse. Therefore, sparse matrix solvers, such as LSQR (Paige and Saunders, 1982) can solve such systems efficiently for h^{est}. The test data sets, described later, yield matrices including up to $\sim \mathrm{320}\phantom{\rule{0.125em}{0ex}}\mathrm{000}\times \mathrm{90}\phantom{\rule{0.125em}{0ex}}\mathrm{000}$ elements. The LSQR algorithm for solving such a matrix required approx. 2 s on a standard PC.
A critical part of the GlaTE inversions includes a proper choice of the weighting parameters λ_{1} to λ_{4}. Parameter λ_{3} is not critical and can be fixed to an appropriate value (e.g., 1.0). The magnitudes of the remaining three parameters must be chosen, such that the system of equations in Eq. (11) is solvable. However, it also needs to be considered that the constraints related to λ_{1} and λ_{2} are subject to significant inaccuracies. It is difficult to predict the accuracy of the modeling constraints, but the accuracy of the GPR data constraints, subsequently denoted as ε^{GPR}, can usually be quantified. Therefore, we have chosen the following strategy.

Initially, we set λ_{1}=1 and choose a low λ_{2} value (i.e., a high λ_{1}∕λ_{2} ratio). Such a ratio indicates a much higher confidence in the GPR data constraints compared with the glaciological modeling constraints. Furthermore, we choose a large value of λ_{4}, which is expected to oversmooth the ice thickness estimates.

With this choice of parameters, a first GlaTE inversion is carried out, and it is checked if a prescribed percentage (e.g., 95 %) of the estimated thicknesses h^{est} matches the GPR data h^{GPR}within their accuracy limits ±ε^{GPR}.

For an overly high value of λ_{4}, it cannot be expected that the prescribed percentage of matching data can be achieved. Therefore, λ_{4} is gradually lowered until the condition, specified in point 2, is met or a prescribed lower threshold of ${\mathit{\lambda}}_{\mathrm{4}}={\mathit{\lambda}}_{\mathrm{4}}^{min}$ is reached. The final smoothing weight, obtained with this procedure, is denoted as ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{4}}$. Since the λ_{1}∕λ_{2} ratio is still large, it is expected that ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{4}}$ is also large because the modeling constraints do not contribute much to the GlaTE inversion. Essentially, a smooth interpolation of the GPR data between the profile lines is performed.

The λ_{1}∕λ_{2} ratio is gradually lowered and step 3 is carried out again (λ_{4} is reset to a high initial value). This iterative procedure is repeated until (i) ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{4}}={\mathit{\lambda}}_{\mathrm{4}}^{min}$ without reaching the prescribed data match, or (ii) the λ_{1}∕λ_{2} ratio has reached a prescribed lower limit.
With decreasing λ_{1}∕λ_{2} ratios, the importance of the glaciological modeling constraints increases, and the contributions of the smoothing constraints need to be lowered to achieve the prescribed data match. Below a certain λ_{1}∕λ_{2} ratio, it will likely no longer be possible to fit a sufficiently large percentage of the data within the limits ±ε^{GPR}, even when ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{4}}={\mathit{\lambda}}_{\mathrm{4}}^{min}$. If this is not the case, the λ_{1}∕λ_{2} ratio could be lowered to an arbitrary low level, but if the confidence in the glaciological modeling constraints is rather limited, it is possible to define a lower threshold where the GlaTE inversion would stop, even when ${\stackrel{\mathrm{\u203e}}{\mathit{\lambda}}}_{\mathrm{4}}>{\mathit{\lambda}}_{\mathrm{4}}^{min}$.
With such a strategy, it is possible to achieve several desirable features of glacier thickness estimations, namely

the GPR data are fitted only within the prescribed accuracy limits and no overfitting is performed,

the contribution of the glaciological constraints are maximized, and

the influence of the (unphysical) smoothing constraints is minimized.
2.2 Performance tests
For testing the GlaTE inversion algorithm, we investigated glacier ice thickness at three sites in the Swiss Alps (Fig. 1 and Table 1). The first site is Vadret da Morteratsch (Fig. 1a). Lying at altitudes between 2050 and 4000 m a.s.l. (Zekollari et al., 2013), the glacier has a typical valleyglacier shape and is located in the Engadin region of Switzerland. In 2015, the tributary glacier Vadret Pers, on the eastern part of the main glacier, detached from the main trunk of Vadret da Morteratsch, but we continue to treat both glaciers as a connected system, since the last available outline of the glaciers in 2015 shows the remnant of the former connection. In 2010, the glacier system covered an area of ≈15 km^{2} and had a length of ≈7.4 km.
The second site, Glacier de la Plaine Morte (2400–3000 m a.s.l., Fig. 1b), is the largest plateau glacier in the European Alps (Huss et al., 2013). The surface slope is shallow with an average slope angle of about 6^{∘} and a short glacier tongue draining towards the north.
The third site is a cluster of small valley flanks and cirquetype glaciers on the eastern flank of the Matter Valley (Fig. 1c) below the Dom peak. From north to south, the glaciers are named Hohbärggletscher, Festigletscher, Kingletscher and Weingartengletscher. The Hohbärggletscher is the largest (2800–4500 m a.s.l.) and longest of the group. The individual glaciers were treated as individual flow sheds during the data analysis, but the α_{GPR} factor was determined for the entire Dom area.
For all sites, the recorded GPR profiles are shown in Fig. 1. Most of the data were recorded with the dual polarization system AIRETH (Langhammer et al., 2018). A grid of profiles was acquired in 2016 on the Glacier de la Plaine Morte and in 2017 on the Vadret da Morteratsch and in the Dom Region. The data were processed as described in Grab et al. (2018), and the bedrock depths and the corresponding ice thicknesses were obtained from the migrated GPR images.
As input data for the glacier models, surface topography and an outline of the individual glaciers was required. As surface topography, we used the swissALTID3D (DTM, Digital Terrain Model Release 2017 © swisstopo, JD100042). The most recent version, covering the individual glaciers, was extracted and downsampled to 10 m resolution. The outline represents the extension of the glacier in 2015–2016. DTM and glacier outlines are displayed in Fig. 1. In accordance with Farinotti et al. (2009), we employed mass balance gradients of 0.05 and 0.09 in the accumulation and ablation zones, respectively.
As an appropriate measure of the accuracy of the GPR data, we considered a relative (depthdependent) quantity ${\mathit{\epsilon}}^{\mathrm{GPR}}=\u2225{\mathit{h}}^{\mathrm{glac}}{\mathit{h}}^{\mathrm{GPR}}\u2225/\left({\mathit{h}}^{\mathrm{GPR}}+{h}^{min}\right)$, where h^{min} is a minimum thickness to avoid unreasonably large relative errors at shallow depths. Values of ε^{GPR}=0.05 and ${h}^{min}=\mathrm{5.0}$ were judged to be adequate. For all three data sets, we employed ${\mathit{\lambda}}_{\mathrm{4}}^{min}=\mathrm{4.0}$ and the prescribed data fit was 95 %. This could be achieved with a minimum λ_{1}∕λ_{2} ratio of 3.0 (initial values for λ_{4} and λ_{1}∕λ_{2} were 50.0 and 5.0).
Figure 2 shows the ice thicknesses distributions, when only glaciological constraints are applied (h^{glac}, Fig. 2a), when only GPR constraints are considered (h^{GPR}, Fig. 2b), results from the GlaTE algorithm (h^{est}, Fig. 2c) and considering the difference between h^{glac} and h^{est} (Fig. 2d). In Fig. 2b, the thicknesses were obtained by natural neighbor interpolation from the GPR data. Since no extrapolation was performed, not all glacierized regions have an ice thickness estimate. h^{glac} and h^{GPR} exhibit increased thicknesses in the western glacier, but only the glaciological constraints indicate an overdeepening in the eastern glacier, thereby indicating that the two models are inconsistent. The results from the GlaTE inversion (h^{est}, Fig. 2c) demonstrate that it is possible to find a smooth model that satisfies both the glaciological and the GPR data constraints.
The corresponding results for the Glacier de la Plaine Morte are shown in Fig. 3. The glaciological model suggests a deep isolated trough slightly east of the center (Fig. 3a). This is not supported by the GPR data, which instead indicate a larger E–Woriented elongated zone of increased thickness (Fig. 3b). Such a feature is also contained in the GlaTE inversion results (Fig. 3c). Furthermore, the glaciological model in Fig. 3a overestimates the ice thickness in the northeastern part of the glacier.
Results from the Dom region show a relatively good match between the glaciological model (Fig. 4a) and the GlaTE inversion result (Fig. 4c). The glaciological model tends to underestimate the maximum thickness in the center of the glacier tongues and overestimate the thickness towards the edges (Fig. 4d). The isolated trough structures (ice thickness >200 m) in the northernmost glacier in the glaciological model (Fig. 4a) are only partially supported by the GPR data (Fig. 4b) and the GlaTE inversion (Fig. 4c). In the southernmost Weingartengletscher, no data constraints exist (Fig. 4b). The nonzero differences in this part (Fig. 4d) are the result of the smoothing constraints. Here, the thickness estimates from the glaciological model are thus more trustworthy.
All the investigations, described in Sect. 2, were based on existing GPR data. Their experimental layouts were designed heuristically using experience from prior surveys. Once a glacier model has been established, one may realize that another GPR survey layout may have provided better information. Therefore, a dense survey grid, as employed for 3D seismic reflection campaigns for hydrocarbon exploration (e.g., Vermeer, 2003) would be the best choice. This, however, would exceed by far the budgets typically available for glacier investigations.
Optimizing the glaciological constraints with only a limited number of GPR data is a chickenandegg problem: identifying the most useful GPR data to be added would require knowledge on where the true ice thickness distribution deviates most from the distribution in the glaciological model, but this would require advanced prior knowledge about the ice thickness that one wants to measure. The problem can be tackled nevertheless by making some specific assumptions (see below).
With our investigations, we address the following questions.

Was the experimental geometry and the amount of data acquired in the three investigation areas adequate?

Do better experimental layouts exist for constraining the ice thicknesses in a costoptimized manner?

Can some general recommendations be made for designing helicopterborne GPR surveys on glaciers?
Due to the lack of knowledge on the true ice thicknesses, we assumed that the GlaTE inversion results, shown in Figs. 2, 3 and 4 are a good proxy for the actual thickness distributions. Without GPR data, the state of knowledge is represented by the glaciological model (Figs. 2a, 3a and 4a). For these models, only 12 % (Morteratsch), 8 % (Glacier de la Plaine Morte) and 14 % (Dom) of the GPR data constraints satisfy the condition $\u2225{\mathit{h}}^{\mathrm{glac}}{\mathit{h}}^{\mathrm{GPR}}\u2225/\left({\mathit{h}}^{\mathrm{GPR}}+{h}^{min}\right)<{\mathit{\epsilon}}^{\mathrm{GPR}}$, and the average ice thickness misfits over the entire glacier area (mean(h^{glac}−h^{true})) (h^{true}= “true” model) are 22, 32 and 23 m for the three data sets, respectively. It should be noted that the glaciological models h^{glac} have been calibrated with α_{GPR}. If no GPR data were to be available, the performance of the glaciological models would have been even worse.
Subsequently, we analyzed which of the profiles j (j=1… nprof, number of profiles) causes the largest discrepancies between h^{GPR} and h^{glac}. For that purpose we define
where index i runs over all n_{j} data points of profile j. ${h}_{ij}^{\mathrm{GPR}}$ and ${h}_{ij}^{\mathrm{glac}}$ represent the measured and modeled ice thickness at data point i of profile j, respectively. The function P is defined as
Since longer profiles would be associated with higher (monetary) data acquisition costs, the discrepancy ${d}_{\mathrm{1}}^{\mathrm{cost}}$ is normalized with a cost factor c_{j}, defined as
where len_{j} represents the length of profile j. This cost function assumes that the acquisition costs increase linearly with profile length, which is realistic because the helicopter costs are typically charged per minute of flight time. To avoid overly short profiles dominating ${d}_{\mathrm{1}}^{\mathrm{cost}}$, the assumption was made that profiles with a length <200 m would incur the same costs (for such short profiles the flight time is typically governed by positioning the helicopter at the starting point of a profile).
The profile associated with the largest discrepancy ${d}_{\mathrm{1}}^{\mathrm{cost}}$ is expected to offer the largest amount of additional information per unit cost. In this virtual experiment, we assumed that one would acquire this profile and subsequently perform a GlaTE inversion, yielding an improved model ${\mathit{h}}^{{\mathrm{est}}_{k}}$. Index k indicates the actual state of the experimental design, that is, k is equal to 1, when adding the first profile. Then, the next profile line to be acquired is identified using
Repeated application of Eq. (15) identifies an optimized sequence for how the profiles should be acquired. Figure 5a, c and e show the evolution of what we call the “data fit curve”, that is, the evolution of
with
For the Vadret da Morteratsch data, there is an approximately linear increase in the data fit curve. Likewise, we observe a corresponding linear decrease in the average model misfit. As discussed in Maurer et al. (2010), cost–benefit curves, such as the d^{fit} graphs in Fig. 5, typically enter into the area of diminishing returns at some stage, that is, the curves exhibit a characteristic kink and flatten out at larger numbers of profiles. This indicates that it becomes increasingly expensive to obtain additional information. The curve in Fig. 5a therefore indicates that the area of diminishing returns was not reached during the Vadret da Morteratsch campaigns and that it would have been useful to acquire more profiles. In contrast, the d^{fit} and average misfit curves for the Glacier de la Plaine Morte and Dom regions (Fig. 5c and e) start flattening out, although we do not observe a characteristic kink in the curves. This indicates that it would have become increasingly expensive to obtain a more accurate ice thickness distribution for the Glacier de la Plaine Morte and Dom field sites.
Figures 6 to 8 show examples of model misfit plots $\left({\mathit{h}}^{{\mathrm{est}}_{k}}{\mathit{h}}^{\mathrm{true}}\right)$ superimposed with the selected profile lines. The corresponding stages of the selection procedure are indicated with dashed black lines in Fig. 5a, c and e. For Vadret da Morteratsch, profiles are selected preferentially in the western part of the glacier because the model fit is already quite good in the eastern region. For the Glacier de la Plaine Morte and Dom regions, no obvious selection patterns can be recognized.
A major limitation of this design experiment is that the true model and the recorded GPR profiles have a strong dependency. When all profiles of a particular region are selected, there is a perfect match between ${\mathit{h}}_{k}^{\mathrm{est}}$ and h^{true}. However, this is the result of our choice of the true model, and thus it does not indicate that this data set is optimal.
To reduce this dependency at least partially, we have generated synthetic data sets that are covering all glacierized areas with a dense grid. We assumed a line spacing of 100 m and an inline sampling interval of 0.5 m, which is representative for the helicopterborne GPR data that we acquired. With such a comprehensive data set, the experimental design procedure should have more flexibility to choose costoptimized suites of profiles.
The resulting cost–benefit curves are shown in Fig. 5b, d and f. As expected, the curves start flattening out after selecting a sufficiently large number of profiles. For Vadret da Morteratsch (Fig. 5b), it seems to be worthwhile to acquire more than the 43 profiles acquired during the actual experiment. After about 70 profiles, the curve starts flattening out. The curves for Glacier de la Plaine Morte (Fig. 5d) indicate clearly that acquiring a larger number of profiles would have been beneficial. After adding about 40 profiles, the d^{fit} curve starts flattening out. For the Dom region, the amount of profiles chosen for the actual survey seems to be adequate (Fig. 5f). After approx. 40 profiles, the curve flattens out.
Using the d^{fit} curves in Fig. 5 seems to be a good option for selecting an appropriate number of profiles, but it is also insightful to consider the associated model misfit curves. Figure 5b, d and f indicate that the average thickness misfit typically approaches a low level before the d^{fit} curves start flattening out.
For the experimental design with the synthetic data, Figs. 9 to 11 show examples of model misfit plots $\left({\mathit{h}}^{{\mathrm{est}}_{k}}{\mathit{h}}^{\mathrm{true}}\right)$ superimposed with the selected profile lines. In contrast to the selection based on observed data from the Vadret da Morteratsch (Fig. 6), the design based on the dense synthetic grid (Fig. 9) yields a better balance of profiles among the eastern and western portions of the glacier. This is the consequence of the larger flexibility of choosing profiles with the dense grid. For Glacier de la Plaine Morte (Fig. 10), it is interesting to note that almost exclusively N–Soriented profiles were chosen. In contrast, predominantly E–Woriented profiles were chosen for the Dom region (Fig. 11). Both observations are governed primarily by the cost factor c_{j} in Eq. (15).
The GlaTE inversion scheme presented in this paper offers numerous beneficial features. A benchmark for its capabilities, compared with other methods, will be evaluated in the framework of the ITMIX2 initiative, which is currently ongoing.
Its main advantage is its versatility, as there are several parameters by which the algorithm can be tuned to the peculiarities of a particular investigation area. However, this is also one of the method's major drawbacks, since the choice of the control parameters may include a considerable amount of subjectivity. This applies primarily to the choice of the weighting factors λ_{1}, λ_{2} and λ_{4}. We consider our strategy for determining these factors to be adequate, but other options may work equally well.
Another potential problem is the determination of the scaling factor α_{GPR} in Eq. (7). It is largely dependent on the available GPR data, and it is assumed that the GPR profiles have a good areal coverage, which might not be always the case. If values for α_{GPR} become available for a large number of glaciers, a statistical analysis might be used to correlate the values with specific features of the glaciers (e.g., average slope, elevation above sea level, size or shape of the glacier, exposure, etc.). This may be helpful in areas where the GPR data coverage is poor or even nonexistent.
In principle, any observations (e.g., boreholes) can be employed as data constraints in Eq. (1), but GPR measurements are typically the main source of information. Migration of the GPR data allows the bedrock reflections to be imaged at the correct position and slope along a profile, but it is possible that the reflections originated from locations away from the profile lines (offplane reflections). This may cause systematic errors affecting the reliability of the results. We note, however, that this is not a problem specific to GlaTE, but rather a general issue affecting GPR data acquired on a sparse grid.
As mentioned in Sect. 2, the system of equations in Eq. (11) can be augmented by any linear constraint. An obvious and, in our view, particularly useful set of constraints could be offered by surface displacement measurements. They can be obtained from differential satellite images and offer full coverage over a glacier. Such constraints could possibly substitute the smoothness constraints in Eq. (11) with a physically more meaningful quantity.
Despite the limitations of our experimental design approach, we judge that our results provided useful insights for designing GPR experiments and that some answers to the questions posed in Sect. 3 can be provided.

Was the experimental geometry and the amount of data acquired in the three investigation areas adequate?
The cost–benefit curves in Fig. 5 indicate that, at least for Vadret da Morteratsch, it would have been useful to acquire more data.

Do better experimental layouts exist for constraining the ice thicknesses in a costoptimized manner?
The experimental layouts in Figs. 6 to 11 do not provide unexpected features but indicate that acquiring a larger number of shorter profiles, instead of recording a few long ones, could be beneficial, but it should be noted that we do not take into account the flight time required to move to the next profile. This could be significant on glaciers with steep mountain flanks.

Can some general recommendations for designing helicopterborne GPR surveys on glaciers be made?
Based on our results, it is difficult to offer general recommendations. For estimating the overall amount of data to be collected, the cost–benefit curves are most indicative. However, in our case studies they do not flatten out clearly, thereby indicating that it would be worthwhile to acquire more data. When highprecision ice thickness maps are required, it is therefore advisable to acquire as much data as can be afforded.
It is common practice to acquire crossing profiles, but from the experimental layouts, shown in Fig. 10, it could be concluded that it is not necessary to acquire a large amount of crossing profiles. From a practical point of view, this recommendation cannot be fully supported. When the signaltonoise ratio of the GPR profiles is poor, it can be difficult to identify the bedrock reflections unambiguously. Due to the importance of crossing profiles, it has been to be judged worthwhile to extend the cost function of the experimental design algorithm, such that crossing profiles are favored.
It is not realistic to adopt a realtime experimental design strategy (i.e., choosing the next profile based on the results of the previously acquired data), as assumed in our virtual experiments in Sect. 3. However, if logistically feasible, it might be useful to employ a twostep acquisition strategy. Initially, only a few profiles could be acquired. After analyzing these data sets, regions where large discrepancies between h^{est} and h^{glac} exist can be identified and thus a suitable set of additional profiles could be acquired with a second campaign.
A MATLAB implementation of GlaTE and the test data sets, shown in this paper, can be downloaded from https://gitlab.com/hmaurer/glate (last access: 29 June 2019).
HM performed the conceptual work, LL and HM performed the coding and carried out the analyses, MG and AB provided data and intellectual input, and HM wrote the paper.
The authors declare that they have no conflict of interest.
We thank Patrick Lathion, Philipp Schaer, and Kevin Délèze from GEOSAT SA; Patrick Fauchère from Air Glacier; Hansueli Bärfuss from HeliBernina; and Lasse Rabenstein and Lino Schmid for acquiring the data. Furthermore, we thank Matthias Huss for fruitful discussions and Daniel Farinotti for an insightful inhouse review, which improved the clarity of the manuscript. The manuscript was further improved by the valuable comments of Ben Pelto, Douglas Brinkerhoff and Fabien Maussion. Finally, we acknowledge the authors of the publicdomain MATLAB TopoToolbox, which proved to be very useful for this project.
This research has been supported by the ETH Zurich (grant no. ETH15 132), the SCCERSOE (Swiss Competence Center for Energy Research, Supply of Electricity) and the Swiss Geophysical Commission.
This paper was edited by Valentina Radic and reviewed by Douglas Brinkerhoff, Fabien Maussion and Ben Pelto.
Clarke, G. K., Anslow, F. S., Jarosch, A. H., Radić, V., Menounos, B., Bolch, T., and Berthier, E.: Ice volume and subglacial topography for western Canadian glaciers from mass balance fields, thinning rates, and a bed stress model, J. Climate, 26, 4282–4303, 2013.
Constable, S. C., Parker, R. L., and Constable, C. G.: Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289–300, 1987.
Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, Amsterdam, the Netherlands, 2010.
Evans, S.: Radio techniques for the measurement of ice thickness, Polar Rec., 11, 406–410, 1963.
Farinotti, D., Huss, M., Bauder, A., Funk, M., and Truffer, M.: A method to estimate the ice volume and icethickness distribution of alpine glaciers, J. Glaciol., 55, 422–430, 2009.
Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., GilletChaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., Ramsankaran, R., Reerink, T. J., Sanchez, O., Stentoft, P. A., Singh Kumari, S., van Pelt, W. J. J., Anderson, B., Benham, T., Binder, D., Dowdeswell, J. A., Fischer, A., Helfricht, K., Kutuzov, S., Lavrentiev, I., McNabb, R., Gudmundsson, G. H., Li, H., and Andreassen, L. M.: How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment, The Cryosphere, 11, 949–970, https://doi.org/10.5194/tc119492017, 2017.
Grab, M., Bauder, A., Ammann, F., Langhammer, L., Hellmann, S., Church, G., Schmid, L., Rabenstein, L., and Maurer, H.: Ice volume estimates of Swiss glaciers using helicopterborne GPR – an example from the Glacier de la Plaine Morte, 2018 17th International Conference on Ground Penetrating Radar (GPR), June 2018, Rapperswil, Switzerland, 1–4, 2018.
Huss, M. and Farinotti, D.: Distributed ice thickness and volume of all glaciers around the globe, J. Geophys. Res.Earth, 117, F04010, https://doi.org/10.1029/2012JF002523, 2012.
Huss, M., Voinesco, A., and Hoelzle, M.: Implications of climate change on Glacier de la Plaine Morte, Switzerland, Geogr. Helv., 68, 227–237, https://doi.org/10.5194/gh682272013, 2013.
Iken, A.: Adaption of the hotwaterdrilling method for drilling to great depth, Mitteilungen der Versuchsanstalt fur Wasserbau, Hydrologie und Glaziologie an der Eidgenossischen Technischen Hochschule Zurich, Zurich, Switzerlnad, 211–229, 1988.
Kamb, B. and Echelmeyer, K. A.: Stressgradient coupling in glacier flow: I. Longitudinal averaging of the influence of ice thickness and surface slope, J. Glaciol., 32, 267–284, 1986.
Langhammer, L., Rabenstein, L., Schmid, L., Bauder, A., Grab, M., Schaer, P., and Maurer, H.: Glacier bed surveying with helicopterborne dualpolarization groundpenetrating radar, J. Glaciol., 65, 1–13, https://doi.org/10.1017/jog.2018.99, 2018.
Linsbauer, A., Paul, F., and Haeberli, W.: Modeling glacier thickness distribution and bed topography over entire mountain ranges with GlabTop: Application of a fast and robust approach, J. Geophys. Res.Earth, 117, F03007, https://doi.org/10.1029/2011JF002313, 2012.
Maurer, H., Curtis, A., and Boerner, D. E.: Recent advances in optimized geophysical survey design, Geophysics, 75, 75A177–175A194, 2010.
Maurer, H., Nuber, A., Martiartu, N. K., Reiser, F., Boehm, C., Manukyan, E., Schmelzbach, C., and Fichtner, A.: Optimized Experimental Design in the Context of Seismic Full Waveform Inversion and Seismic Waveform Imaging, Adv. Geophys., 58, 1–45, 2017.
Menke, W.: Geophysical data analysis : discrete inverse theory, Matlab ed., Academic Press, Waltham, MA, USA, 293 pp., 2012.
Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: A mass conservation approach for mapping glacier ice thickness, Geophys. Res. Lett., 38, L19503, https://doi.org/10.1029/2011GL048659, 2011.
Morlighem, M., Rignot, E., Mouginot, J., Seroussi, H., and Larour, E.: Highresolution icethickness mapping in South Greenland, Ann. Glaciol., 55, 64–70, 2014.
Nye, J.: A method of calculating the thicknesses of the icesheets, Nature, 169, 529–530, 1952.
Paige, C. C. and Saunders, M. A.: Lsqr – an Algorithm for Sparse LinearEquations and Sparse LeastSquares, ACM T. Math. Software, 8, 43–71, 1982.
Rutishauser, A., Maurer, H., and Bauder, A.: Helicopterborne groundpenetrating radar investigations on temperate alpine glaciers: A comparison of different systems and their abilities for bedrock mapping, Geophysics, 81, WA119–WA129, https://doi.org/10.1190/geo20150144.1, 2016.
Schwanghart, W. and Kuhn, N. J.: TopoToolbox: A set of Matlab functions for topographic analysis, Environ. Modell. Softw., 25, 770–781, 2010.
Steinhage, D., Nixdorf, U., Meyer, U., and Miller, H.: New maps of the ice thickness and subglacial topography in Dronning Maud Land, Antarctica, determined by means of airborne radioecho sounding, Ann. Glaciol., 29, 267–272, 1999.
Vermeer, G. J.: 3d seismic survey design optimization, The Leading Edge, 22, 934–941, 2003.
Watts, R. D. and England, A. W.: Radioecho sounding of temperate glaciers: ice properties and sounder design criteria, J. Glaciol., 17, 39–48, 1976.
Zekollari, H., Huybrechts, P., Fürst, J., Rybak, O., and Eisen, O.: Calibration of a higherorder 3D iceflow model of the Morteratsch glacier complex, Engadin, Switzerland, Ann. Glaciol., 54, 343–351, 2013.