Impact of frontal ablation on the ice thickness estimation of marine-terminating glaciers in Alaska

Frontal ablation is a major component of the mass budget of calving glaciers, strongly affecting their dynamics. Most global-scale ice volume estimates to date still suffer from considerable uncertainties related to (i) the implemented frontal ablation parameterization or (ii) not accounting for frontal ablation at all in the glacier model. To improve estimates of the ice thickness distribution of glaciers, it is thus important to identify and test low-cost and robust parameterizations of this process. By implementing such parameterization into the ice thickness estimation module of the Open Global Glacier Model (OGGM v1.1.2), we conduct a first assessment of the impact of accounting for frontal ablation on the estimate of ice stored in glaciers in Alaska. We find that inversion methods based on mass conservation systematically underestimate the mass turnover and, therefore, the thickness of tidewater glaciers when neglecting frontal ablation. This underestimation can amount to up to 19 % on a regional scale and up to 30 % for individual glaciers. The effect is independent of the size of the glacier. Additionally, we perform different sensitivity experiments to study the influence of (i) a constant of proportionality (k) used in the frontal ablation parameterization, (ii) Glen’s temperature-dependent creep parameter (A) and (iii) a sliding velocity parameter (fs) on the regional dynamics of Alaska tidewater glaciers. OGGM is able to reproduce previous regional frontal ablation estimates, applying a number of combinations of values for k, Glen’s A and fs. Our sensitivity studies also show that differences in thickness between accounting for and not accounting for frontal ablation occur mainly at the lower parts of the glacier, both above and below sea level. This indicates that not accounting for frontal ablation will have an impact on the estimate of the glaciers’ potential contribution to sea-level rise. Introducing frontal ablation increases the volume estimate of Alaska marine-terminating glaciers from 9.18±0.62 to 10.61±0.75 mm s.l.e, of which 1.52±0.31 mm s.l.e (0.59±0.08 mm s.l.e when ignoring frontal ablation) are found to be below sea level.

Abstract. Frontal ablation is a major component of the mass budget of calving glaciers, strongly affecting their dynamics. Most global-scale ice volume estimates to date still suffer from considerable uncertainties related to (i) the implemented frontal ablation parameterization or (ii) not accounting for frontal ablation at all in the glacier model. To improve estimates of the ice thickness distribution of glaciers, it is thus important to identify and test low-cost and robust parameterizations of this process. By implementing such parameterization into the ice thickness estimation module of the Open Global Glacier Model (OGGM v1.1.2), we conduct a first assessment of the impact of accounting for frontal ablation on the estimate of ice stored in glaciers in Alaska. We find that inversion methods based on mass conservation systematically underestimate the mass turnover and, therefore, the thickness of tidewater glaciers when neglecting frontal ablation. This underestimation can amount to up to 19 % on a regional scale and up to 30 % for individual glaciers. The effect is independent of the size of the glacier. Additionally, we perform different sensitivity experiments to study the influence of (i) a constant of proportionality (k) used in the frontal ablation parameterization, (ii) Glen's temperature-dependent creep parameter (A) and (iii) a sliding velocity parameter (f s ) on the regional dynamics of Alaska tidewater glaciers. OGGM is able to reproduce previous regional frontal ablation estimates, applying a number of combinations of values for k, Glen's A and f s . Our sensitivity studies also show that differences in thickness between accounting for and not accounting for frontal ablation occur mainly at the lower parts of the glacier, both above and below sea level. This indicates that not accounting for frontal ablation will have an impact on the estimate of the glaciers' potential contribution to sea-level rise. Introducing frontal ablation increases the volume estimate of Alaska marine-terminating glaciers from 9.18±0.62 to 10.61±0.75 mm s.l.e, of which 1.52±0.31 mm s.l.e (0.59±0.08 mm s.l.e when ignoring frontal ablation) are found to be below sea level.

Introduction
Estimates of the spatial distribution of ice thickness are needed as initial conditions for glacier models, for attempting to understand how glaciers respond to climate change and for quantifying their contribution to sea-level rise. Despite this importance, ice thickness measurements around the globe are scarce, performed only on approx. 600 glaciers (Gärtner-Roer et al., 2014) out of more than 200 000 identified in the latest Randolph Glacier Inventory (RGI v6.0, Pfeffer et al., 2014). In order to overcome this under-sampling problem, a number of methods have been developed to infer the total volume and/or the ice thickness distribution of glaciers from characteristics of the glacier surface properties. Some of these methods rely on scaling approaches relating the length, slope and area of the glacier to its total volume (e.g. Bahr et al., 1997;Lüthi, 2009;Radić and Hock, 2011;Grinsted, 2013). Others rely on parameterizations of basal shear stress (e.g. Paul and Linsbauer, 2012;Linsbauer et al., 2012;Frey et al., 2014), observed surface velocities (e.g. Gantayat et al., 2014), or applying the shallow-ice approximation (e.g. Oerlemans, 1997;Cuffey and Paterson, 2010) and/or an integrated form of Glen's flow law (see Farinotti et al., 2017, for a review of all these methods and Farinotti et al., 2019, for a global-scale intercomparison).
One method presented by Farinotti et al. (2009) and successfully applied several times since then (e.g. Morlighem et al., 2011;Huss and Farinotti, 2012;Clarke et al., 2013;Maussion et al., 2019a) combines ice flow dynamics and mass conservation principles to constrain mass fluxes through given glacier cross sections. The method infers ice thickness from estimates of ice fluxes derived from the assumption that ice fluxes balance the surface mass budget (Farinotti et al., 2009). The results are thus sensitive to the spatial distribution of the mass flux and the mass balance. For calving glaciers, the surface mass budget cannot be considered balanced, even assuming equilibrium between glacier and climate. The derived ice thickness estimate for these glaciers hence depends on estimates of frontal ablation.
Frontal ablation (mass loss by calving and frontal melting Pope, 2012) is an efficient process to deliver ice from glaciers and ice sheets into the ocean. It has contributed substantially to sea-level rise in the past and played an important role in the stability of ice sheets and tidewater glaciers during the Pleistocene (Benn et al., 2007). Calving is strongly coupled with dynamical processes inside the glacier. An increase in the ice flux can trigger a calving event and in turn this event can accelerate the movement of the ice. External aspects like ocean temperature, fjord bathymetry and, in polar areas, sea ice concentration along the calving front can also influence the discharge of solid ice to the ocean (Straneo et al., 2013). As a consequence of the diverse nature of calving processes, the development of parameterizations of frontal ablation in numerical ice sheet and glacier models remains an important challenge. There is a wide spectrum of approaches that vary in scale and complexity, justified through the diversity of intended applications of the models (Price et al., 2015).
There have been many successful efforts to represent frontal ablation for individual glaciers (e.g. Ultee and Bassis, 2016;Åström et al., 2014;Todd and Christoffersen, 2014;Oerlemans et al., 2011;Nick et al., 2010). While these achieve encouraging results, it is unlikely that they can be implemented in a global glacier model anytime soon because of the amount and quality of data needed to constrain this type of model. The crevasse-depth criterion proposed by Nick et al. (2010) for example, requires knowledge of surface melt and refreeze rates at the crevasses of the glacier tongue and crevasse depth observations to calibrate and validate these rates. These kinds of observations are hard to obtain for entire glaciated regions: e.g. the 198 calving glaciers in Alaska investigated here or the 3222 glaciers classified as calving (marine-and lake-terminating) glaciers in the RGI v6.0. Other recent calving models that use discrete particles or a full-Stokes model approach (e.g. Åström et al., 2014;Todd and Christoffersen, 2014;Todd et al., 2018) are too computationally expensive to be included in global glacier models that seek to consistently simulate past and future global-scale glacier changes. At the regional and global scale, very few estimates of frontal ablation fluxes of glaciers outside the ice sheets exist (Blaszczyk et al., 2009;Burgess et al., 2013;McNabb et al., 2015;Huss and Hock, 2015). From all the global glacier models published in the last decade, only Huss and Hock (2015) account for frontal ablation of marine-terminating glaciers. However, this model, along with the rest of the ice thickness inversion methods, still suffers from considerable uncertainty associated with the uncertainty of the frontal ablation parameterization.
For improving ice thickness distribution estimates at the global scale, it is thus important to identify and test low-cost and robust parameterizations of frontal ablation that might not resolve all the dynamical processes at the calving front (e.g. subaqueous frontal melting, subaerial frontal melting and sublimation) but that can estimate the amount of ice passing through the terminus of the glacier during a mass balance year. Using the ice thickness estimation module of the Open Global Glacier Model (OGGM v1.1.2), we assess the impact of frontal ablation on the estimate of ice stored in Alaska glaciers classified as marine-terminating in the RGI v6.0 (also referred to tidewater glaciers in this study).
Alaska glaciers cover approximately 12 % of the global glacier area outside of the ice sheets . In the RGI (v6) there are 27 109 glaciers in the region occupying an area of 86 776.6 km 2 , including adjacent glaciers in the Yukon and in British Columbia. From these glaciers, 51 have been classified as marine-terminating (74 km of tidewater margin) and 147 as lake-and river-terminating glaciers (420 km of lake/river margin) occupying an area of 11 962.4 and 16 720.6 km 2 , respectively. Calving glaciers (marine-and lake-terminating) occupy approximately 33 % of the Alaska glacier area ( Fig. 1; Pfeffer et al., 2014;Kienholz et al., 2015).
The glaciers are divided into six subregions in the RGI. Subregions 1 and 3 contain only land-terminating glaciers. Calving glaciers are mostly concentrated in the subregions 4, 5 and 6, along the mountain ranges of the southern Alaska coast (Fig. 1), an area characterized by maritime climate and topography reaching > 5000 m a.s.l. . Glaciers contained in the RGI in this region range in size from a few square kilometres (Ogive Glacier, 2.8 km 2 ) to many thousands of square kilometres (Hubbard Glacier, 3400 km 2 ; McNabb and Hock, 2014).
The subregions 4 and 5 are well-studied glacierized areas of Alaska. McNabb et al. (2015McNabb et al. ( ) presented a 28-year record (1985McNabb et al. ( -2013 of frontal ablation for a subset of marineterminating glaciers that includes the 27 most dominant tidewater glaciers of the region. They represent 96 % of the total tidewater glacier area in the gulf of Alaska. The total mean rate of frontal ablation was estimated to be 15.11 ± 3.63 Gt yr −1 (16.48 ± 3.96 km 3 yr −1 ), over the period 1985-2013. Other studies also reported similar values (e.g. Larsen et al., 2007). Frontal ablation in this region is heavily dominated by two glaciers in particular: Hubbard and Columbia  (McNabb et al., 2015). (b) Regional glacier types and basic statistics of the database (area of glaciers per terminus type, regional contribution to the Alaska area in percent and percentage of the regional area, which cannot be modelled by OGGM).
glaciers (McNabb et al., 2015). Additionally, McNabb et al. (2015) identified 36 actively calving tidewater glaciers in Alaska; 27 of those were used to estimate the total mean rate of frontal ablation presented in McNabb et al. (2015).
We implement a simple parametrization of frontal ablation in OGGM, following the approaches proposed by Oerlemans and Nick (2005) and Huss and Hock (2015). By performing sensitivity studies on the model, we (i) investigate the effect of accounting for frontal ablation on the ice thickness estimation of OGGM and on the ice volume estimate for these glaciers and (ii) study the impact of varying several OGGM parameters (the calving constant of proportionality k, Glen's temperature-dependent creep parameter A, and sliding velocity parameter f s ) on the regional frontal ablation rates of Alaska.
2 Input data and pre-processing

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 Pa-cific, 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. 1 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.

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.

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.

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 elevationdependent 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 −1 ). 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).

Open Global Glacier Model (OGGM) and frontal ablation parameterization
For this study, a simple frontal ablation parametrization is implemented into the Open Global Glacier Model (OGGM v1.1.2). OGGM is developed to provide a global-scale, modular and open-source numerical model framework for consistently simulating past and future global-scale glacier change.
The mathematical framework of the model and its capabilities have been explained in detail by Maussion et al. (2019a).
In this section, we will only describe the modifications done to the mass balance and ice thickness inversion modules, together with the frontal ablation parametrization implemented in order to improve the initialization of the model for marineterminating glaciers. Section 3.3 provides details on the limitation of applying the parameterization to lake-terminating glaciers.

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 3 s −1 ) through a glacier cross section of area S (m 2 ) is defined as follows: with u being the average cross section velocity (m s −1 ). By applying the well-known shallow-ice approximation (Hutter, 1981(Hutter, , 1983Cuffey 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 2.4×10 −24 s −1 Pa −3 ), n the exponent of Glen's flow law (default: n = 3), h 0 the centre-line ice thickness (m) and τ the basal shear stress, defined as follows: with ρ the ice density (900 kg m −3 ), g the gravitational acceleration (9.81 m s −2 ) 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: 5.7 × 10 −20 m −2 s −1 Pa −3 ). We then assume that the centreline velocity is equal to the average section velocity (u ≈ 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 = 2 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 R + , easily computable for each grid point.

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 −2 s −1 ) at elevation z is computed as follows: where P solid i 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 www.the-cryosphere.net/13/2657/2019/ The Cryosphere, 13, 2657-2672, 2019 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 31 (t) is equal to zero, with m 31 (t) being the glacier integrated mass balance computed for a 31-year period centred around the year t (e.g. t * = 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 31 (t)) is then assumed to be equal to the "apparent mass balance" (m =ṁ − ρ ∂h ∂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, m 31 = 0 by construction (a property which is used to calibrate µ * in Eq. 5). q is then obtained by integrating the equilibrium mass balance m 31 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 3 yr −1 ). This flux is then converted to units of specific MB (kg m −2 yr −1 ) by multiplying with the ice density (900 kg m −3 ) 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.

Calving law
To account for frontal ablation of marine-terminating glaciers, we employ a calving law proposed by Oerlemans and Nick (2005) that has already been applied at a large scale by Huss and Hock (2015). The annual frontal ablation flux q calving (km 3 yr −1 ) is computed as a function of the height (h f ), width (w) and estimated water depth (d) of the calving front as follows: k is a calibration parameter (which has a default value of 2.4 yr −1 in this study). The water depth (d) is estimated from freeboard, using elevation and ice thickness (h f ) data obtained from the model output: where E t is the elevation of the glacier surface at the terminus and z w is the elevation of the water body with respect to sea level. The water depth (d) is estimated using the terminus elevation (E t ) obtained by projecting the RGI outline onto the DEM (i.e. the terminus elevation is the top of the cliff). We follow the same definition as Oerlemans and Nick (2005), where d is the bed elevation with respect to sea level. For lake-terminating glaciers, we are not able to estimate a water depth since one would need to know the freeboard of the glacier terminus, i.e. the elevation of the glacier lake surface. For this reason, most of our experiments and results focus on marine-terminating glaciers only (z w is set of 0 m a.s.l.), with the exception of the experiment presented in Sect. 4.2. Unlike Huss and Hock (2015), who estimated the thickness of the calving front (h f ) by scaling approaches, we solve for the ice thickness by prescribing that the amount of ice calved (q calving ) must be equal to the amount of ice delivered to the terminus by OGGM (q, computed from ice deformation and sliding in Sect. 3.1): q calving varies with h f as a polynomial of degree 2. q is a polynomial in h f of degree 5 (with n = 3 in Eq. 2), with an extra term in degree 3 if we account for a sliding velocity (see Eq. 4). Equation (10) is, therefore, a polynomial that can be solved for h f .

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 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.

Implementation
We solve the polynomial in Eq. (10) numerically, via boundconstrained minimization methods (algorithm provided by SciPy, Jones et al., 2001), which leads to a quick convergence. The advantage over an analytical solution is that numerical solvers have the flexibility to be applied to any other formulation of q calving and q, i.e. that this method will still be applicable if a lateral drag parameterization or another formulation for the calving law is added to OGGM in the future.
After finding the solution for the frontal ice thickness (h f ) and the corresponding frontal ablation flux (q calving ), we give this flux back to the mass balance model (Eq. 7), adjust the temperature sensitivity of the glacier µ * and invert for a new ice thickness distribution for the entire glacier (see Sect. 4.1 for results). This always results in a adjustment of µ * towards lower values, in order to lower the ELA, and unbalancing of the steady-state surface mass budget for a frontal ablation flux to exist. Note that this adjustment of mass balance is always necessary (regardless of the choice of model parameters such as k or Glen's A) in order to ensure massconservation and the update of upstream ice thickness.
However, sometimes the flux estimated by the calving law (Eq. 8) is too large to be sustained by the surface mass balance. Even without glacier melt (µ * = 0), the total accumulation over the glacier is too small to close the frontal mass budget. This can be due to several factors: e.g., frontal ablation is overestimated or solid precipitation is underestimated.
The frontal ablation can be overestimated if k and/or the calving law does not represent the dynamics of that particular glacier or if h f is overestimated. In most cases (see Results), it is possible to find a realistic µ * compatible with a frontal ablation flux, but when this is not possible µ * is fixed to zero and the frontal ablation flux q calving is obtained by closing the mass budget instead of using the calving law.

Results
We apply this frontal ablation parameterization to all marineterminating glaciers in Alaska. We study the impact of including this parameterization on the estimated glacier thickness, volume and ice flow velocity. The following sections describe different sensitivity experiments: (i) varying the frontal ablation flux added to the MB model and assessing the impact on glacier volume, (ii) varying several model parameters (Glen's flow law ice creep parameter, A; a sliding parameter, f s ; and the calving constant of proportionality, k) and assessing each parameter's impact to the regional frontal ablation of Alaska, and (iii) showing the impact of different model configurations (obtained from the sensitivity experiments of Sect. 4.4) on the total volume of Alaska marineterminating glaciers. The parameter set-up for each configuration can be found in Table 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 (Mc-Nabb 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 3 yr −1 (2.73 Gt yr −1 ). 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 −1 during 1982-2007, with previous studies estimating 5.5 Gt yr −1 for the same period . 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 3 yr −1 in Eq. 7). If q calving = 0 (Fig. 4a), we estimate the total volume of Columbia Glacier to be 270.40 km 3 , 29.21 % less than the volume calculated if the frontal ablation is added (Fig. 4b), which results in a volume of 349.39 km 3 .
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 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.

Frontal ablation and glacier volume
In this experiment, we assign a frontal ablation flux ranging from 0 to 5 km 3 yr −1 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 2 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 f P Solid i (z)) 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 3 yr −1 . Equations (5) and (7) indicate that a temperature    sensitivity µ * ≤ 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 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.

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. Figure 6. Glacier 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.
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.
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.

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 −1 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 −1 , in the case of excluding sliding (f s = 0), or if k is equal to 0.67 yr −1 in the case of including a sliding velocity (with f s = 5.7 × 10 −20 s −1 Pa −3 ). 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. 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 −1 ). 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 −1 ) 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.

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 con- structed 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.
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 Figure 8. Total 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. 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.

Discussion
We have shown that the model is capable of computing regional frontal ablation estimates by tuning model parameters with published regional-scale estimates of frontal ablation, but the question of model performance for individual glaciers still remains open. In areas with no observational data or previous knowledge of frontal ablation, OGGM could make use of physical constraints (e.g. that µ * must be greater than zero), as well as bathymetry and terminus width estimates to calibrate the model at the glacier scale. In the following section, we will explain such calibrations, discuss other parameters that affect frontal ablation estimates and discuss these estimates for individual glaciers.
In all previous model runs, we used the standard OGGM terminus geometry computation without correcting the width and water depth at the glacier front using potentially known values from other sources. As a result, not all of the glaciers classified as marine-terminating glaciers in the RGI produce a frontal ablation flux in OGGM (six glaciers). This is mainly due to a wrong estimation of the water depth from freeboard. These glaciers typically have a high terminus elevation (e.g. E t = 151.96 m a.s.l. for the Chenega Glacier, RGI60-01.09639), for which the only possible value of q calving that complies with both the calving law and the ice thickness inversion model of OGGM is a q calving = 0, since there is not enough mass turnover to grow a calving front under our mass conservation assumptions (see Sect. 3.3.2). The wrong water depth estimation can thus be best explained by a poor surface altitude estimation at the calving front (E t ). The problematic surface altitude estimation in turn can probably be explained by a mismatch between the acquisition dates of the DEM and the glacier outline. Maussion et al. (2019a) noted that a number of glaciers will suffer from poor topographic information, especially those located in the high latitudes. Most marine-terminating glaciers are located in regions where cloud-free satellite measurements are rare. Therefore, the DEM of these regions might present errors (e.g. a wrong elevation at the terminus and/or a date of data acquisition, which does not match that of the RGI outlines) that will spread to water depth estimations from freeboard (see Eq. 9). The possibility of using higher-resolution DEMs such as the ArcticDEM was explored during this study but was quickly eliminated because of large data voids present on the data, especially for big glaciers (e.g. Hubbard Glacier). However, new datasets such as the TanDEM-X (Wessel et al., 2018) are currently being explored for future versions of OGGM.
For 36 marine-terminating glaciers, we assess the model performance in comparison to the estimates by McNabb et al. (2015), with and without corrections for these errors (Fig. 9a). Calving front widths were corrected with the Alaska Tidewater Glacier Terminus Positions database (Mc-Nabb and Hock, 2014). The database contains terminus positions for 49 marine-terminating glaciers. Since three of these glaciers (Grand Pacific, Hubbard and Sawyer glaciers) were merged with their respective pair branches (Ferris, Valerie and West Sawyer glaciers), we are left with a total of 46 glacier terminus widths. The widths are computed by selecting the terminus positions closest to the glacier's RGI outline and by averaging the widths that resulted from the projection of the vector lines selected. These widths are used to correct OGGM's flow line width at the calving front in the cases where the model is not able to represent the real calving front width. The last flow line width at the front of the glacier is then clipped to the width value estimated from the database. To smooth the transition between the clipped value and OGGM's flow line width, we linearly interpolate between the clipped value and 5 pixels upstream on the flow line. We then correct the modified widths to preserve the same glacier area than the RGI's. By doing this, we slightly modify the altitude-area distribution of the glacier.
Additionally, multi-beam bathymetry data from NOAA National Centers for Environmental Information (2004) was used to estimate the water depth in front of the glacier terminus. These data were used only for glaciers where the DEM resolution would not allow an estimate of the water depth from elevation data and ice thickness (freeboard). The bathymetry data were compiled into a raster format and provided to us by Robert McNabb (personal communication, 2018). Both corrections were used for Fig. 9 only. Figure 9a demonstrates that without calibrating any OGGM parameter (only using the model default values for Glen A, f s and k) but making use of additional data (e.g. terminus positions and bathymetry), we are able to estimate a frontal ablation flux for individual glaciers within the same order of magnitude as those estimated by McNabb et al. (2015). The model root-mean-square error (RMSE) is reduced from 1.08 km 3 yr −1 (mean deviation of 0.28) to 0.53 km 3 yr −1 (mean deviation of 0.11). Even though part of these errors may arise from the fact that glaciers are in a disequilibrium state at the time of the McNabb et al. (2015) estimate, errors in boundary conditions (e.g. topography date that does not coincide with the glacier outline date and uncertainties in the frontal width) and plain model errors also contribute. By using bathymetry and real terminus width estimates we improve the boundary conditions of the parameterization that are highly dependent on the DEM quality.
When these corrections (terminus width and water depth) are not implemented and errors occur while estimating the real terminus geometry, OGGM has to rely on clipping µ * to be larger than or equal to zero, setting a physical limit where the frontal ablation flux for each individual tidewater glacier cannot be larger than its annual accumulation (p f P Solid i (z)). This is not ideal because it implies that all of the glacier's ablation in an equilibrium setting is due to frontal ablation and no surface melt occurs, which is unrealistic in the climate conditions of Alaska. For applications on the global scale, bathymetry data and terminus mapping will be very valuable in regions with poor topographic resolution and where no observations of frontal ablation exist.
Additionally, we compare our final glacier volume estimates with the default configuration and by correcting the terminus geometry. These results are shown in Fig. 9b, together with the consensus volume estimate for each glacier from Farinotti et al. (2019). Figure 9b shows that even if these corrections to the glacier terminus might have a large effect on the frontal ablation flux for some glaciers (e.g. Hubbard glacier), the effect is not as big as if we do not account for frontal ablation at all. This is the case of the consensus estimate from Farinotti et al. (2019) (red bars in Fig. 9b), where the models used do not account (or crudely account) for this extra loss of mass when inverting for the ice thickness.

Conclusions
We have implemented a frontal ablation parameterization into OGGM and shown that inversion methods ignoring frontal ablation systematically underestimate the mass flux and thereby the thickness of calving glaciers. Accounting for frontal ablation in ice thickness inversion methods based www.the-cryosphere.net/13/2657/2019/ The Cryosphere, 13, 2657-2672, 2019 on mass conservation (as listed in Farinotti et al., 2017) increases estimates of the regional ice mass stored in marineterminating glaciers by approximately 11 % to 19 %. While for individual glaciers ice volume may be underestimated by up to 30 % when ignoring the impact of frontal ablation, the effect is independent of the size of the glacier. Implementing a frontal ablation parameterization allows OGGM to represent a non-zero thickness calving front, which is necessarily the case when no ice flux is assumed to cross the glacier terminus. This parameterization is key for initializing the glacier's thickness in the model. The model was able to reproduce previously calculated regional frontal ablation estimates by finding the best combination of values for k, Glen's A and the sliding parameters. Note that this comparison is limited by the equilibrium condition imposed on OGGM during initialization, which is not the case in observations. The best-performing parameter set for transient runs of OGGM may be different.
Our sensitivity studies also show that the differences in thickness, between adding or not adding frontal ablation to the MB model, occur mainly at the lower parts of the glacier but often above sea level. This indicates that not accounting for frontal ablation will have an impact on the estimate of this glacier's potential contribution to sea-level rise.
Additionally, our experiments highlight the need for bathymetry data and terminus mapping, as they may constrain model parameters when the DEM quality is not sufficient to provide a realistic estimate of the terminus geometry.
Code availability. The OGGM software, together with the frontal ablation parameterization module, is coded in the Python language and licensed under the GPLV3 free software license. The latest version of the OGGM code is available on Github (https: //github.com/OGGM/oggm, last access: 8 October 2019, Maussion et al., 2019), the documentation is hosted on ReadTheDocs (https://oggm.readthedocs.io, last access: 8 October 2019), and the project web page for communication and dissemination can be found at https://oggm.org (last access: 8 October 2019). The code and data used to generate all figures and analyses of this paper can be found at https://github.com/bearecinos/cryo_calving_2019 (last access: 8 October 2019) (Recinos, 2019). The OGGM version used for this study is available in a permanent DOI repository (https://doi.org/10.5281/zenodo.2580277, Maussion et al., 2019).
Author contributions. BM and FM are the initiators of the OGGM project and conceived this study. BR is the main developer of the frontal ablation parameterization module and wrote most of the paper. FM is the main OGGM developer and was largely involved with the development of the frontal ablation parameterization module. TR made significant contributions to the OGGM code and the frontal ablation parameterization module and is responsible for the successful deployment of the code in supercomputing environments.
Competing interests. The authors declare that they have no conflict of interest.