2.1 Glacier outlines and local topography

The glacier outlines used in this study are those defined in region 1 of the RGIv6. Four glaciers (Columbia, Grand Pacific, Hubbard and Sawyer glaciers) were merged with their respective pair branches (West Columbia, Ferris, Valerie and West Sawyer glaciers) into a single outline. A local map projection is defined for each glacier in the inventory following the methods described in Maussion et al. (2019a). A transverse Mercator projection is used, centred on the glacier in order to conserve distances, area and angles. Then, topographical data are chosen automatically depending on the glacier's location and interpolated to the local grid. For this study we used a combination of the Shuttle Radar Topography Mission (SRTM) 90 m Digital Elevation Database v4.1 (Jarvis et al., 2008) for all latitudes below 60^∘ N and the Viewfinder Panoramas DEM3 product (90 m, http://viewfinderpanoramas.org/dem3.html, last access: 8 October 2019) for higher latitudes. For Columbia Glacier, we used the digital elevation model (DEM) from the Ice Thickness Models Intercomparison eXperiment (Farinotti et al., 2017, ITMIX) instead.¹ All datasets are re-sampled to a resolution depending on glacier size (Maussion et al., 2019a) and smoothed with a Gaussian filter of 250 m radius.

2.2 Glacier flow lines, catchment areas and widths

The glacier centre lines are computed following an automated method based on the approach of Kienholz et al. (2014). Figure 2a illustrates an example of this geometrical algorithm applied to Columbia Glacier. The centre lines are then filtered and slightly adapted to represent glacier flow lines with a fixed grid spacing (Fig. 2c). The geometrical widths along the flow lines are obtained by intersecting the normals at each grid point with the glacier outlines and the tributaries' catchment areas. Each tributary and the main flow line has a catchment area, which is then used to correct the geometrical widths. This process assures that the flow line representation of the glacier is in close accordance with the actual altitude–area distribution of the glacier. The width of the calving front, therefore, is obtained from a geometric first guess multiplied by a correction factor. This may lead to uncertainties in the frontal ablation computations, as discussed in Sect. 5.

Figure 2Columbia and LeConte glaciers model workflow: (a, b) topographical data preprocessing and computation of the flow lines, (c) width correction according to catchment areas and altitude–area distribution, and (d) thickness distribution before accounting for frontal ablation.

2.3 Regional frontal ablation estimates

Frontal ablation for 27 marine-terminating glaciers presented by McNabb et al. (2015) are used to compare the results of the model and calibrate the calving constant of proportionality k. These estimates were calculated from satellite-derived ice velocities and modelled estimates of glacier ice thickness.

2.4 Climate data and mass balance

The mass balance (MB) model implemented in OGGM uses monthly time series of temperature and precipitation. The current default is to use the gridded dataset Climatic Research Unit time series (CRU TS) v4.01 (Harris et al., 2014), which covers the period of 1901–2015 with a 0.5^∘ resolution. This raw, coarse dataset is downscaled to a higher-resolution grid (CRU CL v2.0 at 10^′ resolution, New et al., 2002), following the anomaly mapping approach described in Maussion et al. (2019a), allowing OGGM to have an elevation-dependent climate dataset from which the temperature and precipitation at each elevation of the glacier are computed and then converted to the local temperature according to a temperature gradient (default: 6.5 K km⁻¹). No vertical gradient is applied to precipitation, but a correction factor p_f=2.5 is applied to the original CRU time series (see Maussion et al., 2019a, Appendix A for more information). The MB model (see Sect. 3.2) is calibrated with direct observations of the annual surface mass balance (SMB). For this, OGGM uses reference mass balance data from the World Glacier Monitoring Service (WGMS, 2017) and the links to the respective RGI polygons assembled by Maussion (2017).

3.1 Ice thickness

The method of estimating ice thickness from mass turnover and principles of ice-flow dynamics in glaciers go back to Budd and Allison (1975) and Rasmussen (1988), whose ideas were further developed by Fastook et al. (1995) and Farinotti et al. (2009). The latter aims to estimate ice thickness distribution from a given glacier surface topography, which can be achieved assuming that the mass balance distribution should be balanced by the ice-flux divergence. This method has been modified in OGGM in order to implement a new ice thickness inversion procedure physically consistent with the flow line representation of the glaciers and taking advantage of the mass balance calibration procedure of OGGM (see below).

The flux of ice q (m³ s⁻¹) through a glacier cross section of area S (m²) is defined as follows:

with $\stackrel{\mathrm{‾}}{u}$ being the average cross section velocity (m s⁻¹). By applying the well-known shallow-ice approximation (Hutter, 1981, 1983; Cuffey and Paterson, 2010; Oerlemans, 1997) and making use of the Glen's ice flow law, we compute the depth-integrated centre-line velocity u of the cross section with the following equation:

with A being the ice creep parameter (which has a default value of $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{24}}$ s⁻¹ Pa⁻³), n the exponent of Glen's flow law (default: n=3), h₀ the centre-line ice thickness (m) and τ the basal shear stress, defined as follows:

with ρ the ice density (900 kg m⁻³), g the gravitational acceleration (9.81 m s⁻²) and α the surface slope (computed along the centre line). Optionally, a sliding velocity u_s can be added to the deformation velocity to account for basal sliding, using the following parametrization (Oerlemans, 1997; Budd et al., 1979):

with f_s a sliding parameter (default: $\mathrm{5.7}\times {\mathrm{10}}^{-\mathrm{20}}$ m⁻² s⁻¹ Pa⁻³). We then assume that the centre-line velocity is equal to the average section velocity ($\stackrel{\mathrm{‾}}{u}\approx u$), which in absence of lateral drag is correct for a rectangular bed shape but is not in the parabolic case, where we neglect the variations in the shear stress (and u) along the parabola. In the parabolic case and with N=3, this results in a section velocity overestimation of a factor of 315∕128 (approx. 2.46) in comparison to the section velocity obtained by integrating the shallow-ice velocity over the parabola. We proceed with this approximation because (i) this factor cannot be computed analytically for any other non-integer value of Glen N or for other bed shapes (e.g. trapezoidal) and (ii) the uncertainties about the true shape of the bed would make the model very sensitive to this choice. The computed flux in OGGM, however, does vary by a factor of 2∕3 depending on whether one assumes a parabolic ($S=\frac{\mathrm{2}}{\mathrm{3}}hw$) or rectangular (S=hw, with w being the glacier width) bed shape. The default in OGGM is to use a parabolic bed shape, unless the section touches a neighbouring catchment or neighbouring glacier (ice divides, computed from the RGI). For the last five grid points of tidewater glaciers, the bed shape is also assumed to be rectangular. Singularities with flat areas are avoided since the constructed flow lines are not allowed to have a local slope α below a certain threshold (default: 1.5^∘; see Maussion et al., 2019a).

Following the approach described in Maussion et al. (2019a), q can be estimated from the mass balance field of a glacier. If u and q are known, S and the local ice thickness h (m) can also be computed by making some assumptions about the geometry of the bed and by solving Eq. (1). This equation becomes a polynomial in h of degree 5 with only one root in ℝ₊, easily computable for each grid point.

3.2 Mass balance and ice flux q

OGGM's mass balance model is an extension of the model proposed by Marzeion et al. (2012) and adapted in Maussion et al. (2019a) to calculate the mass balance of each flow line grid point for every month, using the CRU climatological series as boundary condition. The equation governing the mass balance is that of a traditional temperature index melt model. The monthly mass balance m_i (kg m⁻² s⁻¹) at elevation z is computed as follows:

where $P_i^{\mathrm{solid}}$ is the monthly solid precipitation, p_f a global precipitation correction factor, T_i the monthly temperature and T_melt is the monthly mean air temperature above which ice melt is assumed to occur (default: -1^∘C). Solid precipitation is computed as a fraction of the total precipitation: 100 % solid if T_i≤T_solid (default: 0 ^∘C), 0 % if T_i≥T_liquid (default: 2 ^∘C) and linearly interpolated in between. The parameter μ^* indicates the temperature sensitivity of the glacier, and it needs to be calibrated. In a nutshell, the MB calibration consists of searching a 31-year climate period in the past during which the glacier would have been in equilibrium while keeping its modern-time geometry, implying that the mass balance of the glacier during that period in time m₃₁(t) is equal to zero, with m₃₁(t) being the glacier integrated mass balance computed for a 31-year period centred around the year t (e.g. $t^{*}=\mathrm{1962}$ for most glaciers in Alaska) and for a constant glacier geometry fixed at the RGI outline's date (e.g. 2009 for Columbia Glacier). It should be noted that the mass balance calibration in OGGM excludes MB measurements from tidewater glaciers as reference data, for reasons described below.

This “equilibrium mass balance” (m₃₁(t)) is then assumed to be equal to the “apparent mass balance” ($m=\dot{m}-\mathit{\rho }\frac{\partial h}{\partial t}$) as defined by Farinotti et al. (2009), where the flux of ice q through a glacier catchment area (Ω) is defined as follows:

If the glacier is land-terminating, $\int m_{\mathrm{31}}=\mathrm{0}$ by construction (a property which is used to calibrate μ^* in Eq. 5). q is then obtained by integrating the equilibrium mass balance m₃₁ along the flow line(s). q starts at zero and increases along the major flow line, reaches its maximum at the equilibrium line altitude (ELA) and decreases towards zero at the tongue (Maussion et al., 2019a).

However, this assumption does not hold for tidewater glaciers, where a steady state implies that

where q_calving is the frontal ablation flux of the glacier (m³ yr⁻¹). This flux is then converted to units of specific MB (kg m⁻² yr⁻¹) by multiplying with the ice density (900 kg m⁻³) and dividing by the total glacier area as given by the RGI. A more precise definition would be that q_calving is the average amount of ice that passes through the glacier terminus in a year for a glacier in equilibrium with the climate forcing. This has direct consequences for the calibration of the temperature sensitivity parameter μ^*. With all other things kept equal, two otherwise identical glaciers (one calving, one non-calving) will have to have different temperature sensitivities μ^*: the calving glacier will have a lower μ^*, resulting in a lowered equilibrium line altitude (ELA), a positive surface mass budget and finally to a mass flux through the terminus. The objective here is to allow the model to have a non-zero calving flux, with the goal of improving the glacier thickness inversion computed by OGGM.

3.3 Frontal ablation parameterization

3.3.2 Illustration of the method

We use the LeConte Glacier (see Fig. 2b and d) as a test case to illustrate our solution method. Figure 2d shows the result of the model's default ice thickness inversion procedure, which assumes an ice flux of zero at the terminus (q_calving=0). Note that by default the ice thickness at the glacier front h_f is zero.

First we examine how the frontal ablation flux (q_calving) from the calving law would change if we increase the terminus ice thickness of the glacier, while keeping the freeboard fixed (E_t is the only variable known in Eq. 9 “with certainty”, from the DEM surface elevation at the terminus). Figure 3a shows that the flux remains equal to zero as long as h_f is not thick enough to reach water, after which the water depth is positive and calving occurs. At this point, we are unaware of the real frontal ablation flux for this glacier, but we make some very coarse assumptions:

• Oerlemans and Nick (2005) calving law is perfectly exact,

• the tuning parameter k is known,

• our glacier is in equilibrium with climate (we assume mass-conservation inversion in OGGM),

• ice deformation at the glacier terminus follows Glen's flow law.

Under these assumptions, we set up an experiment where we compute a frontal ablation flux (from the calving law, Eq. 8) for a range of prescribed frontal ice thicknesses (see Fig. 3b, blue line), then give this flux back to the inversion model, which computes a frontal ice thickness according to the physics of ice flow (Fig. 3b, green and orange lines). As shown in Fig. 3b, both curves meet at a frontal thickness value, which complies with both the calving law (q_calving) and the ice thickness inversion model of OGGM (q). Note that changing Glen's deformation parameter A or adding sliding does not change the problem qualitatively: we will still solve a polynomial degree 5 in OGGM, with a new term in degree 3.

Figure 3Idealized experiments applied to the LeConte Glacier. (a) Frontal ablation flux computed by the calving law when prescribing a terminus thickness, with h_f ranging from 0 to 500 m. (b) Terminus ice thickness per frontal ablation flux obtained (i) by the calving law (blue curve, same as a), (ii) by OGGM using ice deformation (orange curve), and (iii) by OGGM using ice deformation and adding a sliding velocity (green curve). (c) Illustration of the ice thickness function from Eq. (10) for a given range of water depth values: (i) without a sliding velocity (orange curve) and (ii) with a sliding velocity (green curve).

Download

Figure 3c displays the same data as Fig. 3b (here as a function of the prescribed water depth), showing more clearly that there are two locations where the zero line is crossed and the condition of Eq. (10) is met. However, only one solution (the larger one) provides a realistic water depth and, therefore, a realistic frontal ablation flux.

4.1 Case study: Columbia Glacier

Columbia Glacier, located in south-central Alaska, is one of the most studied tidewater glaciers in the world. With a detailed record of its retreat since 1976, it is the single largest contributor of the Alaska glaciers to sea-level rise (Berthier et al., 2010; Larsen et al., 2015). The ice flow, ice discharge and tidewater retreat of the glacier are all extensively documented, providing rich insight into the underlying processes that modulate tidewater glacier behaviour and stability (McNabb et al., 2012). These reasons motivated the selection of Columbia Glacier as an exemplary study site to illustrate our results for an individual glacier, while the goal of our approach is the ability to improve the model representation of any calving glacier.

Following the process described in Sect. 3.3, we calculate a virtual frontal ablation for Columbia Glacier of 2.98 km³ yr⁻¹ (2.73 Gt yr⁻¹). This flux represents the estimated amount of ice passing through the terminus of the glacier if the glacier was in equilibrium with the climate for a constant glacier geometry fixed at the RGI outline's date (e.g. 2009 for this glacier). This estimate was obtained using the model's default values for the parameters k, A and f_s. McNabb et al. (2015) estimated a mean frontal ablation of 3.53±0.85 Gt yr⁻¹ during 1982–2007, with previous studies estimating 5.5 Gt yr⁻¹ for the same period (Rasmussen et al., 2011). Figure 4 shows the difference between not accounting for frontal ablation in the mass balance (q_calving=0 in Eq. 7) and accounting for frontal ablation, adding the frontal ablation flux calculated to the MB module (q_calving=2.98 km³ yr⁻¹ in Eq. 7). If q_calving=0 (Fig. 4a), we estimate the total volume of Columbia Glacier to be 270.40 km³, 29.21 % less than the volume calculated if the frontal ablation is added (Fig. 4b), which results in a volume of 349.39 km³.

When computing the ice thickness distribution map of the glacier, the impact of accounting for frontal ablation is mainly reflected in the two adjacent branches of Columbia Glacier (Fig. 4b) and at the glacier terminus (Fig. 4c). An overview of the glacier main centre line profile is shown in Fig. 4c, together with the 2007 thickness map published by McNabb et al. (2012) (dotted green line), a study that provided a reconstructed bed topography and ice thickness based on velocity observations of Columbia Glacier and mass conservation. Figure 4c also includes the result of the “consensus estimate” for Columbia Glacier ice thickness from Farinotti et al. (2019). OGGM's glacier bed estimation without accounting for frontal ablation (grey line), as well as the composite solution from Farinotti et al. (2019) (yellow line), estimate zero thickness at the calving front.

Figure 4Ice thickness inversion results for Columbia Glacier. (a) Thickness distribution before accounting for frontal ablation. (b) Thickness distribution after accounting for frontal ablation, with a frontal ablation flux computed by the model of 2.98 km³ y⁻¹. (c) Columbia Glacier main centre line profile, comparison between the 2007 estimated bed map (green doted line) from McNabb et al. (2012), the consensus estimate from Farinotti et al. (2019) (orange line), and model output before accounting for frontal ablation (grey line) and after accounting for frontal ablation (black line).

By accounting for frontal ablation in OGGM's MB and thickness inversion modules, we can compute a bedrock profile closer to the 2007 bed map, especially close to those points located at the terminus of the glacier. The frontal ablation parameterization allows OGGM to grow a thick calving front at the glacier terminus. Additionally, we observe that both bed estimations from OGGM (grey and black lines, Fig. 4c) diverge primarily below sea level.

4.2 Frontal ablation and glacier volume

In this experiment, we assign a frontal ablation flux ranging from 0 to 5 km³ yr⁻¹ to each glacier classified as potentially calving in the RGI v6.0, keeping the model's default values for the parameters A and f_s. The aim is to calculate the changes in volume for each glacier as a function of the frontal ablation value, while keeping the rest of the aspects that control the volume of the glacier fixed to the default values (e.g. solid precipitation, outline, topography, ice parameters). As a result of the automated workflow of OGGM, we are able to calculate the changes in volume of all 198 calving glaciers in Alaska² for each value in the frontal ablation flux range.

The results of this experiment are shown in Fig. 5, where the frontal ablation fluxes are expressed as a fraction of the annual accumulation ($p_{\mathrm{f}}{P}_i^{\mathrm{Solid}}\left(z\right)$) over each individual glacier. This fraction is de facto normalized to a maximum of 1, since the calving flux cannot exceed the total accumulation. Large glaciers (green and red lines in Fig. 5) will not reach this value in the prescribed calving range of 0–5 km³ yr⁻¹. Equations (5) and (7) indicate that a temperature sensitivity ${\mathit{\mu }}^{*}\le \mathrm{0}$ would imply that the glacier is producing a frontal ablation larger than its annual accumulation. When this happens, OGGM clips the temperature sensitivity μ^* to zero, setting a physical limit to the frontal ablation of each individual calving glacier.

Figure 5(a) Normalized glacier volume and (b) temperature sensitivity (μ^*) of individual glaciers, as a function of the prescribed frontal ablation fluxes normalized by the total accumulation over each glacier. The different colours represent different glacier classes.

Download

Figure 5a shows that the effect of frontal ablation on the glacier volume is systematic, in that accounting for frontal ablation in the MB will always result in an increase in the glacier volume. Even if the frontal ablation fraction is only 0.14 of the total accumulation, a glacier volume can be underestimated by up to 20 % if we ignore this extra source of ablation. However, there is a wide range of sensitivities of the estimated glacier volume to the calving flux, and no simple relation to, e.g. glacier size, was found. Other glacier-specific parameters likely to play a role are the slope, the accumulation area ratio and the total precipitation.

4.3 Effect of frontal ablation on ice velocity

To analyse the effect of frontal ablation on ice velocity, we keep the same model configuration (default values of k, A and f_s) and calculate the average ice velocity along the main flow line for all marine-terminating glaciers that produced a frontal ablation flux. Figure 6 shows the difference between the average velocity output of the model when accounting for frontal ablation and without accounting for frontal ablation. When taking frontal ablation into account, the glaciers experience an increase in ice velocity towards the terminus. This increase in velocity is due to an increase in the mass flux (and, therefore, ice thickness) when we account for frontal ablation.

Figure 6Glacier average velocity differences between the two outputs of the model for a subset of marine-terminating glaciers. The differences are between the model output before accounting for frontal ablation and after accounting for frontal ablation in points along the main flow line. The x axis has been normalized.

Download

These results highlight the importance of applying a frontal ablation parameterization at the initialization stages of the model in order to recreate a realistic tidewater glacier behaviour. Without this extra term on the mass balance, velocities and ice thicknesses go to zero towards the terminus. This is not only a problem for the inversion procedure: these unrealistic features will also affect the dynamical runs realized with the forward model, i.e. any calving parameterization applicable in the future will rely on a realistic bedrock to work properly.

Note that these velocities are not surface velocities but average section velocities. Annual surface velocities would have to be estimated from these values and the vertical profile of velocity in order to be compared to observations. Furthermore, since we run OGGM under an equilibrium assumption, the results presented here will not reflect the transient states that appear in observations. The usefulness of any comparison to other data (observed or modelled) is, therefore, limited. Additionally, some of the velocity maps in Alaska previously published (e.g. Burgess et al., 2013) are computed with many glaciers undergoing significant interannual velocity variability over the observation interval and only one velocity snapshot is included in the maps. Velocities might thus not represent long-term average. However, comparing surface velocities derived from OGGM with observations might be useful when no previous estimates of frontal ablation fluxes exist in a RGI region (e.g. Greenland), providing another way of calibrating OGGM parameters.

4.4 Sensitivity studies in Alaska marine-terminating glaciers

We perform different sensitivity experiments to study the influence of (i) the calving constant of proportionality k, (ii) Glen's temperature-dependent creep parameter A and (iii) sliding velocity parameter, on the regional frontal ablation of Alaska. The results of these experiments are shown in Fig. 7. In the first experiment we vary the calving constant of proportionality k in a range of 0.24–2.52 yr⁻¹ and used the model default values for Glen A and sliding parameter. Figure 7a shows that our estimate for the regional frontal ablation matches the regional estimate by McNabb et al. (2015) if k has an approximated value of 0.63 yr⁻¹, in the case of excluding sliding (f_s=0), or if k is equal to 0.67 yr⁻¹ in the case of including a sliding velocity (with $f_{\mathrm{s}}=\mathrm{5.7}\times {\mathrm{10}}^{-\mathrm{20}}$ s⁻¹ Pa⁻³). It is important to emphasize that the regional frontal ablation from McNabb et al. (2015) is only comprised of 27 glaciers but represents an estimated 96 % of the total frontal ablation of Alaska.

Figure 7Total frontal ablation of Alaska marine-terminating glaciers computed with varying OGGM parameters. The dashed dark line indicates the Alaska regional frontal ablation calculated by McNabb et al. (2015); light grey shading indicates the standard errors as provided in the study. (a) Sensitivity on calving constant of proportionality (k). (b) Sensitivity on Glen's A parameter, the dashed coloured lines represent zero sliding. (c) Sensitivity on sliding parameter (f_s). Crosses in all plots represent the intercepts between OGGM frontal ablation estimates and McNabb et al. (2015). Note the different y axis ranges.

Download

We then keep these two values of k and vary the values of Glen A creep parameter and sliding parameter (f_s). The results are shown in Fig. 7b and c. It is well known that ice flow models are sensitive to the values chosen for parameters describing ice rheology and basal friction (e.g. Enderlin et al., 2013; Brondex et al., 2017). As expected, our frontal ablation estimates are also sensitive to different values of Glen A and sliding parameter but highly dependent on different values of k, at least for the first part of the k values range (0.24–0.80 yr⁻¹). The linear relationship between q_calving and k at the start of the curve in Fig. 7a is mainly a consequence of the calving law used in the parameterization. For larger k values (≥0.8 yr⁻¹) the shape of the curve is due to OGGM's physical constraint of clipping μ^* to zero and calculating the maximum q_calving allowed by the local climate (see Sect. 3.3.3).

Maussion et al. (2019a) showed that both sliding and ice rheology (A) have a strong influence on OGGM's computed ice volume, hence a strong influence on the thickness of the glacier and, in this case, the frontal ablation estimate. Like in Maussion et al. (2019a); Fig. 7 shows that one could always find an optimum combination of Glen A and sliding parameters that lead to (in this case) previously calculated frontal ablation estimates. Enderlin et al. (2013) also showed that when such flow line models are applied to a tidewater glacier, there is a non-unique combination of these parameter values that can produce similar stable glacier configurations, making k, Glen A and f_s parameters highly dependent on observations of either frontal ablation, ice velocity or glacier ice thickness.

4.5 Regional volume of marine-terminating glaciers for different model configurations

Finally, we compute the total volume of marine-terminating glaciers for different “equally good” parameter sets based on the results of Fig. 7a, b and c. Each configuration is constructed by finding the intercepts between the model frontal ablation estimates and the regional estimate from McNabb et al. (2015), including the intercepts to the lower and upper error (see Fig. 7). A summary of the different parameter sets used for each model run can be found in Table 1 and the results of each configuration are shown in Fig. 8. Each configuration was run twice: once setting q_calving=0, then a second time accounting for frontal ablation.

Figure 8Total volume of Alaskan marine-terminating glaciers before (blue) and after (green) accounting for frontal ablation and the total volume below sea level (red) before and after accounting for frontal ablation. The lightly shaded colour bars represent configurations obtained from finding the intercepts between OGGM frontal ablation estimates and the lower and upper error provided by McNabb et al. (2015). The grey bar represents the consensus estimate for these glaciers obtained by Farinotti et al. (2019). The descriptions for each configuration can be found in Table 1.

Download

Similarly to the results shown in Sect. 4.1 and 4.2, Fig. 8 shows that there are significant differences between total volume estimates without and with accounting for frontal ablation. Volume estimates after accounting for frontal ablation are 11.7 % to 19.7 % higher than the volume estimates ignoring frontal ablation, considering all model configurations shown in Table 1, indicating a robust relationship. We find that there are no significant differences between the resulting volumes for different k values and that the differences in volume estimates between configurations are mainly due to adding or ignoring a sliding velocity or varying the value of the Glen A creep parameter.

Additionally, we also calculate the regional ice volume below sea level. The results for Columbia Glacier discussed in Sect. 4.1 might create the impression that the differences in thickness along the main centre line, with and without accounting for frontal ablation, are not relevant for the potential glacier contribution to sea-level rise, since most of the differences in thickness (grey and black line in Fig. 4c) are found below sea level. However, Fig. 8 shows that considering the whole region, a significant fraction of the total volume difference is found above sea level, implying that accounting for frontal ablation will directly impact the estimate of these glacier's potential contribution to sea-level rise. By introducing frontal ablation, the volume estimate of marine-terminating glaciers in Alaska is increased from 9.18±0.62 to an average of 10.61±0.75 mm s.l.e, of which 1.52±0.31 mm s.l.e are found to be below sea level (instead of only 0.59±0.08 without). The uncertainties presented here are the standard deviation of the model configurations shown in Table 1. The consensus estimate from Farinotti et al. (2019) for these glaciers is 7.68 mm s.l.e, 27.58 % lower than our average estimate of 10.61±0.75 mm s.l.e.