1. Introduction
The Greenland Ice Sheet (GrIS) is surrounded by peripheral glaciers (PGs) and ice caps that have different levels of connectivity with the ice sheet: they are either dynamically connected, entirely detached or separated from the ice sheet by well-defined ice divides in their accumulation zone (Rastner and others, Reference Rastner2012; Bjørk and others, Reference Bjørk2018; Liu and others, Reference Liu, Enderlin, Marshall and Khalil2022). All PGs in Greenland cover approximately 12% of the world's glaciated area (excluding the ice sheets) and just 5% of Greenland's ice area, but have lost mass at a rate of ~20 Gt a−1 between 2000 and 2019 (Hugonnet and others, Reference Hugonnet2021) thus playing an important role in Greenland's contributions to sea-level rise. The loss of mass from tidewater glaciers to the ocean through frontal ablation (i.e. calving, subaerial and subaqueous frontal melting) is a major component of the mass budget of the GrIS (Cowton and others, Reference Cowton, Sole, Nienow, Slater and Christoffersen2018a; King and others, Reference King2018; Mankoff and others, Reference Mankoff2019; King and others, Reference King2020) estimated to be nearly 487 ± 49 Gt a−1 (Mankoff and others, Reference Mankoff2019). Therefore, frontal ablation is likely also a major component of the mass budget for most PGs (Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021). The latest global glacier ice thickness inventory of PGs that are dynamically de-coupled from the ice sheet, implies that they can potentially contribute up to 26.8 ± 9.5 mm sea-level equivalent when integrating only the ice volume above flotation (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022). Here we focus only on tidewater glaciers de-coupled from the ice sheet, since their solid ice discharge still remains poorly constrained. We only model those glaciers with a connectivity level of 0 and 1 in the Randolph Glacier Inventory (RGI v6.0, Rastner and others, Reference Rastner2012, i.e. 763 glaciers), which is common practice for mountain glacier change assessments (Hock and Huss, Reference Hock, Huss and Letcher2021).
Constraining and initializing glacier models of tidewater PGs remains a challenging task, due to sparse in situ ice thickness measurements causing unknown ice thicknesses at the calving front. Even with complete data coverage, ice-dynamical models cannot simply use observation-based ice thicknesses or frontal ablation fluxes to initialize model runs: these often heterogeneous data sources need to be assimilated into models in order to ensure physical consistency with the model parameterizations and prevent non-physical shocks when the model is run forward. These types of shocks are caused by the uncertainties that arise from the discretization of continuous data to a numerical environment, data uncertainties and discrepancy between present-day modeled glacier state and observations (see Zekollari and others, Reference Zekollari, Huss, Farinotti and Lhermitte2022, for recommendations on model initialization).
There are few estimates of dynamic mass loss for tidewater PGs (Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021; Bollen and others, Reference Bollen, Enderlin and Muhlheim2022; Kochtitzky and others, Reference Kochtitzky2022), but those still suffer from considerable uncertainty associated with the method selected to estimate the ice thickness at the calving front (e.g. Bollen and others, Reference Bollen, Enderlin and Muhlheim2022), significant uncertainty in the terminus geometry and errors in the true value of the frontal ablation flux (e.g. Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021).
Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) rely on an empirical scaling function to estimate glacier thickness from speed observations across a terminus flux gate. Making use of the relationship between glacier geometry and speed (Enderlin and others, Reference Enderlin2014), Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) construct a non-linear fit between ice thickness observations and satellite surface velocities. From that fit they predict an ice thickness for those tidewater PGs without thickness observations and then calculate a solid ice discharge time series for 585 PGs. However, Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) solid ice discharge estimates do not include mass changes due to terminus retreat or advance during the calculated period (see Enderlin and others, Reference Enderlin2014; Bollen and others, Reference Bollen, Enderlin and Muhlheim2022, for more information). Kochtitzky and others (Reference Kochtitzky2022) address this issue and estimate frontal ablation indirectly, by calculating ice discharge and terminus mass change, while accounting for the climatic mass balance (MB) (i.e. the mass changes due to snow accumulation, surface melt and refreezing) below the flux gate. Yet Kochtitzky and others (Reference Kochtitzky2022) estimates are limited to 537 glaciers and there is no observation prior to 2000.
Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) calibrate the frontal ablation parameterization implemented in the Open Global Glacier Model (OGGM) with three model-independent datasets and two different methods. The first method calibrates the parameterization using surface velocity fields derived from satellite observations. For this method they use and compare surface velocities from the MEaSUREs Multi-year Greenland Ice Sheet Velocity Mosaic (MEaSUREs v1.0 Joughin and others, Reference Joughin, Smith, Howat and Scambos2016) and from the ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities Mosaic (Gardner and others, Reference Gardner, Fahnestock and Scambos2019, ITS_LIVE). The second method uses surface mass balance (SMB) from the monthly output of the polar Regional Atmospheric Climate Model (RACMOv2.3p2; Noël and others, Reference Noël, van de Berg, Lhermitte and van den Broeke2019), statistically downscaled to 1 km resolution following Noël and others (Reference Noël2016). Both calibration methods make equilibrium assumptions for the ice dynamics and MB processes, assuming that the average amount of ice that passes through the glacier terminus in an average MB year must be equal to the amount of ice delivered to the terminus by ice flow if the glacier geometry remains unchanged (see Section 4 in Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021, for details). By assuming such equilibrium conditions, they estimate parameter values (a calving constant of proportionality k in the parameterization) for each individual tidewater PG, compute frontal ablation fluxes and improve the first ice thickness estimation computed by OGGM.
However, the frontal ablation fluxes computed in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) are only an average flux and the ice thickness/volume obtained was not evaluated with independent observations (i.e. data that were not used for calibration). All three studies mentioned above (Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021; Bollen and others, Reference Bollen, Enderlin and Muhlheim2022; Kochtitzky and others, Reference Kochtitzky2022), hence represent independent assessments of frontal ablation for tidewater PGs that can be compared and evaluated using new published region-wide data. Here we make use of the ice thickness estimates from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) to evaluate the ice thickness estimates obtained by OGGM after calibrating the frontal ablation parameterization with the methods described in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021). Additionally, we use flux-gate lengths derived from Kochtitzky and Copland (Reference Kochtitzky and Copland2022) to better constrain the glacier geometry in OGGM.
In this letter, we revise the Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) frontal ablation fluxes and ice thickness estimates in light of the recent publications from Kochtitzky and Copland (Reference Kochtitzky and Copland2022), Kochtitzky and others (Reference Kochtitzky2022), Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) and Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022). Our goal is to improve initialization methods for large-scale glacier models. We identify future research priorities in modeling tidewater PGs, which ultimately may help to improve estimates of future sea-level rise from Arctic glaciers.
2. Data Input
We use the same data input as in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) with a few differences specified in this section.
2.1. Calving front geometry
To improve the glacier geometry at the calving front, we correct the terminus width computed by OGGM. The widths in the model are obtained from a geometric first guess by intersecting normal's from the main flowline at each gridpoint, with the glacier outlines and the tributaries’ catchment areas (for more details see Maussion and others, Reference Maussion2019). However, when the topography has a poor resolution (Fig. 1b) this model task may lead to uncertainties in the frontal ablation flux (Recinos and others, Reference Recinos, Maussion, Rothenpieler and Marzeion2019). To ensure a good representation of the calving front width, we make use of flux-gate lengths derived from Kochtitzky and Copland (Reference Kochtitzky and Copland2022) to correct the calving front width computed by OGGM following the method described in Recinos and others (Reference Recinos, Maussion, Rothenpieler and Marzeion2019).
Fig. 1. Calving front overview for the glacier ID: RGI60-05.05558: (a) Randolph Glacier Inventory outline (RGIv6.2, black solid line) and overview of the glacier area from Sentinel-2 multi-spectral satellite imagery from 2015 (Moon and others, Reference Moon, Fisher, Harden, Simonoko and Stafford2022), (b) OGGM catchment widths and flowline estimation and (c) ITSLIVE surface velocities (Gardner and others, Reference Gardner, Fahnestock and Scambos2019) and OGGM flowlines (red solid line).
A poor estimate of the glacier terminus geometry however, is not the only issue shown in Figure 1: there is a big disconnect between the OGGM flowlines and the ice velocity direction which would likely result in a wrong thickness estimation at that location. Any glacier model using an automatic method to compute center lines (and relying on digital elevation models i.e. Maussion and others, Reference Maussion2019), has the challenge of simulating tidewater glaciers that have more than one terminus. Kochtitzky and Copland (Reference Kochtitzky and Copland2022) have mapped the terminus position for every marine-terminating glacier in the Northern Hemisphere for 2000, 2010 and 2020, including PGs. Such data can be used to identify glaciers with more than one terminus and either manually (or automatically) correct flowlines, and/or subdivide and correct the RGI outlines that need to be separate entities. Here we make use of the Kochtitzky and Copland (Reference Kochtitzky and Copland2022) database to identify glaciers with more than one terminus and remove them from our workflow, as such glaciers cannot yet be modeled by OGGM. Tidewater PGs with more than one terminus represent almost 6% (1909.28 km2) of the glacierized area of interest and can introduce large uncertainties to any regional estimate (see Section 4.1).
2.2. Ice thickness and volume used for evaluation
To evaluate ice thickness and volume estimates obtained from OGGM, we compare them with the results from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022). The ice thickness in Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) is computed on the basis of surface motion and slopes using the shallow-ice approximation (SIA). For the Greenland periphery, they use surface velocities from the MEaSUREs Multi-year Greenland Ice Sheet Velocity Mosaic, v2.0 (Joughin and others, Reference Joughin, Smith, Howat and Scambos2016; Mouginot and others, Reference Mouginot, Rignot, Scheuchl and Millan2017) and slopes were derived from multiple sources of DEMs. The ice-dynamic parameters in the SIA equations are calibrated using ice thickness measurements from the Glacier Thickness Database (Welty and others, Reference Welty2020). The estimates by Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) are independent of the SMB; therefore, unlike most glacier models based on mass conservation (as listed in Farinotti and others, Reference Farinotti2017), there is no need to reconcile upstream mass flux and ice thickness with a frontal ablation estimate. However, their method is strongly dependent on the surface slope, surface velocity and the SIA, and is therefore most uncertain for fast tidewater glaciers, where the surface slope is nearly flat, surface velocities are high and basal processes become increasingly important.
Despite these limitations, ice thickness values at the calving front produced by Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) are likely the truest ice thickness estimate available for tidewater PGs at this time (Recinos and others, Reference Recinos, Maussion, Rothenpieler and Marzeion2019; Welty and others, Reference Welty2020). These modeled ice thickness values are an improvement over OGGM, which shares these limitations in addition to further limitations attributable to the flowline assumption. To compare Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) data with OGGM, we re-project their ice thickness (and uncertainty) products to a local glacier map and into the main glacier flowline (see Fig. 2).
Fig. 2. (a) Ice thickness from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) re-projected to the glacier grid for the glacier with RGI ID 60-05.09915. Note the data gaps over some parts of the glacier. (b) Glacier main centerline profile; comparison between the estimated bed map from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) (red line, the pink shading shows the estimated uncertainty) and the model ice thickness obtained by calibrating OGGM's calving parameterization with two methods: (1) velocity constraint (green) and (2) RACMO SMB constraint (orange). The black solid line represents the glacier surface elevation.
2.3. Frontal ablation fluxes used for evaluation
We evaluate OGGM frontal ablation fluxes with the solid ice discharge from 585 tidewater PGs. Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) computed time series (1985–2018) of solid ice discharge by using a flux-gate method (Enderlin and others, Reference Enderlin2014; King and others, Reference King2018; Mankoff and others, Reference Mankoff2019), remotely sensed glacier speed time series and an empirically estimated ice thickness. However, such empirical estimates of ice thickness suffer from large uncertainty, and the solid ice discharge needed to be adjusted to account for temporal variations in thickness at the flux gate and SMB between the flux gate and terminus. Despite having high uncertainties, Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) represent the only region-wide estimate of frontal ablation fluxes for PGs previous to 2000 and accounts for three different time periods. This allows us to test equilibrium assumptions used in OGGM to compute frontal ablation fluxes (see next section or Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021). We compare total and individual model estimates of frontal ablation to the time-averaged annual glacier discharge from Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) averaged over three different time periods (see Section 4.1). Additionally, we also compare regional frontal ablation fluxes to those from Kochtitzky and others (Reference Kochtitzky2022). A glacier-to-glacier comparison between Kochtitzky and others (Reference Kochtitzky2022) and this study is problematic and could be biased, given that we correct OGGM terminus widths with flux-gate lengths used in Kochtitzky and others (Reference Kochtitzky2022) frontal ablation estimates.
3. Methods
3.1. OGGM: model set up and calibration
The framework of the OGGM (v1.5.3, Maussion and others, Reference Maussion2019) and the frontal ablation parameterization have been explained in detail by Recinos and others (Reference Recinos, Maussion, Rothenpieler and Marzeion2019, Reference Recinos, Maussion, Noël, Müller and Marzeion2021). Here we use the same methods described in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) to calibrate the calving constant of proportionality k in OGGM calving law, simulate frontal ablation fluxes, ice thickness and ice volume of tidewater PGs. The only difference in the OGGM configuration for this study is that we fix the climatological period for the OGGM MB model to the period between 1961 and 1990, in line with the analyses of GrIS mass loss (Rignot and Kanagaratnam, Reference Rignot and Kanagaratnam2006; Fettweis and others, Reference Fettweis2017) and PGs SMB change (Noël and others, Reference Noël2016, Reference Noël, van de Berg, Lhermitte and van den Broeke2019). By making this change and those explained in Section 2.1, we are able to simulate between 75 and 88% of the glacierized area of interest (depending on the data input used for the calibration).
Approximately 0.5% of the glacierized area has data gaps in RACMO2.3p2 (178 small glaciers), $\sim$
0.9% (51 glaciers) in MEaSUREs v1.0 and $\sim$
9% (45 glaciers) in the ITS_LIVE dataset. Another ~6% (44 glaciers) of the study area present errors in the preprocessing stages of OGGM including those glaciers identified with more than one terminus. Approximately 4% (174 glaciers) of the remaining area do not calve in OGGM under any k value. Non-calving glaciers and glaciers with a zero frontal ablation flux estimated by the model have low velocities (<10 m a−1) at the terminus according to ITS_LIVE data. We assume a calving flux of zero for such glaciers. Additionally, $\sim$
13% (126 glaciers) of the study area has a RACMO SMB mean ($\overline{\it{m}}_{(1961-1990)}$
) that is below zero, for such glaciers the steady-state assumption no longer holds, and we are not able to estimate a frontal ablation flux via the RACMO SMB method. For such glaciers (and only for the RACMO SMB method) we also specify a calving constant of proportionality and a frontal ablation flux equal to zero.
4. Results
The total frontal ablation flux is significantly less than that in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021), for both calibration methods (e.g. from 7.38 ± 3.45 decreased to 3.5 ± 1.5 Gt a−1 for the velocity method). This is mainly due to the removal of glaciers with more than one terminus, which consequently introduce large uncertainties in the total frontal ablation flux (see in Section 5 our proposed solution to this problem). However, in general we increase the number of glaciers for which we can find frontal ablation fluxes compared to our previous study. Despite this difference, the main findings in Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021) remain the same: velocity constrained frontal ablation fluxes are significantly higher than those found if the calving parameterization is constrained using RACMO-derived frontal ablation fluxes (see Table 1). In the next sections we compare our frontal ablation fluxes, ice thickness and volume estimates against three region-wide studies (Bollen and others, Reference Bollen, Enderlin and Muhlheim2022; Kochtitzky and others, Reference Kochtitzky2022; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022).
Table 1. Comparison of previous and current frontal ablation fluxes computed by OGGM for the different calibration methods and those from Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) and Kochtitzky and others (Reference Kochtitzky2022) in Gt a−1
In parentheses we note the fraction of the glacierized area of interest (tidewater PGs area = 32 202.54 km2) covered by each method and each study.
4.1. Modeled frontal ablation fluxes evaluation
We compare total frontal ablation fluxes in Table 1 to those estimated by Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021), Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) and Kochtitzky and others (Reference Kochtitzky2022) for different time periods. A glacier-to-glacier comparison among all datasets is challenging, as not all the glaciers have an ice discharge estimate in the earlier period or have a model frontal ablation flux from OGGM (this depends on the calibration method used and data gaps on each data input used for calibration). For each calibration method, we look for the glaciers that have an ice discharge in Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) and a frontal ablation flux estimated by OGGM (on average we can only compare ~280 glaciers or 46% of the glacierized area of interest).
There is a moderate correlation (e.g. r 2 = 0.50 with a 95% confidence interval (CI) [0.42, 0.61], bias = 0.0009 km3 a−1 for ITS_LIVE) between OGGM velocity constrained frontal ablation fluxes and the time-averaged annual glacier discharge from Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) among the median of all time periods (see Figs 3a, b). However, there is a clear disagreement with those estimated using RACMO-derived frontal ablation fluxes under the assumption of a balanced mass budget (see Fig. 3c and Table 1). OGGM is not able to predict individual tidewater glacier dynamics if it relies only on RACMO SMB estimates and the assumption of a closed budget to constrain the calving parameterization. Velocity observations are essential to constrain the model and estimate the dynamic mass loss of tidewater PGs. Frontal ablation fluxes derived using velocity observations to constrain the calving parameterization have a stronger agreement with the median solid ice discharge over the period of 1985–98 (see orange line in Figs 3a, b). Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) refer to this period as a ‘steady-state’ period in line with the analyses of GrIS mass loss and peripheral glacier SMB change (e.g. van den Broeke and others, Reference van den Broeke2016; Noël and others, Reference Noël2017). Figures 3a, b show that OGGM is able to reproduce solid ice discharge for tidewater PGs in an equilibrium setting (time-averaged annual ice discharge over 1960–90), an expected result given our model's ‘steady-state’ assumption to estimate frontal ablation fluxes.
Fig. 3. Frontal ablation fluxes computed by OGGM via the different calibration methods compared with solid ice discharge estimates from Bollen and others (Reference Bollen, Enderlin and Muhlheim2022) for three different periods: long-term median 1985–2018 (blue circles), 1985–98 median (orange stars) and 1999–2018 median (green squares): (a) ITSLIVE, (b) MEaSUREs and (c) RACMO-derived frontal ablation fluxes. Regression lines (solid lines) and statistics are shown in the upper part, we only show statistics for the first period (1985–98), i.e. percent of study area represented in the graph, regression slope, intercept, coefficient of determination (r 2), RMSD and bias. P-values are all smaller than 0.05. Gray solid lines represent slopes equal to 1 and intercepts equal to zero. Error bars are plotted in light gray and they represent 95% CI.
The true frontal ablation flux of each glacier is probably best represented by the velocity calibration method. Satellite velocity measurements provide an insight into ice dynamics and reflect the ocean influence over the movement of the ice at the calving front, whereas the RACMO calibration method is based on climatic MB data and the assumption that this climatic MB is exactly stabilized by frontal ablation. This assumption appears not to be true for most of the PGs, where fluxes derived with the velocity method are significantly larger than those fluxes derived with the RACMO method (see Results section in Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021).
4.2. Modeled ice thickness and ice volume evaluation
We rely on the results from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) to evaluate the estimated ice thickness and volume computed by OGGM via the different calibration methods. Because of the lack of thickness data from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) at the termini of many glaciers, we can only compare between 45 and 53% of the study area. To increase the number of glaciers with thickness data, we average both ice thickness estimates, OGGM and Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022), over the lowest five pixels of each glacier flowline and compare those averages in Figures 4a–c. OGGM ice thickness derived with the velocity calibration method, especially those constrained by ITS_LIVE velocities, correlates well with Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) ice thickness (r 2 = 0.5 with a 95% CI [0.4, 0.57], bias = 28.13 m).
Fig. 4. Model performance. Comparison between model ice thickness (a–c) and ice volume (d–f) computed via the different calibration methods and those from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022). The ice thickness is an average of the thickness at the last five pixels of the calving front. Regression lines (solid lines) and statistics are shown in the upper right corner, i.e. percent of study area represented in the graph, regression slope, intercept, coefficient of determination (r 2), RMSD and bias. P-values are all smaller than 0.05. Gray solid lines represent slopes equal to 1 and intercepts equal to zero and in all scatter plots uncertainty bars are plotted in light gray and they represent 95% CI.
In all three methods and for the majority of the glaciers, the model does not reproduce the ice thickness at the calving front estimated by Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) (RMSD >100 m for all three calibration methods). This is likely due to either OGGM or Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) misrepresents the real ice thickness or OGGM underestimates the frontal ablation flux (see Fig. 3) and is not able to reproduce a thick glacier at the terminus (see orange line in Fig. 2b). Additionally, without calibrating the rest of the ice-dynamic parameters in OGGM (ice stiffness and sliding parameter), the thickness upstream of the terminus is likely overestimated. That could explain the different results regarding the ice volume shown in Figures 4d–f; where the model clearly overestimates the volume of tidewater PGs by ~50% (RMSD slopes 1.54–1.44 for all methods). Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) depend on satellite coverage to estimate the ice thickness; therefore, data gaps in the satellite data mean that the exact volume comparison with a model is problematic, as shown in Figure 2a; parts of the glacier in Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) dataset could be ice-free. Our results highlight the importance of in situ ice thickness observations at the calving front, in order to truly evaluate ice thickness and frontal ablation estimates for tidewater PGs.
5. Future research directions
5.1. Simulating glaciers with more than one terminus
Uncertainties in the calving front geometry illustrated in Figure 1 originate from the DEM-dependent automatized glacier centerline identification method implemented in OGGM (Kienholz and others, Reference Kienholz, Rich, Arendt and Hock2014; Maussion and others, Reference Maussion2019), which misrepresents real glacier flowlines for PGs with more than one terminus. Glaciers with such errors cover $\sim$
6% of our study area and cannot be modeled by OGGM. Outside of Greenland, the glaciated area covered by glaciers with more than one terminus could be more. We propose the use of the terminus positions mapped by Kochtitzky and Copland (Reference Kochtitzky and Copland2022) to identify glaciers with more than one terminus and correct this problem by manually subdividing the RGI outlines so there is only one terminus per RGI entity, or by testing these glaciers with new automatic centerline detection methods that are independent of topographic data. Ultee and Bassis (Reference Ultee and Bassis2020) have successfully applied a new centerline detection algorithm based on tracing ice surface velocities from terminus positions. This approach could be implemented in OGGM, by making use of Kochtitzky and Copland (Reference Kochtitzky and Copland2022) terminus positions and ITS_LIVE surface velocities.
5.2. Implementation of ocean forcing in the model
Tidewater glaciers calve into Greenlandic fjords, a process that can export vast amounts of fresh water and drive an overturning circulation in the fjord, strongly influencing the interactions between the glaciers and the ocean (Davison and others, Reference Davison, Cowton, Cottier and Sole2020). Outlet glaciers of the GrIS interact with the ocean heat, which melts the ice at the edge of the glacier (Cowton and others, Reference Cowton, Sole, Nienow, Slater and Christoffersen2018b; Bevan and others, Reference Bevan, Luckman, Benn, Cowton and Todd2019), changing the glacier's frontal ablation rates and causing retreat and acceleration, thus enhancing the glacier's sea-level contribution (Benn and others, Reference Benn2017; Cowton and others, Reference Cowton, Todd and Benn2019; Slater and others, Reference Slater2022). Therefore, it is reasonable to assume that such interaction also occurs at tidewater PGs. While no global glacier model accounts for oceanic forcing of PG glaciers, it could be an important component of their dynamics. Neglecting these processes could lead to the underestimation of the solid ice discharge computed by the model (as shown in Fig. 3). The implementation of submarine melting is key to represent the dynamical processes at the calving front and necessary to asses sea level rise contributions of tidewater PGs.
5.3. Implementation of basal processes in large-scale glacier models
Modeled glacier dynamics depend critically on the choice of sliding law (Benn and others, Reference Benn, Warren and Mottram2007). Therefore, more complex sliding laws need to be implemented in OGGM and compared to assess the sensitivity of the model to basal processes, and how these could impact regional estimates of ice volume. The current OGGM sliding parameterization is outdated and does not account for effective pressure, which is the difference between ice overburden pressure and the basal water pressure. For a tidewater glacier where basal water flows toward the terminus, a minimal value for the basal water pressure is determined by the depth of the bed below sea level (Benn and others, Reference Benn, Warren and Mottram2007). By using a sliding law which relates the sliding coefficient, and thus the basal velocity, to effective basal pressure, the model could improve the representation of observed ice velocities behind the calving front (Vieli and others, Reference Vieli, Funk and Blatter2000; Downs and Johnson, Reference Downs and Johnson2022).
5.4. Improved observations for tidewater Peripheral Glaciers at the calving front
The NASA Operation Ice Bridge (OIB) mission conducted annual airborne-radar surveys across Greenland from 2009 to 2018 but focused primarily on data acquisition for the ice sheet and outlet glaciers. There are only a handful of ice thickness observations within a few kilometers of PG termini (Bollen and others, Reference Bollen, Enderlin and Muhlheim2022). The lack of ice thickness observations at the calving front hinders a true evaluation of the recent region-wide studies presented here and prevents a full estimate of the sea-level budget for Greenland, with tidewater PGs sea-level rise contribution still uncertain. As with the GrIS, to determine the future glacier stability and mass flux from PGs contributing to sea-level rise, high-resolution coastal ocean bathymetry datasets are required together with near-coastal observations of hydrography and circulation (Schaffer and others, Reference Schaffer2020). Finally, satellite ice velocity products are key to constrain and calibrate models but uncertainty estimates from such data products might not represent the true error magnitude of the velocity measurement (e.g. unrealistically small uncertainties in the ITS_LIVE data, see Gardner and others, Reference Gardner, Fahnestock and Scambos2019, and the product documentation for more information). Underestimates in observational errors give a false confidence in model parameters (Recinos and others, Reference Recinos, Maussion, Noël, Müller and Marzeion2021) and limit model error propagation and quantification.
6. Conclusions
Via a simple calibration of the calving parameterization, OGGM is able to reproduce solid ice discharge for tidewater PGs in an equilibrium setting and to improve the first bedrock topography estimated by the model in this region, which is key to decrease uncertainties in sea-level rise projections for tidewater PGs. However, large uncertainties in the model results remain due to (1) the misrepresentation of the true error in the input used for calibration (e.g. error magnitude from satellite velocities), (2) the non-calibration of ice-dynamic parameters (sliding and ice stiffness) upstream of the flowline and (3) the absence of oceanic forcing in the model. Correlations between OGGM's findings and new region-wide estimates of ice thickness and frontal ablation are limited by the assumption that such region-wide estimates can be representative of the dynamics of each tidewater PGs. The lack of in situ observations in this region makes the validation of model and remote-sensing-derived estimates problematic and an unsolved issue for PGs. Nevertheless, by improving the initial bedrock topography in models (i.e. OGGM) which are based on mass conservation, we decrease uncertainties in sea-level rise projections for tidewater PGs. A good representation of the frontal ablation flux and ice thickness distribution is key to accurately compute the sea-level rise contribution from tidewater PGs.
Acknowledgements
We thank Brice Noël for sharing his downscaled RACMO data which was originally provided for Recinos and others (Reference Recinos, Maussion, Noël, Müller and Marzeion2021), Will Kochtitzky for sharing the flux-gate length and frontal ablation data and Ellyn M. Enderlin for sharing the solid ice discharge time series. BR is supported by NERC standard grant NE/T001607/1 (QUoRUM). The computations were realized on the computing facilities of the University of Edinburgh, Scotland, UK.
Data availability
The OGGM software together with the frontal ablation parameterization module are coded in the Python language and licensed under the BSD-3-Clause license. The latest version of the OGGM code is available on Github (https://github.com/OGGM/oggm), the documentation is hosted on ReadTheDocs (https://docs.oggm.org/) and the project webpage for communication and dissemination can be found at http://oggm.org. The code and data used to generate all figures and analyses of this paper are available in a permanent DOI repository (https://doi.org/10.5281/zenodo.7143947) as well as the OGGM version used for this study (https://doi.org/10.5281/zenodo.6408559).
Open access
