Exploring the impact of a frontal ablation parameterization on projected 21st-century mass change for Northern Hemisphere glaciers

Abstract

Marine-terminating glaciers cover more than one-fourth of the total glacierized area in the Northern Hemisphere outside the Greenland ice sheet. It is therefore crucial to ensure an adequate representation of these glaciers when projecting large-scale glacier mass changes. We investigate how the introduction of marine frontal processes in the modeling chain influences the results of mass change projections, compared to projections neglecting such processes. We find that including frontal processes reduces the projected glacier mass loss, since incorporating frontal ablation in the model's mass-balance calibration results in a decrease in marine-terminating glaciers’ sensitivity to atmospheric temperatures. We also find that retrograde bed slopes lead to increased frontal ablation as the atmosphere warms, while frontal ablation decreases if bed slopes are prograde. These opposing effects have the potential to partly cancel each other when considering large glacier ensembles. Although we do not account for potential future changes in oceanic climate yet, any effect of these would be moderated by around half of today's marine-terminating glaciers becoming land-terminating in the course of the 21st century. While we find a significant influence of ice flow parameters on our results, boundary conditions remain the largest source of uncertainty in our projections.

1. Introduction

Earth's glaciers are rapidly losing mass, causing global mean sea level rise (GMSLR) and changing the hydrology of regions where they are part of the landscape (Huss and Hock, 2018; Zemp and others, 2019; Hugonnet and others, 2021). In regions where glaciers are located close to the ocean, their direct export of fresh water can substantially change the ocean's surface water properties and thus its local and regional circulation as well as ecology (Castro de la Guardia and others, 2015; Meire and others, 2017). Some glaciers are in direct contact with the ocean, often producing icebergs. These marine-terminating glaciers are special insofar as they do not only undergo mass changes due to atmospheric climate perturbations, but also by ice–ocean interaction processes taking place at their front, such as iceberg calving and submarine melt (Straneo and others, 2013). The contact of marine-terminating glaciers’ terminal boundary with ocean water also makes them subject to different ice flow dynamics compared to land-terminating ones. Overall, marine-terminating glaciers comprise more than one-fourth of the glaciated area in the Northern Hemisphere outside the Greenland ice sheet (RGI Consortium, 2017). Hence, they contain such a large amount of ice that their mass loss has the potential to intensify freshwater input to the oceans, thereby increasing GMSLR as well as triggering changes in circulation and ecological patterns. In order to understand and project regional and global glacier mass changes as well as their wider implications, it is therefore necessary to investigate processes occurring at marine-terminating fronts and to incorporate them in numerical models.

One example illustrating the need to include marine frontal processes comprehensively in a numerical glacier model is that neglecting frontal ablation in the surface mass balance model's calibration will result in overestimating marine-terminating glaciers’ sensitivity to atmospheric temperatures. That is because geodetic observations, for instance, do not usually distinguish between surface and frontal ablation. Thus, if such data are used in the surface mass balance model's calibration and frontal ablation is not taken into account, it is implicitly assumed that all of the ablation takes place on the surface. Although the inclusion of frontal ablation in the dynamical model compensates part of the decrease in mass removal due to the lowering of sensitivities to atmospheric temperatures, the question arises whether there is a net effect of actually partitioning these two types of ablation.

Accounting for frontal processes does not only affect model projections, but model estimates of initial ice thicknesses as well (Recinos and others, 2019, 2021). Since ice volume cannot be measured at large scales, models constrained by observations are often used (Farinotti and others, 2009). Among these approaches, those relying on mass-conservation typically do not account for frontal ablation in the mass budget (i.e. not allowing for an ice flux through the glacier's terminal boundary), which results in a systematic underestimation of volumes of marine-terminating glaciers and their thicknesses at the front (Recinos and others, 2019).

So far only one effort to estimate the impact of frontal ablation on global glacier mass change projections has been undertaken, though not focusing on frontal dynamics (Huss and Hock, 2015). However, that work neglected some processes that might be important for modeling marine-terminating glaciers’ volume evolution. Most prominent in this regard is ice dynamics, and particularly sliding, since marine-terminating glaciers’ fronts can be sensitive to dynamic thinning at the front (McFadden and others, 2011). Additionally, Huss and Hock (2015) pooled frontal ablation and surface mass balance in their approach to estimate surface elevation changes, although the mass removal by frontal ablation acts horizontally at the glacier front, as opposed to vertical changes by surface melt.

Here, we use a numerical model capable of simulating ice dynamics in a simplified flowline fashion (Maussion and others, 2019). This approach allows the model to be run on a large number of glaciers using reasonable computing resources. We configure the model in such a way that it is able to capture important features of marine-terminating glaciers’ behavior. For this purpose, we incorporated a frontal ablation parameterization (Oerlemans and Nick, 2005), water-depth dependent basal sliding and the hydrostatic stress balance at the front into both the ice thickness inversion as well as the dynamical core. We then calibrate the model on a glacier-per-glacier basis for all glaciers in the Randolph Glacier Inventory (RGI; Pfeffer and others, 2014), disaggregating the total glacier mass budget into surface mass balance and frontal ablation, where applicable, using independent datasets for each. Thereby we are able to constrain the involved parameters for both mass budget parts separately.

Accounting for frontal processes in the modeling chain allows us to investigate the differences to model runs ignoring it. This is to get an impression of such processes’ relevance for glacier mass change projections and to identify further development prospects for large-scale glacier models. We do so by running the modeling chain once without special treatment of marine-terminating glaciers, as was the standard setup of the model used previous to this work, and once including the aforementioned aspects of marine-terminating glaciers’ dynamics. Furthermore, we examine the uncertainty in such findings caused by the variation of particular unknown parameters.

2. Numerical model

2.1 The Open Global Glacier Model

The Open Global Glacier Model (OGGM) is a flowline model capable of modeling large numbers of glaciers at once (Maussion and others, 2019). Due to the fact that it relies on certain (simplifying) assumptions, this model has a reasonable computational cost. Here, we give a brief overview of the model's functionality.

The RGI (Pfeffer and others, 2014; RGI Consortium, 2017) is the basis of OGGM, similar to several other global glacier models (Marzeion and others, 2020). In its recent version (RGI V6), the coordinates and outlines of ~210 000 glaciers worldwide, which are not connected to the ice sheets, are recorded. The glacier outlines are projected onto a local gridded map for each glacier. Topographical data, based on an appropriate digital elevation model (DEM), is automatically retrieved depending on the glacier's location and interpolated onto the local grid. The grid's spatial resolution is scaled to the square root of the glacier area, with a maximum of 200 m and a minimum of 10 m. Here, we use single, binned elevation-band flowlines, based on the approach described by Werder and others (2020). Dynamical simulations start at the date a glacier was recorded in the RGI. The initial geometry consists of the surface area given by the RGI, and the result of the ice thickness inversion, which will be described in Section 2.4. Simulations before the RGI date are only possible with fixed geometries, since it is generally not possible to find glacier states before the RGI date without large computational effort (Eis and others, 2019, 2021).

The gridded meteorological dataset (monthly temperature and precipitation) is interpolated to the glacier location in a nearest-neighbor manner. For this work the Climatic Research Unit Time-Series dataset version 4.03 (CRU TS 4.03; Harris and others, 2020) is used. The temperature is subsequently corrected using a globally fixed linear lapse rate (6.50°C km⁻¹). For precipitation we do not apply a lapse rate, but a global correction factor (see below). A glacier's monthly surface mass balance for grid point i at elevation z _i is then calculated at each of the flowline's grid points as:

(1)$$m_i( z) = f_{\rm p} P^{\rm solid}_i( z) - \mu \max( T^{\rm m}_i( z) ,\; 0) $$

where f _p is a dimensionless precipitation factor (Giesen and Oerlemans, 2012), $P^{\rm solid}_i( z)$ the solid precipitation (in millimeter water equivalent (mm w.e.)) that is calculated assuming threshold temperatures for solid and liquid precipitation, μ the surface temperature sensitivity (in mm w.e. K⁻¹) and $T^{\rm m}_i( z)$ the temperature above the threshold for ice melt at the glacier surface (in K). In Section 3.2.1 we will further elaborate on the calibration of μ and other involved parameters. Former versions of OGGM relied on an interpolation approach for the surface mass balance calibration, since observational data were sparse and not available on every single glacier, but we are now able to calibrate on a glacier-per-glacier basis (see Section 3). In the following sections we will describe the formulation of ice dynamics and ice thickness inversion in greater detail, since these two aspects of the model were subject to the most changes from Maussion and others (2019) in this work.

2.2 Modulation of the shallow ice approximation for terminal cliffs

In OGGM, the thickness-averaged deformation velocity u _d of a glacier, utilizing the shallow ice approximation (SIA), is computed as follows:

(2)$$u_{\rm d} = {2A\over n + 2} h\tau^n$$

where A is the temperature-dependent ice creep parameter (here we use the default value of 2.4 × 10⁻²⁴ s⁻¹ Pa⁻³), n the exponent of Glen's flow law (here we use n = 3), h the ice thickness (in m) and τ the basal shear stress (in Pascal), which can be approximated as follows:

(3)$$\tau = \rho g h \alpha$$

where ρ is the ice density (here we use 900 kg m⁻³), g the gravitational acceleration (9.81 m s⁻²) and α the surface slope computed numerically along the flowline on a staggered grid:

(4)$$\alpha = {z_{i} - z_{i + 1}\over \Delta x}$$

where z _i is the surface elevation of that gridcell (in m), and Δx the size of the grid that is defined on the glacier (in m).

Equation (4) indicates that α can become arbitrarily large at a glacier's terminus, if there is a discontinuity of the ice thickness. This situation occurs in the presence of a terminal cliff. For a glacier whose thickness decreases smoothly to zero toward the terminus, velocity at the terminus is mostly very small, since the ice thickness (h) and the surface slope (α) are small. Hence, the SIA, conveyed here in Eqns (2) and (3), holds well and a change in Δx will essentially not affect the dynamics at the glacier front. However, for a marine-terminating glacier with a terminal cliff, the grid size dependency can be noticeable, since a halving of Δx will result in an eightfold increase in velocity (see Eqns (2–4)), with ice thickness at the terminus not being negligibly small. In that case, the higher stress would be distributed over a smaller volume of the glacier as well, as it only acts on the last gridcell in the SIA formulation we applied. To tackle these issues of large changes in the stress balance when modeling glaciers with a terminal cliff, we introduce the hydrostatic pressure balance at the terminal boundary of a marine-terminating glacier as an additional force F _H governing frontal dynamics, similar to the approach of Howat and others (2005). This additional force (per unit width) can be calculated as:

(5)$$F_{\rm H} = {1\over 2} g ( \rho_{\rm i} h_{\rm f}^{2} - \rho_{\rm o} d_{\rm f}^{2}) $$

where h _f is the thickness of the glacier front, ρ _i the density of ice (as above), ρ _o the density of ocean water (here we use 1028 kg m⁻³) and d _f the water depth at the glacier front (in m). The additional force is then distributed over a distance (L _F) inland. In order to emulate that this force acts more strongly at the boundary than further upstream, we apply the following weight $w_{i_L}$ for gridcells within L _F:

(6)$$w_{i_L} = {2 i_L\over n_L + 1}$$

where i _L is a gridcell within the n _L gridcells contained in L _F, and for the first gridcell inland i _L = n _L. Then the stress that is thereby added to the driving stress τ (see Eqn (3)) at gridcell i _L is:

(7)$$\tau_{\rm H} = w_{i_L} {F_{\rm H}\over L_{\rm F}}$$

Hence, the mean additional stress is:

(8)$$\overline{\tau_{\rm H}} = {F_{\rm H}\over L_{\rm F}}$$

and it is ensured that:

(9)$$\sum_{i_L = 1}^{n_L} w_{i_L} {F_{\rm H}\over L_{\rm F}} \Delta x = F_{\rm H}$$

This formulation allows for application of the SIA on marine-terminating glaciers with a terminal cliff, but a surface slope is still needed for computing the average velocity through the boundary of the last gridcell (see Eqns (2) and (3)). We approximate it simply as the mean slope over L _F from the front inland. Our approach introduces L _F as a new parameter, similar to a stress coupling length (Enderlin and others, 2016), which is hard to constrain for individual glaciers. Although theoretically there is a dependence on a glacier's ice thickness at the front, we set it to 8 km for all glaciers here. The motivation for this choice is that this value should be higher than the ~4–6 ice thicknesses found by Enderlin and others (2016), ensuring numerical stability in all cases. If a glacier's length (L) is smaller than 8 km, L _F = L. Note that it would be possible in our approach to also include sea-ice or ice melange backpressure in Eqn (5) (Robel, 2017), but we chose to neglect this here for simplicity. Note also that Eqn (6) implies a linearly decreasing additional stress upstream of the glacier front. It would be possible to change this to some kind of non-linear weighting function, if such was found to better represent the physics of the process, which we do not examine closer here. Furthermore, it should be noted that we used a different approach to incorporating the hydrostatic stress imbalance than previous works (e.g. Nick and others, 2009; Enderlin and others, 2013). In our case the hydrostatic stress imbalance is integrated over the glacier–water boundary to get the additional driving force, which is then distributed over a distance upstream and added to the driving stress of the SIA. In the cases of aforementioned previous works, the authors integrated the momentum–conservation equations to find a velocity gradient. Those works used the Nye–Glen rheology to directly calculate the gradients in longitudinal stress needed for their approach from the hydrostatic pressure difference at the calving face.

2.3 Sliding parameterization

At a marine-terminating glaciers’ front, the ice velocity induced by sliding is a relevant part of the dynamics (Benn and others, 2007). Previously, the sliding velocity in OGGM was calculated as:

(10)$$u_{\rm s} = {\,f_{\rm s} \tau^n\over h}$$

where f _s is a sliding parameter (default value 10⁻²⁰ m² s⁻¹ Pa⁻³), based on Oerlemans (1997). This parameterization follows the assumption that sliding is related to basal pressure, which itself is related to the ice thickness h. Basal pressure at the front of marine-terminating glaciers is not only related to the ice thickness though, but to the water depth of the glacier's bed as well. Therefore, we now calculate the sliding velocity as:

(11)$$u_{\rm s} = {\,f_{\rm s} \tau^n\over h^\ast }$$

where h* is the height above buoyancy:

(12)$$h^\ast = h - {\rho_{\rm o}\over \rho_{\rm i}} d$$

Hence, sliding for all gridcells with a bed elevation above the water level will be the same as in Eqn (10). For gridcells close to the front, the sliding velocity can sometimes become too large when using the value for f _s proposed by Oerlemans (1997) (5.7 × 10⁻²⁰ m² s⁻¹ Pa⁻³) in Eqn (11), resulting in numerical instabilities. Therefore, we use that value as an upper bound in the attempt to quantify parameter sensitivity in Appendix B, and apply a value of 10⁻²⁰ m² s⁻¹ Pa⁻³ here. Although this formulation might be an improvement for marine-terminating glaciers, the appropriate sliding parameterization for ice flow is generally still not ascertained (Benn and others, 2007; Stearns and Van der Veen, 2018; Zoet and Iverson, 2020).

2.4 Ice thickness inversion

For consistency, the changes to OGGM explained in the previous sections do not only need to be incorporated into the dynamical model core, but into the ice thickness inversion as well. That means we now numerically solve for the ice thickness by the following polynomial:

(13)$$q = {2A\over n + 2} ( f_{\tau} \tau ) ^n h^2 + f_{\rm s} ( f_{\tau} \tau) ^n r_h$$

where q is the ice mass flux per unit width (in m² a⁻¹), f _τ an amplification factor related to the additional driving stress caused by the hydrostatic stress balance at the front (see Section 2.2, Eqn (7)):

(14)$$f_{\tau} = {\tau_{\rm H}\over \tau} + 1$$

and r _h the inverse of the relative height above buoyancy:

(15)$$r_h = {h\over h^\ast }$$

For gridcells that are more than L _F (see Section 2.2) upstream of the front and that do not have a bed elevation below the water level, this equals the ice thickness inversion approach usually applied in OGGM, because f _τ and r _h equal one.

A further peculiarity of marine-terminating glaciers that needs to be taken into account here is the occurrence of frontal ablation. We parameterize the latter following Oerlemans and Nick (2005):

(16)$$Q_{\rm f} = k d_{\rm f} h_{\rm f} w_{\rm f}$$

where Q _f is the frontal ablation flux (m³ a⁻¹), k the water-depth sensitivity parameter (a⁻¹; hereafter named frontal ablation parameter) and d _f, h _f and w _f are the water depth, ice thickness and width at the glacier front (all in m). As introduced to OGGM by Recinos and others (2019), the mass budget closure requires Q _f to be balanced by ice discharge, which is the dynamical ice flux through the terminal boundary of the inverted glacier. That, in turn, implies a larger ice thickness at the glacier's terminus compared to an inversion neglecting frontal ablation, assuming the same ice flow parameters (see Eqn (13)). Although simple and thus readily applicable, this frontal ablation parameterization is limited insofar as it does not explicitly capture physical processes, but lumps them into one parameter (k). The two main processes relevant in that regard are brittle fracturing and submarine melt, which is related to ocean/water temperatures and subglacial discharge. Note that our framework is also not capable of capturing sediment dynamics and hence the time evolution of proglacial submarine moraines that might influence frontal dynamics/ablation (see, e.g. Oerlemans and Nick, 2006; Brinkerhoff and others, 2017).

One further change compared to previous studies on the ice thickness inversion of marine-terminating glaciers is that we neither assume the water level to necessarily be at 0 m a.s.l., nor prescribe the freeboard to be within a range of 10–50 m (Recinos and others, 2019, 2021). This change is motivated by several factors: (1) the DEMs used can be erroneous, (2) the RGI outlines can be erroneous and the incorrect geometry (i.e. width) at the front derived from these outlines for the elevation-band flowlines can deteriorate the result of the ice thickness inversion, (3) the assumed values for the flow parameters (A, f _s) in Eqn (13) and the frontal ablation parameter (k in Eqn (16)) are uncertain and thus the water level may have to be shifted in order to satisfy Eqn (16) and find a solution for Eqn (13) and (4) a maximum freeboard of 50 m would mean that ice thicknesses at the front could not exceed ~400 m without going into flotation. However, dealing with floating tongues would involve shelf dynamics, which we cannot model using the SIA without special treatment and a refined grid at the front (Vieli and Payne, 2005). Therefore, we inhibit glacier states that feature a floating tongue. Thus, we do not directly seek for a water depth at the front, which ensures that the ice flux given by the frontal ablation parameterization equals that of the apparent mass balance and ice dynamics, but for a value of r _h (Eqn (15)) as explained below. By doing so, we can make sure that the ice thickness inversion never results in a floating tongue, and with a value for the freeboard this can be translated into the frontal ice thickness:

(17)$$h_{\rm f} = {( {\rho_{\rm o}}/{\rho_{\rm i}}) r_h ( z_{\rm f} - z_{\rm w}) \over r_h ( ( {\rho_{\rm o}}/{\rho_{\rm i}}) - 1) + 1}$$

where z _f is the surface elevation of the last gridcell according to the DEM used, and z _w the water level (all in m). If the initial guess for the freeboard (i.e. z _w = 0 m above sea level) results in an error of the numerical solver, we shift the water level until the algorithm successfully finds a value for r _h and thereby for h _f. Tuning of the initial freeboard estimate means that we could calculate different water levels for the same glacier for different values of the flow parameters (A, f _s), and the frontal ablation parameter (k). Here it should be noted that a shift of the water level does not imply a shift of the surface elevations recorded in the DEM. It is merely a numerical attempt to allow for the compensation of inconsistencies, which may arise from model approximations and errors/uncertainties in the observational data. Thereby we are able to ensure a consistent solution of the ice thickness inversion for every glacier with any given set of parameters. For the majority of glaciers shifting the water level is not necessary and if it is, the shift is often relatively small (<100 m in 90% of the cases).

2.5 Frontal ablation in the dynamical model

We apply the same frontal ablation parameterization in the dynamical model as in the ice thickness inversion procedure (see the above section). Since frontal ablation does not act vertically, as surface melt does, one has to decide how to remove the volume calculated with Eqn (16) from the gridded glacier in a time-stepping scheme. For that, we use two reservoirs: one is the temporally accumulated frontal ablation flux (Q _f) and the other one is the temporally accumulated ice flux through the terminus cross section (Q _t). Note that here, unlike in Eqn (16), Q _f ≠ Q _t, for in the dynamical model we do not assume a steady state situation. Then, in every time step we remove the accumulated Q _f from the accumulated Q _t. If the remaining Q _f is still large enough, entire gridcells can be removed from the front. Vice versa, if Q _f < Q _t over a certain time interval, the accumulation of Q _t can lead to an advance of the glacier. Furthermore, if the thickness of one or more gridcell(s) falls below flotation in a certain time step, the part of this volume which is contained in gridcells beyond the one adjacent to the last gridcell above flotation is removed and added to the frontal ablation output variable. That is in order to prevent shelf dynamics (see above). Because most marine-terminating glaciers outside of the ice sheets do not have a floating tongue anymore (Copland and Mueller, 2017), we do not anticipate that neglecting shelf dynamics will have significant influence on our results. For future model developments it might be considered to incorporate stress-related criteria, linked to the findings of Bassis and Walker (2011), to confine the height above buoyancy in the ice thickness inversion and dynamical model.

3. Data and calibration

3.1 Data

3.1.1 Mass change above sea level

We use the mass changes estimates for each glacier in the RGI over 2010–20 provided by Hugonnet and others (2021). Since reanalysis data of meteorological conditions are available over the same time period, it is possible to calibrate surface mass balance models for each glacier individually. The data of Hugonnet and others (2021) are based on computing differences in the surface elevations of glaciers, derived from DEMs, between different points in time. However, it does not include any mass changes occurring below sea level, since the satellite data it is based upon can only detect changes above the water line. This is problematic for estimating total mass changes and when calibrating models of marine-terminating glaciers, because part of the mass budget would be disregarded, if it is not corrected for.

3.1.2 Frontal ablation and mass change below water level

We use the satellite-derived data from Kochtitzky and others (2022) for estimating mass changes below sea level and to obtain frontal ablation estimates. Doing so is necessary to prevent an erroneous calibration of the surface mass balance model described above (Section 2.1), which will be further elaborated on in Section 3.2.1. The necessity of using mass changes below sea level in addition to frontal ablation estimates arises from the fact that frontal ablation estimates include the mass change below sea level, while the Hugonnet and others (2021) data are not able to capture it, which would lead to inconsistencies.

The dataset given by Kochtitzky and others (2022) is largely derived from satellite observations of velocity, glacier area and observed and modeled ice thickness. It does not include estimates for the Southern Hemisphere and lake-terminating glaciers. We thus solely focus on marine-terminating glaciers (outside the Greenland ice sheet) in the Northern Hemisphere. Moreover, we exclude Flade Isblink Ice Cap from our work here, because it possesses problematic outlines in the RGI and a floating tongue, which we are currently not able to model with OGGM (Recinos and others, 2021; Müller and others, 2022). As it also contains a lot of ice, we rather neglect it here than largely distort our (regional) results by modeling it erroneously. We acknowledge that there are some smaller glaciers with a floating tongue in far northern Canada and Greenland, but their impact on our results should be minor, as their ice volume is small (Copland and Mueller, 2017) and many floating tongues have collapsed in the last few decades (White and Copland, 2019; Kochtitzky and Copland, 2022).

3.1.3 Atmospheric forcing data

The temperature and precipitation data needed as boundary conditions for the surface mass balance model are taken from the CRU TS 4.03 dataset for historical runs (from the RGI record date of a glacier until 2020) and from different general circulation models (GCMs) from the CMIP6 archive for projections (Eyring and others, 2016). In order to avoid potential step changes in the atmospheric forcing when switching from the reanalysis data to GCM output data, OGGM features a function that adjusts the GCM data based on its bias to a climate period in the historical data (here we used 1981–2018). A list of used GCMs, and for which scenarios they provide data, is given in Table A1. The CRU TS dataset has a spatial resolution of 0.5°, and the range of the GCM ensemble's spatial resolution is 0.75–2.0°. As calibrating the surface mass balance parameterization is done for 2010–20 (see Section 3.2.1), the CRU data are used for that purpose as well.

3.2 Calibration

3.2.1 Calibration of surface mass balance

The arrival of abundant satellite-derived data for glaciers worldwide offers the new possibility to calibrate glacier-specific parameters without the need to interpolate for glaciers with no observational data, or use regionally aggregated values (as, e.g. in Marzeion and others, 2012, Radić and others, 2014 or Huss and Hock, 2015). Since we have observational values as well as meteorological data from reanalysis now for each individual glacier, the sensitivity to atmospheric temperatures in the temperature index surface mass balance model of OGGM (see Eqn (1)) can be approximated as follows:

(18)$$\mu = \left(\,f_{\rm p} P_{\rm solid} - {\Delta M_{\rm awl} + C + f_{\rm bwl} \Delta M_{\rm f}\over A_{\rm RGI}}\right){1\over T_{\rm m}}$$

where:

• • ΔM _awl is the observed annual mass balance above sea level of a glacier (mm³ a⁻¹) as given by Hugonnet and others (2021),

• • C is the observed annual frontal ablation rate of a glacier as given by Kochtitzky and others (2022) (mm³ a⁻¹),

• • ΔM _f is the observed annual rate of mass change due to area changes in the terminus region of a glacier (mm³ a⁻¹; hereafter named retreat volume) as given by Kochtitzky and others (2022),

• • f _bwl is the assumed fraction of ΔM _f occurring below the waterline

• • A _RGI is the glacier surface area of a glacier as given by the RGI (mm²),

• • T _m is the annually accumulated temperature above the threshold for ice melt at the glacier surface (K a⁻¹).

The second term in the parentheses consequently represents the observed specific surface mass balance, and f _bwlΔM _f the mass balance below sea level (ΔM _bwl). Note that our convention assigns positive frontal ablation to the case of mass removal, and we neglect the case of a positive frontal mass budget in the calibration as well as in the dynamical model.

As indicated in Section 2.1, the value of μ in Eqn (18) depends on the global parameter values chosen for f _p (2.5), the threshold temperatures for liquid/solid precipitation (2/0°C), and the threshold temperature for ice melt (−1°C). Here, we adopted parameter values that were previously derived from a leave-one-glacier-out cross-validation procedure with annual in situ mass-balance measurements (Maussion and others, 2019). Because the data of Hugonnet and others (2021) have large uncertainties associated with it for annual values of individual glaciers, it makes most sense to use an average over longer time intervals, although it is not possible to constrain interannual variability in that way. Here we use the time period 2010–20. The motivation for this is the assumption that most of the recording dates in the RGI lie before that interval and thus potential spin-up effects caused by assumptions in the ice thickness inversion procedure are attenuated (see Section 5.3). We also ignore the fact here that the uncertainty for longer time intervals given by Hugonnet and others (2021) can nevertheless be quite large for individual glaciers, since we focus on large spatial scales and thus assume that uncertainties of individual glaciers will, at least partially, cancel each other out.

We use estimates of mass changes due to changes in the terminus area as well as of mean frontal ablation rates (ΔM _f, and C in Eqn (18)), given by Kochtitzky and others (2022) for the time interval 2010–20, for each marine-terminating glacier in the Northern Hemisphere outside the Greenland ice sheet (see Section 3.1.2). Hence, we can adjust the Hugonnet and others (2021) data by assuming that 75% of the mass change in the terminus area happens below sea level and adding it to the mass balance (f _bwl = 0.75 in Eqn (18)). Although the assumption of 75% is arbitrary, this part of the total Northern Hemisphere's glacier mass budget is only about one-fifth of the total frontal ablation (see Table 1), and thus a change from 75% to, for instance, 87.5% (flotation), will presumably not make a large difference. We investigate the implications of that assumption further in Appendix B.

Table 1. Mass budget components of marine-terminating glaciers as annual mean values over 2010–20 estimated for different RGI regions and the Northern Hemisphere given by Kochtitzky and others (2022) (see Section 3.1.2)

ΔM _f is the mass change due to area changes at the glaciers’ fronts (retreat volume), C the estimated annual frontal ablation used for model calibration, A _mt the marine-terminating glacier area in the region given by the RGI and n the number of marine-terminating glaciers in that region. The number of glaciers with the largest observed frontal ablation rates, together contributing 90% of the regions’ total frontal ablation, are given in parentheses after n.

3.2.2 Calibration of frontal ablation parameterization

In order to disaggregate the total mass change of marine-terminating glaciers, we do not only calibrate the surface mass balance, but the frontal ablation parameter (k in Eqn (16)) as well. For this, we take satellite-derived estimates of the average frontal ablation rates over 2010–20 (see Section 3.1.2). We use an iterative procedure which seeks for a value of k that produces a frontal ablation estimate over this time interval within the uncertainty bounds of the frontal ablation data used. During that procedure, we use the same value for k in the ice thickness inversion and a subsequent (historical) dynamical run, which is forced with CRU TS 4.03 meteorological reanalysis data.

3.3 Setup of model runs

Here we briefly describe how we set up the different types of projections compared in the next section. For the projections including frontal processes of marine-terminating glaciers, we first calibrate the glaciers’ sensitivity to atmospheric temperatures (μ in Eqns (1) and (18)) and then the frontal ablation parameter (k in Eqn (16)) as explained in the previous section. Following that we apply the ice thickness inversion procedure and run the model for a historical period starting at the individual glacier's RGI recording date and ending in 2020 with CRU TS 4.03 data as atmospheric boundary conditions. From there we switch to the individual members of the GCM ensemble given in Table A1 as the source of atmospheric boundary conditions and run the model for each member until 2100. The projections we compare these to are conducted in a similar manner, but exclude C and f _bwlΔM _f in Eqn (18) in the μ calibration. Furthermore, in these runs we neither include frontal ablation in the ice thickness inversion nor in the forward runs. This means that in those projections marine-terminating glaciers are treated as if they were land-terminating, and thus like in previously published OGGM projections. For results labeled Northern Hemisphere, the model is run for the RGI regions 1–15.

4. Results

4.1 Calibration/ice thickness inversion

Inspecting the observational estimates of the frontal mass budget used for the calibration given in Table 1, it is visible that most of the frontal ablation occurred in Svalbard and the Russian Arctic. Although ~66% of the estimated frontal ablation between 2010 and 2020 occurred there, Svalbard and the Russian Arctic account for only 39% of the marine-terminating glacier area considered here and for 38% of the number of marine-terminating glaciers in the Northern Hemisphere (n in Table 1). In contrast, Arctic Canada North covers ~35% of the marine-terminating glacier area, but accounts for only 16% of the number of marine-terminating glaciers and for only 8% of the annual frontal ablation. About 20% of the Northern Hemisphere's frontal ablation was observed in Alaska, but that region only accounts for 3% of the marine-terminating glaciers and covers only 7% of the marine-terminating glacier area. In contrast to Alaska, the Greenland periphery shows a high count of marine-terminating glaciers, and percentage of the region's volume contained in these (36 and 37%), while only accommodating 6% of the Northern Hemisphere's frontal ablation. Despite marine-terminating glaciers containing 8%, and 14% of the ice volume in Arctic Canada South, and Iceland, respectively, the total amount of frontal ablation in these two regions is so small that we neglect them in the further analysis. Table 1 additionally shows that in most regions, between 23 and 33% of all marine-terminating glaciers contribute 90% of the regions’ total frontal ablation. Only in Alaska is this number notably lower at 15%. For the whole Northern Hemisphere this value is also <20%, illustrating that most of the findings that will be discussed in the following sections are caused by a rather small number of glaciers.

Table 2 displays the results of the frontal ablation parameterization's calibration and the ice thickness inversion. The highest percentage of estimated regional ice volume stored in marine-terminating glaciers is found in Svalbard and the Russian Arctic, also where most of the Northern Hemisphere frontal ablation was observed. In most regions the area-weighted average of the frontal ablation parameter (k) lies between 0.15 and 0.63; the only region (besides Iceland) that markedly falls out of this range is Alaska. This might point to the high turnover of Alaskan glaciers, since in Svalbard and the Russian Arctic we found much smaller values for k, but the reduction of the temperature sensitivity is quite similar; another hint at this is the rather high frontal ablation but low retreat volume per glacier given in Table 1. In some individual regions (e.g. Greenland periphery) and the whole Northern Hemisphere, the area-weighted std dev. of k is larger than the average, which indicates rather heavy-tailed distributions, since k cannot be negative. The same can be seen for the reduction of glaciers’ melt sensitivities to atmospheric temperatures (Δμ _mt).

Table 2. Results of the frontal ablation parameterization's calibration and ice thickness inversion for different RGI regions and the Northern Hemisphere

A _mt is the percentage of regional area covered by marine-terminating glaciers according to the RGI, V _mt the estimated percentage of volume contained in marine-terminating glaciers, ΔV _mt the change in volume of marine-terminating glaciers due to including frontal ablation in the ice thickness inversion, Δμ _mt the average and 1 − σ std dev. (weighted by surface area) of the change in the glaciers’ melt sensitivities to atmospheric temperatures due to including frontal ablation and retreat volume in the surface mass balance calibration, k the average and 1 − σ std dev. (weighted by surface area) of the frontal ablation parameter, bias the mean difference to the observational estimates of annual frontal ablation and RMSE the root mean squared error of the modeled frontal ablation with respect to the calibration dataset. Note that values were computed using the default parameter values in Table B1 and exclude Flade Isblink IceCap.

Furthermore, Tables 1 and 2 show that the regions with the largest observed frontal ablation rates also exhibit the largest (1) increases in computed volume and (2) reductions in sensitivities to atmospheric temperatures when the inversion and calibration account for frontal ablation. These changes are attenuated by the frontal volume lost during retreat. As we assume that 75% of ΔM _f happens below sea level for marine-terminating glaciers and is hence not captured by the data of Hugonnet and others (2021), we add a corresponding amount to the total mass budget (Section 3.2.1 as in Eqn (18)). Therefore, if the difference between the absolute value of retreat volume and total frontal ablation (per glacier) is low, the reduction of the surface temperature sensitivity will be rather low. Accordingly, a higher reduction in sensitivity to atmospheric temperatures will allow for more additional ice flux (due to less surface melt) and thus result in more ice volume in the ice thickness inversion.

Our calibration routine succeeds for most of the glaciers (n = 1322 out of 1367). Table 2 also shows that in most regions the bias and root mean squared error (RMSE) are low. The bias is smaller than 5% and the RMSE below 1% in most regions. An unsuccessful calibration can be either caused by model errors on a glacier (n = 3), or an irreducible difference between the modeled and the observed frontal ablation being higher than the estimated uncertainty of the calibration dataset. The latter case can be caused by erroneous RGI outlines and/or DEMs yielding erroneous frontal geometries. Otherwise, erroneous calibration data, or the application of inadequate parameter values (e.g. for the precipitation scaling factor or flow parameters) could inhibit a realistic modeling of the affected glaciers, leading to a failing calibration. Seeing that on the hemispheric scale the bias is 2% and the RMSE is close to that of the calibration dataset provides some confidence in the reliability of the results.

4.2 Projected frontal ablation

Figure 1 displays the projected frontal ablation accumulated over the 21st century (using the conversion of 1 mm sea level equivalent (SLE) ≈362.5 gigatons (Gt)). It also shows these values excluding the amount of volume lost below the water level that has to be replaced by fresh water in order to cause GMSLR, which indicates the regionally aggregated stability of front positions. Moreover, it can be seen that the slope of the projected accumulated frontal ablation for the Northern Hemisphere decreases over the course of the 21st century, which means that the annual rate of frontal ablation decreases. It is noteworthy that the projected amount of frontal ablation is almost independent of the applied climate scenario/atmospheric forcing. We interpret that this is due to two different effects of the atmospheric climate on frontal ablation that cancel each other out to a certain degree in the Northern Hemisphere. Those are shown for two example glaciers in Figure 2. The left-hand panels (a) and (c) show the exemplary behavior of a glacier which experiences less frontal ablation in higher emission scenarios, more typical for the Greenland periphery and Alaska. Such a behavior is determined by a mostly continuous prograde slope of the bed topography, which means that with increased surface melt the glacier will experience less frontal ablation. That is because its thinner front is situated in shallower water in the same years and our frontal ablation parameterization explicitly depends on these two variables (h _f and d _f in Eqn (16)). The glacier shown on the right, panels (b) and (d), represents a behavior more typical to the region Arctic Canada North, meaning more frontal ablation in a higher emission scenario (see Fig. 1). This behavior is characterized by glaciers being close to flotation and a retrograde bed slope. Higher surface melt thus triggers a rapid retreat when such glaciers are forced beyond a bathymetric pinning point.

Fig. 1. Projected accumulated frontal ablation for the Northern Hemisphere and different RGI regions for the remainder of the 21st century. Dashed lines exclude the amount of volume loss below the water level that has to be replaced by fresh water and thus does not contribute to GMSLR (ρ _fw is the density of fresh water: 1000 kg m⁻³). N is the number of GCMs used for the respective emission scenario. The 1 − σ std dev. of the GCM ensemble for scenarios SSP1 2.6 and SSP5 8.5 is displayed as shading. The right y-axes display gigatons (Gt) converted to mm SLE. Note the different scales for each region.

Fig. 2. Lateral view of two example glaciers’ states at different points in time (5 year increments from 2020 to 2100), simulated using the BCC-CSM2-MR climate model data for atmospheric boundary conditions. The top panels (a and b) show the glacier states in simulations forced with the emission scenario SSP1 2.6, and the bottom panels those in simulations forced with SSP5 8.5. Note the different scales for both glaciers.

Although we find these two counteractive effects of atmospheric warming on projected frontal ablation, the total projected amount is slightly larger in the higher emission scenarios for the whole Northern Hemisphere. This is probably due to the flotation criterion we implemented being triggered more often in higher emission scenarios owing to enhanced surface thinning. It is likely that there is a positive feedback of surface melt and thinning at marine-terminating glacier fronts, as the melt-induced surface thinning decreases the height above buoyancy and thereby increases the sliding velocity (see Eqns (11) and (12)). If this increases the total frontal velocity in the absence of increased mass transport to the terminus, it will inevitably lead to further thinning and retreat.

Concerning regional frontal ablation estimates, it is interesting that there is little projected mass change below the water level in Alaska. This indicates that front positions of glaciers contributing to the largest amounts of frontal ablation in this region are mostly stable in our simulations, while the total number of marine-terminating glaciers is decreasing significantly (see Fig. 3). Consistent with what we described above, Arctic Canada North is the only region for which we project slightly higher frontal ablation rates toward the end of the century than in 2020 (increasing slope in Fig. 1). Furthermore, we find rapidly increasing frontal ablation rates for this region already in the last few years of the calibration period (2010–2020; not shown in Fig. 1). This leads to the cumulative total frontal ablation estimates being about twice what we found for Alaska by 2100, where frontal ablation rates decrease throughout the first half of the 21st century. Additionally, we project more accumulated frontal ablation by 2100 in the higher than in the lower emission scenarios for this region. Cook and others (2019) found that marine-terminating glacier retreat patterns in the Canadian Arctic Archipelago showed no significant correlation with subsurface ocean temperature changes until 2015. Assuming the oceanic influence in that region remains small in the future and acknowledging that we cannot model the effect of oceanic melt in our current modeling framework, an explanation would be that a significant number of marine-terminating glaciers there are close to flotation today (and have undulating/retrograde bed slopes). Such a geometric configuration makes them prone to rapid retreat and thereby enhanced frontal ablation when experiencing stronger surface melt/thinning. Our simulations result in more frontal ablation throughout the 21st century under the higher emission scenarios in Svalbard and the Russian Arctic as well. On the other hand, we project less accumulated frontal ablation by 2100 in the higher than in the lower emission scenarios for Alaska and the Greenland periphery.

Fig. 3. Number of glaciers with volume below the water level for the Northern Hemisphere and different RGI regions. The uncertainty in the GCM ensemble is displayed as dotted black lines. The upper line is the positive 1 − σ std dev. for the scenario SSP1 2.6, and the lower line the negative 1 − σ for SSP5 8.5. Note the different scales for each region.

The number of glaciers in each region containing volume below the water level in a certain year of our simulations is displayed in Figure 3. It shows that this number depends on the atmospheric forcing applied, with higher emission scenarios leading to less marine-terminating glaciers in general. The most drastic decrease is simulated for the Greenland periphery, consistent with the relatively large retreat volume given in Table 1. The least (relative) decrease is simulated for Svalbard. Interestingly, the GCM ensemble's std dev. and the spread between emission scenarios differ among the regions shown in Figure 3, indicating differences in the sensitivity of marine-terminating glaciers’ dynamic response to atmospheric forcing.

4.3 Northern Hemisphere and regional mass change

The main aim of this work is to compare projections of glacier mass change including marine frontal processes with projections disregarding these. When comparing such results concerning GMSLR, three effects will play a role: (1) the surface temperature sensitivity of marine-terminating glaciers is reduced and therefore their response to temperature changes is dampened, (2) frontal ablation takes place and compensates a certain amount of the dampened surface ablation and (3) retreating marine-terminating glaciers lose volume below the water level, which does not contribute to GMSLR because it already displaced ocean water. Hence, the amount of volume loss below the water level has to first be translated to freshwater volume and then be subtracted from total glacier mass loss when calculating GMSLR contributions. Figure 4 shows the total projected accumulated glacier contribution to GMSLR, including marine frontal processes and accounting for volume changes below water level. The values are relatively similar to what has been previously presented by Marzeion and others (2020) for the regions displayed. Table 3 shows the estimated cumulative GMSLR contribution at the end of the century given by our two types of projections as well as by the OGGM projections published by Marzeion and others (2020). It is visible that for most of the regions considered here the glacier-by-glacier surface mass balance calibration does not significantly alter the estimates. Only for Svalbard larger differences are recognizable. For the entire Northern Hemisphere our estimates are consistently lower than the ones previously derived from OGGM, but it is out of the scope of this work to investigate whether this comes from the different calibration method or differences in the applied atmospheric forcing ensemble.

Fig. 4. Projected accumulated glacier contribution to GMSLR for the Northern Hemisphere and different RGI regions taking frontal ablation and mass changes below water level into account. N is the number of GCMs used for the respective emission scenario. The 1 − σ std dev. of the GCM ensemble is displayed as shading for the different scenarios. Note the different scales for each region.

Table 3. Northern Hemisphere glaciers’ contribution to GMSLR by 2100 estimated by projections including frontal processes and by those not doing so as well as by OGGM projections published by Marzeion and others (2020), for different RGI regions and all glaciers in the Northern Hemisphere

Columns represent different emission scenarios. Note that the scenarios given here correspond to RCP2.6, RCP4.5 and RCP8.5 used with a different GCM ensemble in Marzeion and others (2020).

Figures 1 and 4 indicate that frontal ablation is a large part of the projected total mass budget for glaciers in Svalbard and the Russian Arctic, while it is a small part in the Greenland periphery and Alaska (see Fig. 1). Overall, our results imply that atmospheric forcing will play an increasingly dominant role throughout the 21st century for the Northern Hemisphere's glaciers (outside the Greenland ice sheet) in higher emission scenarios, as the rate of frontal ablation decreases while that of the total mass loss increases in most cases. This is also suspected to be the case for the Greenland ice sheet (Goelzer and others, 2020).

Figure 5 shows the accumulated difference in Northern Hemisphere glaciers’ GMSLR contribution, relative to 2020, between the two different types of projections. It is visible that less glacier mass loss above the water level is projected for the Northern Hemisphere when taking frontal ablation into account. This is due to the lowering of the sensitivity to atmospheric temperatures, which is not compensated by the amount of projected volume above the water level removed by frontal ablation. Understandably, such a reduction has less of an absolute effect in lower emission scenarios, since atmospheric temperatures increase less in those during the 21st century. The absolute difference in reduction of glacier mass loss contribution to GMSLR between the lowest and highest emission scenarios considered here is ~5 mm SLE (Fig. 5). Of this amount, ~2 mm SLE can be attributed to differences in mass loss below the water level; in the lower emission scenarios, glaciers retreat less and thus less mass below the waterline is lost. In the highest emission scenario a total of ~7 mm SLE are lost below the water level. Evidently, the absolute difference in projected glacier contribution to GMSLR between different emission scenarios is smaller when applying our marine-terminating glacier framework than in projections neglecting marine frontal processes.

Fig. 5. Difference in GMSLR contribution between projections including marine frontal processes and those neglecting these processes for the Northern Hemisphere and different RGI regions. N is the number of GCMs used for the respective emission scenario. The 1 − σ std dev. of the GCM ensemble is displayed as shading. Note the different scales for each region.

The only region shown in Figure 5 for which the opposite is true is the Greenland periphery, where the temperature sensitivity is only reduced by an (area-weighted) average of 8%. Regarding that region, it is also intriguing that our projections result in up to 0.5 mm SLE more mass loss above sea level when including marine frontal processes in the model, but accumulated frontal ablation that can actually contribute to GMSLR is only ~0.2 mm SLE (see Fig. 1). This suggests that other dynamical effects play a role here. Another outstanding feature of Figure 5 is the difference in projected contribution to GMSLR applying the highest emission scenario (SSP5 8.5) in the region Arctic Canada North. In that situation the difference flattens out markedly during the last third of the 21st century, which could be caused by an increase in frontal ablation due to more glaciers rapidly retreating in this region and scenario (see the section above), and/or by the disappearance of glaciers for which the decreases in sensitivity to atmospheric temperature were rather large.

Regarding percentages, the contribution of glacier mass loss in the Northern Hemisphere to GMSLR by the end of the 21st century is reduced by ~9% for the different emission scenarios (see Table 4). If we only consider Northern Hemispheric marine-terminating glaciers, the reduction is ~30%. Concerning individual regions, it can be seen in Table 4 that the reduction of mass loss above sea level reaches up to 30% in Svalbard and the Russian Arctic. The marine-terminating glaciers of Alaska lose ~60% less mass above sea level in our simulations. The other extreme is the Greenland periphery, where we estimate that marine-terminating glaciers contribute up to 8% more to GMSLR when accounting for frontal ablation.

Table 4. Percent difference in glaciers’ contribution to GMSLR by 2100 between projections including frontal processes in the modeling chain and those not doing so, for different RGI regions and all glaciers in the Northern Hemisphere

Lower rows for each region display the values only regarding marine-terminating glaciers. Columns represent different emission scenarios.

5. Discussion

5.1 Calibration/ice thickness inversion

Although a comparison with previously estimated values for k is not straightforward, because these studies either used one value for a whole region (Recinos and others, 2019) or a different calibration strategy (Recinos and others, 2021), it appears that our calibrated values are consistently higher. Additionally, the mentioned studies calibrated k only using the ice thickness inversion routine and not a dynamical model, further impeding a meaningful comparison. Concerning our ice thickness inversion, the volume increases due to including frontal ablation we found for the regions Alaska and Greenland periphery are close to previously published ones (see Table 2). Recinos and others (2019) found a 16% increase in ice volume for Alaska, and our estimate for the Greenland periphery lies between results from the two calibration methods (1.2 and 9.5%) applied by Recinos and others (2021). The larger volume increase in Alaska is related to the larger reduction in the temperature sensitivity (Δμ _mt in Table 2), which itself is related to the difference between retreat volume and total frontal ablation (ΔM _f and C in Table 1). This difference indicates how destabilized frontal positions in a region are, as it indicates whether the ice volume flux arriving at the glacier fronts can balance frontal ablation. An even larger difference (per glacier) than in the Greenland periphery was observed in Arctic Canada North, suggesting that dynamic effects are at play there, since it was found by Cook and others (2019) that glacier mass changes in that region were mostly forced by changes in atmospheric conditions until recently. It could therefore indicate that, as described above, the retreat (and a potential increase in frontal ablation) is caused by the influence of surface thinning and bathymetry on frontal dynamics and is mainly driven by an imbalance in the surface mass budget.

The interrelation of surface mass balance and frontal dynamics furthermore points to the larger problem of how to partition marine-terminating glaciers’ mass changes between surface and frontal mass budget, because such glaciers could be in disequilibrium and retreat even with a positive surface mass balance. Such an assessment of surface versus dynamic (or frontal) disequilibrium is complicated by the fact that it is a transient problem and both mass budget parts are dependent on geometric changes and hence on topography/bathymetry. Similar problems apply to the ice thickness inversion procedure we used, which generates the initial ice thicknesses at the RGI date. It has the caveat that it is a transient problem as well for glaciers not in a steady state (Rabatel and others, 2018). For simplicity, we assume steady state at the RGI date (Maussion and others, 2019), which is inaccurate, given that most glaciers had a glacier-wide negative mass balance, and were thus not in a steady state, around their RGI recording date. This problem is difficult to solve, especially in the presence of frontal ablation, but may be tackled with inverse methods relying on velocity observations as in Millan and others (2022). The equilibrium assumption might result in too much ice volume at the beginning of our simulations (the RGI recording date), thereby potentially inducing spin-up effects in the first years of simulation. Since most RGI dates lie before 2010 (~98%), the effect on our calibration should be small. On the other hand, we use the area at the RGI date in Eqn (18), which leads to incongruities, since we use the 2010–2020 mass change and frontal ablation data.

5.2 Northern Hemisphere and regional mass change

For most of the examined regions we find good agreement between our projections that do not include marine frontal processes and those previously generated using OGGM (published by Marzeion and others, 2020); only for Svalbard our estimates are considerably lower. This is despite the utilization of a glacier-by-glacier calibration in this work compared to an interpolation-based approach in the earlier work. Projections including marine frontal processes result in lower GMSLR estimates by 2100 in all regions except for the Greenland periphery, as described in Section 4.3. The total Northern Hemisphere's GMSLR contribution's reduction of ~9% could be considered rather small. Still, the projected mass changes for individual regions, and for the entirety of Northern Hemisphere marine-terminating glaciers, can be altered quite strongly by changes to the model described in the sections above.

One caveat here is that we are not able to account for increasing rates of frontal ablation due to increasing ocean temperatures in the current framework. Thus, frontal ablation rates might be underestimated for higher emission scenarios in our projections toward the end of the 21st century. Another process neglected in our simulations is the influence of subglacial discharge on frontal ablation (Slater and others, 2015). An increase in subglacial discharge by stronger surface melt over the course of the 21st century could enhance frontal ablation rates in relation to increasing atmospheric temperatures. As we do not include the influence of increases in thermal forcing and/or subglacial discharge on frontal ablation, the ratio of frontal ablation to total mass change decreases more strongly over time for higher emission scenarios. Explicitly including submarine melt, and thus oceanic forcing more comprehensively, might increase the dependence of frontal ablation projections on the climatic forcing. Such a potential increase in frontal ablation would be compensated by the fact that fewer glaciers will be marine-terminating in 2100 than in 2020, however. Kochtitzky and Copland (2022) have shown that several marine-terminating glaciers already retreated to land between 2000 and 2020.

A further limitation is that the dataset we use for calibrating the frontal ablation parameterization does not cover lake-terminating glaciers. Alaska, for example, contains quite large glaciers lake-terminating glaciers, and including them could thus further enhance the influence of frontal ablation on glacier mass change projections. The number of lake-terminating glaciers is likely to increase as glaciers retreat. Finally, we suppose that including the Southern Hemisphere would strongly enhance the impact of incorporating frontal ablation on large-scale glacier mass change projections. That is because glaciers of the Antarctic periphery store by far the most ice below sea level of all RGI regions (Farinotti and others, 2019).

5.3 Further sources of uncertainty

Glacier mass change and its interrelation with glacier dynamics is a complex issue. Examining marine-terminating glaciers adds a level of complexity, as additional processes and variables are involved compared to land-terminating glaciers. Hence, the numerical modeling of such glaciers requires additional parameters and adjustments, which are subject to uncertainties (Hourdin and others, 2017). In order to test the influence of certain parameters on the results presented in Section 4, we conducted additional experiments described in Appendix B.

External sources of uncertainties include the datasets we use as initial (RGI, DEMs) and atmospheric boundary conditions (CRU, GCMs), because they transmit their uncertainties to the model during calibration, ice thickness inversion and dynamical projections. Parameterizations also add to the uncertainty, as shown exemplarily for the flow parameters in Appendix B. Moreover, parameterizations always simplify and approximate physical processes that are not explicitly incorporated in the model formulation. Applying inappropriate parameterizations will thus introduce systematic errors to the model. For instance, a sliding parameterization that includes lubrication of the glacier bed by percolating surface melt, and therefore produces more sliding in lower parts of a glacier than in upper parts (with the same values of h, d and α), might be more appropriate than the one we applied. Also, applying the same value for certain parameters to every glacier (i.e. global parameters) and keeping parameters constant throughout the simulations are simplifications that lead to further inaccuracies. Additional assumptions and simplifications in our modeling framework that we acknowledge are: (1) negligence of floating tongues/ice shelves, lateral drag, and sediment dynamics/proglacial moraines, (2) the equilibrium assumption in the ice thickness inversion, (3) the usage of elevation-band flowlines and (4) the omission of explicit submarine melt (i.e. oceanic thermal forcing as well as subglacial discharge) in the frontal ablation parameterization.

Ultimately, it is not our aim in this work to give as accurate projections as possible of future glacier mass changes, but to get a grasp of the influence of including frontal processes on these. In that regard, the abovementioned uncertainties do not challenge our main conclusions, but future efforts to project mass changes of marine-terminating glaciers should take them into account. For obtaining an increased accuracy, one could devise a calibration strategy that simultaneously constrains additional variables with observational estimates by, for instance, implementing a multi-objective optimization of both the flow parameters and the frontal ablation parameter. Variables that might be included in this are: ice thickness and velocity as well as (frontal) area changes. Such a procedure would likely improve the ability of the model to simulate glacier dynamics properly.

6. Conclusion

Since a large portion of the Northern Hemisphere's glacier mass is contained in marine-terminating glaciers, taking processes occurring at their fronts into account more rigorously is a step toward more robust glacier mass change projections. For the lowest and highest emission scenarios, we project glacier GMSLR contributions by 2100 to be 74 ± 26 and 131 ± 44 mm SLE, respectively. We find that in projections accounting for marine frontal processes, the GMSLR contribution at the end of this century is reduced by ~9%, and by up to ~13 mm SLE, compared to projections neglecting such processes. Though this might be a small impact regarding the whole Northern Hemisphere (excluding the Greenland ice sheet), the effect of this is an ~30% difference when only regarding marine-terminating glaciers or individual regions. Such substantial impacts on regional results would have important consequences for potential changes in regional ocean circulation and ecology caused by increased freshwater input due to glacier mass loss. For the Greenland periphery we find a converse effect; up to ~2% more total and 8% more marine-terminating glacier mass is lost when including frontal processes in the modeling chain. Another interesting finding is that the difference in glacier mass change projections between the emission scenarios is smaller in 2100 when applying our marine-terminating glacier framework compared to projections neglecting marine frontal processes.

Because numerical modeling of glacier dynamics requires boundary conditions and parameterizations, calibration is a crucial step. Here, we calibrated the parameters responsible for the removal of ice mass from the glaciers: the sensitivity to atmospheric temperatures, and the frontal ablation parameter. In addition, we tested the variance in our results caused by other relevant parameters: the flow parameters, and the assumed fraction of mass changes in the terminus area below the water level. We find that varying the flow parameters, especially Glen's A, has a considerable effect on the projected amount of frontal ablation and total glacier mass change. The variance in results induced by this is smaller than that caused by differences in the atmospheric forcing data, namely the GCM ensemble. Still, for future studies that focus on most accurately projecting global glacier mass change, or on regions with much ice mass contained in marine-terminating glaciers, it might be considered to simultaneously calibrate the flow parameters and the frontal ablation parameter by constraining additional variables (i.e. velocity, thickness, area change) using observational estimates. Ultimately, repeating the exercise laid out in this work with calibration data for the Southern Hemisphere (excluding the Antarctic ice sheet), and lake-terminating glaciers, will probably reveal even more strongly that robustly modeling water-terminating glaciers has a significant effect on glacier mass projections in comparison to not doing so. One further crucial future step will be to simulate the effect of changing ocean temperatures and subglacial discharge on frontal ablation. Finally, glacier mass change reconstructions could benefit from accounting for frontal ablation and mass changes below sea level as well, especially when used to reconstruct sea level changes.