Modeling a Century of Change: Kangerlussuaq Glacier's Mass Loss From 1933 to 2021

Kangerlussuaq Glacier (KG) is a major contributor to central‐eastern Greenland mass loss, but existing estimates of its mass balance over the last century are inconsistent, and specific drivers of change remain poorly understood. We present a novel approach that combines numerical modeling and a 1933–2021 climate forcing to reconstruct its mass balance over the past century. The model's final state aligns remarkably well with present‐day observations. The model reveals a total ice mass loss of 285 Gt over the past century, equivalent to 0.68 mm global sea level rise, 51% of which occurred since 2003 alone. Dynamic thinning from ice front retreat is responsible for 88% of mass change since 1933, with short‐term ice front variations having minimal impact on centennial mass loss. Compared to earlier studies, our findings suggest that KG lost 59% (or 301 Gt) less mass over the century than previously thought.


Introduction
Various studies have shown that the net ice loss of the Greenland Ice Sheet (GrIS) has accelerated over the past two decades.Overall, ice loss increased from 35 ± 29 Gt yr 1 between 1992 and 1996 to 257 ± 42 Gt yr 1 between 2017 and 2020, the latter corresponding to 17.7% of the global sea level rise driven by the ice sheets (Otosaka et al., 2023).A few major marine-terminating outlet glaciers, Jakobshavn Isbrae, Zachariae Isstrøm, and Kangerlussuaq Glacier (KG), have been the main contributors to the increase in mass loss over the past recent decades and jointly drained approximately 13% of the GrIS between 1972 and 2018 (Mouginot et al., 2019).For example, the ice loss from Jakobshavn Isbrae increased from 12.5 Gt yr 1 during 2000-2009 to 23.1 Gt yr 1 during 2010-2018 (Mouginot et al., 2019).Most studies mainly focused on the past two to three decades and tried to understand the most crucial dynamics in the current climatological context.While the dynamic changes to the GrIS of the past few decades have had a significant impact on the short-term, it is unlikely that the present-day mass loss rates are representative of the GrIS on a longer time scale (Howat et al., 2011;Nick et al., 2009).Mass loss from glaciers is typically partitioned between changes in surface mass balance (SMB) and ice dynamics (dynamic thinning).The former involves either a reduction in accumulation, an increase in ablation, or both, while the latter encompasses dynamic thinning generally triggered by the retreat of the ice front (B.M. Csatho et al., 2014;King et al., 2020).To comprehensively assess the drivers of long-term mass loss, it is essential to dissect the contributions from SMB and ice dynamics over many decades.Extending the analysis beyond the last three decades and considering the drivers of the past century ice loss is essential to improving our understanding of the current mass loss trends, and will put future projections of mass loss into historical context.
In this work, we focus on KG, the third largest marine-terminating glacier in Greenland in terms of discharge, that has experienced a ∼20 km retreat since the beginning of the last century (Khan et al., 2020).KG was responsible for approximately 3% of the mass balance of the Greenland Ice Sheet (GrIS) and approximately 5% of the ice discharge between 1972 and 2018 (Mouginot et al., 2019).Observational studies have shown that KG lost 80 Gt of ice between 2004 and 2008, a threefold increase in mass loss relative to pre-2004 rates (Howat et al., 2007;Luckman et al., 2006), confirmed in Mouginot et al. (2019), who found a cumulative mass loss between 1972 and 2004 of 2.7 Gt and 159.8 Gt after 2004.During its observational record, KG has at least experienced four retreat events of 5 km or more: between 1932 and 1933 (Vermassen et al., 2020), in 1996(Bevan et al., 2019), in the winter of 2004-2005(Luckman et al., 2006), and again between 2016and 2018(Brough et al., 2019), the latter three characterized by an absence of the usual sustained winter advance (Bevan et al., 2019).Only a few observational studies have set out to investigate present-day events in a long-term transient context, going beyond the first satellite observations in 1972.Khan et al. (2014) suggests that KG experienced substantial thinning before 1930 and minor thinning between 1981 and 1998.Kjeldsen et al. (2015) finds that KG was responsible for 14% of the total GrIS mass loss between 1900 and 1983 and claims that SMB is the primary driver of mass loss on centennial time scales.On the contrary, Khan et al. (2020) found that ice dynamics explained 90% of the centennial mass loss and calculated a total mass loss for KG between 1880 and 2012 of 1,381 ± 178 Gt.We note that Khan et al. (2020) and Kjeldsen et al. (2015) did not consider the same area.Kjeldsen et al. (2015) estimated mass loss from the 2012 terminus and upstream, while Khan et al. (2020) estimated mass loss for the area between the Little Ice Age (LIA) terminus and upstream.Furthermore, the mass loss estimates of Khan et al. (2020) and Mouginot et al. (2019) disagree, when using the 1972 ice mass as a reference for the mass balance, that is, the initial mass of Mouginot et al. (2019).So far, the methods used to study KG have relied on correlating present-day retreat with thinning and assuming that this correlation was the same a century ago.The inconsistencies between these studies underline that the centennial mass balance of KG and the specific drivers of change remain poorly understood, highlighting the need for numerical modeling to constrain the sparse historical data with known physics.

Observations and SMB Reconstruction
Here, we use the Ice-sheet and Sea-level System Model (ISSM) (Larour et al., 2012) to simulate KG between 1933 and 2021.To initialize the model, we use a digital elevation model (DEM) of the ice surface in 1933 from Khan et al. (2014); Kjeldsen et al. (2015) together with estimates of the basal condition.The initial 1933 DEM is based on surface elevation derived from a trim line-based method adopted from Khan et al. (2020), which extrapolates the cross-sectional and longitudinal profiles of KG 1978KG -1987 to the position of the lateral moraine and ice margin moraine of the LIA maximum, assuming the profiles to be constant.Near the ice margin, the 1933 DEM has an uncertainty of 12 m (B.Csatho et al., 2008;Kjeldsen et al., 2015).Figure 1b shows the bed topography from BedMachine Greenland v4.10 (Morlighem et al., 2017).Assuming that bed topography has not changed significantly since 1933, we can compute the ice thickness in 1933, Figure 1e, based on the 1933 DEM shown in Figure 1d.We compute the initial grounding line position using a flotation criterion, leading to an initial area of 2,300 m 2 of floating ice at the front.We use the depth-averaged temperature from the ISSM contribution to ISMIP6 (Goelzer et al., 2018) to obtain a spatially variable rheology parameter in Glen's flow law (Glen, 1958) with n = 3 and temperature dependence given in Cuffey and Paterson (2010, p. 75).We assume that the depthaveraged temperature is approximately constant on centennial time scales (Seroussi et al., 2013).We weaken the ice rheology in the shear margins by a Glen enhancement factor of 0.58, found by searching in the vicinity of commonly chosen values to obtain more realistic transitions from slow to fast flow velocities (Cook et al., 2022;Graham et al., 2018;Greve & Blatter, 2009, p. 58).The sensitivity of this choice is shown in Figure S1 in Supporting Information S1.After initialization, the model is exclusively driven by external forcing from a SMB model (SMB) and the ice front positions.
We found that the Box SMB (Box et al., 2013), which spans 1840-2012, has a bias toward mass loss when compared to RACMO 2.3p2 SMB (Noël et al., 2018(Noël et al., ), spanning 1958(Noël et al., -2021, and the Regional Atmosphere Model (MAR) (Tedesco & Fettweis, 2020), spanning 1950-2021, when evaluated after 1958.Therefore, we employ a reconstruction of the SMB before 1958 by anomaly correcting the Box SMB.We compute SMB anomalies based on the 1972 and 1990 SMB trend, where the ice sheet in Central Eastern Greenland was in mass balance (Mouginot et al., 2019): where SMB R is the RACMO SMB, SMB B is the Box SMB, 18 is the number of years between 1972 and 1990, and we evaluate t ∈ [1933; 1958].The reconstruction effectively adds the Box SMB anomaly from the period before 1958 to the average RACMO SMB, computed in a period of mass balance, to reconstruct a RACMO-like SMB between 1933 and 1958.After 1958, RACMO was used without any modifications; in addition to the temporal correction, we also need to address the reconstruction spatially, as RACMO SMB does not provide SMB values beyond the spatial limits of the 1958 ice domain.Therefore, we used the spatial average of the first term of Equation 1 to fill in missing SMB values outside the 1958 extent.This reconstructed RACMO SMB spans 1933 to 2021 both in space and time, and we denote it SMB * R and in a similar way, we obtain an anomaly corrected MAR SMB denoted SMB * M , for later comparison.
The historical observations of the ice front position from 1933, 1966, 1972, and 1981, shown in Figure 1f, originate from various sources described in more detail in Khan et al. (2014) and Kjeldsen et al. (2015).Additionally, we utilize the two recent automatic terminus delineation data sets called Calving Front Machine (CALFIN) (D.Cheng et al., 2021) and AutoTerm (Zhang et al., 2022) respectively, with observations spanning 1972-2021.The data sets yield 1188 front observations, with a temporal distribution summarized in the histogram shown in Figure S2 in Supporting Information S1.The front observations force the modeled front migration in time, using a level-set method (Bondzio et al., 2016).We use an adaptive time-stepping scheme following the Courant-Friedrichs-Lewy condition (Courant et al., 1928) and obtain front positions for all time steps by linear interpolating between the observed front positions.As both MAR and RACMO are similar in quality, we used a model forced by SMB * R and all available front observations as a reference model.We explore the sensitivity to these choices in the results section.

Ice Sheet Model and Basal Friction
We use ISSM to employ the Shelfy Stream Approximation, defined on a non-uniform triangular finite element mesh, locally refined based on present-day observed surface velocities.The element sizes range from 350 m in the fast flow region to 12 km in the interior, resulting in a model with 26,784 elements.The model includes grounding line dynamics based on hydrostatic equilibrium, and applies a basal melting rate of 20 m yr 1 under floating ice, when it forms.This choice of basal melting rate did not affect the final results, as the floating extensions are small throughout the simulation.
We estimate basal friction, based on the Budd friction law (Budd et al., 1979) formulated as in G. Cheng et al. (2022), by solving an inverse problem (Morlighem et al., 2010) based on a mosaic of observed surface velocity, with the nominal year of 2007 (Joughin et al., 2016(Joughin et al., , 2018) ) and present-day ice geometry.The resulting velocity misfit is shown in Figure S3 in Supporting Information S1.The inland boundary conditions of velocity and thickness are prescribed from these observations and assumed constant in time.We initialized the model with surface elevation and bedrock topography with the nominal year 2007; see Figures 1b and 1c.The inversion procedure only estimates the friction field beneath the present-day ice, but the ice sheet covered a larger surface area in 1933.We extrapolate the friction fields to the extent of the 1933 front position by correlating the friction coefficient to the bed topography, see Figure S4 in Supporting Information S1, and assume that the inverted friction coefficient remains constant in time (Bondzio et al., 2017).We correlate bed topography and friction coefficient using a fourth-degree polynomial, which had a better present-day trade-off between velocity, flowline, and surface elevation misfits compared to other simple extrapolation strategies, for example, 2d linear, nearest neighbor and constant value extrapolation.We added a constant to the extrapolated friction coefficient to control the magnitude of the initial ice velocities, the sensitivity to this choice is seen in Figure S5 in Supporting Information S1.We have summarized the methodology in the flowchart in the Figure S6 in Supporting Information S1.Using the 1933 ice geometry, the estimated basal conditions, and the external forcing in the form of the reconstructed SMB and the front observations, we run the transient model forward in time, completely unconstrained by any additional data from 1933 to 2021.

Experiments
We conducted various experiments to investigate the main contributing factor to the mass loss of KG during the last century.By fixing the 1933 front in time, such that mass loss is only forced by SMB; we investigate how much the SMB and the retreat of the ice front contribute individually relative to the overall modeled mass loss.Additionally, we investigate the sensitivity to the choice of SMB model.We also present four different front-forcing experiments, using a subset of the total available front observations, denoted as s A, B, C, and D. The experiments use front observations from the following years: A: 1933, 2021, B: 1933, 1981, 2021, C: 1933, 1966, 1981, 1999, 2021, D: 1933, 1966, 2021.The front-forcing experiments explore how much the short-term variations of the front position affect the mass loss on a centennial time scale.Experiment A only uses the first and last front observations, such that the front retreats steadily throughout the transient run, as intermediary front positions are obtained by linear interpolation.Experiment B introduces one more front, as close to the middle of the time series as possible, and C uses five front observations spaced approximately equally in time.Finally, D only uses the first, the last, and the 1966 front, which was a relatively retreated position compared to the 1980s; see Figure 1f.

Model Misfit
We compare our modeled velocities to observations between 2007 and 2018 to assess the performance of our model when run freely forward from 1933.The average velocity misfit in Figure 2a, for MEaSUREs ITS_LIVE v1 mosaics (Gardner et al., 2022;Shepherd et al., 2020) between 2007 and 2018 shows that the misfit is highest close to the shear margin in the fast flow region with maximum magnitudes of 2 km yr 1 which is ∼22% of maximum front velocities of KG.We adopt the SMB error estimate for RACMO 2.3p2 from Mouginot et al. (2019), who estimate the error for KG based on a mass balance reference period by assuming that the bias between observations and RACMO 2.3p2 expresses the model error and that it is constant in time.This yields a constant error for KG of ±1.1 Gt yr 1 , equivalent to 0.024 m yr 1 for the KG domain.Figure 2b shows the difference between modeled and observed surface elevation from ICESat-2, h model h obs (Khan et al., 2022), averaging to 39.4 m in total, corresponding to a bias of 0.45 m yr 1 over the century.We further explore how the model compares to observations in Figure S7 in Supporting Information S1, which shows the modeled surface elevation change and the corresponding misfit, compared to ICESat-2 estimates from Khan et al. (2022), between 2003 and 2021 in different basin sections.Our model predicts a basin-wide average change in surface elevation of 2.56 m over 18 years, with more thinning in the fast flow area and slight thickening in some interior areas.The area-averaged change in surface elevation and corresponding root mean squared error is 9.97 ± 3.66 m yr 1 below 1,000 m surface elevation, 4.48 ± 1.13 m yr 1 between 1,000 and 1,500 m, 0.87 ± +0.32 m yr 1 between 1,500 and 2,000 m, 0.06 ± 0.24 cm yr 1 between 2,000-2,500 m and 3.33 ± 5.33 cm yr 1 above 2,500 m. Figure 2c shows that the model has a bias toward too low velocities in the years 2012 to the summer of 2014 and too high velocities between 2018 and 2021 in the first few kilometers upstream of the front.The average error along the flow lines for all times is 5.1 × 10 3 km yr 1 .

Mass Balance of the Past Century
In Figure 3a we show the reconstructed mass balance of KG together with discharge and SMB, computed as Mass Balance = SMB Discharge.Figure 3a shows that it is an increase in discharge which has driven the mass loss of KG since the middle of the 90s until 2021.The integrated discharge and SMB can be found in Figure S8 in Supporting Information S1.We investigate the cumulative mass balance over the last century in Figure 3b, where the black curve, based on RACMO SMB and all available forcings, serves as a relative reference for the other results.Our model predicts a glacier-wide mass loss between 1933 and 2021 of 285 Gt.In Figure S9 in Supporting Information S1, we present a 5-year and a 20-year moving average of the yearly mass loss rates.Prior to year 2000, the mass loss rate fluctuated around a mean of 2 Gt yr 1 , with some local departures, for example, between 1968 and 1970, where the mass loss rate briefly dropped to 7.4 Gt yr 1 .Since the early 2000s, the mass loss rate underwent a significant increase, and for the past 20 years, the average mass loss rate was 7.2 Gt yr 1 .Notably, the rate of change of the mass loss rates increased 22 times from 0.01 Gt yr 2 pre-2000 to 0.22 Gt yr 2 post-2000.
Our mass balance can be compared to the results of Mouginot et al. (2019) and Khan et al. (2020), referencing their results to our modeled mass loss in 1972 and 2004, respectively.The trend of the reference mass balance shows good relative agreement with results from Mouginot et al. (2019) between 1972and 2019and Khan et al. (2020) 1972-2012, displaying similar patterns of equilibrium between 1972 to approximately 1994 and significant acceleration of the mass loss since the early 2000s.Our results lie within the uncertainties of Mouginot et al. (2019) and predict a lower absolute mass loss than both Mouginot et al. (2019) (23 Gt difference in 2019) and Khan et al. (2020) (60 Gt difference in 2012) in relative terms, as seen in Figure 3b. Figure 3b also shows the result of models using alternative SMB models, the cyan results from using MAR and yellow from Box SMB (Box et al., 2013), resulting in a final 20% and 24% relative difference to the reference model.The Box SMB forced model does not replicate the period of mass balance in the 70 and 80s.The brown curve, the control run, is the most significant relative difference of 88% in Figure 3b. Figure S10 and Table S1 in Supporting Information S1 illustrate how the discrepancies between a given model run and the reference model develop throughout the simulation period.The control run starts with a minor discrepancy of 0.4 Gt and increases to 252.3 Gt at its maximum with a median difference of 55.4 Gt.The Box-based mass loss has a relatively consistent offset in the whole period with a median difference of 62.44 Gt, a minimum of 0.2 Gt, and a maximum of 94.9 Gt.In summary, control underestimates the mass loss significantly compared to reference, while Box-based mass loss overestimates the mass loss by a similar magnitude.
Figure 3c shows various front-forcing experiments, with A, B, C, and D resulting in 21%, 24%, 29%, and 1% relative mass balance differences to the reference model, respectively.The experiments represent various simplified front migration patterns, all characterized by consistent retreat over the century.The experiments show that centennial mass loss can be reproduced with an accuracy of 20-30% without accurate frontal dynamics.Table S1 in Supporting Information S1 shows the relative difference between the four front-forcing experiments and the reference model.Table S1 in Supporting Information S1 shows that front-forcing experiments A and C have similar magnitude statistics over the simulation period, with a median difference of 0.13 and 8.7 × 10 11 Gt, respectively, together with similar minimum values of 27.4 and 12.3 Gt and maximum values of 53.42 and 81.89 Gt, respectively, which can be interpreted as the average front retreat having a similar efficacy as the more detailed front migration in C. Likewise, front-forcing experiments B and D also share similar magnitude in the statistics, but with opposite signs.

Centennial Mass Balance Reconstruction
As seen in Figure 2, the centennial reconstruction of KG exhibits an excellent level of agreement with present-day observations of ice velocity (Gardner et al., 2022;Shepherd et al., 2020) and surface elevation (Khan et al., 2022).The fact that the model agrees so well with present-day observations is an achievement, given that the model ran freely for 88 years, starting from an initial ice surface and only forced by climate data and ice front positions.The model yields a mass balance time series, which shows similar patterns as other studies, as shown in Figure 3b.The median velocity misfit level of 47.457 m yr 1 and 12.19 m in surface elevation highlights the model's efficacy in capturing the essential dynamics of the glacier and its response to climate forcing.
Prior assessments of mass loss for the KG basin have primarily relied on remote sensing techniques and data interpolation, yielding absolute estimates that significantly differ from our current modeling results.While previous estimates by Khan et al. (2020) suggested a total mass loss of 510 Gt between 1933 and 2012, our approach, which leverages known ice physics through numerical modeling, estimates a more conservative 209 Gt mass loss during the same period.We find a mass loss rate between 1933 and 2021 of 3.2 Gt yr 1 and 2.7 Gt yr 1 between 1933 and 2012 which is less than half than that reported in Khan et al. (2020).As seen in Figure S8 in Supporting Information S1, the moderate and stable mass-loss rates before the early 2000s underline that prior Geophysical Research Letters 10.1029/2023GL106286 studies have most likely overestimated the mass loss before the satellite age, as the present mass-loss rates are significantly higher than those of the last century.
The model displays sensitivity to the choice of SMB, changing the absolute mass loss by approximately ±20%.While Box SMB is defined back to 1840, it seems that it smooths the mass loss curve out in Figure 3a, compared to the RACMO and MAR mass loss curves, which have many short-term perturbations throughout the transient run.The quality difference between the SMB models is most likely because RACMO and MAR are comprehensive regional climate models that simulate various atmospheric processes, in contrast to Box SMB.

Basal Friction
We implemented various other friction laws, such as the Schoof & Gagliardini (Gagliardini et al., 2007;Schoof, 2007) and a plastic variation of the Budd law (Cornford et al., 2020).However, they either led to unrealistic initial velocities or unexpected mass accumulation.While we could get reasonable inverse solutions with the different friction laws with present-day data, only the Budd friction-based model produced realistic initial velocities and a meaningful centennial mass balance trajectory.The inadequacy of the other friction laws arose from a combination of the extrapolation method and uncertainties in initial ice geometry.Like the findings in Berends et al. (2023), which show how the friction coefficient effectively compensates for uncertainties in unmodeled physics, the extrapolation compensates for errors in the initial state.In all other cases than Budd, the constant offset had to be set at a very high level to produce realistic initial velocities, which in turn resulted in slowing down the ice too much, causing a significant build-up of ice on the long transient runs.In Text S1 in Supporting Information S1, we show that for the Schoof & Gagliardini law, the initial velocity will grow cubically with effective pressure, making this law more sensitive to uncertainties in initial ice thickness, explaining why it was not suitable for the overall approach explained here.

Impact of Ice Front Variations
The control experiment reveals that the migration of the ice front over the last century explains 88% of the total mass loss.The total mass loss might be reproduced accurately with less information, that is, three front observations in B and D. However, it might also diverge, depending on whether the observation in the middle of the time series favors a higher or lower discharge.Despite experiment D showing that it is possible to reproduce the final absolute value of the mass balance with only a few front observations, the overall behavior of the mass balance curve is off.The final 1% difference is likely a coincidence, emphasizing that one should be careful when interpreting these results.More front positions are necessary to recreate intermediary states of the glacier.As shown in Figure 1f, the 1966 front is more retreated than in the 1970s and 1980s, and the resulting mass loss is almost linear throughout the century.To reproduce the mass balance for 1972-1990 requires the front to readvance from the 1966 position.Regarding the final and overall discrepancies between the various frontforcing experiments and the reference mass loss over the whole century, we find that all the front-forcing experiments have a bias below 30% compared to the reference model.While our front-forcing experiments do not directly consider seasonal variability from year to year, they explore how front migration accuracy impacts mass variations on centennial time scales.In this sense, our results align with Felikson et al. (2022), who showed that disregarding seasonal oscillations in the front position leads to a bias in mass loss of the same order of magnitude.

Conclusions
Using a combination of historical records and present-day observations, we presented a transient ice flow model of KG between 1933 and 2021 that agrees excellently with present-day ice velocity and surface elevation observations.The model reveals that KG was losing mass at a steady rate until the early 2000s.Since 2003, KG lost 51% of the total centennial mass loss, corresponding to an average mass loss rate of 7.2 Gt yr 1 .Control experiments highlight the role of ice front retreat as a primary driver of the modeled mass loss of KG, explaining 88% of the total mass loss.The short-term variability of the ice-front position is responsible for up to 29% of the mass loss on longer time scales.These findings have important implications for predictions of future mass loss and, consequently, sea level rise projections.Our study emphasizes that KG's recent acceleration constitutes a recent shift in behavior and underscores the pivotal role of sustained ice front retreat in steering the overall mass loss.

Figure 1 .
Figure 1.(a) Observed mosaic surface velocity (Joughin et al., 2016, 2018), with the Kangerlussuaq Glacier (KG) basin shown in red in the small Greenland map.(b) Bed topography from BedMachine Greenland v4.10.(c) Ice surface elevation from BedMachine Greenland v4.10 (Howat et al., 2014, 2015).Data in panels (a-c) pertains to the nominal year 2007.(d) 1933 ice surface elevation from Khan et al. (2014); Kjeldsen et al. (2015), (e) 1933 Ice thickness resulting from the difference between the 1933 surface and the bed topography.(f) Front observations 1933-2021 and location of KG from Khan et al. (2014), Kjeldsen et al. (2015), D. Cheng et al. (2021) and Zhang et al. (2022).Panels (a-e) show the same area with the same scale, while (f) is zoomed in.Background satellite image mosaic in all panels is from MacGregor et al. (2020), and the black dashed line marks the model domain.

Figure 3 .
Figure 3. (a) The mass change rate (red), computed as Mass Balance = SMB Discharge, for the reference model decomposed into SMB (blue) and ice discharge (green).The smallest filled symbols are monthly mean values; the empty symbols are yearly mean values, and the solid lines are fitted with 4th-degree polynomials to the monthly values, illustrating the trend.(b) The full plot is a zoomed in version of the shaded area in the small panel and shows cumulative mass change, with negative values denoting mass loss and positive values denoting mass gain.The black curve is the output from the reference model, which uses SMB * R , the cyan is the same model but uses SMB * M , the yellow curve uses SMB B , and the brown curve shows the control experiment, which also uses SMB * R , keeping 1933 front position fixed.The gray dots and the shaded area are mass balance and error estimates from Mouginot et al. (2019) (referenced to 1972), and the purple diamonds are estimates from Khan et al. (2020) (referenced to 2004).(c) Cumulative mass change from experiments with different front forcings.The blue curve is Experiment A, the green curve is Experiment B, the orange curve is Experiment C, and the brown curve is Experiment D. The final relative differences in 2021, compared to the reference model, are shown as bars on the right side of both panels (b, c).