2.1 Glacier model and idealized geometry

We use a 1-D model of marine-terminating glacier dynamics that assumes flow is dominated by longitudinal stretching and Weertman-style sliding (see Supplement for equations). Ice velocities are calculated using the shallow shelf/stream approximation (SSA), ice thickness is calculated through mass conservation, and the terminus position is calculated using an implicit flotation condition. Following the numerical approach laid out in Schoof (2007), the governing equations are solved on a grid stretched to the glacier extent at each time step. That is, if the glacier advances or retreats, the grid nodes shift with respect to a fixed spatial coordinate. The model does not simulate any floating ice, so the last grid point always represents the grounding line (i.e., a terminus at flotation). The grid spacing is several kilometers through most of the interior but is significantly refined (∼100 m) near the grounding line. Ice flux across the grounding line is not parameterized but rather calculated directly from the local stress balance. Marine-terminating glacier models using this numerical approach have been benchmarked and used in many prior studies (e.g., Schoof, 2007; Robel et al., 2018).

The majority of the simulations in this study investigate an idealized glacier whose terminus position is near a small bed peak (Fig. 1). Our standard geometry has an inland prograde bed slope of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}$ and a sharp peak 350 km from the ice divide. On its landward face, the peak rises 88–100 m (depending on experiment) in 4 km, and its seaward face has a slope of approximately $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ (Fig. 1). We chose this idealized peak in order to provide a very distinct boundary between stable and unstable terminus positions. However, we introduce more realistic bed geometries with many bed peaks of varying geometries in Sect. 5. The time-averaged surface mass balance (SMB) is +0.7 m yr⁻¹ through most of the interior and smoothly decreases over the last ∼50 km to approximately −1.4 m yr⁻¹ at the terminus (Fig. 1a, top).

We simulate ocean forcing very generally by adding a frontal ablation term to the mass conservation equation at the grid point closest to the terminus (see Supplement for numerical implementation). This frontal ablation term amounts to an additional output flux beyond that of the local ice velocity. Compared to the case without the frontal ablation term, this results in a retracted terminus with a steeper surface slope such that ice flow balances the additional output flux. For real glaciers, flux anomalies could be driven by variable calving, submarine melt, or a combination; here we simply interpret these as flux anomalies at the terminus driven by variable ocean conditions. The most salient aspect is that frontal ablation is a localized flux term that we can force to vary in time. We note that because the frontal ablation is implemented at a single grid point, its effect on surface slopes and the steady-state terminus position depends on the grid size (see Supplement, Fig. S2) and also creates numerical convergence issues for very fine grids (<30 m). However, we use a consistent grid scheme whenever comparing simulations, so this does not affect our overall conclusions. Our standard glacier has a mean frontal ablation rate of 30 m yr⁻¹. This is within the range of observed submarine melt rates, which vary by orders of magnitude (cf. Motyka et al., 2011; Mouginot et al., 2015).

2.2 Climate and glacier variability

Internal climate variability, which arises from chaotic fluctuations in atmospheric and ocean circulation, directly affects accumulation, surface melt, and ocean forcing on marine-terminating glaciers. We simulate the effects of internal climate variability on both SMB and frontal ablation as a first-order autoregressive process (AR-1), which is a common model for realistic climate variability with temporal persistence (i.e., “memory”; Hasselmann, 1976). We generate time series of annual climate anomalies (zero-mean) using a Fourier-transform method (see Percival et al., 2001; Roe and Baker, 2016; Christian et al., 2020). This allows us to generate random frontal ablation and SMB time series that are correlated with each other (i.e., reflecting the same regional climate anomalies) but which have different levels of prescribed persistence. At each model time step, we add the corresponding anomaly to the mean SMB or frontal ablation term in the model's continuity equation, averaging over the time step if it is greater than 1 year. For most simulations, we use a model time step of 5 years; we found little effect on terminus fluctuations with time steps of 1–10 years as the high frequencies of climate variability are strongly damped by the ice dynamics either way (see Supplement, Fig. S3).

We set a decorrelation timescale of 10 years for frontal ablation anomalies, which emulates the decadal variability known to be an important aspect of ocean forcing in both Greenland and Antarctica (e.g., Andresen et al., 2012; Straneo et al., 2013; Jenkins et al., 2016). We explore a range in the magnitude of frontal ablation variability across our experiments, with standard deviations from σ_FA=12 to 27 m yr⁻¹ (when sampled annually). We set a lower bound of zero on absolute frontal ablation values (i.e., we truncate anomalies more negative than 30 m yr⁻¹). This avoids numerical issues that can follow strongly negative ablation anomalies, which add ice to the terminus. We note that this lower bound introduces an asymmetry in the variability for large σ_FA, which affects the mean state as variability is increased. However, based on tests without a zero bound, this appears secondary to the main effects of increasing overall variability (see Supplement, Fig. S6). It should also be noted that other nonlinearities in the ice dynamics affect the mean state when variability is imposed (e.g., Robel et al., 2018).

For SMB, we set a decorrelation timescale of 1.5 years, consistent with higher-frequency atmospheric variability. We set σ=0.4 m yr⁻¹ near the terminus, tapering to σ=0.1 m yr⁻¹ inland, reflecting the greater variability in snowfall and melt near ice-sheet margins (e.g., Fyke et al., 2014).

Figure 2A synthetic glacier's response to stochastic variability, when the terminus is near a peak in bed topography. Normalized climate anomalies are shown on the upper axes, and grounding-line fluctuations are below. Dashed line indicates the location of the bed peak. A long period of stable, stationary fluctuations precedes unstable retreat off of the bed peak. The histogram (right) shows the probability distribution of stable fluctuations.

Download

Figure 2 shows a single realization of stochastic SMB and frontal ablation forcing (the latter with σ_FA=12 m yr⁻¹) and the response from the idealized glacier shown in Fig. 1. The grounding line fluctuates in a roughly 4 km range on the seaward side of the bed peak for over 14 kyr before retreating permanently from the peak due to one particularly persistent (but random) anomaly in both frontal ablation and SMB. Retreat continues over several centuries until the terminus reaches a new stable mean state several tens of kilometers inland (not shown), consistent with the marine-ice-sheet instability. This simulation illustrates the basic null hypothesis for attribution: a sustained retreat from a bed peak initiated purely by natural variability. It also highlights the stochastic nature of such retreats: the glacier fluctuates for a long time in close proximity to the topographic threshold, before permanently retreating. The frontal ablation and SMB anomalies that triggered retreat occurred purely by chance in the sequence of stochastic forcing.

2.3 Ensemble framework

The stochastic nature of noise-driven retreat from a bed peak (Fig. 2) suggests that for attribution, it may be more productive to consider such retreats as discrete events rather than trends, as well as to focus on their likelihood rather than their total magnitude, which is set largely by bed topography. Our method is to run large ensembles (i.e., collections) of simulations with the glacier model, in which each simulation is forced with an independent realization of random climate variability drawn from the same distribution. This assumes that the statistics of the forcing variability are related to constant physical characteristics within the climate system, but the particular sequence of anomalies is fundamentally random. Because any realization of such climate variability is equally plausible, the ensemble thus provides a statistical population from which we can estimate the likelihood of a particular event, such as a sustained retreat from the bed peak. We refer to these ensembles as “aleatory” ensembles because the only difference between members is the realization of random forcing.

Unless otherwise indicated, we assess glacier retreats within an experimental interval of 150 years, corresponding to the approximate time frame of significant anthropogenic forcing (e.g., IPCC, 2021). Within this interval, we count ensemble members with termini that start seaward of the bed peak and the retreat more than 4 km behind the bed peak. This threshold corresponds to the inland extent of the retrograde bed (Fig. 1), so it is a reasonable indicator that the retreat is due to a loss of terminus stability and that the terminus will continue retreating (indeed, this is observed in longer simulations). However, since stability thresholds can be ambiguous for transient glacier states (Sergienko and Wingham, 2021; Robel et al., 2022), we refer to retreats exceeding this threshold as rapid and sustained retreats. The threshold is within the range of many observations of retreats from bed peaks (e.g., Catania et al., 2018; Wood et al., 2018). For attribution assessments on real glaciers, the threshold for what counts as a retreat within ensemble simulations would be informed by the observed retreat extent. We return to the question of defining retreats in the discussion section.

Rather than starting all simulations with a strictly steady-state glacier, we initialize simulations with a 250-year period of stochastic forcing. This is necessary because of noise-induced drift that occurs at the onset of stochastic forcing due to nonlinearities in ice dynamics (e.g., Robel et al., 2018). Indeed, we find that the steady-state grounding line position is closer to the bed peak than the long-term mean under noisy forcing, slightly enhancing the likelihood at the beginning of simulations initiated from a steady state (see Supplement, Figs. S4–S5).

From each aleatory ensemble, we estimate the probability of sustained retreat as the number of retreats within the 150-year experimental interval (as defined above) divided by the total number of simulations with termini seaward of the peak at the beginning of the interval. That is, we are fundamentally focusing on a conditional probability of industrial-era retreat (i.e., conditioned on the glacier not having already retreated). We ran several long ensemble simulations to assess how this conditional probability varies in time and find it to be fairly stable after the noise-induced drift decays, though with some additional caveats at millennial timescales (see Supplement, Fig. S5). We assess the role of anthropogenic forcing by comparing the probability of retreat between an ensemble with constant mean climate and another ensemble with an anthropogenic trend in the mean climate added to all members (Sect. 4).

The aleatory ensembles are the key tools of the attribution framework we propose and are necessary because of the fundamentally chaotic nature of climate variability. A few previous studies have analyzed aleatory ensembles for ice-sheet dynamics, focusing on uncertainty in future response (Tsai et al., 2017; Robel et al., 2019; Tsai et al., 2020), but the present study is the first application to the question of attribution. However, a different category of ensemble methods also exists for quantifying “epistemic” uncertainty in the design or parameter choices within the model. In principle, such epistemic uncertainty can be reduced by adding observational or theoretical constraints, whereas aleatoric uncertainty is irreducible due to the lack of predictability of climate variability. These ensemble methods have recently seen wider use for ice-sheet models. Examples include intercomparison projects targeting structural differences between models (e.g., ISMIP6; Nowicki et al., 2016), perturbed-parameter ensembles (e.g., Applegate et al., 2012; Aschwanden et al., 2019; Nias et al., 2019; DeConto et al., 2021), and ensembles of synthetic subglacial topography (MacKie et al., 2021). We also employ ensemble methods in this broader epistemic category. In the next section, we assess the sensitivity of the probability of retreat to several key parameters by running aleatory ensembles over a range of parameter perturbations (note that this means two layers of ensemble design over two different types of uncertainty). Finally, we also consider an ensemble of different bed topographies (Sect. 5) in order to investigate attribution in the context of a regional population of heterogeneous glaciers.

We note again that the synthetic geometries and climate parameters we use are closer in scale to marine-terminating glaciers in Greenland than Antarctica. However, the principles of the ensemble attribution framework constitute a general approach for assessing glacier variability and retreats near topographic thresholds. These methods could thus be applied to other marine-terminating glaciers in a wide range of contexts.

3.1 Which model parameters affect the probability of rapid retreat?

With multiple uncertain parameters that must be prescribed in any model, it is important to consider how sensitive retreat probabilities estimated from an ensemble will be to such uncertainties. We assess this by running groups of aleatory ensembles, where each ensemble has slightly different model parameter values (model parameters are still identical between the 1000 simulations within a single aleatory ensemble). We perturb the bed peak height (from 88 to 100 m), the mean interior SMB (64 to 71 cm yr⁻¹), and the friction coefficient (reductions of 2 % to 22 % from the default value). We vary only one parameter at a time, but for each group of perturbations, we also apply three levels of frontal melt variability (σ_FA of 15, 21, and 27 m yr⁻¹).

Figure 3c shows the retreat probabilities of 51 ensembles, each with 1000 members, across this parameter space. Retreat probabilities are plotted against the distance between the initial steady-state terminus position (before being forced by climate variability) and the bed peak. Although the steady-state terminus position is not strictly the same as the time mean of fluctuations (see Supplement, Fig. S4; Robel et al., 2018), it still provides a metric for comparing how these parameter perturbations affect how close the system is to the topographic threshold. From Fig. 3c, we can see that the probability of rapid retreat varies primarily along two dimensions; the proximity of the glacier terminus to the bed peak and the magnitude of noisy forcing. Regardless of which model parameter is changed (besides σ_FA), the probability of retreat ultimately collapses onto curves that are solely a function of the steady-state distance of the terminus from the bed peak (i.e., the different colored markers in Fig. 3c). It is possible that parameter perturbations can affect dynamical stability in other ways (e.g., Parizek et al., 2013), but the similarity of these curves for qualitatively different parameters (sliding, mass balance, and bed geometry) indicates that in this case, proximity to the bed peak is the main effect. Similarly, the effect of increasing the climate variability is fairly consistent for a given initial steady state regardless of which glaciological parameters are perturbed to achieve that state.

Setting up attribution experiments for real glaciers will require prescribing multiple uncertain parameters, and the probability of retreat can change quickly across a relatively small range (Fig. 3c). It will thus be important to evaluate how sensitive estimates of retreat probability are to these uncertainties. We note this may also include choices in numerical methods, such as the manner in which frontal ablation is prescribed (see Supplement, Fig. S2). Although there are many parameters to consider, Fig. 3c shows that the priority should be for sensitivity tests to adequately sample the plausible range of mean terminus position prior to retreat and of the magnitude of forcing variability. Both of these key constraints will likely require better observations of glacier state and climate from before the satellite era.

3.2 What causes individual glacier retreats?

Beyond providing a means to estimate the probability of glacier retreat, aleatory ensembles also provide large synthetic datasets that can offer further insight into the process of sustained retreat in the context of climate variability. In Fig. 4 we analyze a 1000-member ensemble, run for 500 years with frontal ablation variability that has σ_FA=18 m yr⁻¹. We use the standard idealized geometry with a bed peak of 90 m. From this ensemble, we find 375 simulations that retreat more than 4 km, which we will refer to as the “retreat group”. These retreats begin at random times within each simulation. In order to analyze them together, in Fig. 4a we synchronize the time series such that the year 0 corresponds to the last time that the terminus is forward of the bed peak. We conduct the same synchronization for each ensemble member's corresponding time series of frontal ablation forcing so that the anomalies causing retreat are lined up (Fig. 4c). Then, for comparison to the retreat group, we select the simulations that temporarily retreat at least 1 km behind the peak but then recover instead of retreating entirely. These are also synchronized to the onset of retreat (Fig. 4b and d). Comparing the ensemble-mean forcings for both synchronized groups demonstrates that the frontal ablation anomalies that trigger sustained retreat are not systematically larger than those that allow recovery (Fig. 4e). However, they are on average more persistent into the decades following the onset of retreat. Note that in some individual cases visible in panel (c), this persistence is not a continuous excursion but multiple positive anomalies in quick succession.

Figure 4Unstable retreats are driven by persistent anomalies in climate forcing. (a) Unstable retreats from a 1000-member ensemble, synchronized to the time that retreat begins (only 50 members shown for clarity). (b) Synchronized members from the same ensemble that retreated >1 km from the bed peak but then recovered. (c) Synchronized frontal ablation anomalies for the retreat group in (a), along with the group mean (red). (d) As for (c) but for the recovery group. (e) Group-mean frontal ablation anomalies plotted together. On average, the key difference between unstable retreats and recoveries is the persistence of the anomalies rather than the maximum at year 0. (f) Autocorrelation function of the AR-1 noise used to simulate frontal ablation anomalies (green) and autocorrelation of the same AR-1 noise plus a 150-year trend of similar magnitude to the variability. A trend will enhance persistence on the multidecadal time frames that set apart the retreat and recovery groups in (e).

Download

The importance of persistent climate anomalies for triggering sustained glacier retreats is related to the timescales of transient ice dynamics. Consider an initial terminus fluctuation driven by anomalous frontal ablation. If the terminus retreats past the bed peak and into deeper water, discharge will increase due to the strong dependence of ice flux on grounding-line thickness (Schoof, 2007). Independent of the initial forcing, this drives dynamic thinning and further retreat (i.e., the marine-ice-sheet instability mechanism begins). These changes are not instantaneous; ice flow near the terminus evolves on multidecadal timescales (Robel et al., 2018), so retreat is reversible if the climate forcing anomaly recovers before significant changes in the inland ice flow occur (Fig. 4b). However, the longer the terminus persists behind the bed peak, the more interior ice is lost to the increased discharge. At some point, dynamic thinning and retreat will proceed to the point where the terminus cannot recover even if the initial forcing reverses, and thus the marine-ice-sheet instability takes over in driving the retreat.

The groups of simulations in Fig. 4 reflect only one aleatory ensemble, in which the time series of forcing variability are drawn from the same statistical distribution. That is, the average persistence of the variability is prescribed, but individual excursions from the mean vary in duration, purely by chance, as would be expected in natural internal climate variability. However, an external forcing trend systematically adds persistence to any time series (Fig. 4f), with a positive trend increasing the time that positive excursions are above the pre-trend mean. We can anticipate that such extra persistence will affect the probability of rapid glacier retreats, as we explore in the next section.

6.1 Uncertainty in the climate forcing

Our results from large aleatory ensembles show how anthropogenic trends in climate forcing affect the probability of rapid marine-terminating glacier retreat in the presence of bed variations, even when the trend is weak compared to natural climate variability (Fig. 5). The comparison between early- and late-onset trends highlights that a glacier's long-term integration of climate forcing is fundamental to the increased probability of retreat. The early stages of anthropogenic forcing may thus play an important role in an overall attribution assessment. These effects are very clear in the synthetic experiments, in which the difference between ensembles – that is, an anthropogenic climate trend – is simply imposed. However, this trend must ultimately be inferred from observations and models of climate, which is an attribution task of its own. When targeting real glaciers, it will be important to evaluate assumptions about the onset and magnitude of anthropogenic trends built into the model simulations.

What do observations tell us about the onset of climate trends in polar regions? Instrumental records spanning the industrial era are scarce or nonexistent in the vicinity of many glaciers, but a few air-temperature records in Greenland extend over two centuries (e.g., Vinther et al., 2006) and show warming trends since the early 19th century. At a regional scale, the multi-proxy reconstructions of Abram et al. (2016) also show an early (∼1830) onset of surface warming in the Arctic but not the Antarctic. The few long-term ocean records off of Greenland show near-surface (0–40 m) temperature variations roughly correlated with air temperatures (Ribergaard and Buch, 2008), though strong variability obscures the onset of local trends. Additionally, mid-latitude Atlantic waters, a key source of thermal forcing for Greenland's glaciers (e.g., Straneo et al., 2011), have been accumulating heat since the late 19th century (Roemmich et al., 2012). These records do not alone partition natural vs. anthropogenic effects but do provide evidence of early-onset trends.

Comparisons of global climate model simulations with and without anthropogenic greenhouse gas emissions show a significant anthropogenic warming signal beginning in the late 19th century and accelerating in the late 20th century (e.g., Haustein et al., 2017; IPCC, 2021). However, the timing of global-mean temperature changes may not translate directly to the local changes relevant to glaciers. Melt thresholds can add nonlinearity to forcing trends, changing the time of emergence (e.g., Trusel et al., 2018). Transient warming is also slower for ocean variables due to the thermal inertia of the ocean, and localized radiative feedbacks produce further regional variations in surface warming, especially at high latitudes (e.g., Armour et al., 2013). Additionally, the spatiotemporal signature of some forcing agents is relevant at high latitudes. For example, anthropogenic aerosols had an outsized effect in the Arctic and North Atlantic in the 20th century, and several studies have concluded that they effectively canceled the greenhouse-gas forcing in the mid-20th century in the Arctic (e.g., Fyfe et al., 2013; England et al., 2021). This would imply a more step-wise anthropogenic forcing trend, which would be important to account for in attribution experiments. Note that time variations in aerosol forcing also complicate interpretations of Atlantic multidecadal variability (Mann et al., 2020). For glaciers sensitive to this variability, these effects would be important to consider when defining the timescales of stochastic forcing in model experiments.

Despite these challenges, it is still possible to get a sense of the external forcing component using output from coupled climate model ensembles. Single-model large ensembles (in which ensemble members differ only due to internal variability) can be used to estimate externally forced trends at any model grid point (e.g., Deser et al., 2012), and unforced simulations are available for comparison. Statistical methods for estimating the forced and internal components in model output and observations also continue to improve (e.g., Mckinnon et al., 2017; Wills et al., 2020).

Our results suggest that assessing uncertainty in the anthropogenic component of local climate forcing will be very important for understanding the robustness of attribution statements because of the way glacier dynamics integrate progressive forcing. Ultimately, attribution statements are always made with reference to a world without anthropogenic forcing, which only exists in models. In order for an attribution analysis to provide useful insight, the assumptions about this counterfactual world have to be clearly presented and their limits recognized.

6.2 Initial glacier conditions

Uncertainty in the preindustrial glacier state is also important to consider in attribution studies. We have shown that a glacier's initial proximity to a topographic threshold is a key factor in its probability of retreat (Fig. 3), but few terminus records are available before the satellite era (e.g., Goliber et al., 2021). Moreover, the terminus will have fluctuated due to variability in the preindustrial climate (e.g., Fig. 2), so even if early terminus positions are known, discrete observations still leave some uncertainty in the mean state.

It is worth considering limiting cases and what they would imply for attribution. One limit is a mean preindustrial terminus position far enough from a topographic peak that it is virtually impossible for natural terminus variability to breach the threshold. In this case, an observed retreat from the peak could likely be strongly attributed to external forcing. Note that this would also require a large forced response to drive the terminus toward the peak prior to rapid retreat. The other limit is a terminus position perched right at the peak so that natural variability may quickly drive retreat. We expect it is unlikely to encounter this limit because such a glacier would have already retreated due to past variability in a stationary climate. The likely targets for attribution studies – that is, those glaciers that have retreated from bed peaks during the observational era but for which the anthropogenic role is ambiguous – will fall somewhere between these limits. It will be necessary to test the sensitivity of conclusions in an attribution study to the range of preindustrial conditions compatible with available observations. Fortunately, assumptions about the mean state affect the probability of retreat in the same direction for ensembles with and without trends in forcing (unlike assumptions about the forcing, discussed above). That is, initial conditions closer to a bed peak make retreat more likely whether there is a trend in the forcing or not, but the difference between the two scenarios is still a way to assess the anthropogenic role despite uncertainties in initial conditions. This would also apply to other parameters common to the scenarios with and without trends, such as the magnitude of stochastic climate variability. Nevertheless, improved constraints on pre- or early-industrial-era glacier state remain an important goal for providing context for recent changes. To extend observations of glacier state before satellite or aerial campaigns, trimlines, marine sediments, and bathymetry may offer useful constraints (Andresen et al., 2012; Csatho et al., 2008; Smith et al., 2017).

Finally, it is also worth considering what an ice sheet's long memory of Holocene climate variations implies for the stability of preindustrial termini with respect to bed peaks. An important and broad question remains as to why glacier termini should be near topographic thresholds to begin with. On one hand, it may seem unlikely to find any glaciers near thresholds if there is a non-trivial probability that random variability will eventually drive retreat. On the other hand, glacier termini may be transiently persistent at bed peaks even without long-term stability (Robel et al., 2022). As an ice sheet responds to natural climate variations (at all timescales) and its margin evolves over variable topography, one might expect to find most termini near peaks at any given time. Indeed, we observe a similar condition in our model spinup with synthetic random topographies (Fig. 6), although this was not strictly intended to emulate past climate. It is also important to consider here that glaciers affect their landscape via erosion, sedimentation, and lithostatic deformation, and so topography may not be random with respect to the terminus. For example, proglacial sedimentation or the character of over-deepened topography immediately upstream of the terminus might yield clues as to how persistent termini have been, which could help bound the probability of noise-driven retreat. Investigating these geomorphic processes in the context of climate variability could help refine our view of terminus stability in both the preindustrial and modern climate.

6.3 Defining glacier “retreat”

Our numerical experiments focus on the general phenomenology of rapid marine-terminating glacier retreats in the context of climate variability. Accordingly, we assume certain thresholds for retreat in the simulations. However, for analyses of specific glaciers, it will be necessary to define what aspects of the observed change are the target for attribution, and therefore what counts as reproducing the glacier retreat “event” in simulations. For probabilistic attribution, it is generally the case that the more narrowly an event is defined, the harder it is to attribute to external forcing because fewer simulations will resemble the observation in detail (van Oldenborgh et al., 2021). For example, if bed topography permits sustained retreat for many decades once initiated (i.e., lacking stabilizing topographic features), then the cumulative retreat observed so far depends partly on when retreat was triggered.

The varied timing of terminus retreats in Greenland illustrates this issue. Before the widespread retreats of the last few decades, many marine-terminating glaciers also retreated following rapid Arctic warming in the 1920s and 1930s (Bjørk et al., 2012; Vermassen et al., 2019; Andresen et al., 2012). Some, such as Upernavik and Kangerdlussuaq glaciers, have retreated nearly monotonically through deeper bathymetry since then, although with some variation in the rate of retreat (Khan et al., 2014; Vermassen et al., 2019). For attributing early-onset retreats, one might focus on the full observed change, assessing the number of stochastic simulations that reproduced the full magnitude of observed retreat. This would likely select for stochastic anomalies that trigger retreats early in the simulation, reducing the relative importance of long-term trends. Alternatively, one could take a binary approach, focusing on any sustained retreat over some threshold initiated in the simulation period. This would admit later-onset retreats into the analysis, allowing anthropogenic trends in the forced ensemble to have more effect. However, in such a case, the full magnitude of observed retreat would not necessarily be part of the formal attribution assessment, and this would have to be clearly stated. It is reasonable to hypothesize that retreats triggered early in the industrial era were more contingent on natural climate anomalies compared to more recently initiated retreats. Some glaciers in West Greenland retreated abruptly in the early 2000s and have since stabilized on new bed peaks (Catania et al., 2018), offering very discrete changes that would be potentially easier to define as targets for attribution experiments. Future work might compare these recent retreats against early-onset, more continuous retreats.

These issues illustrate possible tradeoffs associated with defining the target glacier retreat event. Once again, attribution is inherently contingent on the assumptions and framing of the analysis. Choices on how to define the observed change will affect the final quantitative elements of the attribution statement (e.g., how much more likely was a given retreat made by anthropogenic changes). Our main point here is that such choices have to be carefully considered and clearly communicated in future attribution studies.