Committed retreat of Smith, Pope, and Kohler Glaciers over the next 30 years inferred by transient model calibration

Abstract. A glacial flow model of Smith, Pope and Kohler Glaciers is calibrated by means of control methods against time varying, annually resolved observations of ice height and velocities, covering the period 2002 to 2011. The inversion – termed "transient calibration" – produces an optimal set of time-mean, spatially varying parameters together with a time-evolving state that accounts for the transient nature of observations and the model dynamics. Serving as an optimal initial condition, the estimated state for 2011 is used, with no additional forcing, for predicting grounded ice volume loss and grounding line retreat over the ensuing 30 years. The transiently calibrated model predicts a near-steady loss of grounded ice volume of approximately 21 km3 a−1 over this period, as well as loss of 33 km2 a−1 grounded area. We contrast this prediction with one obtained following a commonly used "snapshot" or steady-state inversion, which does not consider time dependence and assumes all observations to be contemporaneous. Transient calibration is shown to achieve a better fit with observations of thinning and grounding line retreat histories, and yields a quantitatively different projection with respect to ice volume loss and ungrounding. Sensitivity studies suggest large near-future levels of unforced, i.e., committed sea level contribution from these ice streams under reasonable assumptions regarding uncertainties of the unknown parameters.


Introduction
Smith, Pope, and Kohler Glaciers, three narrow (∼ 10 km wide), interconnected West Antarctic ice streams, have ex-hibited substantial thinning and speedup in recent years.As these ice streams are smaller than neighboring Thwaites and Pine Island Glaciers -the contribution of Smith Glacier to total Amundsen Embayment grounding-line flux is ∼ 7-8 times smaller than that of Pine Island or Thwaites (Shepherd et al., 2002) -focus is often placed upon these larger ice streams, with regard to both modeling and observations of the ice shelves and sub-shelf environments (e.g., Jenkins et al., 2010;Tinto and Bell, 2011;Favier et al., 2014;Joughin et al., 2014).However, high thinning rates have been observed near the Smith terminus, even larger than that of Pine Island and Thwaites (Shepherd et al., 2002;McMillan et al., 2014).Additionally, substantial retreat of the Smith grounding line has been observed (Rignot et al., 2014), suggesting that the ice stream may be subject to the same instability thought to be underway on Thwaites (Joughin et al., 2014).As such, there is a need to develop a quantitative dynamical understanding of the causes of this retreat; and if possible, to determine whether it will continue at similar rates.
The problem of projecting ice sheet behavior is challenging, in part due to incomplete understanding of physical processes (Vaughan and Arthern, 2007), but also due to difficulties in estimating the state of an ice sheet at any given time.Unlike other components of the climate (Taylor et al., 2012), ice sheet models cannot be "spun up" to the present state, as the required historic forcing fields are not available.Rather, the models must be initialized from observations, which are mostly limited to surface properties such as surface elevation and velocity.A widely used methodology is one to which we will refer as "snapshot" calibration, first introduced by MacAyeal (1992), and which solves an inverse or optimal Published by Copernicus Publications on behalf of the European Geosciences Union.D. N. Goldberg et al.: Ice stream forecast through transient calibration control problem.In this technique, an optimal set of parameters relating to sliding stress (and possibly ice-shelf viscosity) is found through a least-squares fit of the ice model's nonlinear momentum balance to a given velocity field.Time dependence is not considered, since the momentum balance (or rather stress balance) is non-inertial.We choose the term "snapshot" because it applies to ice velocity and geometry at a single instant, assumed to be the same for both data sets.
A number of studies have employed snapshot calibrations to make near-future projections of the behavior of Pine Island and Thwaites Glaciers in response to varying forcing scenarios (Payne et al., 2004;Joughin et al., 2010;Favier et al., 2014;Joughin et al., 2014;Seroussi et al., 2014).These studies have deepened our understanding of the behavior of these ice streams.However, the use of snapshot calibrations in ice sheet projections is potentially problematic: any temporal inconsistencies among data sets can lead to nonphysical transients which persist for decades, which is not ideal if the goal is projection on a similar timescale.Inconsistency between the data and model discretization can have a similar effect.For instance, co-located gridded velocity and thickness data require interpolation for application to a model whose discretization staggers these fields, potentially leading to transient nonphysical artifacts (Seroussi et al., 2011).An oft-used approach is to allow the model to adjust to these inconsistencies before conducting experiments.The model may then have drifted to a state far from contemporaneous observations, with potentially different sensitivities.
As the observational record grows, so does the availability of data for the same geographic areas at multiple points in time.It is sensible, then, to make use of this temporal resolution for the purpose of constraining the time-evolving state of an ice stream, with the significant benefit of producing initial conditions for forecasting from a realistic past trajectory.Such an approach, which we term "transient calibration", is well developed in other areas of geophysics, e.g., in oceanography where it is known as "state and parameter estimation" (Wunsch and Heimbach, 2007), or reservoir modeling where it is known as "history matching" (Oliver et al., 2008).Here we present the results of such a calibration, applied to Pope, Smith and Kohler Glaciers.Calibrated model parameters are the result of an inversion in which a time-evolving model produces an optimal fit to a 10-year time series of surface elevation and velocity observations.The model is then run for an additional 30 years.In the transiently calibrated run, rapid grounding-line retreat continues for another decade, but then slows, while loss of grounded ice remains near constant at ∼ 21 km 3 a −1 (or ∼ 0.06 mm a −1 sea level contribution).We show that the predicted high levels of ice loss are relatively insensitive to any future changes in forcing, and to any systematic errors in our calibration.
Transient calibration of a model of an Antarctic ice stream with temporally resolved plan-view data has not previously been carried out, though we point out that Larour et al. (2014) used methods similar to those used in this study to infer sur-face mass balance over the Northeast Greenland Ice Stream over a 6-year period from laser altimetry.No future projections were made in their study.
We proceed with detailing what we mean by "snapshot" vs. "transient" calibration of an ice flow model, and show how ice sheet observations are used in this process (Sect.2).We then describe the observational data (Sect.3), as well as the model and the details of the calibration used in this study (Sect. 4).Results of the calibration and projection are presented in Sect.5, followed by an investigation of the sensitivity of these results to plausible uncertainties in the parameter estimates (Sect.6).

Snapshot calibration
A widely used approach for single-time observations is to invert for uncertain control variables, using a stress balance model, via the adjoint or Lagrange multiplier method.MacAyeal (1992) applied such an optimal control method, in which the misfit between model velocity, u, and observed velocity, u * , is minimized with respect to unknown (or uncertain) variables λ (often referred to as a control variables), subject to the constraint that the velocity satisfies the nonlinear stress balance, written in the generic form L(u, λ) = 0.The misfit (or cost) function is expressed as where u i and u * i are at location i (grid cell or node), and η(u i ) 2 the uncertainty of the observation.The constrained optimization problem may be turned into an unconstrained one by introducing Lagrange multipliers µ i : where L i (u, β) = 0 is the discretized form of the stress balance at node i.By finding a saddle point of J with respect to the parameters and Lagrange multipliers, an extremal point of J snap is found in parameter space with the stress balance enforced exactly.The coefficient β of the linear sliding law is often used as the control variable λ.J snap is sometimes extended with an additional "smoothing" term that penalizes small-scale variations in the control parameters (e.g., Morlighem et al., 2010).The ice geometry (i.e., surface and bed elevation) is assumed to be known exactly.
In MacAyeal (1992), the model considered is the shallow shelf approximation (SSA) (Morland and Shoemaker, 1982;MacAyeal, 1989) and the control variable is β as above.Development of sophisticated glacial flow codes and the consideration of ice-shelf physics have led to the use of alternative or augmented control spaces (e.g., Larour et al., 2004;Joughin et al., 2010;Lee et al., 2015) and the use of higherorder stress balances (e.g., Morlighem et al., 2010;Goldberg and Sergienko, 2011;Petra et al., 2012).
The Lagrange multipliers µ i are then used to calculate the gradient of J snap with respect to the control variables λ, which can in turn be used to carry out the minimization of J snap via gradient descent or quasi-Newton optimization methods.The µ i are found by solving the adjoint of L , the linearization of the operator L. The adjoint method is popular for snapshot calibrations in glaciology due to the fact that L is self-adjoint, i.e., the adjoint operator, can be solved by the same code used to solve L if the dependence of ice viscosity on strain rates is ignored.

Transient calibration
When observations distributed in time are available together with a time-evolving model, the "snapshot" calibration can be extended to what we term "transient" calibration, which consists of optimizing agreement of the model with observational data at multiple time levels, with both the nonlinear stress balance and ice thickness evolution enforced as model equations.This is equivalent to the following constrained cost function, which should be compared against J snap : where s is ice surface elevation, the superscript k is the time index, and the asterisk indicates observational values.χ (u) ki and χ (s) ki are equal to 1 if there is an observation at cell i and time step k, 0 otherwise.ω u and ω s are weights to impose relative importance of observations.The Lagrangian J now extends to one with time-evolving Lagrange multipliers, i.e., where the model equations are written in generic form x (k+1) = F (x (k) ), and x represents the model state, i.e., the minimal set of variables needed to step forward the model and to evaluate J trans .
Minimization of J trans can be carried out in a similar manner, by use of its gradient with respect to the control vector.However, gradient calculation is more complicated, now requiring a time-dependent adjoint model, which can be derived via the continuous-form adjoint of the model equations, as has been done for simplified ocean models (Tziperman and Thacker, 1989), or by means of algorithmic differentiation (AD; Griewank and Walther, 2008).Used extensively in ocean modeling (e.g., Heimbach et al., 2005;Wunsch and Heimbach, 2013), the use of AD tools in land ice modeling is becoming increasingly common (Heimbach and Bugnion, 2009;Goldberg and Heimbach, 2013;Larour et al., 2014).
In this framework, the control parameters may now be chosen to be time-dependent.However, doing so is meaningful only if physically justified and if sufficient information is available to constrain the larger control space.In the following, unless stated otherwise, time-independent parameters are used.

Observations
The time-dependent observations of velocity and surface elevation in Eq. ( 4) come from two recently generated data sets.One contains InSAR-derived surface velocities of the Smith Glacier region, binned annually to a 500 m grid for the years 2006-2010 (Joughin et al., 2009;Medley et al., 2014).Velocities are available for floating and grounded ice.Coverage is not spatially uniform, but greater in later years.
The other data set is a series of annual surface digital elevation maps (DEMs) from 2001 to 2011 on a 1 km grid.Coverage is consistent between years, but data are not available seaward of the 1996 grounding line (Rignot et al., 2014), or on slow inter-stream ridges.Figure 1 shows the geographic region of study along with the acceleration and thinning recorded by the transient data sets.The 2001 surface is not from 2001 measurements, but is simply an extrapolation backward in time from later years.Further details of this data set are given in Appendix A.
In addition to these time-dependent data sets, we use the BEDMAP2 bed topography (Fretwell et al., 2013) and the MEaSUREs (Making Earth System Data Records for Use in Research Environments) (450 m grid) data set (Rignot et al., 2011).We also use the Arthern et al. (2006) accumulation data set to estimate ice temperatures in the region, as explained in Appendix B1.

Model and calibration setup
The land ice model used in this study is that described in Goldberg and Heimbach (2013).The model's stress balance is depth-integrated, similarly to the shallow shelf equations, but the effects of vertical shearing are represented (Goldberg, 2011).Grounding-line migration is implemented through a hydrostatic floatation condition.As described in Goldberg and Heimbach (2013), the model has been successfully differentiated using the AD software tool Transformations of Algorithms in Fortran (TAF; Giering and Kaminski, 1998).We solve the land ice equations in the domain shown in other fields are interpolated to this grid.This allows for resolution of the relatively narrow ice streams.However, the domain does not include ice shelf seaward of the 1996 grounding line, as we explain below.We do not account for the effects of firn on ice dynamics.
The observations used in our transient calibration are those described in Sect.3. The initial ice thickness in each model run is from the 2001 DEM.Subsequent DEMs are applied to the cost function J trans at the end of each model year, as are velocity constraints in the years and locations available.For the snapshot calibration we use ice geometry from the 2002 DEM; as we do not have 2002 velocities, MEaSUREs velocities are used as constraints.As discussed below, in transient calibrations the domain excludes ice shelves.We carry out snapshot calibrations in the same domain to enable comparison, and the resulting parameters become initial guesses in our transient calibration.Similar to other ice model calibrations, the basal sliding parameter β 2 is a control parameter.Our other control parameters, less common in glaciological inversions, arise from the nature of the transient calibration and the data sets used, as explained below.
Our results in Sect. 5 are generated assuming timeinvariant control parameters.In Sect.6.2 we allow for timedependent parameters, and consider the implications of the results.

Boundary stresses as control parameters
Our transient surface observations only give values inland of the 1996 grounding line.Time-resolved annual velocity observations are provided for the ice shelves, but only from 2007 to 2010.Including ice shelves in our domain, then, would require estimation of transient ice-shelf thickness from 2001-2011.Such an estimate would be very poorly constrained (see Sect. 7 for a discussion on this topic).We overcome this problem by formulating an open boundary estimation problem (see, e.g., Gebbie et al. (2006) for an oceanographic analogue), with the 1996 grounding line as the downstream boundary of the domain (see Fig. 1).Stresses at the grounding line, which would otherwise be part of the stress balance solution, must now be imposed along this boundary.The action of the membrane stress tensor (Hindmarsh, 2006) along a horizontal boundary has two components: normal membrane stress σ and shear membrane stress τ , as explained in more detail in Appendix B2.In the model, σ and τ can be defined along any horizontal boundary, floating or grounded.These boundary stresses are not known a priori, and we treat them as unknown spatially varying (along boundaries) control parameters to be estimated via calibration, with two unknowns (σ j and τ j ) for each rectangular cell boundary j .Where the domain borders a slowmoving ridge velocities are set to zero, and boundary stresses are not applied.

Boundary volume flux as a control parameter
In our transient calibrations, the ice flux into the domain must be estimated.This is due to the incomplete coverage of the time-dependent velocities, which leaves the upstream regions poorly constrained, leading to anomalously high thinning.To address this we consider boundary fluxes q x and q y as control parameters at x and y facing boundaries, respectively.These boundary fluxes enter the model through the continuity equation, which is solved via a finite-volume scheme, and are treated as constant over a cell boundary.Boundary fluxes are not imposed along the internal boundaries with slow-moving ridges, or where boundary stresses are imposed.Note that q x and q y are only used in transient calibration; for snapshot calibration, MEaSUREs velocities do not lead to high thinning rates in these regions, despite no-flow conditions at the upstream boundary.

Calibration results
Our snapshot calibration recovers MEaSUREs velocities to high accuracy.The RMSE with observations is reduced from 140 m a −1 for the initial guess down to 50 m a −1 -but error is actually much lower any most areas outside the margins of the narrow western branch of Kohler entering Dotson Shelf (Fig. 4a).The control parameters adjusted in the snapshot calibration, β, σ and τ , are then used in a transient (but noncalibrated) run from 2002 to 2011.The degree to which this run agrees with the transient observations is demonstrated in Fig. 2a, c, and the top row of Fig. 3.For velocities in the snapshot-calibrated run, the misfit for 2010 -the last year in which velocity observations are available -is largest in Kohler and Smith glaciers, and is up to at the boundary with the slow-moving ridge, which may be because the no-flow condition imposed there by the model is not accurate.By 2011, modeled surface elevation within 20-30 km of the grounding line is ∼ 100 m higher than observed, a misfit that is larger than the impact of the thinning signal itself over the period of integration.The misfits grow with time, and so only the final years are shown at this level of detail.Figure 3 gives surface error along the flow transects from Fig. 1c.
Relative to the time integration with initial state and parameters obtained from the snapshot inversion, the transient calibration gives good agreement, especially with respect to surface elevation (Fig. 2d).The 2011 surface elevation misfit field looks very different to the one inferred from the snapshot calibration, with uniformly small misfits.Figure 3 (bottom row panels) shows the reduction in transient surface elevation misfit along the transects.On Smith and Kohler, misfit in 2010 velocity has decreased, though it is still substantial (Fig. 2b).The relatively low decrease in velocity misfit between snapshot and transient calibration can be explained by our choices of ω u and ω s , which favor surface elevation.
Also the grounding-line behavior is very different between the two simulations.In the snapshot-calibrated run there is almost no retreat, while in the transiently calibrated run the 2011 grounding line has retreated considerably.The modeled 2011 grounding line is not completely coincident with the observed grounding line of Rignot et al. (2014) (digitized and plotted for comparison), particularly in the western part of the Smith/Kohler grounding region.The cause for this discrepancy is unclear; but in any event, the ice in this region does unground in our simulation, it is simply delayed by 5-10 years (see Sect. 5.3).

Adjustment of control parameters
Aside from the boundary volume fluxes q x and q y , the control parameters in the snapshot and transient calibrations have a one-to-one relationship.Thus it is interesting to examine how the parameters are adjusted for transient calibration.In both the snapshot and transient calibration, we infer an area of very weak bed (basal stresses of ∼ 10 Pa or less) in the fastest moving parts of the glaciers (Fig. 4b).The most striking adjustment of basal stress parameters is a strengthening of the bed under the trunks of Pope, Smith and Kohler Glaciers (Fig. 4c).
This strengthening is offset by a decrease in backstress along the grounding line (Fig. 4d).It is possible that our snapshot calibration is equifinal, i.e., that there is more than one combination of boundary stresses and bed parameters to reproduce imposed velocity and elevation observations.In this case our snapshot calibration does not correctly estimate the dynamic state of the system, while the additional information provided by the transient observations allows us to find a better balance between boundary stresses and basal strength.Alternatively, it may be that the temporal mismatch between velocity and altimetry in the snapshot calibration demands a more extensive weak-bedded region than is realistic, with additional buttressing required to match velocities at the grounding line.
A noticeable feature of the transiently calibrated solution is that of "negative buttressing"; i.e., the normal membrane stress in some locations is larger than what would be felt without any ice shelf.This could be for a number of reasons.It is possible the model, and the fit to observations, is insensitive to small-scale oscillations in the boundary stress field.However, it could also be due to errors in the bed topography data: as detailed in Appendix B2, boundary stresses are expressed as a fraction of unconfined membrane stress, which depends on bed depth.Negative buttressing could be compensating for an assumed bed that is too shallow.Finally, the "negative buttressing" may be very real features of the ice sheet.Schoof (2006) demonstrated that even in the absence of an embayed ice shelf, alternating patterns of ridges and ice streams can lead to ice shelf buttressing, whereby the fast-moving streams essentially "pull forward" the ice on the slow-moving ridges.Such a situation could yield negative buttressing at the ridges.Inspection of Fig. 4d shows that the negative values occur at the slow-moving regions in between the narrow, fast-flowing streams.
Another noticeable feature is the "ribbed" pattern that appears in the β 2 field, but not in the snapshot-calibrated field.The cause of this discrepancy is uncertain.It is possible that the observed velocities could be well-represented in the snapshot calibration without these features, but that they are necessary to fit to surface observations.However, it may be related to the "smoothing" term mentioned in Sect.2.1.In both models, a Tikhonov regularization term (i.e., the squareintegrated gradient of β) is added to the cost function -this is done because ice model velocities are insensitive to highwavenumber variations in the basal sliding coefficient, and these scales are poorly constrained (Morlighem et al., 2010).The term is multiplied by a weighting coefficient -and this coefficient is chosen on the basis that β 2 should not vary by a considerable amount on scales smaller than the membrane stress scale (∼ 5 km).Importantly, this weighting is the same for both transient and snapshot calibrations.In making this choice, we may have implicitly introduced a degree of subjectivity to our estimations (Arthern, 2015).Introducing prior information in a more objective manner is beyond the scope of this study, but it is an important goal and should not be overlooked in the future.

Projected ice loss and behavior
The or the portion of a grounded column that would be supported by ocean pressure, and thus is an indicator of sea level contribution.To calculate VAF from the observational data, thickness h obs and height above floatation h af must be inferred from surface and bed data as follows: where ρ i = 918 kg m −3 and ρ w = 1028 kg m −3 are ice and ocean densities, respectively, s obs is surface elevation from the transient DEM set, and R is BEDMAP2 bed elevation.h af is then integrated for VAF.Both snapshot and transient calibrations predict continued contribution to sea level rise.The transiently calibrated model projects ∼ 21 km 3 a −1 grounded ice volume loss from 2011 to 2041 (∼ 0.06 mm sea level equivalent), while the snapshot calibrated model suggests ∼ 25 % more.Thus there is a quantitative impact of the initial state, and therefore of the type of calibration used, on projected sea level contribution from the region.There is an even more pronounced impact on projected grounding line retreat: in the snapshot-calibrated run, almost no ungrounding takes place, while in the transiently calibrated run ungrounding is significant (Fig. 5b).Given the much closer fit of the transiently calibrated simulation to surface observations in a least-square sense, we accept this simulation as a better estimate of the dynamic state of the glaciers in the region.
Spatial patterns of projected grounding-line position for the transiently calibrated run show significant retreat from 2011-2021 (Fig. 6), followed by a slight slowdown in retreat.In contrast, thinning rates remain high throughout the 30-year integration.Grounding-line retreat does not proceed down the deep troughs incised by Smith and Kohler Glaciers, suggesting the retreat predicted by Rignot et al. (2014) might not happen in the near-term.We argue that this is because the troughs are quite narrow, and lateral stresses from areas of shallower bed limit grounding-line retreat.However, other studies suggest grounding line retreat in Amundsen and Bellingshausen ice streams can be episodic rather than sustained due to details of bed geometry (Joughin et al., 2010;Jamieson et al., 2012).Furthermore, while melting under parts of the ice shelf external to the model are implicitly accounted for through boundary stresses and their sensitivities (Sect.6.1), we do not apply melt rates to ice that goes afloat within the domain.Such effects could lead to stronger retreat than what is shown.Thus, we cannot discount further rapid grounding line retreat in the future (i.e., beyond 2041), particularly since thinning rates remain high throughout our simulation.Spatial patterns for the snapshot-calibrated simulation actually show slight thickening in some areas downstream of the observed 2011 grounding line, and otherwise show a more even pattern of thinning (i.e., it is not skewed downstream).
The imposed mass fluxes at the inland boundary are not expected to influence the results: the time scale (30 years) is less than the diffusive time scale for grounding line changes to propagate across the domain (e.g., Payne et al., 2004), which we calculate to be ∼ 150 years based on a nominal surface slope of 0.01, thickness of 1400 m, and velocity scale in the upstream regions of 100 m a −1 (Cuffey and Paterson, 2010).
Finally, it is important to realize that these projections are unforced: the estimated parameters and boundary conditions β, σ and τ (and q x , q y where applicable) are held constant over this time period, and no submarine melt is applied to any areas which unground.This is the basis for referring to the projected grounded ice loss as committed (Price et al., 2011).

Uncertainty of sea level contribution projection
The projection of committed grounded volume loss of 21 km 3 a −1 over the next 3 decades from 2011 onward is subject to uncertainty due its implicit dependence on model parameters.The adjoint capabilities of the model allow us to estimate reasonable bounds on this uncertainty through cal- culation of sensitivities to these parameters, which can be integrated against parameter field perturbations.For instance, Fig. 7a shows the adjoint sensitivity of transiently calibrated VAF loss to the basal sliding parameter β 2 .We refer to this quantity as δ * (β 2 ), and it can be interpreted as follows: assume the β 2 is subject to a perturbation P (x, y).Then the perturbation to VAF loss that follows from this parameter perturbation is given by where D is the model domain.δ * (R), the sensitivity of VAF loss to topography, is plotted in Fig. 7b.Note that the influence of β is sign-definite, i.e., decreasing β anywhere increases ice loss, while lowering the bed only increases ice loss upstream of the projected 2041 grounding line.
If we assume an error of 100 % for each basal sliding parameter -an unlikely scenario, as this would affect the fit to observations -ice loss projections would change by at most 57 %.Other parameters have lower influence, assuming reasonable uncertainties.100 % error in the boundary stress parameters would change the ice loss projection by at most 1 %.The influence of input fluxes q x and q y is similarly small.The full range of bed elevation errors associated with the BEDMAP2 data set would change the projection by at most 30 %.These values are based on linear sensitivities, while our model is nonlinear -but the results are borne out by experiments with finite perturbations.Of course, these fields would not vary independently -but based on these relatively low sensitivities we anticipate that the projected mass loss value is not overwhelmed by its uncertainty.Thus, our conservative uncertainty analysis suggests a level of committed sea level contribution from the region.
The above estimation of uncertainty bounds is tentative.Our inverted parameters have no a priori estimates or uncertainties, and our minimization does not provide a posteriori uncertainties or covariances.Thus we are unable to provide accurate confidence intervals on ice loss based on observational uncertainty.Estimation of a posteriori uncertainties based on observational uncertainties may be possible e.g., through methods that infer the Hessian of the cost function (Kalmikov and Heimbach, 2014;Isaac et al., 2014).Enabling such calculations within our estimation framework is a future research goal.

Time dependence of control parameters
Our adjoint-based calibration framework allows for the estimation/adjustment of control parameters that vary not only in space, but also in time (e.g., Wunsch and Heimbach, 2007 2013).Justification for doing so derives from the physical interpretation of these parameters, e.g., boundary stresses representing far-field stresses in the ice shelves, which could change due to crevassing or ocean melting.We investigate whether such time dependence can be inferred from the observations.In our framework, parameters vary piecewiselinearly over predefined time intervals of uniform length.For instance, with intervals of 5 years, and over the interval from t = 5 years to t = 10 years, σ j (the normal stress at face j ) takes on the values The parameters σ (5) j and σ (10) j (and σ (0) j ) are distinct for each cell face, and constitute additional parameters for the system.Thus, the greater the temporal resolution, the more calibration parameters are involved.Considering the increase in size of the parameter space, the additional information is only meaningful if it improves the fit of the calibration.
To facilitate the discussion we define an annual cost function, i.e., a breakdown of J trans by year.That is, for each year k we define In Fig. 7c this value is plotted by year for different experiments.The annual cost functions resulting from the snapshot and transient calibrations are plotted (although recall that the snapshot calibration is not designed nor intended to explicitly reduce the transient misfit reflected by J trans ).Results from two additional calibrations are shown as well.In the first, the β 2 parameter is assumed time-invariant, but boundary stresses are allowed to vary linearly over the 2001-2011 period as described above.In the second, boundary stresses are constant while β 2 is allowed to vary linearly in time.In each case, the number of degrees of freedom which describe the time-variant control doubles.The cost function J trans is reduced, but the reduction is relatively small (∼ 20 %).We display the estimated parameters for the linear-in-time boundary stresses experiment in Fig. 7d, by plotting buttressing at the beginning and the end of the simulation -or, more accurately, −γ (x,0) σ and −γ (x,10) σ , where γ σ determines buttressing level (cf.Eq.B3) and the number in the superscript has the same meaning as in Eq. ( 9) above.Results are displayed relative to the time-independent parameters found above.The pattern corresponds to a slight loss in buttressing from 2001-2011, albeit of a smaller magnitude than its temporal average.(Note that the loss of basal stress due to grounding-line retreat, found to be an important mechanism by Joughin et al. (2014), is resolved by our model and therefore not implicit in inferred boundary stresses.)Also, the pattern is slightly different at Smith as compared to Kohler and Pope.The corresponding inferred pattern of time-dependent β 2 (not shown) is somewhat noisy but contains a clear signal of bed weakening under fast-flowing regions just upstream of the 2011 grounding line.
We emphasize that the above results should be regarded with caution due to the relatively small reduction in J trans resulting from additional degrees of freedom.However, we are not aware of a quantitative measure to determine whether the improvement is significant, i.e., whether the inferred timedependent adjustment of the parameters can be regarded as real, or just "noise".It is also possible that the small reduction of the cost function is due to the shortness of the estimation period, over which the distinction between timevarying vs. time-mean controls does not influence the solution significantly.However, the pattern of temporal buttressing change is at least plausible given observed submarine melt rates (Pritchard et al., 2012) and loss of ice rumples and pinning points (Rignot et al., 2014).Thus the information presented may be of use in future studies of the region that include ice shelves, as it could be used to accept or reject various ice shelf forcing scenarios on the basis of resulting changes in buttressing.Questions regarding the level of temporal data resolution required to constrain time-varying parameters, and of appropriate criteria to identify overfitting of such parameters, are targets for future work.

Discussion
We do not hold our snapshot calibration to be the best possible in the sense of reproducing spatiotemporally resolved observations.For this calibration we used MEaSUREs velocities, which have a much later time stamp than the ice geometry used.This choice was made because no 2002 velocity data were available.Nevertheless, our results demonstrate that a snapshot calibration with non-contemporaneous data, or data sets that might be inconsistent with each other if used at face value in a dynamical framework, cannot be expected to reproduce time-dependent behavior, whereas transient calibration can take account of time-varying data in order to better reproduce observations, thereby giving more confidence in near-future projections of ice sheet behavior.The nonlinear least-squares framework ensures that mutually incompatible data sets can be properly weighted, i.e., interpolated by the model dynamics, instead of having to be simultaneously fulfilled exactly.Importantly, within such a framework increased care must be taken to provide useful error estimates for each observational element (the η entries in Eq. 4).This requires understanding of measurement errors, potential systematic biases, and representation errors.
While transient calibration can potentially constrain timevarying behavior of poorly known control parameters, care must be taken that the increase in dimension of the parameter set yields an improved fit with observations.Otherwise, the additional information provided (relative to time-invariant parameters) may be of limited use.For our calibration, we see that allowing for time-varying control parameters only provides a small improvement of fit, and thus we do not reject the null hypothesis that far-field buttressing (and bed strength) did not change from 2001-2011.While it is possible that buttressing did decrease over this time, it is also possible that some perturbation to the system occurred long before observations began, and the 2001-2011 retreat is just a continued response to this perturbation.More investigation is needed regarding the details of how temporal observational sampling is able to constrain temporal structure of poorly known parameters.
As explained in Sect.4, the decision was made to remove the ice shelves seaward of the 1996 grounding line from the domain in favor of boundary stresses.It is worth briefly considering the implication of this decision.BEDMAP2 draws ice shelf thickness data for the region from Griggs and Bamber (2011), who give an effective timestamp of January 1995.It is likely that the change in thickness from this date to 2001 was both non-negligible and roughly on order with the change in thickness over the 2001-2011 window (Paolo et al., 2015).Apart from BEDMAP2 our only available ice shelf data are velocities in 2007-2010 (2006 had little ice shelf coverage); we do not possess any data regarding ice shelf thickness change over time.In order to model the evolution of the ice shelves, then, it would be required to estimate 2001 ice shelf thickness, as well as the spatially and temporally varying melt rates and effective Glen's Law ice stiffness parameter (A) from 2001 to 2011.Data from grounded ice (such as velocities and surface elevation) are not sufficient to infer such detailed information about ice shelves, as modeling studies indicate that grounded ice evolution might be insensitive to melt rates and ice stiffness over large parts of the ice shelves (Goldberg and Heimbach, 2013).Thus estimates of the above parameters would need to be made, with only velocity information at the end of the decade as a basis of improving the estimates.Such a strategy would be ill-posed owing not only to the limited temporal coverage, but also to the fact that both ice shelf thickness and Glen's Law parameter determine velocities.Thus the approach, while not impossible, would require very careful quantification of a posteriori parameter uncertainty -which, as stated previously, requires more sophisticated computational tools than those used for this study.However, incorporation of ice shelf data and sim- In addition to the control parameters discussed above (boundary stresses, upstream fluxes, and sliding parameters), two others were initially investigated: adjustments to initial (2001) surface elevation, and adjustments to bed elevation.These fields were considered to be potentially important for observational agreement, as the 2001 DEM from which the initial condition is derived is a backward-in-time extrapolation of later measurements, and bed topography is considered to be a source of uncertainty for ice flow (Durand et al., 2011;Morlighem et al., 2011;Sun et al., 2014).However, significant adjustments were not found for either (the inversion adjusted initial surface on the order of millimeters, and the bed on the order of meters), and their inclusion did not improve the fit to observations.Thus these control variables were not considered further.We point out these results may depend somewhat on our prior assumptions of their variability, which is implicitly imposed by the scaling of cost function gradients (see Appendix B3), and we stress the importance of choosing conservative and unbiased prior information in future transient ice sheet calibrations.
We briefly consider potential reasons for the discrepancy between our modeled 2011 grounding line and that of Rignot et al. (2014).As mentioned in Sect.4, we do not account for the effects of firn density in our model, neither have our transient surface data been corrected for firn.As the depth of the firn layer can affect the floatation condition (e.g, Griggs and Bamber, 2011), it is reasonable to ask whether these omissions can explain the disagreement between our modeled 2011 grounding line and observations.Figure 8 gives a detailed comparison between the modeled and observed grounding lines, as well as the 2011 grounding line inferred from the 2011 DEM and the BEDMAP2 data via Eq.( 6).There is slight disagreement between the latter two grounding line estimates, but it does not explain the erroneously grounded region in our model.Rather, we suggest this region is anomalously thick (and therefore grounded) due to buttressing from the small grounded "island" at the Smith Glacier grounding line, which is not visible in the Rignot et al. (2014) data.Furthermore, we point out that grounding line agreement is not explicitly accounted for in our transient cost function.Still, future studies should account for firn effects in order to achieve better agreement with grounding line observations.

Conclusions
Generalizing optimal control methods based on steady-state adjoint models well-known in glaciology to those using a transient forward and adjoint model, enables us to perform model calibration based on simultaneous state and parameter estimation through a nonlinear least-squares fit of a model to time-resolved observations.We perform such a transient calibration for the grounded portion of the Smith, Pope and Kohler Glacier region based on velocity and surface observations covering the years 2001-2011.This transient calibration is compared with a "snapshot" calibration of the same region based on instantaneous (and assumed contemporaneous) observations.The transient calibration agrees far better with spatially and temporally resolved observations, giving increased confidence in near-future behavior predicted by the model.
Extending the simulations beyond the 2001-2011 calibration period, both snapshot-and transiently calibrated models are run in "predictive mode" from 2011 to 2041, without any changes in boundary conditions or external forcing.Both show a significant sea level contribution.That of the transiently calibrated model is nearly 20 % smaller, but with significant grounding line retreat and grounding lineconcentrated thinning.
Sensitivity calculations suggest that, under reasonable assumptions regarding parameter uncertainties, a committed grounded ice loss of ∼ 21 km 3 yr −1 can be expected from the region, even in the absence of external forcing or climateinduced feedbacks.Our sensitivity analysis does not replace a comprehensive uncertainty quantification of projected ice volume loss, and a more complete end-to-end uncertainty propagation chain is needed for transient ice model calibration.
As the catchment of Smith, Pope and Kohler Glaciers is relatively small, the potential for sea level contribution is not as large as that of Thwaites and Pine Island (Joughin et al., 2010(Joughin et al., , 2014)).Nevertheless, the volume loss from these glaciers is quite high given their size, and our projection shows no indication of it slowing in the next few decades.Furthermore, significant thinning of the region could affect flow of nearby ice streams by changing surface gradients.The methodology of transient calibration introduced in this study -which has not previously been applied to a marinebased Antarctic ice stream -could be applied to other regions of Antarctica to better constrain near-future behavior.To do this, better availability of spatially and temporally resolved observations, for both grounded and floating ice, along with credible error estimates for each observational element will be essential.stress balance were again solved, but only over the grounded part of the domain, with S • n imposed to be equal to the same (σ (θ ), τ (θ )), then velocities and stresses within the ice would be the same.(This is mathematically true for the depth-averaged hydrostatic stress balance used in this study; while it does not hold for the general Stokes balance, any nonhydrostatic effects will likely be limited to the vicinity of the grounding line.)In other words, the effect of the ice shelf on grounded velocities (and thickness evolution) is imposed solely through σ (θ ) and τ (θ ).Thus, in our runs, the boundary of the computational domain is internal to the ice body (and initially coincides with the grounding line).As our model has a rectangular grid, this boundary is not a continuous line but a collection of cell faces, some directed in (i.e., normal to) the x direction and some in the y direction (Fig. B1b).We implement σ and τ as a set of parameters, with a separate value for each cell face.Effectively, we implement a Neumann boundary condition, albeit one that does not depend uniquely on the ice thickness and bed depth, as is the case for a calving cliff.Rather, the boundary condition is a forcing that needs to be estimated.These parameters are expressed not as stresses but as an excess fraction of the unconstrained membrane stress.Thus, σ = (1 + γ σ ) σ cf , τ = γ τ σ cf (B3) and γ σ , γ τ are the actual parameters.Notice that in this formulation σ and τ depend on bed depth at the cell face according to the topographic data set (in this case BEDMAP2, Fretwell et al., 2013).
In some of our simulations, the boundary of the domain does not remain coincident with the grounding line, as there is grounding-line retreat.The grid cell faces along which stresses are imposed do not follow the grounding line in this case; rather, they remain fixed and we effectively impose the stresses on a portion of the shelf.However, they are still imposed far from the calving front, and σ and τ (s) are still representative of buttressing within the ice shelf.
In Fig. 4d we distinguish between σ (x) and τ (x) , boundary stresses along faces normal to the x direction (and likewise γ (x) σ , τ (x) σ ), and σ (y) and τ (y) .Note than σ (x) and τ (y) enter into the x momentum balance (and are therefore more relevant to flow predominantly in the x direction).

B3 Normalization of gradient information
When carrying out adjoint-based inversions or state estimations with heterogeneous control fields, the units of the different control variables must be accounted for.For instance, the boundary stress parameters as described above nominally vary between 0 and 1 (dimensionless), while values on the order 10 4 m 2 a −1 were found for the input flux parameters.Thus for a given stress parameter σ i and a given flux parameter q j , one might expect ∂J trans ∂σ i to be several orders of magnitude larger than ∂J trans ∂q j .The gradient with respect to the parameter set, and thus the search direction in parameter space, would be overwhelmed by the gradient with respect to input fluxes.This issue is addressed by normalizing the cost function gradient by nominal "unit" values, where the unit value corresponds to the type of parameter.In our inversion, values of 0.1, 5 × 10 4 m 2 a −1 , and 10 Pa (m a −1 ) −1 were used for boundary stresses, input fluxes, and basal sliding parameters, respectively.Additionally, values of 1 and 10 m were used for adjustments to the initial surface and the bed elevation, respectively (see Sect. 7 of main text).The normalization factor for the initial condition was chosen since this value was in line with the errors applied to the surface observations.The factor for the bed was chosen due to the relatively small bed adjustments required by mass continuity considerations for this region (Morlighem et al., 2011;Rignot et al., 2014).
www.the-cryosphere.net/9/2429/2015/The Cryosphere, 9, 2429-2446, 2015 Figure 1.(a) Ice speed in the Pope/Smith/Kohler and Crosson/Dotson system.The white contour is the grounding line as given by BEDMAP2, and the magenta contour represents the limits of the transient surface elevation data set.The rectangular box shows the subdomain used for our state estimate simulations -boundary stresses are imposed along the black contour and boundary fluxes are imposed along the light blue boundaries.(b) Norm of velocity change between 2006 and 2010 within the model domain, excluding the areas of no coverage in either 2006 or 2010.(c) Cumulative surface thinning, 2001-2011 in the surface elevation data set.The shaded region shows where data are available.(d, e, f) Hövmoller plots of cumulative thinning along transects in (c) in descending order.

Figure 2 .
Figure 2. Difference between (top panels) modeled and observed velocities in 2010 (the last year available) and (bottom panels) modeled and observed surface elevation in 2011.Left panels: snapshot calibration.Right panels: transient calibration.The magenta contours represent modeled grounding lines in 2011.In (d), the green hatches give the 2011 grounding line position reported by Rignot et al. (2014) (digitized from the publication).

∼Figure 3 .Figure 4 .
Figure 3.Comparison of transient misfit of modeled surface elevation between snapshot and transient calibration along different flow lines.From left to right, panels correspond to flow lines in Fig. 1c in descending order.Top row panels: snapshot calibration.Bottom row panels: transient calibration.
Figure 5. (a) Sea level contribution from the region and (b) total ungrounded area in domain through 2041 based on snapshot and transient calibrations (solid curves) and inferred from the DEM data, BEDMAP2, and Eq.(6) (red hatches).

Figure 6 .
Figure 6.Cumulative thinning since 2001 (shading) and grounding-line position (red contours) in 40-year run from transient calibration.The 2021 and 2031 grounding lines are shown in successive plots with green and brown contours, respectively.

Figure 7 .
Figure 7. Sensitivity of grounded volume (Volume above floatation, or VAF) loss from the domain over the 40-year integration to (a) the sliding parameter β 2 and (b) bed topography R (see Sect. 6.1 above for explanation).(c) Annual cost functions (cf.Eq. 10) for various calibrated model runs.(d) −γ (x) σ (cf.Eq.B3) at the initial and final times in the "time-dep.boundary stress" estimation referred to in (c).Values are plotted relative to the red curve in Fig. 4c.(Note the difference in scale from Fig. 4c.)

D
. N. Goldberg et al.: Ice stream forecast through transient calibration ulation into transient calibration procedure is an important goal, and future efforts should try to achieve this goal with the above limitations in mind.

Figure 8 .
Figure 8.A detailed comparison of modeled grounding lines, the grounding line implied by the data used in the modeling study, and directly observed grounding-line position.The red shaded area represents the portion of the domain which is ungrounded in 2011, inferred from floatation with the 2011 surface DEM and BEDMAP2, and assuming ice and ocean densities of 918 and 1028 kg m −3 , respectively.The blue contour is the modeled 2011 grounding line, and green hatches give the 2011 grounding line position from Rignot et al. (2014).The thin black contour is the computational boundary, and the thick black contour the 1996 grounding line.Note that theRignot et al. (2014) data do not extend to Pope Glacier.
Figure B1.(a)Planform visualization of depth-integrated normal and shear stress along a vertical front or grounding line.In the case of the ice shelf, the stress balance must be solved within the glacier and ice shelf, and stresses along the grounding line depend on this solution.If these grounding-line stresses were imposed along the calving cliff, velocities in the glacier would be the same in both cases.(b) Schematic of representation of boundary stresses through parameters.Shaded cells represent computational domain, and white cells represent area where an ice shelf would be, were it included in the domain.Separate degrees of freedom describe normal and shear stress at each cell face.