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, Reference Huss and Hock2018; Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). 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, Reference Castro de la Guardia, Hu and Myers2015; Meire and others, Reference Meire2017). 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, Reference Straneo2013). 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, Reference Recinos, Maussion, Rothenpieler and Marzeion2019, Reference Recinos, Maussion, Noël, Möller and Marzeion2021). Since ice volume cannot be measured at large scales, models constrained by observations are often used (Farinotti and others, Reference Farinotti, Huss, Bauder, Funk and Truffer2009). 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, Reference Recinos, Maussion, Rothenpieler and Marzeion2019).
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, Reference Huss and Hock2015). 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, Reference McFadden, Howat, Joughin, Smith and Ahn2011). Additionally, Huss and Hock (Reference Huss and Hock2015) 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, Reference Maussion2019). 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, Reference Oerlemans and Nick2005), 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, Reference Pfeffer2014), 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, Reference Maussion2019). 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, Reference Pfeffer2014; RGI Consortium, 2017) is the basis of OGGM, similar to several other global glacier models (Marzeion and others, Reference Marzeion2020). 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 (Reference Werder, Huss, Paul, Dehecq and Farinotti2020). 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, Reference Eis, Maussion and Marzeion2019, Reference Eis, van der Laan, Maussion and Marzeion2021).
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, Reference Harris, Osborn, Jones and Lister2020) is used. The temperature is subsequently corrected using a globally fixed linear lapse rate (6.50°C km−1). 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:
where f p is a dimensionless precipitation factor (Giesen and Oerlemans, Reference Giesen and Oerlemans2012), $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−1) 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 (Reference Maussion2019) 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:
where A is the temperature-dependent ice creep parameter (here we use the default value of 2.4 × 10−24 s−1 Pa−3), 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:
where ρ is the ice density (here we use 900 kg m−3), g the gravitational acceleration (9.81 m s−2) and α the surface slope computed numerically along the flowline on a staggered grid:
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 (Reference Howat, Joughin, Tulaczyk and Gogineni2005). This additional force (per unit width) can be calculated as:
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−3) 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:
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:
Hence, the mean additional stress is:
and it is ensured that:
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, Reference Enderlin2016), 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 (Reference Enderlin2016), 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, Reference Robel2017), 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, Reference Nick, Vieli, Howat and Joughin2009; Enderlin and others, Reference Enderlin, Howat and Vieli2013). 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, Reference Benn, Hulton and Mottram2007). Previously, the sliding velocity in OGGM was calculated as:
where f s is a sliding parameter (default value 10−20 m2 s−1 Pa−3), based on Oerlemans (Reference Oerlemans1997). 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:
where h* is the height above buoyancy:
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 (Reference Oerlemans1997) (5.7 × 10−20 m2 s−1 Pa−3) 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−20 m2 s−1 Pa−3 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, Reference Benn, Hulton and Mottram2007; Stearns and Van der Veen, Reference Stearns and van der Veen2018; Zoet and Iverson, Reference Zoet and Iverson2020).
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:
where q is the ice mass flux per unit width (in m2 a−1), 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)):
and r h the inverse of the relative height above buoyancy:
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 (Reference Oerlemans and Nick2005):
where Q f is the frontal ablation flux (m3 a−1), k the water-depth sensitivity parameter (a−1; 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 (Reference Recinos, Maussion, Rothenpieler and Marzeion2019), 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, Reference Oerlemans and Nick2006; Brinkerhoff and others, Reference Brinkerhoff, Truffer and Aschwanden2017).
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, Reference Recinos, Maussion, Rothenpieler and Marzeion2019, Reference Recinos, Maussion, Noël, Möller and Marzeion2021). 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, Reference Vieli and Payne2005). 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:
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, Reference Copland and Mueller2017), 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 (Reference Bassis and Walker2011), 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 (Reference Hugonnet2021). 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 (Reference Hugonnet2021) 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 (Reference Kochtitzky2022) 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 (Reference Hugonnet2021) data are not able to capture it, which would lead to inconsistencies.
The dataset given by Kochtitzky and others (Reference Kochtitzky2022) 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, Reference Recinos, Maussion, Noël, Möller and Marzeion2021; Müller and others, Reference Müller, Friedl, Palmer and Marzeion2022). 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, Reference Copland and Mueller2017) and many floating tongues have collapsed in the last few decades (White and Copland, Reference White and Copland2019; Kochtitzky and Copland, Reference Kochtitzky and Copland2022).
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, Reference Eyring2016). 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, Reference Marzeion, Jarosch and Hofer2012, Radić and others, Reference Radić2014 or Huss and Hock, Reference Huss and Hock2015). 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:
where:
• ΔM awl is the observed annual mass balance above sea level of a glacier (mm3 a−1) as given by Hugonnet and others (Reference Hugonnet2021),
• C is the observed annual frontal ablation rate of a glacier as given by Kochtitzky and others (Reference Kochtitzky2022) (mm3 a−1),
• ΔM f is the observed annual rate of mass change due to area changes in the terminus region of a glacier (mm3 a−1; hereafter named retreat volume) as given by Kochtitzky and others (Reference Kochtitzky2022),
• 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 (mm2),
• T m is the annually accumulated temperature above the threshold for ice melt at the glacier surface (K a−1).
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, Reference Maussion2019). Because the data of Hugonnet and others (Reference Hugonnet2021) 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 (Reference Hugonnet2021) 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 (Reference Kochtitzky2022) 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 (Reference Hugonnet2021) 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 (Reference Kochtitzky2022) (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 (Reference Hugonnet2021), 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−3). 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 (Reference Cook2019) 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 (Reference Marzeion2020) 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 (Reference Marzeion2020). 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 (Reference Marzeion2020), 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 (Reference Marzeion2020).
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, Reference Goelzer2020).
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, Reference Recinos, Maussion, Rothenpieler and Marzeion2019) or a different calibration strategy (Recinos and others, Reference Recinos, Maussion, Noël, Möller and Marzeion2021), 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 (Reference Recinos, Maussion, Rothenpieler and Marzeion2019) 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 (Reference Recinos, Maussion, Noël, Möller and Marzeion2021). 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 (Reference Cook2019) 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, Reference Rabatel, Sanchez, Vincent and Six2018). For simplicity, we assume steady state at the RGI date (Maussion and others, Reference Maussion2019), 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 (Reference Millan, Mouginot, Rabatel and Morlighem2022). 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, Reference Marzeion2020); 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, Reference Slater, Nienow, Cowton, Goldberg and Sole2015). 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 (Reference Kochtitzky and Copland2022) 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, Reference Farinotti2019).
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, Reference Hourdin2017). 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.
Data
The data of the projections including frontal ablation can be accessed at Zenodo (https://doi.org/10.5281/zenodo.7550643). Further data inquiries can be made to the first author. The model code repository used for this work is archived at Zenodo (doi: 10.5281/zenodo.7547966) and its main part will be added to the OGGM core (github.com/OGGM/oggm) in the foreseeable future.
Acknowledgements
We thank the two anonymous reviewers for helping to improve the manuscript with detailed and thoughtful comments as well as the editor Jonathan Kingslake for handling the revision process. We also thank Beatriz Recinos for the groundwork in OGGM and fruitful discussions on the matters presented in this work. Additional thanks go to Timo Rothenpieler for his steady technical support.
Author contributions
J.-H.M. designed the research together with F.M., developed the changes to the model code, conducted the simulations and statistical evaluation, and wrote the manuscript. F.M., L.U. and B.M. contributed to the manuscript. W.K. and L.C. supplied and assisted with the data, and made contributions to the manuscript.
Financial support
This research has been supported by the Deutsche Forschungsgemeinschaft (DFG, grant nos. IRTG 1904/3 ArcTrain and MA6966/4-1). The article processing charges for this open-access publication were covered by the University of Bremen.
Appendix A
Table A1. GCMs used for projections and their availability for the emission scenarios considered in this work
Table B1. Parameter values tested for the estimation of parameter sensitivity
Appendix B
B.1 Parameter sensitivity test setup
The two flow parameters (A and f s in Eqns (2), (11) and (13)) play a significant role in our simulations, since they are changing the outcomes of the ice thickness inversion, and the behavior of ice dynamics. We conducted additional simulations to get a notion of these parameters’ influence on the results. For those experiments we calibrated the frontal ablation parameter (k) with different values of the flow parameters and then reran the projections with one GCM and one emission scenario (BCC-CSM2-MR and SSP3 7.0). We additionally varied the assumed fraction of retreat volume below sea level (f bwl) to get an impression of how strong the effect of this is. Table B1 shows the parameter space we investigated. The modeling chain was run once for every possible parameter combination; hence 27 times in total. Moreover, we conducted runs in which we halved the stress coupling length (see Section 2.2) applying the different values for the flow parameters given in Table B1 while keeping f bwl at 0.75.
B.2 Parameter sensitivity results
Here we discuss the uncertainty in our results due to the influence of the flow parameters (Glen's A, and the sliding parameter f s; see the previous section, and Eqn (13)). Though such uncertainties apply to land-terminating glaciers to some degree as well, marine-terminating glaciers can be especially sensitive to ice thickness changes, which are a compound result of ice thickness inversion, ice dynamics, frontal ablation and surface mass balance. This is because they can be prone to acceleration and rapid retreat when approaching flotation. Moreover, our frontal ablation parameterization is directly dependent on ice thickness, and the water depth estimated via the ice thickness inversion (see Eqn (16)). Also, with low values for these parameters, the ice thickness inversion will result in very thick glaciers and thus more ice that can potentially be removed and (if situated above the water level) add to GMSLR. Additionally, the flow parameters impact the simulation of ice dynamics and thus the geometric adjustment of glaciers to a certain perturbation in the forcing. In that sense, two dynamical feedbacks, with opposite signs, play a role: the elevation change of a glacier's terminus which is a negative feedback, and the change of a glaciers’ surface elevation which is a positive feedback.
All the simulation results presented above were computed with certain parameter values (standard values named default in Table B1). As explained in Section B.1, we therefore conducted additional experiments to get an impression of these parameters’ influence. In Figure 6 it is visible that varying Glen's A has the largest impact on the results concerning the difference between taking frontal processes into account and not doing so (relating σ A to σ tot). Only the GCM ensemble's variance ($\sigma _{\rm GCM}^{\rm ssp}$
) has an even stronger effect throughout the 21st century, though the decreasing variance in the high emission scenario during the second half of the displayed period is intriguing. Plausibly, varying the assumed fraction of mass changes in the terminus area below the water level has a linearly growing influence on the results ($\sigma _{f_{\rm bwl}}$
), while it is the overall weakest of the three investigated parameters. It can also be seen that there are positive covariances between the parameters, which is probably because varying the two flow parameters (A, f s) has a somewhat similar effect. Regarding temporally accumulated frontal ablation, varying Glen's A has an even stronger relative effect, while varying f bwl has nearly none (see Fig. 7). The large influence of A is mostly based on the fact that it has a strong influence on the initial ice thickness found by the inversion procedure (see Section 2.4) and this, in turn, strongly influences frontal ablation. Interestingly, the large variance in frontal ablation projections caused by varying A (σ A) shown in Figure 7 is not directly translated to the effect of including frontal processes in the projections (Fig. 6).
Fig. 6. Northern Hemisphere std dev. (1 − σ) of differences between projections accounting for marine frontal processes and those that do not, caused by varying different parameters. The yellow line displays covariances of the three tested parameters, and $\sigma _{\rm GCM}^{\rm ssp}$
the variance of the GCM ensemble in the respective scenario with default parameter set. The subscripts in the legend refer to the respective parameter varied (see Table B1), and σ tot to the ensemble of model runs’ total variance.
Fig. 7. Same as Figure 6, but for accumulated frontal ablation estimates. $\sigma _{\rm GCM}^{126}$
is excluded, due to its similarity with $\sigma _{\rm GCM}^{370}$
and thus better visibility.
There are two effects of varying A attenuating each other regarding glacier mass loss over the course of the century. At first, more frontal ablation takes place with a lower A value due to the larger initial (frontal) ice thicknesses. The high frontal ablation overcompensates the reduced sensitivity to atmospheric temperatures, resulting in more glacier mass loss compared to projections not including marine frontal processes. Over the course of the century, as atmospheric temperatures increase, the effect of reduced sensitivities to these takes hold and the projections including frontal ablation result in less mass loss. As the ice thickness inversion produces less ice volume with a high A value, less ice will be available for further melt and frontal ablation toward the end of the century, which means that the difference in ice melt cannot grow more, while it does so in the simulations with a low A value (see Fig. 9). Such competing effects of varying A probably explain the peculiar shape of the std dev. in Figure 6.
The influence of changing the sliding parameter (f s) on simulated frontal ablation is similar (see $\sigma _{f_{\rm s}}$
in Fig. 7), although applying no sliding (i.e. f s = 0) results in slightly less frontal ablation than applying the default value (f s = 10−20; see Fig. 10). That is despite including sliding results in thinner glaciers at the front and thus less mass to be removed by frontal ablation (see Eqn (16)). Though sliding also enhances dynamical thinning at the front, more mass is apparently removed from the front by the flotation criterion. Applying the high value for f s results in the least frontal ablation though, probably because of the (frontal) ice thicknesses resulting from the inversion procedure being much smaller.
Concerning the uncertainty in total accumulated contribution to GMSLR by Northern Hemisphere glaciers, it is visible in Figure 8 that varying Glen's A has again the strongest influence and produces a 1 − σ std dev. of close to 20 mm SLE. Varying the sliding parameter (f s) results in a std dev. of ~10 mm SLE, while varying f bwl has an effect of <1 mm SLE. The non-linearity of the variance induced by the flow parameters can most likely be explained by their influence on the simulation of the dynamical/geometric feedbacks mentioned above in combination with the initial glacier volumes and frontal ablation. Figure 11 shows that simulations with lower values for the flow parameters generally result in higher GMSLR contribution estimates, which presumably is mostly due to the higher amounts of initial ice volume available for melt over the remainder of the century. Moreover, we find that the variance caused by the differences in the GCM ensemble ($\sigma _{\rm GCM}^{\rm ssp}$
) is larger than that caused by varying any of the parameters considered here.
Fig. 8. Same as Figure 6, but for estimated accumulated GMSLR contribution.
Fig. 9. Estimated differences between projections accounting for frontal processes and those that do not, averaged over different parameter values. Panels (a–c) display averages over the nine possible different parameter sets with the respective parameter value corresponding to the three different columns in Table B1. Panel (d) displays the results of all 27 parameter sets possible in Table B1 in gray shading with the mean/median as well as the result for the parameter set chosen for the computation of results in Section 4 (default in Table B1) highlighted.
Fig. 10. Same as Figure 9, but for accumulated frontal ablation estimates.
Fig. 11. Same as Figure 9, but for accumulated GMSLR contribution.
We additionally conducted runs in which we held f bwl constant at 0.75, but halved the stress coupling length (L F = 4 km in Eqns (7–9)) and applied the different parameter combinations for the flow parameters as in the experiments described above. The resulting accumulated frontal ablation estimates are 13% higher in 2100 on average, which results in slightly higher average GMSLR contribution estimates (0.7 mm SLE) compared to results applying the higher value for L F. The increase in frontal ablation produced by lowering the stress coupling length makes sense, since the additional force due to the hydrostatic pressure imbalance at the glacier front (see Eqn (5)) is distributed over a shorter distance. This causes a stronger thinning close to the front, which in turn causes more ice removal by the flotation criterion.
The analysis presented in this subsection demonstrates that further work on finding appropriate parameter values for individual glaciers is required in order to improve the accuracy of mass change projections for (marine-terminating) glaciers. We would like to emphasize that the default parameter set we used for computing our main results in Section 4 was not chosen for its suitability for all glaciers worldwide, but for its applicability to most glaciers without causing numerical errors in the dynamical model.
Figures 9–11 show that results obtained with the parameter set chosen for the simulations presented in the sections above are reasonably close to the parameter ensemble's mean and median though. In Figures 6–8 the GCM ensemble's std dev. with the default parameter set is given as well. This is the largest source of uncertainty in accumulated mass changes, though it differs between emission scenarios. The difference between GCMs is relatively small concerning accumulated frontal ablation in the whole Northern Hemisphere. This hints that the Northern Hemisphere's accumulated frontal ablation is somewhat independent of the atmospheric forcing in our simulations, although annual frontal ablation rates may differ between different GCMs.
Open access
