SERMeQ Model Produces a Realistic Upper Bound on Calving Retreat for 155 Greenland Outlet Glaciers

The rate of land ice loss due to iceberg calving is a key source of variability among model projections of the 21st century sea level rise. It is especially challenging to account for mass loss due to iceberg calving in Greenland, where ice drains to the ocean through hundreds of outlet glaciers, many smaller than typical model grid scale. Here, we apply a numerically efficient network flowline model (SERMeQ) forced by surface mass balance to simulate an upper bound on decadal calving retreat of 155 grounded outlet glaciers of the Greenland Ice Sheet—resolving five times as many outlets as was previously possible. We show that the upper bound holds for 91% of glaciers examined and that simulated changes in terminus position correlate with observed changes. SERMeQ can provide a physically consistent constraint on forward projections of the dynamic mass loss from the Greenland Ice Sheet associated with different climate projections.


Introduction
The Greenland Ice Sheet is currently the largest single contributor to global mean sea level rise (van den Broeke et al., 2017). It discharges ice mass to the ocean through three main processes: release of surface meltwater, submarine melting where ice is in contact with the ocean, and the detachment (calving) of icebergs. The ice mass lost to submarine melting has only recently been directly observed (Sutherland et al., 2019) and remains difficult to estimate for the whole ice sheet (Beckmann et al., 2018), but it is clear that enhanced surface melting and calving processes have resulted in increased mass discharge since the late 1990s (Enderlin et al., 2014;Khan et al., 2014;van den Broeke et al., 2016).
Processes that control surface melt are increasingly resolved in regional models (Mottram et al., 2017;Noël et al., 2018). Iceberg calving, by contrast, remains poorly understood, with multiple contradictory parameterizations incorporated into ice sheet/glacier models Levermann et al., 2012;Morlighem et al., 2016). Furthermore, iceberg calving can remove mass more rapidly than is possible through melting alone, contributing to rapid tidewater glacier retreat through mechanisms like tidewater glacier instability Meier and Post (1987) and the recently described Marine Ice Cliff Instability (Bassis & Walker, 2012;Pollard et al., 2015).
Simulating discharge from the Greenland Ice Sheet is further complicated by the local factors affecting ice discharge at the nearly 200 outlet glaciers that connect the ice sheet to the ocean (e.g., Catania et al., 2018;Enderlin et al., 2018). For all but the largest outlets, iceberg calving occurs at smaller scales than are captured in continental-scale ice sheet models. Existing estimates of dynamic mass loss from Greenland outlets have come from extrapolating perturbations on the largest outlets Price et al., 2011), simulating the sea level contribution from only selected outlets (Choi et al., 2017;Morlighem et al., 2019) or simulating the entire ice sheet at a spatial resolution of 500 m (Aschwanden et al., 2016(Aschwanden et al., , 2019 to resolve about 30 of the nearly 200 glaciers that drain the Greenland Ice Sheet. Despite these achievements, more than 100 outlet glaciers, responsible for approximately one third of current Greenland Ice Sheet discharge (Enderlin et al., 2014), are not routinely simulated, and their dynamics cannot necessarily be inferred from the dynamics of larger outlets. Another layer of spatial complexity arises in that many outlet glaciers collect ice from several interacting tributary branches that are themselves also smaller than typical ice sheet model grid scale. The small scale of tributary glacier networks feeding outlets makes them especially challenging to simulate in continental ice sheet models, requiring model resolution of hundreds to tens of meters to adequately resolve.
A more fundamental challenge in projecting mass loss due to calving is the incompatibility of fracture-driven iceberg calving with the assumption of continuum deformation inherent in most ice sheet models (e.g., Greve, 2000;Price et al., 2015;Winkelmann et al., 2011). Simple empirical parameterizations can relate calving rate to continuous variables, such as proglacial water depth (Brown et al., 1982;Hanson & Hooke, 2000), but may not hold into the future as climate forcing enters a new statistical regime. Physically based calving laws, such as the fracture field approach developed by Albrecht and Levermann (2012) or von Mises calving law developed for Greenland by Morlighem et al. (2016), often impose an empirically adjustable calving rate parameter. Recent work has sought to simulate ice failure using continuum damage mechanics, with some success in a variety of case studies (Borstad et al., 2012;Duddu et al., 2013;Krug et al., 2014;Mercenier et al., 2019;Sun et al., 2017). However, at present the evolution of the damage field through a damage production function is also empirical, with multiple tuned parameters that are poorly constrained by laboratory or field measurements (Emetc et al., 2018). Another recent approach couples a granular model that allows true fracture and calving to a finite-element model that solves an approximation to the Stokes equations for viscous deformation, offering a very promising basis for process-scale simulation of fully dynamic calving . Unfortunately, the coupled approach remains too computationally expensive for century-scale projections. Despite their promise, neither continuum damage models nor granular calving models have been able to reproduce observed multiannual evolution of calving front positions in Greenland.
Improving projections of 21st century sea level rise requires models that can (i) reproduce observed patterns of glacier advance and retreat that vary among Greenland outlets and (ii) efficiently simulate dynamic discharge and iceberg calving from individual outlet glaciers for a spectrum of climate scenarios. To address this, we have developed a simple model to simulate advance, retreat, and dynamic mass loss due to calving on networks of marine-terminating glaciers (Bassis & Ultee, 2019;Ultee & Bassis, 2016. Our model framework, called SERMeQ, is able to directly simulate decade-to-century-scale evolution of hundreds of outlet glaciers in response to surface mass balance forcing across multiple climate scenarios. This explicit simulation capability, together with recent observations of more than 200 Greenland outlet glacier termini (Joughin et al., 2017a), makes it possible to evaluate our model's performance in a wide range of glacier environments. Here, we show that SERMeQ bounds retreat rates and reproduces patterns of present-day observed changes in terminus position of 155 Greenland outlet glaciers, providing one of the largest validations of any calving parameterization. Supported by this validation, our model physics can be incorporated into global glacier and ice sheet models to constrain iceberg calving. This validation also provides a basis for physically consistent, SERMeQ-derived constraints on future dynamic mass loss from the Greenland Ice Sheet.

SERMeQ Ice Dynamics Model
SERMeQ-the Simple Estimator of Retreat Magnitude and ice flux (Q), sermeq meaning "glacier" in Greenlandic-is a width-averaged, vertically integrated model that determines centerline glacier surface elevation corresponding to a given terminus position. The ice dynamics are based on a perfectly plastic limiting case of a viscoplastic rheology (Bassis & Ultee, 2019;Nye, 1951), with modifications to allow calving at a grounded ice water interface (Ultee & Bassis, 2016) and interaction between multiple tributary glaciers (Ultee & Bassis, 2017). Our flowline-modeling approach is compatible with other flowline-based models such as the Open Global Glacier Model (Maussion et al., 2019), but SERMeQ focuses specifically on near-terminus dynamics of marine glaciers to simulate the calving process.
Rather than imposing an empirical calving rate, SERMeQ self-consistently calculates the maximum rate of terminus advance or retreat at each time step for a given climate forcing. Terminus position evolves in response to near-terminus stretching, bedrock topography, and changes in catchment-wide surface mass balance as described in Ultee (2018) and Bassis and Ultee (2019), (1) is the ice thickness, U = U(x, t) the ice velocity,ȧ =ȧ(x, t) the net ice accumulation rate, H y the thickness at which effective stress within the ice reaches its yield strength (Equation S1 . Glacier ID 1, which is in the Disko Bay region, appears in the lower left; glacier IDs increase clockwise around the map border. Blue diamonds mark the map location of each outlet we simulated, and every 10th glacier ID is labeled and connected to its outlet location in black. A table of MEaSUREs glacier IDs and names appears in the supporting information. Border spaces with no bar correspond to outlets where data were not sufficient to initialize a SERMeQ simulation or where our analysis indicated SERMeQ would not be applicable (see section 2). Yellow bars and map stars show the case study glaciers highlighted in Figure 3. Colored overlay on the satellite map is ice velocity derived from Sentinel-1 observations (ENVEO, 2017), shown on a logarithmic scale such that fast-moving outlet networks appear brighter than slow-moving inland ice.
in the supporting information), and all terms are evaluated at the instantaneous terminus position, x = L(t) (see supporting information Text S1 and S2). For a change in terminus position determined from Equation 1, SERMeQ calculates a new steady-state glacier surface elevation profile and calculates change in glacier volume above buoyancy (supporting information Figure S1). The latter produces a net contribution to global mean sea level (example in supporting information Text S1, not evaluated in this validation exercise).
The only adjustable model parameters are ice temperature T, which is used to calculate the horizontal stretching rate U/ x at the terminus, and yield strength y , which is used to calculate the yield thickness H y (supporting information Text S1-S3). Both are material quantities that can be independently constrained by laboratory and field measurements. Crucially, we do not tune either of our parameters to match changes in terminus position. Comparison of simulated with observed terminus position thus provides a completely independent validation.
Here, we extend the physical realism and applicability of our model to demonstrate that it can simulate calving retreat of a wide variety of marine-terminating glaciers. Novel elements of SERMeQ specific to this application include upstream forcing with surface mass balance from a regional climate model

Identification of Flowline Networks
We first identified 181 Greenland outlet glaciers that have multiple terminus positions recorded in Joughin et al. (2017a). For each glacier, we then defined a network of interacting flowlines with spatially variable width by tracing ice surface velocity from Joughin, Smith, Howat, and Scambos (2017b, and see supporting information Text S5). We extracted ice surface and bed elevation from BedMachine Version 3  and applied a Gaussian filter to produce width-averaged topography. Where the data suggested the presence of transient ice tongues, we removed the floating portion from consideration and simulated the grounding line as the "terminus." We removed three glaciers with long, persistent ice tongues, as SERMeQ is unable to simulate ice tongue evolution. Thirteen of the 181 outlets had initial termini grounded above sea level, and iceberg calving is thus unlikely to dominate dynamic mass changes there. We removed those 13 glaciers from consideration as well. Noisy or missing velocity data that produced flowline networks with unphysical bed topography caused us to remove 10 additional outlets, leaving 155 glaciers for our analysis.
For the remaining 155 outlet glaciers, we defined the initial terminus as the grounded ice point along our central flowline that lies closest to the centroid of the 2006 terminus reported in Joughin et al. (2017a). We optimized a single parameter, the yield strength of ice, to best fit the 2006 observed surface profile, as described in Ultee and Bassis (2017). We used a best guess ice temperature T of −10 • C for all glaciers. We then found the catchment-wide, annual mean surface mass balance forcing for each outlet,ȧ in Equation 1, from HIRHAM regional climate model reanalysis (Lucas-Picher et al., 2012;Mottram et al., 2018;Rae et al., 2012), and simulated resulting changes between 2006 and 2014 in ice extent (Figures 1-3) and volume above buoyancy (supporting information Figure 1). Finally, we compared the simulated changes in terminus position with observed changes reported in Joughin et al. (2017a) for the same period. Because our optimization of y considers only the initial observed surface profile, and the changes in terminus position are an independent response to changes in surface mass balance, this comparison examines an independent model prediction that is not tuned to match observations.

Comparison With Observations
We extracted all available terminus position records from Joughin et al. (2017a) for each year within our simulated period: 2006, 2007, 2009, 2013, and 2014. Each terminus position record consists of one or more points; records with multiple points trace across-flow variation in terminus position. We projected all available points from a given record onto the central flowline of the corresponding glacier network, and we identified the space between the most seaward and most landward points of that projection as the "observational range." We also tracked the change over time in the position of the terminus centroid projected on the flowline, which we identified as the "observed (terminus centroid) retreat rate." Finally, we compared the simulated retreat rates with the observed terminus centroid retreat rates ( Figure 2) and the simulated terminus positions with the observational range (Figures 3 and 4a). Figure 1 shows the total retreat we simulated for each glacier between 2006 and 2014, arranged by approximate outlet position. SERMeQ simulates less than 5 km of length change during the observed period on most outlets. There is no relationship between outlet glacier latitude and magnitude of upper-bound retreat: Simulated glacier response to downscaled climate reanalysis forcing is not a simple function of annual average temperature. Dynamic glacier response depends on glacier geometry, as previous studies have also highlighted Catania et al., 2018;Felikson et al., 2017).

An Upper Bound on Calving Retreat for 155 Greenland Outlets
Equation 1 includes an assumption that the glacier calving front is a yield surface, which produces a theoretical upper bound on calving retreat for a given glacier geometry and surface mass balance (see Bassis & Ultee, 2019). Thus, in the absence of errors in the bed geometry and surface mass balance used, we anticipate that SERMeQ-simulated rates of retreat will generally overestimate observed rates. Figure 2 shows that SERMeQ satisfies this expectation and overestimates retreat for 91% (108/119) of glaciers for which more than two terminus position observations are available to constrain the observed retreat rate.
The bulk model results shown in Figures 1 and 2 summarize multiannual change in terminus position simulated across Greenland. Figure 3 compares observed and simulated terminus position change, for example glaciers where SERMeQ underestimates, overestimates, or correctly captures the observed rate of retreat. Apuseeq Anittangasikkaajuk, which is 2 km wide at the terminus and has a small floating ice tongue, is one of a handful of outlets where SERMeQ underestimates observed retreat (Figure 3a). The simulated terminus positions are still within the (small) observational range in that case. SERMeQ strongly overestimates retreat of Helheim Glacier, a large and high-flux glacier on Greenland's east coast whose terminus approaches flotation (Figure 3b). On Sermeq Kujalleq (Danish: Jakobshavn Isbrae), a very large and well-studied outlet glacier on the southwest coast of Greenland, the simulated retreat of 6 km is comparable to observed retreat ( Figure 3c).

Upper-Bound Retreat Rates Are Realistic
A useful upper bound on calving retreat would consistently overestimate the rate of retreat (Figure 2), simulate terminus positions relatively close to observed termini, and correlate with observed changes. We quantify SERMeQ's performance on the latter indicators in Figure 4.
The histogram in Figure 4a summarizes 404 comparisons of simulated versus observed terminus positions, normalized by each glacier's observational range for each year, such that values within ±1 indicate simulated terminus positions within the observed range. About 40% of simulated terminus positions fall within that range, and 55% of simulated terminus positions are within twice the range of the observed-that is, the simulations are relatively close to the observations. Most simulated terminus positions are more retreated than the observed (positive x axis values in Figure 4a), as expected for an upper bound.
Because we present an upper bound on retreat rather than a calibrated model fit, we do not expect a linear relationship between simulated and observed retreat. Instead, we assess Spearman's rank correlation coefficient for each glacier's terminus positions over time. The coefficient ranges from −1 to 1, where positive indicates that retreat is observed when the model simulates retreat, advance is observed when the model simulates advance, and larger magnitudes of observed and simulated change correspond. Of the 155 glaciers we simulate, is positive for 103, as shown in Figure 4b. For 62 glaciers simulated, ≥ 0.5 and significant at the p < 0.1 level, which indicates a moderately strong and statistically significant relationship between simulated and observed terminus position over time. Only two glaciers have negative significant at the same level. The mean over all 155 glaciers is 0.5.

Discussion
Our simulated upper-bound rate of terminus retreat/advance emerges as a dynamic glacier response to climate forcing and glacier geometry (Equation 1) and does not rely on any tuning to match observations. The two model parameters, yield strength of glacier ice y and ice temperature T, are physical quantities constrained by laboratory and field observations, and neither is optimized against observed retreat rates. The yield strengths we use for most Greenland outlet glaciers simulated here range from 50-250 kPa (supporting information Text S3), well within the range of 50-500 kPa suggested by previous works (Cuffey & Paterson, 2010;Nimmo, 2004;O'Neel et al., 2005). We use an ice temperature of −10 • C, which is also within the range expected from simple physical scaling (van der Veen, 2013), observations (Clow et al., 1996), and modeling (Greuell & Konzelmann, 1994). It is possible that an improved match to observed retreat rates could be found if we did allow parameters ( y , T) to vary within and between glacier catchments or over time (Text S3). However, that would sacrifice the physical upper bound in favor of empirical tuning that cannot be independently constrained by laboratory or field observations. The upper-bound retreat rate computed from Equation 1 can far exceed the observed rate, as shown in Figures 2 and 3b. There are three notable sources of discrepancy between modeled and observed retreat rates: (1) quality of available model input data, (2) performance of automated flowline selection algorithm, and (3) presence of floating ice. First, on small outlets that are rarely visited or studied in detail, the bed topography and climate reanalysis data used as input for SERMeQ may be poorly constrained. As a result, the simulated glacier evolves in response to conditions that do not accurately reflect the local environment, and the simulated change in terminus position is more likely to be inaccurate (Text S6 and Figure S6). Second, on small or slow-moving outlets, or where there are gaps in Sentinel-1 velocity data, our method for tracing flowlines (Text S5) is prone to error. As a result, the simulated glacier has unrealistic geometry and may flow over bedrock features that are not present in a true central flowline of the outlet. Finally, where floating tongues are present, we remove them and simulate the first grounded grid point as the "terminus." This can change the near-terminus stress state, in some cases exposing an unstable wall of thick ice and initiating rapid retreat. Effects (1) and (2) are likely responsible for the underestimated retreat of Apuseeq Anittangasikkaajuk; effect (3) is likely responsible for the overestimated retreat of Helheim Glacier (see Text S6). The first two effects can be mitigated with improved observational data and manual data processing where possible. The third effect reflects upper-bound retreat dynamics that are currently held in check by floating ice but which we speculate could be activated if that floating ice were removed.
The 91% satisfaction of the intended upper bound on retreat rate (Figure 2) supports the utility of our model for producing upper bounds on calving retreat and dynamic mass loss. In contrast to existing estimates of the 21st century calving loss, our approach does not impose a uniform calving rate or outlet glacier speedup factor (DeConto & Pollard, 2016;Goelzer et al., 2013Goelzer et al., , 2020Graversen et al., 2011;Pfeffer et al., 2008); instead, we calculate a theoretical maximum rate of calving retreat that can vary by glacier (Bassis & Ultee, 2019). The result is a physically consistent bound on terminus position change that correlates with observed changes for most glaciers (Figure 4b). Further, our model can track terminus retreat and mass loss from multiple interacting branches of a glacier tributary network (Ultee, 2018;Ultee & Bassis, 2017), ensuring that potentially important contributions to sea level are not overlooked. Within ice sheet scale models, our method could be implemented as a calving criterion at grounded ice-ocean interface cells or used as a module to enhance resolution of outlet glacier networks.
The current version of SERMeQ does not explicitly simulate frontal ablation by submarine melting, which can be a large component of mass loss from both floating tongues and grounded glacier fronts Rignot et al., 2010;Wood et al., 2018). Our derivation of Equation 1, which we emphasize is an upper bound on retreat rate, is consistent with high submarine melt that prevents the glacier terminus from advancing (see Text S4 and Ma, 2018;Ma & Bassis, 2019). However, changes in ocean conditions over time can affect glacier terminus dynamics such that the rate of terminus position change becomes closer to or farther from the theoretical maximum. For example, a decrease in submarine melt rate has been implicated in the recent slowing of Sermeq Kujalleq's retreat (Khazendar et al., 2019). Future implementations Geophysical Research Letters 10.1029/2020GL090213 of our method in larger-scale models may therefore benefit from modifications to account for time-varying submarine melt rates.

Conclusions
We have applied a flowline network model of ice dynamics, SERMeQ, to evaluate an upper bound on annualto decadal-scale calving retreat of 155 Greenland outlet glaciers in response to variable climate forcing. Comparison with nearly a decade of terminus position records from MEaSUREs (Joughin et al., 2017a) shows that the model bounds retreat rate for 91% of glaciers examined and that 55% of simulated terminus positions are within twice the observed range. SERMeQ can also evolve upstream surface elevation with each change in terminus position and compute the resultant loss of ice mass above buoyancy (supporting information Text S1; Ultee, 2018). The upper bound on retreat rate that we construct with SERMeQ will produce a corresponding high-end estimate of the loss of grounded ice mass, consistent with efforts to find an upper bound on the ice dynamics contribution to the 21st century sea level rise. Our approach is especially promising in constraining the dynamic sea level contribution from smaller outlet glaciers that are difficult to resolve in larger-scale continental ice sheet models.

Conflicts of Interest
The authors have declared that no conflict of interest exists.

Data Availability Statement
Data on Greenland outlet glacier terminus position and surface ice velocity come from the MEaSUREs project (Joughin et al., 2010(Joughin et al., , 2017a(Joughin et al., , 2017b, available from the National Snow and Ice Data Center. Surface mass balance forcing comes from the HIRHAM regional climate model for Greenland, maintained by the Danish Meteorological Institute and available online (https://prudence.dmi.dk/data/temp/RUM/ HIRHAM/GREENLAND/). Python code for data processing, simulation, and analysis is maintained in a public GitHub repository, which can be inspected online (https://doi.org/10.5281/zenodo.4075607). The core model code of SERMeQ is available online (https://doi.org/10.5281/zenodo.4075612).