1. Introduction
1.1. Ice-flow models
Numerical ice-flow models are widely used to solve problems in glaciology that cannot be solved analytically (e.g. Reference Van der VeenVan der Veen, 1999;Reference HookeHooke, 2005, p. 288). We can use an ice-flow model in combination with any available information about past and present ice-sheet geometry, ice-sheet internal structure and climate variables that, for example, can be determined from ice-penetrating radar (e.g. Reference Conway, Hall, Denton, Gades and WaddingtonConway and others, 1999; Reference Vaughan, Corr, Doake and WaddingtonVaughan and others, 1999), from ice cores (e.g. NorthGRIP Members, 2004) or from glacial- geological reconstructions (e.g. Reference Denton and HughesDenton and Hughes, 2002; Reference StoneStone and others, 2003). A variety of these geophysical and paleoclimatic data have often been collected together near ice-core sites, and ice-flow models can assist in ice-core interpretation by revealing strain history and origin sites for ice in a core. This motivates us to investigate ice-sheet and climate histories in the vicinity of an ice divide. While an ice-core site is chosen because the history of ice flow there is usually simpler to decipher than at other sites on an ice sheet, the ice-divide thickness and the ice-divide location can change due to temporal variations in accumulation rate and ice flow. In addition, spatial variations alter particle trajectories within the ice, and the spatial dimension should be included in any ice-flow model that is used to interpret an ice-core record where ice-divide migration may have been significant. Since climate information that is recorded in the ice has been affected by the history of ice flow (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 650), ice-flow models that capture ice-sheet transients are needed in combination with laboratory analyses of ice-core samples to infer the correct climate history; this is a motivation for finding a solution to the problem of modeling transient ice flow with a limited domain.
Some of the information about ice-sheet history that is sought from data in the vicinity of an ice divide can be inferred by solving an inverse problem. However, computational efficiency is required when solving inverse problems, which can require many iterations of the ice-flow model. Computational efficiency is also necessary when using a higher-resolution model or when running the model over long timescales. Limiting the model domain to include only the relevant portions of the ice sheet is a way to reduce computation time. This approach is also simple, which is an objective of our ice-sheet model, and has a similar motivation to minimal glacier models (e.g. Reference OerlemansOerlemans, 2011). In addition, when the model domain is limited we do not need to make estimates of observable quantities in regions where model-parameter values and boundary conditions are unconstrained, or where data are unavailable. However, limiting the domain of a transient ice-flow model can lead to an ill-posed problem because accurate calculation of the boundary values in the limited-domain model requires additional information to ensure that the limited domain evolves consistently with the full domain within which it exists; this is also a motivation for finding a solution to the problem of modeling transient ice-flow with a limited domain.
In practice, this additional information can be provided to the limited-domain model by embedding the limited- domain model (also known as a regional model) in a full- domain model (also known as a global model). There are at least two approaches to embedding a limited domain in a full domain. In a commonly used approach, the limited- model boundary values are provided directly from calculations performed within a full-domain model where some interval grid interfaces coincide with boundaries of the limited-domain model;this is referred to as a 'nested' model. In our new approach, the limited-domain boundary values are provided from calculations performed within the limited-domain model that rely on the behavior characteristics of the full-domain model, rather than specifically on its detailed evolution.
1.2. Nested ice-flow model
Numerical models of physical processes often use nesting schemes. An example is a regional climate model that is driven by the lower-resolution output of a global climate model (e.g. Reference Christensen and SolomonChristensen and others, 2007, ch. 11; Reference Salathé, Mote and WileySalathe and others, 2007), or a nested ocean-circulation model (e.g. Reference Blayo, Debreu, Chassignet and VerronBlayo and Debreu, 2006). In ice-sheet modeling, higher resolution and/or higher-order physics that are important in specific regions of an ice sheet have been nested in a full ice-flow model that otherwise has coarser resolution and/or simplified physics. For example, several models of ice-sheet evolution have incorporated regions with higher spatial resolution, which are nested in a threedimensional (3-D) thermomechanically coupled model of the entire ice sheet (e.g. Reference Marshall and ClarkeMarshall and Clarke, 1997; Reference FastookFastook, 2005;Reference Huybrechts, Rybak, Pattyn, Ruth and SteinhageHuybrechts and others, 2007, Reference Huybrechts, Rybak, Steinhage and Pattyn2009). These modifications for regions of fast ice flow, for ice-shelf flow, or for regions requiring higher-order physics, result in full ice-sheet evolution that is more realistic, while remaining computationally tenable compared to a full-domain model at the resolution or sophistication of the nested component. In another example, Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (2002) used a full ice-flow model to calculate lateral boundary conditions for their local anisotropic flow model. The local model and the full-domain model were calculated synchronously, so the calculated quantities within both models were always consistent.
1.3. Limited-domain ice-flow model
We define a limited-domain ice-flow model as a model whose spatial calculation domain includes only a limited portion of an ice sheet. In our case, our model domain is limited to the vicinity of an ice divide (Fig. 1). In transient flow, the ice-surface profile is specified as an initial condition. A full-domain model always has a natural zero-flux boundary condition or else a calving condition at its terminus. However, because a limited domain has no natural boundary condition, evolution within the limited domain is an ill-posed problem. The ice flux crossing the arbitrarily defined boundaries must be specified or calculated. Calculating that boundary flux in a way that is inconsistent with the preferred external ice sheet can lead to numerically driven ice-sheet transients and to ice-sheet behavior that is incompatible with that preferred external domain. We make the limited-domain problem well-posed by using impulse-response functions derived from a full- domain model to calculate variations in ice flux due to variations in accumulation rate and ice flow that originate within the limited domain, and we prescribe variations in ice flux due to variations that originate external to the limited domain.
A limited-domain ice sheet is embedded in a full-domain ice sheet. Surface elevation S(x,t) is a function of horizontal position x and time t. In this case, bed elevation B(x) is a function only of x. Ice fluxes q(x,t) at both the left and the right boundaries of the limited domain depend on variations in accumulation rate and in ice flow that originate internal to the limited domain and external to the limited domain.
In this paper, we discuss the general problem of modeling transient ice flow with a limited domain. Our primary focus is on treatment of the boundary conditions. In order to illustrate our points, we must use some specific ice-flow model. We chose a flowband model that uses the shallow- ice approximation (SIA;e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 322) because it is a simple representation of ice flow in the ice-sheet interior. In the Appendix we present the details of our solution for ice-sheet evolution within a limited- domain flowband. We emphasize that our general approach could also be applied to transient ice-flow problems using other ice-flow models of different complexity. Another case is a 3-D region that has natural (no-flux) lateral boundaries (e.g. valley walls or ice-sheet catchment boundaries).
1.4. Ice-flux variations within the limited domain
In the absence of physics that allow bifurcation of solutions (e.g. elevation-precipitation feedback (Reference BodvarssonBodvarsson, 1955) or the tidewater-glacier cycle (Reference PostPost, 1975;Reference Meier and PostMeier and Post, 1987)), a full-domain model should return to its original steady configuration when the mass-balance forcing returns to its original configuration. We also expect a limited- domain model to show the same behavior. We use behavior characteristics of a full-domain model to ensure that ice- sheet evolution in the limited-domain model is consistent with evolution in the full-domain model. First, we embed the limited-domain model in a full-domain model that includes an ice-sheet terminus, so that the boundary condition is known. If the limited domain includes an ice divide, the full domain should include both termini. Second, we characterize the behavior of this full ice sheet by calculating its response to an impulsive perturbation in accumulation rate. Third, we use the impulse-response functions at the limited-domain boundaries to control the flux entering or leaving the limited domain at each time-step in response to actual variations in accumulation rate. This enables the limited-domain model to adjust to any climate change at a glaciologically realistic rate that is compatible with the full ice-sheet model experiencing the same climate changes. We illustrate our general approach to establishing well-posed boundary conditions using the specific limited- domain model and the specific full-domain model described in the Appendix.
1.5. Ice-flux variations external to the limited domain
Ice-flow variations within the limited domain can also be driven by variations in ice flow and climate that originate external to its arbitrary boundaries. For example, ice-shelf loss, or variations in sea level or ice-stream activity can ultimately affect the flow of interior ice (e.g. Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). The need for additional information from outside the limited domain in order to correctly model ice- sheet evolution on a limited domain is a drawback. However, it may be possible to infer those externally forced fluxes by solving an inverse problem using internal-layer architecture within the limited domain; this is a future research direction. In this paper, we show that we can correctly calculate the evolution on the limited domain in the presence of external forcing if we know its impact on the ice flux crossing the limited-domain boundaries.
2. Boundary Conditions for a Limited- Domain Model
The key question that we address in this paper is how to assign appropriate boundary conditions on a limited- domain model. We illustrate the concepts with the simplest form of limited-domain model, i.e. a 1.5-dimensional (1.5-D) flowband.
2.1. Mass conservation
When the limited domain is a flowband, we can calculate ice-thickness evolution by solving the mass-conservation equation (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010, p. 333),
where x follows the center of the flowband, h(x,t) is the ice thickness, q(x,t) is the total volumetric ice flux in the flowband of width W(x), and b(x, t) is the accumulation rate (or ablation rate). Ice thickness h(x,t) is the difference between the surface elevation S(x, t) and the bed elevation B(x). As discussed in the Appendix, we illustrate the concepts using the SIA, and we solve this conservation equation numerically using an implicit approach (e.g. Reference PatankarPatankar, 1980). To solve for ice-thickness evolution, we need to calculate the ice flux q(x,t). The accumulation-rate history ?(x, t) is prescribed, or can be inferred as part of an inverse problem if internal-layer data are available (e.g. Reference Waddington, Neumann, Koutnik, Marshall and MorseWaddington and others, 2007;Reference KoutnikKoutnik, 2009). Without information from the external ice sheet, thickness evolution in a limited domain is an ill-posed problem; using our specific ice-flow model to illustrate our points, we show how to make this a well-posed problem.
2.2. Flux boundary conditions uninformed by external domain
A steady-state mass-balance pattern produces a well- defined steady-state ice-flux pattern. However, flux depends on both ice thickness and slope, so in a limited domain there can always be a trade-off between thickness and slope. A thinner steady solution with a steeper slope (i.e. a smaller external ice sheet) can have the same flux distribution as a thicker solution with gentler slope (i.e. a larger external ice sheet). In a transient model, no extrapolation scheme based on values from inside the limited domain, regardless of its order, can calculate the correct ice flux q(x,t) crossing the limited-domain boundaries in any particular full ice sheet that we might want it to represent. While the calculated flux may adequately represent behavior of some whole ice sheet, we will show that which whole ice sheet it represents is determined by numerical truncation errors from calculations on a discrete grid. To represent a particular whole ice sheet, a correct flux treatment on each limited-domain boundary requires information from both sides of that boundary.
In a limited-domain model, the ice thickness is known only at the initial time-step;as the ice sheet evolves, there is no further explicit constraint on ice thickness at the limited- domain boundaries. Therefore, the boundary treatment must incorporate some information about the span of the full ice sheet. In addition, we know that the ice flux q(x,t) crossing the limited-domain boundaries through time depends on accumulation rate and ice flow both within the limited domain and external to the limited domain; therefore any extrapolation scheme based solely on information from within the limited domain will be inadequate.
Figure 2a shows four points on an initial steady-state ice surface calculated using Eqn (A9). The ice flux crossing the left and right boundaries was calculated at each time-step by extrapolating the thickness and slope from interior points using spatial grids with two different horizontal-grid resolutions, i.e. 300 and 600m. The ice-sheet surface including these four points was then tracked over time under steady- state forcing. Instead of holding the original steady-state profile, the ice sheet drifted toward different steady states (dashed lines in Fig. 2a), depending on the grid spacing used. Figure 2b shows how the ice surface at the four locations changes over time in response to these erroneous changes in the calculated boundary flux. Different hori- zontal-grid resolutions can give different extrapolated values of the boundary flux that lead to different ice-sheet evolution. This evolution is physically possible because, for an ice sheet with an unknown span, there are infinitely many surface profiles that have the same flux profile across the limited domain, and there is no ice-thickness information provided after the initial time;a full-domain model with a known span and ice thickness equal to zero at the terminus does not have this problem.
(a) Solid line shows the initial steady-state surface profile on a limited domain. With boundary conditions that incorporate no information about the ice sheet outside the limited domain, the ice sheet evolves to a different steady-state profile under steady-state climate forcing. The thick and thin dashed lines show the new steady-state surfaces obtained with horizontal grid resolution of 600 and 300m respectively. (b) Elevations of points 1–4 on the initial ice surface in (a) are tracked through time. (c) Evolution of ice flux at the four locations calculated with a grid resolution of 600 m.
Instead of using an extrapolation to calculate the boundary flux, we could estimate the flux crossing the last downstream finite-volume edge with a kinematic calculation. The kinematic flux (Eqn (A1)) is calculated by integrating the continuity equation (Eqn (1)). The accumulation rate ?(x, t) is known, and the rate of ice-thickness change at the current time-step ?(x, ti ) can be estimated from the known value of ?(x, t i) at the previous time-step. Calculating the boundary flux kinematically using Eqn (A1) does allow the model in Figure 2 to hold steady state. However, when real transient forcing occurs, a limited- domain model with kinematic boundary conditions fails to return to its original steady-state profile after the original steady-state forcing is re-established; this is a critical point because, in the absence of any catastrophe or cusp-related behavior, a real ice sheet should return to its original state after the original forcing was re-established. In order for our numerical limited-domain model to be able to return to its original steady state and behave like a real ice sheet, the model must keep track of the original steady-state volume. Equation (A1) allows the volume to change through time without any memory; in the model there was no accounting of how the ice volume changed from the original state. Modifying Eqn (A1) by dropping the term in ?(x, t) allows the model to retain information about the initial steady state. By instantly exporting any additional volume that enters the domain in a given time-step through transient accumulation, its volume never changes. However, instantaneous export is not glaciologically realistic. Therefore, kinematic-flux calculations are also inadequate to represent flux crossing the limited-domain boundaries. These simple tests highlight the challenge of solving for transient ice flow on a limited domain and demonstrate that extrapolating boundary conditions to limited-domain boundaries without regard for properties of the full-domain ice sheet in which it exists, can produce ice-sheet evolution that is incompatible with that external ice sheet.
2.3. Flux boundary conditions informed by external domain
In a well-posed ice-sheet evolution problem, an ice sheet returns to an initial steady state when forcing terminates. By controlling the ice thickness, slope and velocity at the limited-domain boundaries, a full-domain ice sheet controls the volume of ice leaving the limited domain at any time. This suggests that in order to have a well-posed limited- domain problem, we must also control the flux that leaves the domain in a glaciologically reasonable way, driven by climate variations and ice-flow changes both within the domain and exterior to the domain. In this paper, we focus on calculating the flux change at the limited-domain boundaries due to variations in accumulation rate and ice flow that originate within a limited domain.
2.3.1. Temporarily embed the limited-domain model in a full-domain model
First, we temporarily embed our limited-domain model in a full-domain model that provides a more generic context than traditional nesting schemes because it is used in order to calculate impulse-response functions that characterize this full domain. The shape of the full-domain model is matched to the shape of the limited-domain model over the horizontal extent of the limited domain.
To illustrate these necessary boundary-condition concepts, we use a specific numerical solution for a steady-state surface profile to extend our limited-domain model (Eqn (A9));details are provided in the Appendix. To calculate the full ice-sheet shape, we must prescribe a mass-balance pattern outside the limited domain. In Figure 3a, the prescribed accumulation-rate pattern across the limited domain was uniform. We will use a realistic estimate of mass balance outside the limited domain if it is available. Here we apply a simple mass-balance pattern with a zone of uniform accumulation rate c and a zone of uniform ablation rate a (see Appendix), which follows the steady-state analytical model by Reference PatersonPaterson (1972) so that we can validate our steady-state numerical model. The uniform ablation rate over the new portion of the domain has a rate determined by a ratio of accumulation rate c to ablation rate a, which we chose to be c/a = 0.2. In this case, the average accumulation rate is 20 cm a-1 and the average ablation rate near the terminus is 1 m a-1.
(a) The prescribed accumulation-rate and ablation-rate pattern used to generate a full-domain surface profile by solving Eqn (A9). The solid black line shows the known mass-balance rate for the limited domain, and the dashed line shows the estimated values across the extended domain. (b) Ice-surface profile for a limited domain that crosses an ice divide (solid line), and for the corresponding full domain (dashed line) within which it is embedded. The bed topography is flat. The dot at the right-side boundary of the limited domain marks the spatial position used in subsequent figures.
2.3.2. Find response functions for the full-domain model
Second, we numerically calculate the nonlinear evolution of the full-domain ice sheet in response to an impulsive perturbation in accumulation rate. The ice sheet is initially in a steady state described by thickness h
0(x) and mass-balance pattern 
. Reference NyeNye (1965) linearized Eqn (1) and discussed the idea of an impulse-response function h
1 (x, t – t
1) for ice thickness that described the subsequent deviations in ice thickness along the full-domain surface, from the application of the impulsive perturbation at time t1
, until the surface had (effectively) returned to steady state. Reference NyeNye (1960) showed that ice-thickness response h1 (x, t – t1
) and ice-flux response q1(x, t – t1
) could be related using kinematic-wave theory. In the Nye linear perturbation theory, the impulsive perturbation to the accumulation rate was spatially uniform. However, actual variations in accumulation rate may be non-uniform; below, we use the full nonlinear evolution equation (Eqn (1)) to explore the sensitivity of ice-sheet evolution at the limited-domain boundaries to spatial changes in accumulation rate that may be non-uniform. The accumulation-rate pattern for the impulse-response calculation is given by the sum of the steady-state ('datum') pattern and a deviation 
in accumulation rate from the datum state at time t
We numerically find the impulse-response function for the impulse described by
at a time t 1, where the spatial function bδ is the ice- equivalent thickness (m) added by the impulse δ(t-t 1) (a Dirac delta function; e.g. Reference Arfken and WeberArfken and Weber, 1995, p. 81). The total volume added to the left side of the domain at time t1 by this impulse in Eqn (3) is
Equivalent expressions exist for the right side of the domain. After receiving this impulsive addition of volume, the ice sheet must evacuate an equal volume of ice over time in order to return to its original steady state. Time is discretized in units of Δt, so the average rate of addition of new volume to the left side of the limited domain in the time-step Δt that includes time t
1 is 
given by
In general, our impulse-response functions F
L(t-t
1) for ice flux at the left side of the domain and F
R(t-t
1) at the right side of the domain are causal time-delay filters (e.g. Reference GubbinsGubbins, 2004) determined by the ice-sheet response from the time t
1 when the perturbation was applied to the time t
1 + T
L or t
1 + T
R when the ice thickness at the boundary returns to within a specified threshold of the original steady- state value. For example, our chosen thickness threshold was 10-4bδ. Each impulse-response function is then normalized so that it integrates to unity over its length 
In Section 3.1.2, we apply these concepts to two types of mass-balance perturbations and each type will have its own impulse-response functions. The first type evaluates response to spatially distributed accumulation-rate changes over the domain, while the second type evaluates flux changes due to divide migration. There is a third term that is prescribed to account for changes in ice flux through the limited-domain boundary that are caused by processes external to the limited domain, such as ice-shelf loss, sea- level change or ice-stream activity.
2.3.3. Boundary values for the limited-domain model
Third, we calculate the time-dependent ice fluxes Q L(t) and Q R(t) at the left and right limited-domain boundaries in response to the impulse in Eqn (3) where
Steady-state ice flux q 0(x) in a flowband at time t 0 is given by
where q 0(x 0) is the initial ice flux at any position x 0. When x in Eqn (7) is at a limited-domain boundary, we calculate ice- flux deviations from this known steady-state value at later times.
The time-dependent ice flux Q L(t) at the left limited- domain boundary (Eqn (6a)) comprises this steady-state flux q 0(x L) and three flux-variation terms, which are introduced below:
Again, we present only equations for the left side of the domain at x
L. However, equivalent equations apply on the right boundary at x
R. The first flux variation 
accounts for variations in accumulation rate and ice flow that originate within the limited domain, in response to a time series of volume changes B
L(t) in Eqn (4), filtered through an impulse-response function 
. At the left boundary, the flux response to a generic history of volume changes B
L(t) can be expressed by the convolution
where the notation (g*f)(t) denotes the convolution of a function g(t) with a filter f(t). This impulse-response function characterizes the response, at the left limited- domain boundary, to an accumulation-rate impulse affecting the full domain; however, its role is to evacuate from the limited domain just the volume B
L(t) (Eqn (4)) that accumulates within the limited domain itself at each time t
1. It evacuates each volume B
L(t
1) through the boundary x
L over time, starting at time t
1, when volume B
L(t
1) was deposited, and finishing at time 
later. Considering the contributions at time t from all previous times, the integral in Eqn (9) starts at 
because volume perturbations B
L(t
1) from all previous times 
have been completely expelled by time t. The surface in the limited domain can return to its initial steady-state configuration when the climate returns to its original steady state, because all additional volumes due to mass-balance variations, and only those additional volumes, have been exported. A similar expression holds for the right-side boundary. The numerical solution of the convolution in Eqn (9) at time-step k can be written as
where 
is the number of points in the impulse-response filter.
The second flux-variation term 
is nonzero if at some time t > t
0 the divide has moved to a new position x
div(t). If that new position is to the left of the steady-state divide position x
div(t
0), then some of the steady-state accumulation rate that initially contributed to ice flux through the left boundary is now directed toward the right boundary of the limited domain. In a manner similar to Eqn (9), this loss of ice input to the left of the divide can be expressed as a volume contribution and represented as an average rate Ḋ
L(t) during the time-step Δt at time t by
An equal and opposite ice volume per unit time Ḋ R(t) is now directed toward the right boundary, i.e. since flux is taken to be negative on the left side of the domain and positive on the right side of the domain Ḋ R(t) = Ḋ L(t). In a manner similar to Eqn (10), the flux response at the left boundary at time t following a history of ice-divide variations can be expressed by the discrete convolution
where 
is the appropriate impulse-response function.
The third flux-variation term 
accounts for ice- sheet evolution that is driven by factors outside the limited domain. This time series must be prescribed, or else inferred as part of an inverse problem.
3. Results
The flow model within a limited domain is in general fully nonlinear, so the accuracy of the solution depends on the suitability of the impulse-response functions that provide the boundary conditions. These boundary conditions are linear in the net accumulation-rate variations within the limited domain. The full-domain solution and the limited-domain solution should be equivalent, and we show that this can be achieved. In the tests presented here, to illustrate our concepts, we use a limited-domain ice sheet with an initial ice thickness of 1000 m at the left-side boundary (~1125 m initial maximum ice thickness), a flowband width that is uniform, and a domain that crosses an ice divide with a limited-domain length of 25 km;this is the surface profile shown in Figure 3b. The mean accumulation rate is 20cma-1. The initial steady ice flux leaving the left-side boundary is -2500 m3 a-1 per meter width. The ice divide is defined as the location with the highest surface elevation.
In an impulse-response calculation, an accumulation-rate perturbation, as in Eqn (2), has a small magnitude and is typically assumed to have a spatially uniform distribution (e.g. Reference NyeNye, 1960). While an impulse-response function can well characterize the response of an ice sheet to spatially uniform variations in climate (e.g. Reference HookeHooke, 2005, p. 373-375), actual accumulation-rate variations may not be spatially uniform. Furthermore, an impulse-response function is associated with a specific ice-sheet geometry, and this geometry may change over time. We must quantify the ability of a given impulse-response function to characterize ice-sheet behavior as the limited-domain model experiences climate variations and evolves to different ice-sheet geometries.
In our formulation of boundary conditions for limited- domain models, four different impulse-response functions 
are required to control ice- flux variations across the boundaries of the limited domain from accumulation-rate variations (Eqn (9)) and from ice- divide migration (Eqn (11)). We illustrate the sensitivity of the solution for ice-sheet evolution to different representations of the four impulse-response functions in this problem, using the particular limited-domain model described in the Appendix. For the calculations of ice-sheet evolution below, we set the external-flux forcing term ΔQ
ext(t) = 0 (Eqn (7)) on both sides of the limited domain.
3.1. Sensitivity to spatial extent of an accumulation- rate perturbation
3.1.1. Full domain vs limited domain only
The spatial extent over which the accumulation-rate variations are distributed will affect the response of the ice sheet. For example, accumulation rate may change uniformly over the entire ice sheet, or only over the limited domain. Figure 4a shows the limited-domain boundary response functions associated with a uniform perturbation that extends only over the limited domain, while Figure 4b shows the boundary response functions associated with a uniform perturbation that extends over the full domain; all response functions correspond to the right-side boundary of the limited-domain ice sheet in Figure 3b. While an accumulation-rate perturbation that extends only over the limited domain is improbable, Figure 4 demonstrates the sensitivity of the impulse-response function to the spatial extent of the perturbation.
(a) Ice-flux impulse-response function for an impulsive uniform accumulation-rate perturbation that extends spatially only over the limited domain. (b) Ice-flux response function for an impulsive uniform accumulation-rate perturbation that extends over the full domain. Response functions in (a) and (b) correspond to the right-side boundary of the ice sheet in Figure 3b. (c) Ramped step changes in accumulation rate of 5% from the steady-state value of 20 cma–1. (d) Divide-thickness response following accumulationrate changes in (c) in limited-domain model, using filters in (a) and (b) compared with full-domain model.
The boundary ice-flux response function in Figure 4a differs from the response function in Figure 4b because the response function depends on the spatial extent of the impulsive perturbation. Figure 4d shows that correct ice- sheet evolution can be achieved if the correct response function is used in the limited-domain model calculations in response to the step change in accumulation rate shown in Figure 4c. If actual accumulation-rate perturbations extend only over the limited domain, the response functions in Figure 4a will be appropriate. However, in reality, we expect that actual spatial and temporal variations in accumulation rate will occur over the full span (or a significant portion of the full span) of the ice sheet, and will not be restricted to the arbitrary extent of our limited domain. In the subsequent sensitivity tests and results shown here, we use impulse-response functions associated with accumulation-rate perturbations that extend over the entire span of the full-domain model.
3.1.2. Sensitivity to narrow vs broad perturbations
The spatial pattern of the impulsive perturbation in accumulation rate over an ice sheet will affect the response time. While Reference Johannesson, Raymond and Waddingtonjohannesson and others (1989) showed that spatial variations in the accumulation-rate perturbation have an insignificant effect on the total change in volume, we show that spatial variations in the accumulation-rate perturbation can significantly affect the impulse-response function at the limited-domain boundaries, and actual variations in accumulation rate across the limited domain may not be spatially uniform. We compare the impulse-response functions associated with accumulation-rate perturbations across the full domain that are (1) spatially uniform across the full domain, (2) linearly varying across the domain, and (3) strongly peaked at the divide. A delta function can be represented as the limit of a Gaussian distribution (e.g. Reference Arfken and WeberArfken and Weber, 1995, p. 81),
where α is the width of the Gaussian distribution, and a ! 0 in the definition of the delta function.
We evaluate Gaussian functions over the span of the full- domain ice sheet from the divide to the terminus L (m) with α = 0.5L and α = 0.1 L (by definition, α = 0.1 L gives a narrower peak). Figure 5a shows the shapes of the four accumulation-rate perturbations, Figure 5b shows the impulse-response functions and Figure 5c shows the step response at the limited-domain boundary. Previous glacio- logical analyses use the step response (e.g. Reference NyeNye, 1960), where integration by parts shows that the step response is the temporal integral of the impulse-response function (e.g. Reference GubbinsGubbins, 2004). Each accumulation-rate perturbation in Figure 5a adds the same total ice volume impulsively over the full domain. By definition, the impulse-response function integrates to unity over the full response time. The response time τir is an e-folding time, which is the time to reach (1 -e-1) of the new equilibrium value. Figure 5c shows that the re-equilibration of the ice sheet to a narrow spike in accumulation rate is very different from the re-equilibration of the ice sheet to a uniformly distributed accumulation-rate perturbation. In Figure 5d we show the evolution of surface elevation at the limited-domain boundary in response to the accumulation-rate history in Figure 4c. Ice-thickness response is influenced by the response function used to calculate the ice-sheet evolution. When the actual accumulation-rate perturbation is spatially uniform, but an impulse- response function for a spatially restricted spike perturbation is used to characterize ice-sheet behavior, ice added to the domain is evacuated too quickly.
(a) Accumulation-rate perturbations that are spatially uniform (dashed line), linearly varying (bold line), a Gaussian function (Eqn (13)) with α = 0.5L (thin line), and with α = 0.1L (gray line) where L is the half-span of the full domain, from the divide to the terminus; these are examples of the perturbation term in Eqn (2), and all perturbations add the same ice volume to the full domain at one time-step. The vertical line marks the limited-domain boundary on the right side of the limited domain (as in Fig. 3b) where the impulse-response functions are evaluated. (b) Boundary impulse-response functions corresponding to the perturbations in (a). The impulse-response time τIR is the e-folding time for the step-response curves in (c), given by the integrals of the impulse-response functions. (d) Ice-sheet evolution at the limited-domain boundary (as in Fig. 3b) from the limited-domain model in response to uniformly distributed accumulation rate in Figure 4c using the four response functions in (b).
In reality, it is unlikely that the spatial pattern of accumulation-rate perturbations across the limited domain will be as spatially restricted as a spike function. In addition, ice-sheet evolution is determined primarily by long-term spatially averaged changes in accumulation rate, and not by localized, even if large in magnitude, excursions in the accumulation rate. Therefore, to calculate the impulse- response function , in Eqns (9) and (10), we use a spatially uniform accumulation-rate perturbation, i.e., in Eqn (3), 
However, changes in ice input directed toward one boundary or the other due to ice-divide migration are better characterized by a localized function 
in Eqn (3). Therefore, to calculate the impulse-response functions 
, in Eqn (12), we use an accumulation-rate perturbation that is a peaked function with α = 0.1 L (Eqn (13)) centered at the initial divide.
3.2. Sensitivity to ice-sheet geometry
Reference Johannesson, Raymond and WaddingtonJóhannesson and others (1989) showed that the volume response time τV for an ice sheet could be estimated as the ratio of the maximum ice thickness to the ablation rate at the terminus. Therefore, the impulse-response function is dependent on the ice thickness and is associated with specific ice-sheet geometry. For example, our impulse-response functions should produce longer response times for thicker ice sheets. Since the ice thickness and the ice-divide position can change in a transient problem, we must take into account how this will affect the impulse-response functions used to characterize ice-sheet behavior.
Figure 6b shows the impulse-response functions, and Figure 6c shows the step response at the right-side boundary of the limited domain (marked in Fig. 6a), for the ice sheets in Figure 6a with initial ice-surface elevation at the left-side limited-domain boundary of 1000, 800 and 1200 m. We calculate the evolution of our standard ~1000 m ice sheet in response to the accumulation-rate history in Figure 4c, but we use impulse-response functions for an ice sheet that is 200 m thicker (~1200 m), and for an ice sheet that is 200 m thinner (~800m), than the standard ice sheet shown in Figure 6d. All solutions are similar.
(a) Three steady-state surfaces with elevations of 800m (bold line), 1000m (dashed line) and 1200m (thin line) at the left edge of the domain (at –20km along flowband). The dots mark the right boundary, where the time-varying values in the lower panels are obtained. (b) Impulse-response functions and (c) step response to a uniform accumulation-rate perturbation across a full-domain model for the three ice-surface elevations from (a). The impulse-response time τIR is the e-folding time. (d) Surface evolution at the right boundary of a limiteddomain model with surface elevation of ~1000 m, in response to uniformly distributed accumulation-rate changes in Figure 4c, using the three response functions in (b). The dashed line is the correct solution.
3.3. Efficient transient calculations
To accurately calculate ice-sheet evolution with a limited- domain model, we must prescribe the accumulation-rate history ?(x, t) over the limited domain, and we must prescribe the histories ΔQeL ext(t) and ΔQR ext(t) of externally forced changes in ice flux on the limited-domain boundaries. The accumulation-rate and external-forcing histories, together with the impulse-response functions, contain all of the information needed to calculate ice-sheet evolution within a limited domain. If the correct impulse-response functions are used at the limited-domain boundaries, the ice sheet can thicken and thin in exactly the same way as a full- domain model.
For computational efficiency, we use a single set of impulse-response functions throughout the calculation for ice-sheet evolution. While these simple limited-domain calculations using response functions are at least five times faster to compute than the same full-domain calculation, there would be less of a computational advantage if the response functions were recalculated more frequently. Given other uncertainties in the problem, and given our objective of an efficient calculation, this is a reasonable approach. Therefore, we assume that an impulse-response function generated with a spatially uniform accumulation- rate perturbation is appropriate to control ice-flux variations due to accumulation-rate variations and we use an impulse- response function generated with a narrow Gaussian- function accumulation-rate perturbation to control ice-flux variations at the boundaries associated with ice-divide migration. With these assumptions, a limited-domain model may not evolve in exactly the same way as a full-domain model, but it will evolve in a physically similar way compared to the full-domain model. In general, deciding when to recalculate impulse-response functions so that they suitably represent ice-sheet response as the ice sheet evolves will be a problem-specific consideration.
If an impulse-response function associated with a spatially uniform perturbation is used to obtain boundary conditions for ice-volume variations that are not strictly uniform, the ice sheet will thicken and thin in a way that is similar, but not identical, to a full-domain model. In Figure 7a, the accumulation-rate history varies in space and time across the limited domain, but remains at the steady- state values elsewhere. We chose this accumulation-rate history in order to isolate how changes in the spatial pattern of accumulation rate influence a limited-domain solution found with impulse-response functions that correspond to a spatially uniform accumulation-rate perturbation; in this case, the impulse-response functions were not updated to account for changes in the actual pattern of accumulation rate or ice thickness over time. Figure 7b shows the ice-surface evolution, and Figure 7c shows the ice-divide migration found with the limited-domain model and found with a full-domain model. The solutions from the limited-domain model and from the full-domain model are not exactly the same because the impulse-response functions used in the limited-model solution correspond only to a uniform perturbation in accumulation rate, whereas the actual accumulation-rate perturbations are spatially variable. This error is minor compared to uncertainties in the actual accumulation rate and exter- nal-flux forcing histories.
(a) Accumulation-rate history that varies in space and time across the limited domain but remains constant outside the limited domain. (b) Ice-surface evolution at the right-side boundary of the limited domain (as in Fig. 3b) from the full-domain model (dashed line), and from the limited-domain model (solid line) using impulseresponse functions that assumed uniform changes outside as well as inside the limited domain. (c) Ice-divide position from the limiteddomain model (solid line) and from the full-domain model (dashed line).
4. Discussion
4.1. Nonlinear approach
Impulse-response functions have traditionally been associated with studies in which the ice-flow equations were linearized (e.g. Reference NyeNye, 1960, Reference Nye1965;Reference HindmarshHindmarsh, 1996; Reference Nereson, Hindmarsh and RaymondNereson and others, 1998). Results were therefore applicable only to small changes in surface elevation or ice flux. Those impulse-response functions described the linearized evolution of the entire glacier profile following a small but uniform mass-balance impulse.
Our approach is different; the full nonlinear form of the equations is always used, both for a full-domain model when determining the impulse responses at the limited- domain boundaries, and for the limited-domain model when solving for the surface and ice-flow evolution within it. We assume linearity only in how and when the boundary conditions release the variations in ice volume that collect inside the limited domain (e.g. if twice as much extra ice collects in a given year, it will be released through the boundaries on the same timetable, but at twice the amount in each subsequent time-step). While the impulse-response functions are found numerically as the response of a nonlinear model to small perturbations, the (nonlinear) model itself is not restricted to small changes in accumulation rate, ice thickness or flux.
4.2. Estimating external forcing
A realistic ice sheet will experience variations in ice flow and climate across its entire extent (not just in a limited portion near the ice divide), and those changes that originate external to the limited domain may dominate the long-term evolution of the ice sheet. For example, on glacial- interglacial timescales, ice-sheet interiors respond to externally forced changes in global ice volume. In particular, ice-sheet margins can respond directly to changes at the ice/ ocean boundary by advancing or retreating, and this affects ice thickness and ice-divide position in the interior (e.g. Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). Since the climate history prescribed on the limited domain has no information about externally forced changes in ice flux at the limited-domain boundaries (beyond the impact of a spatially uniform b1), variations in external-flux forcing ΔQL ext(t) and ΔQR ext(t) must be prescribed.
Figure 8a shows an accumulation-rate history over the full extent of the ice sheet that changes through time from a spatially uniform pattern, to a pattern that has a strong spatial gradient across the divide; the variation across the limited domain is the same as in Figure 7a. In this case, the accumulation-rate variations outside the limited domain are an external forcing that also drives divide migration. Figure 8b shows the ice-thickness evolution at the right-side boundary (as in Fig. 3b) and Figure 8c shows the ice-divide migration for a full-domain model, for a limited-domain model using the incorrect external forcing, ΔQext(t) = 0, and for a limited-domain model using the correct external forcing. Figure 8b and c confirm our expectation that the impulse-response functions alone cannot provide enough information to facilitate accurate ice-sheet evolution when there are externally forced ice-flux variations. Even in this simple case with relatively small changes in accumulation rate outside the limited domain, additional information that can describe the correct variation in external forcing is required for the limited-domain model to produce the correct ice-sheet history.
(a) Accumulation-rate history that varies in space and time across the full domain; the non-uniform variations outside the limited domain produce external forcing. This history is used to calculate ice-sheet evolution in the full-domain model; the portion covering the limited domain (bounds shown with white lines) is used to calculate ice-sheet evolution in the limited-domain model, with impulse-response functions that assumed uniform variations outside the limited domain. (b) Ice-surface evolution over time at the right-side boundary of the limited domain (as in Fig. 3b) calculated with the full-domain model (dashed line), the limiteddomain model with no information about external forcing (solid line) and the limited-domain model with the correct external-flux forcing (gray line). (c) Ice-divide position. The solutions from the limited-domain model and from the full-domain model are the same when the correct external-flux forcing ΔQ ext (t) is prescribed.
The correct value of ΔQ ext(t) at the limited-domain boundaries over time may come directly from a separate full-model calculation (e.g. output from a 3-D model), or may be estimated using a proxy for externally forced perturbations (e.g. changes in sea level). However, it may be that the history of external forcing is largely unknown. We suggest that in future work, it may be possible to infer the change in ice flux due to external forcing at the limited- domain boundaries as part of an inverse problem, if internal- layer architecture has retained that information. Specifically, an efficient limited-domain model for transient ice flow, with well-posed boundary conditions following the procedures described here, could be used to infer histories of ice thickness, of ice-divide position, of accumulation rate and of external forcing that are consistent with internal-layer architecture and the modern ice sheet (Reference KoutnikKoutnik, 2009).
5. Conclusions
A limited-domain model can efficiently and realistically calculate transient ice flow near an ice divide; we have demonstrated that this can be a well-posed problem. There are two key insights that promote efficiency and accuracy in this problem. First, rather than calculating ice-sheet evolution using a limited-domain model that is nested in a full- domain model that must be run for the entire duration of the experiment, we embed the limited-domain model only at specific times in the calculation. Second, we use the full- domain model only to provide information about characteristic behavior of the full ice sheet, not to provide the exact value of the required boundary conditions. We characterize the behavior of a full ice sheet using impulse-response functions calculated from a full-domain model. While our results directly support modeling ice flow in the vicinity of ice-core sites near ice divides, our approach using response functions to set appropriate boundary conditions on a limited-domain model may also apply to other problems where the response at the boundaries can be characterized using either a simple function, or a separate model, or available data.
We have shown how the impulse response at the boundaries of a limited-domain ice sheet depends on the spatial extent and on the spatial pattern of the accumulation- rate perturbation, in addition to the ice-sheet geometry. Therefore, in the calculation of ice-sheet evolution within a limited domain, different boundary impulse-response functions are needed to distinguish ice-sheet response to variations in accumulation rate from response to ice-divide migration. For variations in accumulation rate we calculate the impulse response of a nonlinear full ice-sheet model to a spatially uniform accumulation perturbation across the full domain. For ice-divide migration, we calculate the response of the same model to an impulsive and spatially restricted perturbation in accumulation rate centered near the divide. The response histories at the limited-domain boundaries are extracted to form the left and right boundary impulse- response functions. In this paper, we have illustrated how these specific response functions can well characterize ice- sheet evolution in a limited-domain model. Future work related to a specific portion of Greenland or Antarctica could explore the sensitivity of the limited-domain model results to different calculations of the impulse-response functions depending on the spatial variation in the accumulation perturbation and changes in ice-sheet response over time.
The response functions can produce realistic ice-sheet evolution only for variations in accumulation rate that originate within the arbitrary bounds of the limited domain, or are uniform across the entire ice sheet. Because long-term evolution near an ice-sheet center can be strongly influenced by variations in climate and ice flow over the entire ice sheet, accurate ice-sheet evolution with a limited- domain model requires an accurate estimation of externally forced changes in ice volume. An estimate of this external forcing could come from an independent 3-D model (similar to a nesting scheme), or we could solve an inverse problem to find histories of accumulation rate, ice thickness, ice- divide position and external forcing that are consistent with internal-layer architecture and the present-day ice sheet (Reference KoutnikKoutnik, 2009);this is a focus of our future work. By making a transient limited-domain ice-flow model well-posed, we have created an efficient approach to an inverse problem that can use data across a limited portion of the ice sheet, and would otherwise be computationally expensive because the forward problem must be run many times to solve the inverse problem.
Acknowledgements
This work was supported by US National Science Foundation grant NSF-ANT0440666, and by NSF International Research Fellowship Program (IRFP) award 0853407 to M.R.K. We thank Jim Fastook and an anonymous reviewer for comments, and Scientific Editor Ralf Greve for effective handling of our manuscript.
Appendix: Limited-Domain Model
To illustrate our points about appropriate boundary conditions, we choose to use a simple flowband model, which is a 1.5-D representation of the ice surface and a 2.5-D representation of ice flow, where the boundaries are defined by two nearby flowlines and the volume is bounded vertically beneath these flowlines. Figure 9 illustrates the geometry of a typical flowband. Variations in the bed topography B(x) and in the flowband width W(x) along the flowband are specified. We assume that glacialisostatic, tectonic or erosional processes do not change the bed topography over time. The ice thickness h(x,t) is related to the ice-surface elevation by h(x,t) = S(x,t)–B(x). When the initial ice surface at time t
0 is estimated using a steady-state surface calculation, the required boundary conditions are (1) the ice-surface elevation S
0(x
0,t
0) at one point x
0 at the first time-step, (2) the ice flux q
0(x
0,t
0) entering or leaving the domain at one boundary x
0 at time t
0, (3) the spatial and temporal accumulation rate, , and (4) externally forced changes in ice flux; all of these values must be known, estimated or solved for as a part of an inverse problem.
By integrating Eqn (1) from x 0 where ice flux is specified, the ice flux at the end of the domain x end can be represented kinematically by
where q 0(x 0 ,t) is the history of ice flux entering or leaving the flowband domain, ḣ(x,t) is the rate of change in ice thickness, and ṁ(x,t) is the basal melt rate.
Dynamically, the flux of ice passing through a crosssectional area W(x) ×h(x,t) at any point x and at any time t, is related to the depth-averaged horizontal velocity 
in that cross section by
We can calculate using the SIA (e.g. Cuffey and Paterson, 2010, p. 322). The SIA is a simplifying assumption that can be applied in cases where the ice thickness h(x,t) is much smaller than any horizontal length scale L over which ice-sheet properties such as bed topography, ice thickness, velocity and sliding fraction change significantly; the SIA is applicable where derivatives of velocities and stresses with respect to x are generally much smaller than derivatives with respect to z. The constitutive relationship for ice flow (Glen, 1955) using the SIA is
where 
the simple-shear strain rate along a horizontal plane, A(T(x,z,t)) is the temperature-dependent ice-softness parameter, τ
xz
is the component of the shear-stress tensor acting horizontally along a horizontal plane, and we choose the flow law exponent n = 3 (e.g. Cuffey and Paterson, 2010, p. 58). The shear component of the strain-rate tensor along a horizontal plane is
Following the SIA, the derivatives of velocities with respect to x are negligible, giving 
(e.g. Cuffey and Paterson, 2010, p. 303). Using the flow law given by Eqn (A3) for the SIA, and assuming that the temperature is uniform with depth for each position in x, T(x,z,t) = T(x,t), the depth-averaged horizontal velocity can be found by integrating ∂u/∂z twice over depth z,
where ρ is density, g is gravitational acceleration, S(x, t) is the ice-surface elevation, and h(x, t) is the ice thickness.
Geometry of an ice-sheet flowband with a limited domain. The change in flux from steady state must be calculated on the left side q(x, t) and on the right side q(x R, t) of the domain (as in Fig. 1), and the spatial and temporal history of accumulation rate ḃ x,t is prescribed as a boundary condition. For this model, the bed topography B(x) and the width function W(x) do not change in time.
If the ice temperature is not uniform with depth at each location, we can solve for an effective isothermal softness parameter à (x, t) that gives the same depth-averaged velocity and ice flux as a depth-varying temperature dependent softness parameter A(x, z, t). We calculate the effective isothermal value by equating the depthaveraged ice velocity with the actual softness profile A(T(x, z, t)),
to the depth-averaged velocity from Eqn (A5), and solving for A(T(x,t)), given as here as à (x, t),
where 
is a nondimensional height above bedrock, defined by
Using Eqn (A5) for the depth-averaged velocity, but with the effective isothermal softness parameter from Eqn (A7), and then by representing average velocity in terms of ice flux q(x, t) and ice thickness h(x, t) using Eqn (A2) and h(x, t) = S(x, t) – B(x), we can rearrange this equation to formulate a nonlinear ordinary differential equation for a steady-state ice surface,
We use this ice-surface profile S0(x,t0) as the initial condition to solve for the ice-thickness evolution h(x,t) that we find by solving Eqn (1); therefore, all the values used to solve Eqn (A9) apply at the initial time, t = t0. While our study is in the vicinity of an ice divide, and while the SIA may inadequately represent the velocity field near the divide, the SIA can adequately describe ice-surface evolution and ice-divide migration (e.g. Hindmarsh, 1996). However, a different representation of the velocity field could also be used. For example, to capture ice-divide flow we could use a full Stokes model, or we could use velocity shape functions based on full-stress flow solutions (Nereson and Waddington, 2002). In addition, we could couple our mechanical model to a thermal model and investigate the response to temperature forcing; because that is not the focus of this work, we choose to use a simple SIA model to demonstrate application of the well-posed boundary conditions.
We use the finite-volume method (FVM; e.g. Patankar, 1980; Versteeg and Malalasekera, 2007), also known as the control-volume method, to discretize our domain in order to find a numerical solution to Eqn (1). In the FVM solution, ice thickness is evaluated at each volume center, and ice flux is evaluated across each volume edge.
Following Patankar (1980, p. 57), we use a fully implicit scheme, and we have verified that our fully implicit solution matches an appropriately time-stepped fully explicit solution. In addition, we invoke under-relaxation (e.g. Patankar, 1980, p. 67) to stabilize our procedure to solve a nonlinear problem that has been linearized between iterative calculations of updates to the solution values.
The transient problem is nonlinear because the ice flux q (x,t) in Eqn (1), calculated dynamically using Eqn (A2) with Eqn (A5), is a nonlinear function of ice thickness h (x,t) and surface slope dS/dx; the ice thickness and surface slope are the values we are trying to find. To address this nonlinearity, we use an iterative procedure and we stop iterations when changes in the solution are smaller than a threshold value (e.g. Patankar, 1980, p. 47; Waddington, 1982, p. 239; Van der Veen, 1999, p. 226); here we use a threshold value of 10–6 m.
Full-domain model
As sketched in Figures 1 and 3b, the limited-domain ice sheet is embedded in a full-domain ice sheet. Our limited-domain model crosses an ice divide, and therefore must be embedded in a full domain that extends off both sides of the divide. Any ice-flow model could be used to embed (extend) the limited domain, but here we use the numerical calculation for a steady-state ice surface given by the solution to Eqn (A9). In this model we can specify variations in bed topography and in mass-balance rate across the entire domain. However, these values may be unknown. If we do not have a good estimate of the bed topography and the flow patterns, we can uniformly extend the bed elevation and the flowband width across the full domain using their values at the limited-domain boundaries. If we have a realistic estimate of the mass balance outside the limited domain then we can use these values. Here we chose a mass-balance distribution given by a zone of uniform accumulation rate c over part of the domain, and a zone of uniform ablation rate a over the rest of the domain, with the two zones separated at the equilibrium line R; this follows the analytical solution for an ice-sheet surface given by Paterson (1972) and allowed us to compare our numerical model to this analytical model. We prescribe a ratio of accumulation rate to ablation rate of c/a = 0.2. However, this could be any ratio, and in different applications a large ablation rate may be most realistic (e.g. a calving margin in Antarctica). For this mass-balance pattern, the solid curve in Figure 3b shows a limited-domain surface that crosses an ice divide with a maximum thickness of ~1125m and a limited extent of ~25 km. The dashed curve in Figure 3b shows the extension of this limited surface to a full domain that includes both ice-sheet termini and has a full extent of ~70 km.
In a model that includes an ice-sheet terminus, the ice flux and the ice thickness approach zero at the terminus. This means that the velocity at the terminus will also approach zero if we follow Eqn (A5). As pointed out by Nye (1960, 1963a,b), this results in the non-physical situation of an immobile terminus; in reality the terminus should be able to advance and retreat. To address this problem, the region near the terminus can be replaced by a wedge with a defined angle to the surface (e.g. Nye, 1963a,b; Waddington, 1982, p. 247). Mass conservation and the flow law are satisfied separately in the wedge, which we represent by one gridcell at the end of the domain. While Lam and Dowdeswell (1996) suggested that an adaptive-grid scheme should be used at the terminus in order to best represent ice-sheet behavior in this type of numerical model, we are solving for the full-domain response to small perturbations only. As a result, the length changes are small, and our simple terminus treatment adequately captures terminus variations.
