A framework for estimating the anthropogenic part of Antarctica’s sea level contribution in a synthetic setting

The relative contributions of anthropogenic climate change and internal variability in sea level rise from the West Antarctic Ice Sheet are yet to be determined. Even the way to address this question is not yet clear, since these two are linked through ice-ocean feedbacks and probed using ice sheet models with substantial uncertainty. Here we demonstrate how their relative contributions can be assessed by simulating the retreat of a synthetic ice sheet setup using an ice sheet model. Using a Bayesian approach, we construct distributions of sea level rise associated with this retreat. We demonstrate that it is necessary to account for both uncertainties arising from both a poorly-constrained model parameter and stochastic variations in climatic forcing, and our distributions of sea level rise include these two. These sources of uncertainty have only previously been considered in isolation. We identify characteristic eﬀects of climate change on sea level rise distributions in this setup, most notably that climate change increases both the median and the weight in tails of distributions. From these ﬁndings, we construct metrics quantifying the role of climate change on both past and future sea level rise, suggesting that its attribution is possible even for unstable marine ice sheets.


Introduction
The West Antarctic Ice Sheet (WAIS) has undergone dramatic changes over the satellite era, characterized by ice acceleration [1], thinning [2], retreat [3], and ice loss [4].The WAIS currently contributes approximately 10% of global sea level rise (SLR) [5,6] and could add tens of centimeters over the coming decades, possibly dominating by the end of the century [7].However, despite being key symbols of anthropogenic climate change [8,9], Antarctic ice loss, and thus associated SLR contributions, are yet to be formally attributed to anthropogenic climate change [10].
A robust causal relationship between WAIS ice loss and anthropogenic climate change is yet to be established because of strong internal variability in the region's climate as well as ice-ocean feedbacks which perpetuate ice loss [10].There are several lines of evidence highlighting their complex interplay.While WAIS retreat was initiated in the 1940s [11][12][13], after an approximately 10,000-year quiescent period [14], anthropogenic influence on key climatological drivers in the region only became significant in the 1960s [15].This suggests that the trigger for retreat would have occurred even without anthropogenic forcing.Following its initiation, WAIS retreat was likely sustained by ice-ocean feedbacks [16][17][18][19][20][21] (figure 1).Most notably, retreat of this marine ice sheet across a retrograde bed (upward sloping in the flow direction) is associated with increased ice flux across the grounding line (where the ice transitions from sitting on bedrock to a floating ice shelf), which promotes further retreat [22,23] (figure 1).Thus, one possibility is that the ongoing ice loss was triggered naturally in the 1940s and retreat is dominated by self-perpetuating feedbacks, playing out on the long timescales on which ice-sheets evolve [11,13,15,24].However, this retreat cannot be purely self-sustaining, independent of external forcing, because ice discharge remains responsive to ocean variability [25][26][27].This picture is further complicated by a proposed centennial scale warming of the Amundsen Sea [24,28], which is partly attributed to anthropogenic changes in large-scale climate systems [15,[28][29][30].While all of these processes may contribute to the ongoing ice loss, the relative contributions of a historical trigger, ice-ocean feedbacks, and changes in climatic forcing are still unknown.
Determining the role of anthropogenic climate change in SLR from the WAIS is important for providing causal evidence to support recourse for the myriad social [e.g.31], economic [e.g.32], and ecological [e.g.33] impacts of SLR, which are borne primarily by poorer and low-lying island nations [34].This is particularly pertinent in light of the recent outcomes of the COP27 conference, in which a loss and damage fund was established to compensate countries for the harm inflicted by anthropogenic climate change.In addition, attribution (or lack thereof) has implications for the future of the WAIS: if the observed ice loss is due solely to internal variability and ice-ocean feedbacks, SLR is likely already committed and irreversible; whereas, a significant anthropogenic component might suggest that ongoing contributions strongly depend on future greenhouse gas emissions.
Despite the importance of this question, an outline of how to address it is not yet clear.Progress has been made towards such by ref. [35], who considered how ice sheet retreat from a local topographic high under variable forcing may be attributed, using a one dimensional ice sheet model.Using a set retreat threshold as the event to be detected, they showed that while an observation of large retreat under a single realization of stochastic climatic forcing does not necessarily indicate that anthropogenic climate change was present in the forcing (figure 1), even modest anthropogenic trends in forcing make retreat more likely when averaged over multiple realizations.They conclude that a probabilistic approach, with multiple realizations of forcing, must be taken if robust attribution statements are to be made.Additionally, they showed that model parameter choices have a large impact on the likelihood of retreat, and thus the attribution statement; this suggests that multiple model parameters should be considered simultaneously in the attribution assessment, particularly when these are poorly constrained.
Here, we consider how the anthropogenic component of SLR contributions from WAIS may be determined, which uses a Bayesian approach integrating multiple realizations of forcing; we build upon ref. [35] in two main ways: firstly, we consider SLR contributions, rather than retreat, as the metric to be attributed.By using SLR as the attribution metric, we are able to quantify the role of anthropogenic climate change for observed SLR within any interval, rather than only exceedance of a single, predefined retreat threshold.This alleviates the common event definition problem which commonly impacts attribution studies [36].Secondly, we explicitly account for the role of variable model parameters in the attribution assessment.Bayesian approaches naturally permit the joint probability density of multiple model parameters, which may be poorly constrained in general, to be represented within a projection of SLR [37].This avoids the need to specify the precise values of model parameters at the outset, which yield very different attribution results depending on the particular choice of parameters in the framework of ref. [35].
More specifically, we consider parameter variability in the parametrization of ice shelf basal melting, which is calibrated by comparing the resulting ice shelf basal melt rate fields with output from a more detailed ocean model.This procedure represents a hybrid approach that sits between parametrizations of basal melting and coupled ice-ocean models, and calibrates melting directly, rather than only indirectly via its effect on ice flow.We demonstrate how the anthropogenic component of SLR contributions may be determined by considering the retreat of a synthetic marineterminating ice sheet, which is highly susceptible to ice-ocean feedbacks and subject to forcing with strong internal variability, the characteristic features that are thought to obscure signals of anthropogenic climate change in SLR contributions from the WAIS.We demonstrate how uncertainties associated with poorly constrained model parameters interact with uncertainties associated with stochastic climate forcing, identifying that it is necessary to consider both, a feature that is lacking in current SLR projections.To the best of our knowledge, this is the first time such uncertainties have been considered simultaneously in an ice sheet modelling exercise.
We explicitly construct distributions of SLR which simultaneously account for parametric uncertainty (that arising from poorly constrained model parameters) and aleatory uncertainty (that arising from an ice sheet's variable response to different realizations of stochastic forcing).These distributions also reveal characteristic signatures of anthropogenic forcing on distributions of SLR from marine ice sheets, which we describe, and allow us to construct a metric describing the influence of anthropogenic forcing on SLR in this system.We conclude that even in highly unstable marine ice sheets, the impact of anthropogenic forcing is detectible in principle, given sufficiently large simulation ensembles as well as a full treatment of model parameter uncertainty.We finish with a brief discussion of the challenges associated with determining the role of anthropogenic forcing on SLR contributions from the WAIS, which are avoided in our use of a synthetic configuration.These include uncertainty in other model parameters, uncertainty in the initial state, and uncertainties in climatic forcing.

Interactions between aleatory and parametric uncertainties in sea level rise projections
We adopt a Bayesian approach in which parametric and aleatory uncertainties are simultaneously accounted for.As is standard, parametric uncertainty is accounted for by performing multiple simulations with different model parameters spanning the parameter space (for each realization of forcing), with the resulting SLR contributions weighted according to the level of agreement between a simulated quantity and its ground truth [e.g.[38][39][40][41][42][43].It is straightforward to incorporate aleatory uncertainty into such an approach (see methods) by placing no preference on the specific realization of forcing.Although accounting for parametric uncertainty in this way is now standard, no study has yet probed the interaction between parametric and aleatory forcing uncertainties, primarily because of the computational expense of doing so [40], since multiple simulations with different model parameters must be run for each additional realization of forcing.
To illustrate the approach, we focus on parametric uncertainty arising from the use of a parametrisation of ice shelf basal melting.Parameterisations of basal melting are often used instead of coupled-ice ocean models to reduce computational expense (in coupled ice-ocean models, the ocean component typically represents the vast majority of the expense [44]).Coupled ice-ocean models remain computationally intractable for the large ensembles of simulations [44] required to incorporate both aleatory and parametric uncertainty.However, parameterisations of melting neglect processes that have been shown to be important in determining basal melting [e.g.16,45,46], and simulations employing parameterisations have been shown to yield basal melt rates which result in poor skill at reproducing observed grounding line retreat [47] and ice loss [48][49][50], compared to coupled ice-ocean models.Our approach can be considered a hybrid between a parametrisation of melting and a coupled ice-ocean model: we use a parameterisation of basal melting for computational efficiency and adopt a Bayesian approach to the model parameters within: simulations are weighted by comparing their predictions of basal melt rates with those from an offline ocean model at different snapshot times throughout a simulation (methods); the ocean model thus plays a role analogous to a ground-truth in a traditional Bayesian update, i.e. it is the information assimilated into the model.It should be noted that this is a slightly different philosophy to a typical Bayesian update in ice sheet modelling, in which agreement with satellite observations, rather than with results of more detailed models, are typically used to update probabilities.We employ a common melt rate parameterisation in which melting has a quadratic dependence on ocean temperature and scales linearly with a dimensionless parameter M , which is independent of the ocean temperature (methods).The melt rate calibration procedure is only capable of calibrating the melting aspects of the flow model; other parameters, such as those related to basal sliding and ice viscosity, which are important in determining ice flow (and thus SLR) are not calibrated.Other studies [e.g.[38][39][40][41][42][43] have established procedures for calibrating many such aspects of ice-sheet models using observational data; the novelty of our calibration method is that it permits precise calibration of basal melt rates, which have, to the best of our knowledge, only previously been indirectly calibrated via the effect of melting on ice flow.In practice, all parameters with an important effect on ice dynamics should be calibrated (see 'Discussion'), but our use of a generic ice sheet configuration (described below) allows us to neglect them, and focus on errors arising purely from poor melt rate parametrisation skill.
Our example configuration features a prominent seabed ridge (figure 2a) on which the ice shelf is stably grounded (figure 2b) during an initialization stage with temporally constant ocean forcing, corresponding to typical conditions in the Amundsen Sea offshore of the WAIS (methods).This grounding line position, located at a topographic high, is reminiscent of the WAIS configuration prior to the 1940s [11] and renders the system highly sensitive to ice-ocean feedbacks once grounding line retreat has been initiated [49].We consider evolution from this steady state under variable ocean forcing, which is imposed by varying the depth of the pycnocline in the ambient ocean conditions (figure 2c-d).The ocean forcing includes a stochastic internal variability component, which mimics the observed amplitude [51,52] and persistence [35] of internal variability in ocean conditions in the Amundsen Sea on decadal and interdecadal timescales.Superimposed on this forcing is either an anthropogenic trenda 100 m/century linear shallowing of the pycnocline, illustrating a plausible historical anthropogenically driven trend in Amundsen Sea conditions [28,53] -or no trend, representing the counterfactual scenario in which no anthropogenic climate change has taken place (figure 2g).For both of these scenarios (referred to as anthropogenic and counterfactual, respectively), we perform simulations with 40 independent realizations of forcing (the realizations in each of the two ensembles are also independent).Although accumulation rates also feature notable interdecadal internal variability, and are projected to display an anthropogenic trend in the future [7], this variability is smaller than in the ocean forcing for WAIS.In addition, changes in melting, rather than accumulation, are understood to have been the dominant driver of recent WAIS retreat [54,55], and having multiple forcings, each with a unique anthropogenic trend complicates the attribution task somewhat.
For each realization of forcing, we perform simulations sampling the parameter space of M .Requiring that the ice shelf remains stably grounded at the ridge crest during the initialization phase, and retreats under forcing corresponding to the warmest observed conditions applied constantly, restricts us to considering the range 0.5 < M < 1.5 (methods); we sample this range by taking M ∈ {0.5, 0.75, 1.0, 1.25, 1.5}.Thus, the total number of simulations is 400 (2 ensembles × 40 members × 5 M values).
Examining the response to a single illustrative realization of forcing (figure 2e), for different melt parameters M , highlights the interplay between stochastic forcing and parameter variability, elucidating the inextricable relationship between aleatory and parametric uncertainty.On the centennial scale, this realization of forcing features two prominent warm periods (figure 2e).During the first warm period (between approximately t = 20 and t = 40 years), retreat is triggered in those simulations with the largest values of M (M = 1, 1.25, 1.5; figure 2f).These retreats are initiated towards the end of the first warm period (figure 2f), when the time-integrated melt anomaly has caused enough ice shelf thinning to reduce ice shelf buttressing to the level at which retreat is initiated.Accordingly, retreat is initiated soonest in the simulation with the largest melt parameter M (figure 2f), which has the highest melt rates and thus accumulates the time-integrated melt anomaly most rapidly.Once initiated, retreat proceeds at a rate approximately independent of forcing (figure 2f), suggesting that, once triggered, retreat is set primarily by ice-ocean feedbacks, although it remains weakly responsive to changes in forcing.Simulations with smaller M (lower melting) remain grounded at the ridge crest during the first warm period.Retreat is initiated in the M = 0.75 simulation during the second warm period, again towards the end of the period.A simulation with the same realization of forcing but with the anthropogenic trend removed, and M = 0.75, does not retreat during this period (note that this simulation is outside the ensemble structure outlined above, for which anthropogenic and counterfactual ensembles are independent): the integrated melt anomaly required to initiate retreat is achieved more easily during a given time period if there is an anthropogenic trend in the forcing, than if not.
Under the same realization of forcing, SLR may be highly non-linear in M (figure 2h).For example, SLR contributions in the highlighted curve in figure 2h increase by 1800% (from 0.15 mm to 2.91 mm after 100 years) when the melt rate parameter is increased from M = 1 to M = 1.25.This strong sensitivity demonstrates the necessity of considering a range of parameter values in determining SLR contributions, particularly when the system is susceptible to ice-ocean feedbacks, or so-called tipping points may be passed.Furthermore, there are simulations in the anthropogenic ensemble which yield lower SLR than simulations in the counterfactual ensemble (figure 2h), and this behavior is strongly influenced by the value of M .Thus, an observation of high SLR under a single realization of forcing is not necessarily an indicator of strong anthropogenic influence (figure 1).Taken together, these results -a strong sensitivity to the parameter M and to the specific realization of forcing -demonstrate that parametric and aleatory uncertainty must be simultaneously accounted for in SLR distributions, and thus any framework attempting to determine the role of anthropogenic trends in forcing in them.
The non-linearity of SLR in M also demonstrates how single-point parameter calibration (where the set of model parameters are specified based on agreement with a single metric, say the total melt flux out of an ice shelf cavity) may be problematic.Such single-point calibrations are often applied when tuning melt rate parameterisations [e.g.50, 56, 57].In the example presented here, the mean melt rate at the start of the simulation (at the end of the initialization stage, which is performed separately for different values of M ) is only weakly dependent on the melt rate parameter M (supplementary figure 3d), owing to a feedback between melting and ice geometry (methods).As a result, a small change in the single target metric to be matched would result in a large change in the selected value of M (supplementary figure 3d), which would ultimately result in a large change in the simulated SLR at the end of the simulation (figure 2h).In other cases where the target metric is more sensitive to parameters, a small change in the target metric would be expected to result in a small change in the selected parameter, but this may also ultimately result in a large change in the simulated SLR at the end of the simulation, owing to the non-linearity of SLR in M .

Influence of anthropogenic forcing on sea level rise probability distributions
Applying the Bayesian melt rate calibration procedure (methods), yields, for each time in each simulation, a distribution of SLR associated with the particular realization of forcing applied (supplementary figure 4).Then, by marginalizing over the realizations of forcing (methods), we obtain calibrated probability distributions of SLR for both anthropogenic and counterfactual ensembles, at each time (figure 3a).
The time evolution of both ensembles display qualitatively similar behavior.The evolution of the distributions can be categorized into two temporal parts: 'tail emergence' and 'shift towards tails' (figure 3c).At early times, the distributions are symmetric (figure 3a), with low skewness (figure 3c) reflecting retreat having not been triggered in any simulations.As retreat begins to be triggered in individual simulations, the 'tail emergence' period begins: a tail emerges (skewness increases, figure 3c), supported by increasing SLR contributions from those already retreating simulations, and kurtosis increases (figure 3d), indicating that the relative weight in the tails is reducing (kurtosis quantifies the proportion of weight placed in the tails, with low kurtosis corresponding to heavy tails).The timescale on which the tails emerge depends on the forcing (see below).Median SLR remains small in the tail emergence period (figure 3b).
As retreat is triggered in an increasing number of ensemble members, weight begins to shift to the tails; the 'shift towards tails' period begins when skewness and kurtosis reach a maximum (figure 3c-d).Beyond this maximum, weight moves towards the tails (kurtosis reduces, figure 3d) and, in response to this, the median increases (figure 3a), continuing to the end of the simulation.(The median is a more appropriate metric than the mean given the skewed data.)Both medians display a non-linear evolution, reflecting non-linear SLR contributions in individual simulations once retreat has been initiated.Although the precise details of the evolution of the distributions depends on both the system and the forcing (see below), we expect that this qualitative behavior is generic in marine ice sheets with tipping points under high variability stochastic forcing.
Despite these qualitative similarities between the anthropogenic and counterfactual distributions, there are clear quantitative differences, which highlight the importance of the anthropogenic trend in forcing.Firstly, the tail emerges sooner in the anthropogenic ensemble (figure 3c), because retreats are initiated sooner when a trend in forcing is imposed (supplementary figure 1).This is despite the anthropogenic additional forcing being zero at the start of the simulation (figure 2g), highlighting the role played by increases in forcing during the time period in which the destabilizing integrated melt anomaly is accumulating: if forcing did not change over this period (or, if the changes did not matter), the first retreats would take place at approximately the same time in both ensembles.This is consistent with [15], who suggest that the current retreat of WAIS was triggered naturally in the 1940s, but may have subsequently failed to recover due to increasing influence of anthropogenic forcing towards the start of the 1960s.Secondly, the maximum skewness is lower, and achieved sooner, in the anthropogenic case (figure 3c).In a given time period, retreat is triggered in a greater proportion of simulations in the anthropogenic ensemble than in the counterfactual ensemble (supplementary figure 1), resulting in probability distributions shifting more quickly towards the heavy-tailed regime.This difference in retreat rate triggering is because, as time proceeds, melt anomalies under anthropogenic forcing become increasingly large, so a shorter positive anomaly duration is required to initiate retreat.More specifically, with a linear anthropogenic trend, the time-integrated melt anomaly scales with the square of time, which rapidly outweighs any time-integrated negative internal component: the system is more vulnerable to long-lasting trends in melting than to short term variability.Finally, and most importantly, on the centennial timescale, both the median is larger, and the kurtosis smaller, in the anthropogenic ensemble than in the counterfactual ensemble; i.e. not only does anthropogenic forcing increase the median of the distribution, it also results in greater weight in the tails: extreme events, with high SLR contributions, have relatively large probabilities in the anthropogenic ensemble.This emphasizes the need to consider the shape, as well as the spread (e.g. the variance), when communicating how emissions pathways affect future SLR scenarios with policymakers.
Figures 3b-d also indicate how summary statistics differ between the calibrated and uncalibrated distributions, with the latter obtained by setting the posterior probability equal to the prior probability (methods), i.e. all values of M are weighted equally.In both ensembles, parametric calibration of M has an important effect on the median, evidencing the need to apply parametric calibration in projections of SLR from ice sheets.Reduced uncertainty in projections is often (perhaps implicitly) cited as a key benefit of parametric calibration [e.g.38,40]; whilst our simulations provide evidence to support this, displaying increased kurtosis (reduced weight in the tails; figure 3d) in the calibrated case, there remain large uncertainties in calibrated distributions (figure 3a).This suggests that aleatory uncertainty is an unavoidably large part of uncertainty in projections of SLR from ice sheets, particularly those highly susceptible to ice-ocean feedbacks, and cannot be neglected: parametric calibration alone is not sufficient, and there is irreducible uncertainty in SLR from marine ice sheets.

Quantifying signals of anthropogenic trends in forcing
The role of anthropogenic climate change in individual weather events is often framed as an anthropogenic enhancement [36]: how many times more (or less) likely was the event made by anthropogenic climate change?Having constructed distributions of SLR in both anthropogenic and counterfactual cases, the ratio of these -the anthropogenic enhancement ratio (AER) -naturally emerges as a metric to quantify how many times more likely an observed SLR was made by the presence of an anthropogenic trend in forcing, and go beyond the qualitative comparisons of the previous section.An AER of 2, for example, indicates that anthropogenic forcing made a given SLR contribution 100% more likely (or, equivalently, twice as likely).The AER for our ensembles is shown in figure 4a, where values along each line of constant time represent the ratio between the anthropogenic and counterfactual probability distributions (as shown for specific times in figure 3a).Note that, because the AER can be constructed for any time throughout the simulation, past and future SLR are equally applicable -the present has no special status.Therefore, attribution statements may be made for either past or future SLR contributions (or both).
There is a band in which the AER is infinite, which is caused by the tails of the anthropogenic distribution extending to higher SLR values than those in the counterfactual distribution (figure 4a).An observation of SLR in this band would have been impossible without anthropogenic climate change-no counterfactual simulations produce this value.The band spreads out in time from an area close to the origin (recall that the tail of the anthropogenic distribution emerges soon after the start of the simulation) at a rate that is set by the retreat of the individual simulation with the highest SLR.
The AER is generally increasing in SLR, indicating that a higher SLR over many realizations of forcing is a stronger indicator of anthropogenic climate change.This demonstrates the importance, and value, of accounting for aleatory uncertainty: under a single realization of forcing, higher SLR does not necessarily indicate a strong influence of anthropogenic climate change (figure 1), but it does when appropriately averaged over many realizations of forcing.This also highlights the shift from a binary yes-no question, to a probabilistic approach, that necessarily takes place when accounting for aleatory uncertainty [35].The AER has a slightly banded structure (figure 4a), which results from the finite size of our ensembles (in the limit of infinite ensemble members, the proportion of retreats initiated would be smooth, whereas because of the finite size of our ensemble, the proportion of retreats initiated oscillates around the trend in this quantity, with periods when relatively more, and periods when relatively few, retreats are initiated compared to the background trend, see figure 1).While we expect that the banding would disappear as the number of realizations of forcing goes to infinity, we note that increasing this number is particularly computationally expensive when accounting for aleatory and parametric uncertainty simultaneously.
In practice, observed SLR follows a single trajectory through this AER space, such as the selected simulations shown in figure 4b-d, in which retreat is triggered after approximately 20, 40, and 60 years, respectively (figure 4a).Their values are indicative of the clear signal of anthropogenic climate change: at the end of the simulation, AER is approximately 2.5, 3.9 and 2.2, respectively, corresponding to increases in probability of 150%, 290%, and 120%.Once retreat has been triggered, the AER remains fairly constant.It is worth noting that these values are perhaps modest compared to glaciological attribution studies applied to mountain glaciers [e.g.58,59].This is a direct consequence of our choice of setup: we consider a scenario in which internal variability is relatively large compared to the anthropogenic trend and is highly susceptible to ice-ocean feedbacks (and these selected trajectories don't enter the tail band, for which AER → ∞).
From a policy perspective, a third useful question, beyond how to address and how to quantify the role of anthropogenic trends in forcing, is: what is the uncertainty in this quantification?Having constructed distributions associated with each realization of forcing (which the distributions shown in figure 3a are the mean over), such uncertainties can be probed.To do so, we bootstrap values of the distributions from individual realizations of forcing to determine a confidence interval (methods)a measure of the likely spread in AER -around our central estimates (figure 4b-d).Uncertainty in AER is generally smaller along contours corresponding to later retreat (figure 4b-d).This is commensurate with relatively few simulation trajectories entering the region in and around the tail band, leading to increased uncertainty: although the central estimate of anthropogenic enhancement is itself largest in the tails, that is where the uncertainty in the value is greatest.We expect that this error bound would reduce with increasing numbers of realizations of forcing.Thus, we expect that real world attribution studies will have to grapple with the limitation that increasing ensemble size is required to reduce uncertainty in the role of anthropogenic forcing, but to do so requires substantial additional computational resources.

Discussion
The example presented here provides a path towards assessing the role of anthropogenic climate change in SLR contributions from the West Antarctic Ice Sheet, including both quantifying the strength of the anthropogenic signal and its uncertainty.Our use of a Bayesian framework allows us to treat parametric uncertainty within attribution assessments and avoids the need to specify a single event to be detected.By abstracting and considering a generic ice sheet, we are able to focus on errors in melting, with the hope that the melt calibration approach may help to bridge the considerable gap in fidelity to observations between parameterisations of melting and coupled ice-ocean simulations.
Determining the precise influence of anthropogenic climate change on SLR contributions from the WAIS requires simulations to be performed using geometries and parameters that represent real world conditions.Here we identify three key classes of problems which must be overcome in doing so: computational challenges, initial state challenges, and challenges arising from uncertainty in climatic forcing.Computational challenges arise from the large number of simulations required to appropriately account for parametric and aleatory uncertainty.In considering a generic marine ice sheet, we are able to neglect uncertainty arising from model parameters governing basal sliding and ice viscosity, as well as processes such as damage [e.g 19, 60], calving [e.g.46,61,62], sliding law [e.g.63], and ice rheology [e.g.64] which might obscure (or amplify) long-term climatic trends in the forcing, but should be included in assessments of SLR and thus its attribution to anthropogenic climate change.Additional parametric uncertainties can be succinctly integrated into the Bayesian approach taken here [37], and should be calibrated with observations.The computational challenge is particularly pertinent given that a high spatial resolution must be used to ensure correct representation of ice sheet key processes [e.g.65].In addition, the effect of parameters which control the strength of Bayesian updates must be explored; although we find that varying these parameters within reasonable ranges does not qualitatively change the results (methods), they may influence the precise values of anthropogenic enhancement.It should also be noted that, ideally, multiple different ice sheet models should be used in order to assess structural uncertainties arising from those processes not represented in some ice sheet models [37], further adding to the computational challenge.
Determining the initial state -the configuration of the ice sheet prior to the era of anthropogenic influence -also represents an crucial challenge.Projections of ice sheet evolution are sensitive to their initial states, similar to numerical weather forecasts [66], but relatively little is known about the configuration of the WAIS prior to the satellite record beyond broad bounds on grounding line locations [11].One particular challenge in this regard is determining the ice front position, which typically remains fixed in ice sheet models, but may have a strong impact on ice shelf buttressing and thus retreat potential.Additionally, ice sheet memory of Holocene conditions must be considered: here, we have assumed that the ice sheet is in steady state at the onset of a trend in forcing; in practice, however, there is evidence of a slow retreat of the WAIS over the Holocene [14].Given the long timescales on which ice sheets fully respond to changes in forcing, knowledge of this state may be retained by the ice sheet, and thus affect the likelihood of retreat.
Finally, challenges associated with uncertainties in climatic forcing must be overcome.Here, we assumed that the anthropogenic trend is known and well characteristed, but in practice this must itself be inferred from observations and models of climate, representing an attribution challenge in itself.For WAIS, this is complicated by the compound drivers of changes: ocean warming drives retreat, but trends in ocean warming are primarily driven by trends in winds [24].Additionally, anthropogenic trends in accumulation, not considered in this study, must be considered simultaneously with trends in ocean warming; in the future, trends in accumulation are expected to partly offset ice loss from WAIS [7], potentially obscuring the anthropogenic signal.
The work presented here can be considered as a framework for producing calibrated distributions of SLR, in addition to their application to attribution statements.We have demonstrated that both aleatory and parametric uncertainty are important components of ice sheet SLR projections, and suggest that future assessments of SLR from ice sheets must account for these sources of uncertainty.As we have shown, parametric calibration reduces uncertainty, but the susceptibility to ice-ocean feedbacks renders broad distributions inevitable [67]: much like other aspects of the climate system [68], ice sheets have irreducible uncertainty.The glaciological community must become more comfortable with these fundamental aspects of uncertainty and appropriately communicate them to policy-makers and stakeholders.
By constructing calibrated distributions of SLR contributions, we showed that anthropogenic climate change increases both the median of distributions, and the relative weight of their tails: much like many other weather events [69], even modest anthropogenic climate change can make extreme scenarios many times more likely.
Using these distributions, we constructed a metric to quantify the role of anthropogenic forcing, concluding that even in highly unstable marine ice sheets, the impact of anthropogenic forcing is detectable in principle, given sufficiently large simulation ensembles forced by profiles with and without an anthropogenic trend, as well as a full treatment of model parameter uncertainty.In other words, attribution studies are tractable for the WAIS.The implications of attributing ice loss from the WAIS, both for the harms caused by SLR, and the implications for the future of the WAIS, provide strong motivation to pursue such studies.

Sea level rise contributions accounting for parametric and aleatory uncertainty
For a given trend in forcing, denoted F, (i.e. after specifying whether the trend is anthropogenic or counterfactual), the probability of a given SLR, ∆SLR, accounting for aleatory and parametric uncertainty may be expressed as [37] P (∆SLR|F, Here, N is the space of model parameters, n is the total number of realizations of forcing, R i is the specific realization of forcing (with i a dummy index), and I 0 represents the initial conditions.The expression (1) follows from a first-principles probabilistic expression of SLR [37], after assuming that each specific realization of forcing has equal probability, P (R i ) = 1/n, and that the initial state I 0 is known.For our specific application of (1), N is the space of melt rate parameters, 0.5 < M < 1.5.Note that the expression (1) does not include any account of model structural uncertainty, which arises from the approximations that ice sheet models make, as well as their incomplete representation or omission of physical processes [37].Such uncertainties can only be accurately probed by performing the same numerical experiments with an ensemble of different ice sheet models, typically in a model intercomparison exercise [e.g.70] and is therefore beyond the scope of this work.(It should be noted that the WAVI ice sheet model used herein demonstrates good agreement with other state-of-the-art ice sheet models in the most recent ice sheet model intercomparison exercise [70].)Note that constructing distributions of SLR using the calibration procedure outlined below requires values of SLR to be known for all parameter values, but simulations provide only a finite amount of observations.Here, we obtain SLR as a function of M by linearly interpolating between individual M (see figure 4e).

Melt rate calibration
The calibration of model parameters M enters distributions of SLR through the probability P (M |R i , F, I 0 ), which appears in (1) (here we use the specific parameter name M , rather than the generic name N ).Following a standard Bayesian approach, we assume a prior distribution on the parameters M (with hyperparameter µ), which is then updated as new information is assimilated through the likelihood.In our case, this assimilated information is melt rates from an offline ocean model (see below); denoting this information by O, Bayes' rule states that The first term in the numerator on the right-hand side of (2) represents a likelihood function, describing how the prior distribution (second term in the numerator on the right-hand side) is updated to assimilate ocean model results.The prior distribution describes the state of belief in model parameters N prior to comparison with the ocean model.The left-hand side of (2) represents the posterior distribution -the distribution of parameters M following assimilation of ocean model information.The denominator of the right-hand side of (2) simply acts to normalize the probability distribution.
Here, we assume a Gaussian prior, which maximizes the relative entropy when only estimates of the prior mean µ and standard deviation σ P are available [71,72]: Here α is a normalization constant, which ensures that the distribution (3) integrates to unity when initialization bounds on M are imposed (see 'Ice Sheet Model Initialization' below).σ P can be thought of as describing the strength of confidence in the initial estimate of M , which is centered about the hyperparameter µ: a low (high, respectively) σ P corresponds to high (low) confidence that the hyperparameter µ represents the true value of M .In the results contained herein, we use µ = 1.25, based on agreement in the mean melt rate after the initialization stage (in this case, a mean melt rate of 23 m year -1 , which can be thought of as an arbitrary piece of prior information).We use σ P = 0.2 (supplementary figure 4), representing somewhat weak confidence that the value M = µ represents the true value of M .Supplementary figure 6b shows a plot of the Gaussian prior (3) as a function of M for different values of σ P with µ = 1.25.
To determine the likelihood P (O|M, µ), we first specify calibration timeslices τ = {τ 1 , . . .τ n } and, for each timeslice, run the ocean model in the geometry set by the ice-only model.After doing so, we have two melt rate fields, from the parameterisation of melting and from the ocean model, respectively, and for each timeslice k = 1, . . ., n. (Note that the ocean model depends on the melt rate parameter M via the ice-shelf cavity geometry.)A melt error functional D j is determined by comparing these two fields.The particular choice of the form of the D j is subjective, reflecting how melting should be penalized.Here, we take D j to be the mean absolute error in the two melt fields on grid cells below 500 m depth.This reflects the fact that deep areas, typically close to grounding lines, have disproportionately large impacts on the dynamics of the grounded ice [73][74][75].
From the timeslice errors D j , we determine an average error D = (1/n) n j=1 D j .The likelihood is then determined from an exponential error model, Here σ L is a melt error covariance, which describes how harshly errors in the melt rate from the parameterisation are penalized (with respect to the ocean model): for low σ L , errors are penalized more harshly, whereas for high σ L , errors are penalized less harshly.In the limit σ L → ∞, each parameter value M is assigned equal weight, and the posterior distribution is identical to the prior (supplementary figure 4f).Supplementary figure 6b shows a plot of the exponential error model ( 6) as a function of D for different values of σ L .In the results presented here, we use σ L = 10 m/year.In general, the σ L should be on the same order of magnitude as errors in melting; in our simulations, melt errors are typically on the order of 10s of meters per year (supplementary figure 4a).
To assimilate timeslice errors into the Bayesian update, we require P (O|M, µ) as a function of M , but the simulation only provides sparse points (figure 4b).To overcome this, we interpolate between the data points using a smoothing spline fit, via the fit function in MATLAB.
Supplementary figure 6c shows the AER as a function of SLR and time (i.e. as in figure 4 of the main text) for different values of the prior parameter σ P and melt error covariance σ L within reasonable ranges.We see that, while varying these parameters adjusts the precise value of the AER, the overall picture -that higher observed SLR are concomitant with stronger anthropogenic influence -remains.The small exception to this is for large σ P and small σ L , for which the anthropogenic signal is most obscured (see below) and a band of AER < 1 emerges close to the tail.This is a finite size effect, and would disappear in the limit of a large number of simulations, emphasizing the need for large ensembles of simulations.
For smaller σ L , errors in melting are penalized more harshly; in this study, smaller σ L tends to shift weight towards smaller M , which typically display smaller errors in melting (see supplementary figure 4a-b for an example from one realization of forcing).Simulations using a smaller value of M require a larger time-integrated forcing anomaly to achieve the same integrated melt anomaly required to initiate retreat.Simulations in which this is achieved in the anthropogenic case and not in the counterfactual case, tend, therefore, to be observed later on average, when the ensemble mean difference in forcing is greater.Thus, for a given time, the ratio of ensemble members which have retreated in the anthropogenic ensemble to those which have retreated in the counterfactual ensemble is closer to unity for smaller M , leading to a dampened anthropogenic effect.Conversely, smaller σ P shifts weight towards M = µ = 1.25 (in this case), which is at the higher end of the M range considered here, enhancing the anthropogenic effect.

Details of ice sheet configuration
The setup of the generic marine ice sheet configuration is very similar to that of [49], who interrogated how ice-ocean feedbacks perpetuate retreat of an ice sheet from a seabed ridge using a coupled ice-ocean model under constant forcing scenarios.In this setup, the bathymetry (figure 2a) can be expressed as the sum of along-flow and cross-flow components: where B y (y) = − 500 + 600 sin π 2 + πy 5 × 10 4 m.
Here, x and y are co-ordinates in the along-and cross-flow directions, respectively (the ridge is aligned along the cross-flow direction, see figure 2a).The cross-flow bathymetry contribution, B y (y), corresponds to a symmetric valley-like configuration, whose margins are located 500 m below sea level and whose center is 1100 m below sea level; the along-flow bathymetry contribution, B x (x), corresponds to a Gaussian ridge with height 400 m and lengthscale σ b = 1.1 × 10 4 m, which is superimposed on the valley at a position centered on x = 265 km.Following [49], ice rheology is described by Glen's law with flow exponent n = 3.A constant rate factor A = 2.94 × 10 −9 a −1 kPa −3 is applied everywhere, except for within 5 km of the ice margins (i.e. for y < −20 km and y > 20 km), where the rate factor is set to A = 5.04 × 10 −9 a −1 kPa −3 ; this is to mimic the narrow, low viscosity, shear margins which are characteristic of WAIS outlet glaciers, particularly Pine Island Glacier [76].The sliding coefficient is set to 20 m a -1 kPa -1 everywhere.Surface accumulation varies linearly from 15 m a −1 at the ice divide (x = 0 km) to 1 m a −1 at x = 150 km and is set to a constant value of 1 m a −1 between x = 150 km and the ice front (x = 300) km.The resulting total surface accumulation, 67.5 Gt a −1 , closely matches observations [77], while the spatial pattern respects reduced accumulation with reducing altitude.

WAVI ice sheet model
SLR contributions are determined from simulations using the Wavelet-based Adaptivegrid Vertically-integrated Ice-sheet model (WAVI) [72,78], a finite volume ice sheet model including a treatment of both membrane and simplified vertical shear stresses [79].WAVI uses a regular solution grid (here 1 km in both directions), which is refined dynamically during the solution procedure to facilitate solution speed and accuracy.WAVI assumes a fixed ice front position, which is set to x = 300 km (this is equivalent to prescribing a calving law that the calving flux is equal to the normal ice velocity at the ice front).

Melt rate parameterisation
Melting in the ice sheet model is parameterized according to a quadratic temperature law [80], (10) Here, M is a (variable) dimensionless melt rate parameter, T a is the ambient temperature far from the ice shelf base (see below), T f is the local freezing temperature and Γ = 0.56 m yr -1 • C -2 plays the role of an exchange coefficient between temperature and melt rate.(Using the nomenclature of [50,81], Γ = γ T [ρ w c p /(ρ i L)] 2 , where γ T is an exchange velocity, ρ w is water density, ρ i is the ice density, c p is the specific heat capacity of water, L is the latent heat of fusion).The formulation (10) essentially encodes two mechanisms which strongly affect ice shelf basal melting: (1) ice shelf melting is governed by the turbulent heat flux from the ocean to the ice, which varies like the product of ocean temperature and velocity; (2) ocean velocity increases with the local thermal forcing (T a − T f ) as meltwater is released, increasing the buoyancy forcing and thus circulation strength.This parameterisation has been used in numerous ice sheet modelling studies [see 44, and references therein], including the latest ISMIP assessments [81].
As is standard, we assume that the local freezing point depends linearly on pressure and salinity, T f = λ 1 S a + λ 2 + λ 3 z b , where λ 1 = −5.73× 10 −2 • C is the liquidus salinity slope, λ 2 = 8.32 × 10 −2 • C is the liquidus intercept, λ 3 = 7.61 × 10 −4 • C m -1 is the liquidus depth slope, S a the ambient salinity (see below), and z b is the depth of the ice shelf base.
We take a layered structure for the ambient temperature and salinity (figure 2cd), parameterized solely via the depth of the pycnocline centre, P c (which is in general time-dependent), and the pycnocline half-width w: The profiles (11) and ( 12) are piecewise linear functions of depth (figure 2b): they are constant in both an upper (temperature −1 • C, salinity 34 PSU, corresponding to Winter Water) and lower layer (temperature 1.2 • C, salinity 34.6 PSU, corresponding to Circumpolar Deep Water), which are separated by a pycnocline of 2w m thickness, across which the temperature and salinity vary linearly.These piecewise linear profiles are approximations to typical conditions in the Amundsen Sea [26,52].Here, we take w = 200 m, corresponding to a pycnocline width of 400 m, which is consistent with observations [51,52].Time varying stochastic forcing is applied by varying the pycnocline center (see 'Stochastic Forcing' below).

MITgcm ocean model
Ocean model melt rates used as calibration data are calculated by resolving the ice shelf cavity circulation using the Massachusetts Institute of Technology General Circulation Model (MITgcm) [82].The procedure applied to determine ocean model melt rates at timeslices τ 1 , . . ., τ n under a given forcing P c (t) is as follows: (1) run the ice sheet model (with parameterized melting) under this forcing profile; (2) use the output of this to determine ice shelf geometries at timeslices t = τ 1 , . . ., τ n ; (3) for each of these geometries, run the ocean model in this geometry, with forcing applied via a restoring boundary condition corresponding to the profiles P c (τ k ).The restoring boundary condition is applied at the downstream end of the domain at x = 360 km (figure 2a), where the temperature and salinity are restored to vertical profiles T a and S a over a distance of five horizontal grid cells with a restoring timescale of 12 hours.An example of melt rates fields ṁk param and ṁk ocean-model produced by this procedure is shown in figure 2.
The ocean model grid has 55 layers with a vertical spacing of dz = 20 m, and a horizontal resolution of dx = 1 km.We use the MITgcm in hydrostatic mode with an implicit nonlinear free surface scheme, a third-order direct space-time flux limited advection scheme, and a non-linear equation of state [83].The Pacanowski-Philander [84] scheme parameterizes vertical mixing.Constant values of 15 and 2.5 m 2 s -1 are used for the horizontal Laplacian viscosity and horizontal diffusivity, respectively.The equations are solved on an f -plane with f = −1.4× 10 −4 s -1 .For each geometry, the MITgcm is run for three months, using a timestep of 30 seconds, after which the configuration is in quasi-steady state.The ocean model melt rate is taken as the melt rate after three months of the simulation.The drag coefficient in the three-equation formulation of melting [85] used in the MITgcm is taken to be 9 × 10 -3 ; this value ensures that the ocean model melt rate in the post-initialization geometries (see 'Ice Sheet Model Initialization') closely matches observed total meltwater flux values [e.g.52] from Pine Island Glacier.

Ice sheet model initialization
Following [49], we apply a two-stage initialization procedure, outlined in figure 3a.In the first initialization stage, the ice geometry is timestepped from an initial configuration in which the ice-surface is 150 m above sea level for 50 years (note that WAVI uses a hydrostatic flotation condition, so specifying the ice surface and bed elevation prescribes the ice thickness everywhere).Following this, the ice is approximately in steady state, with ice shelf geometry shown in figure 3c.
In the second stage of the initialization procedure, melting is turned on (figure 3).The ice geometry is then timestepped from that at the end of the first initialization stage for fifty years using a constant ocean forcing with P c = −500 m.This pycnocline depth corresponds to typical conditions offshore of the WAIS (i.e.neither warm not cold) [51,52].In the following, we refer to warm forcing as constant forcing with P c = −400 m, corresponding approximately to the shallowest recorded pycnocline depth [51].Similarly, we refer to cold forcing as constant forcing with P c = −600 m, corresponding approximately to the deepest recorded pycnocline depth [51].The second initialization stage is performed independently for each value of M .The (Mdependent) state at the end of the second initialization stage (figure 3c) is then used as the initial condition in the following retreat simulations (figure 3).
Note that for a consistent estimate of SLR contributions from simulations with different values of M , we require similar initial conditions, chosen to be a grounding line at or near the seabed ridge crest.For M 1.5, the ice retreats irreversibly down the ridge during the second initialization stage.We therefore consider only M values smaller than this.In addition, we impose that a constant warm forcing applied to the shelf should initiate retreat (WAIS retreat was, in practice, hypothesised to be initiated with forcing oscillating between warm and cold [11]); we found that for M 0.5, no ice sheet retreat was initiated under warm forcing.Therefore, we restrict ourselves to the range 0.5 ≤ M ≤ 1.5.Note that this restriction is consistent with our Bayesian framework: it is equivalent to setting the prior density to zero outside the range 0.5 ≤ M ≤ 1.5, based on observational constraints.
During the second initialization stage, the ice shelf thins in response to applied melting, but the grounding line does not retreat (figure 3c).The mean melt rate after the second initialization stage is only weakly dependent on M (figure 3b).If the geometries at the end of the second initialization were identical for different values of M , the mean melt rate in the simulation with M = 1.5 would be 3 times as large as that with M = 0.5 (black dashed line in figure 3b); however, owing to temperature-depth effects, this value is only approximately 1.1 times (approximately 23.5 m year -1 in the M = 1.5 case versus approximately 21.3 m year -1 in the M = 0.5 case, see figure 3b).As the ice shelf thins in response to melting, it shallows, exposing it to colder ocean conditions, reducing melt rates sharply and restricting further thinning (the melt rate is proportional to (T a − T f ) 2 , which varies sharply with depth, particularly in the depth range occupied by the ice shelf in the second calibration phase, see figure 3d).

Stochastic forcing
Following the two stage initialization proceedure outlined above, stochastic forcing is applied via ambient ocean conditions: where P c,0 = −500 m is the pycnocline depth in the second stage of the initialization procedure, T (F) is a forcing-scenario-dependent (i.e.anthropogenic or counterfactual) trend (see below), A is the amplitude of random forcing, and R(t) is a first-order autoregressive process, containing the stochastic part of the forcing.In the results shown here, we use A = 100 m, which agrees with observed internal variability in the Amundsen Sea [52].In a first-order autoregressive time-series, the following value is decomposed into a component proportional to the current entry, whose constant of proportionality describes the persistence timescale of the variability, and an additive white-noise term.We take the same autocorrelation function as [35], with interdecadalto-decadal timescales well represented.Anthropogenic and counterfactual ensembles are distinguished via the trend T (F): realizations of forcing from the counterfactual ensemble have no trend added to them, T = 0 m; realizations of forcing in the anthropogenic ensemble have a linear trend, T = A 0 (t/100 yrs), where A 0 = 100 m is the per-century shallowing trend of the pycnocline (figure 2g).

Bootstrapping distributions of sea level rise
Each of the n realizations of forcing yields a parametrically-calibrated distribution of SLR for each time in the simulation.Thus, for any time and any SLR, we have n values of the distributions from both anthropogenic and counterfactual ensembles (supplementary figure 5).An uncertainty estimate in the anthropogenic enhancement ratio is constructed by bootstrapping these values -resampling from these n values with replacement (here, we sample 1000 times); the resulting set yields a standard deviation λ = λ(SLR, t) for both anthropogenic and counterfactual ensembles (supplementary figure 5).Using subscripts to denote the ensemble (that is, counterfactual or anthropogenic), the upper bound shown in figure 4b-d

Fig. 1
Fig.1Sea level rise from marine ice sheets is not necessarily an indicator of climate change.Schematic diagram demonstrating how an ice sheet configuration that remains stable under a realization of forcing including anthropogenic climate change (orange) may experience runaway retreat under a different, counterfactual realization of forcing with no anthropogenic climate change (green).As a result, grounding line retreat (filled dots in ice shelf configurations) and SLR are much higher in the counterfactual case.Once initiated (say, at the star), retreat from a topographic high is sustained by ice-ocean feedbacks.

Fig. 2
Fig. 2 Strong dependence of simulated marine ice sheet sea level rise on both forcing and model parameters.(a) Bathymetry (given by equation 7) of the marine ice sheet configuration.(b) Initial ice thickness along the dashed centerline in (a) for M = 1.The gray line indicates sea level.(c)-(d) Ambient temperature Ta (c) and salinity Sa (d) used in the parameterisation of melting and as restoring boundary conditions in the ocean model (methods).Pc denotes the pycnocline center, which parameterizes the piecewise linear forcing profiles and is oscillated to mimic variability.(e) Time evolution of a single realization of forcing and (f) corresponding SLR contributions for different values of M ∈ {0.5, 0.75, 1.0, 1.25, 1.5} (the arrow indicates the direction of increasing M ).Blue and red regions in (e) indicate whether the forcing is warmer (shallower pycnocline) or colder (deeper pycnocline) than during the calibration phase, where Pc = -500 m (black horizontal line), and shaded red regions indicate two prominent warm periods.The black dashed line indicates the 100 m/century anthropogenic trend in the pycnocline depth.(g) Time evolution of pycnocline centres Pc in all realizations of forcing.Here, orange curves correspond to forcing scenarios with an anthropogenic trend of a 100 m/century shallowing of the pycnocline, while green curves correspond to a counterfactual scenario, with no trend in the forcing (methods).In both cases, faint curves indicate individual ensemble members, while solid curves indicate ensemble means, and dashed lines indicate the externally imposed trend (the T (F ) term in equation (13)), i.e 100m/century and 0m/century shallowing of the pycnocline in the anthropogenic and counterfactual cases, respectively.(h) SLR after 100 years as a function of M for a subset of the different realizations of forcing.Each line corresponds to an individual realization of forcing, and colors indicate whether the forcing is drawn from the anthropogenic (orange) or counterfactual (green) ensemble.Blue hue points correspond to the points shown in panel f.The arrow indicates the curve referred to as the 'highlighted' curve in the main text.

Fig. 3 Fig. 4
Fig. 3 Influence of anthropogenic forcing on distributions of sea level rise.(a) Time evolution (running bottom to top) of distributions of SLR from ensembles with an anthropogenic trend in forcing (orange) and with a counterfactual trend (i.e.no-trend) in forcing (green).Filled markers indicate the median of the distributions at the corresponding time.(b)-(d) Summary statistics of the distributions in (a) as follows: (b) median, (c) skewness and (d) kurtosis.In each, the dashed lines indicate the corresponding summary statistics for distributions obtained without parametric calibration, obtained by assigning equal likelihood to each value of M .