Boundary layer models for calving marine outlet glaciers

We consider the flow of marine-terminating outlet glaciers that are laterally confined in a channel of prescribed width. In that case, the drag exerted by the channel side walls on the glacier affects extensional stress at the grounding line, and as a result, the flux through the grounding line. Lateral drag in turn is affected by the length of the floating ice shelf when the latter is present, and therefore by calving. Using two calving laws, one due to Nick et al based on a model for crevasse propagation due to hydrofracture, and the other simply asserting that calving occurs where the glacier ice becomes afloat, we 5 pose and analyse a flowline model by two methods: direct numerical solution and matched asymptotic expansions. The latter leads to a boundary layer formulation that predicts flux through the grounding line as a function of depth to bedrock, channel width, basal drag coefficient, and a calving parameter. By contrast with unbuttressed marine ice sheets, we find that flux can decrease with increasing depth to bedrock at the grounding line, reversing the usual stability criterion for steady grounding line location. We show how this anomalous behaviour relates to the strength of lateral versus basal drag on the grounded portion of 10 the glacier, and to the specifics of the calving law used.


Introduction
In the theory of laterally unconfined marine ice sheet flow, a standard result is that flux through the grounding line is an increasing function of bedrock depth (Weertman, 1974;Thomas and Bentley, 1978;Schoof, 2007a). This leads to the conclusion that grounding lines can have stable steady states only when the ice sheet bed has sufficiently steep down-flow slope (Fowler, 2011;15 Schoof, 2012): a slight advance in grounding line position into deeper water leads to an increase in flux through the grounding line, causing the ice sheet retreat back to its original position. Analogously, a slight retreat leads to a reduction in flux through the grounding line and a re-advance of the ice sheet.
There are a number of mechanisms that can alter the flux-to-bedrock-depth relationship, including the appearance of 'hoop stresses' in an ice shelf fringing the ice sheet (see Pegler and Worster, 2012;Pegler, 2016, though these may require unrealis- 20 tically large ice shelves), the fact that bedrock elevation can actually change due to loading and unloading of the lithosphere (Gomez et al., 2010), thermomechanically mediated changes in basal friction (Robel et al., 2014(Robel et al., , 2016, and lateral drag due to geometrical confinement of the flow into a channel (Dupont and Alley, 2005;Jamieson et al., 2012). The latter is probably the most significant mechanism; when ice flows through a channel, drag can be generated by the side walls of the channel. Drag at We consider a flowline model for a rapidly sliding, channelized outlet glacier with a parameterized representation of lateral drag. The model has the same essential ingredients as those in Dupont and Alley (2005) Jamieson et al. (2012) Nick et al. (2010), Hindmarsh (2012) Pegler et al. (2013) and Pegler (2016). Figure 1 shows the physical domain. Mass accumulates over the glacier and is advected seaward by ice flow. Mass is ultimately lost by flow across the grounding line and eventual 5 calving of icebergs. Let x be the along-flow coordinate and t time, while u(x, t) and h(x, t) are width-averaged ice velocity and thickness, respectively. If b(x) is bed elevation and w(x) the width of the outlet channel, each assumed constant time, then we model force balance and mass conservation as where subscripts x and t denote partial derivatives. Here, ρ i and ρ w are the densities of ice and water, respectively, and g is acceleration due to gravity, while a is surface mass balance. To keep our analysis easier to follow, we ignore any basal melting occurring at the base of the floating portion of the ice; the effect of melting will be addressed in a separate paper. The indicator function θ is given by in other words, θ = 1 if and only if the ice thickness is above flotation and the glacier is grounded.
The parametersB and n are the usual parameters in the Glen's law rheology for ice (Paterson, 1994). We neglect any complications associated with the dependence of ice viscosity on temperature or moisture content, and treatB as well as n as constant. C is a drag coefficient in a power-law basal friction law, with m the corresponding exponent (e.g. Budd et al., 1979;Fowler, 1981); on theoretical grounds, it is often assumed that m = 1/n. Note that other friction laws have also been 20 considered in boundary layer models for unconfined ice sheets (Tsai et al., 2015). The basal drag term only applies where ice is grounded, corresponding to θ = 1. For simplicity, we neglect the possibility that C may depend on additional degrees of freedom such as basal water pressure or temperature.
The second termB w −1/n−1 h|u| 1/n−1 u is a parameterization of lateral drag, withB another constant. We assume that lateral drag is exerted on both, grounded and floating ice. A more sophisticated treatment of lateral drag would require a 25 domain with two horizontal dimensions, and an additional equation representing force balance in the direction perpendicular to the channel axis (e.g. Wearing et al., 2015).
It is worth noting however that (1a) becomes accurate in one of two mutually incompatible limits: (i) a wide channel where lateral drag vanishes (this is the one-dimensional flow case previously studied in e.g. Schoof (2007a)), or (ii) a long and narrow glacier where extensional stress is insignificant and there is no significant flow transverse to the channel axis. By 'extensional 30 stress', we mean the non-cryostatic part of normal stress on a vertical surface placed across the flow, that is, 2B|u x | 1/n−1 u x . In the narrow-channel case, assuming no slip at the channel side walls, (1) predicts the correct width-averaged velocity if we put (see also e.g. Raymond, 1996, for details on flows dominated by lateral shear) B = (n + 2) 1/n 2 (n+1)/nB ; smaller values ofB can be justified if there is actual sliding at the lateral margins of the ice. We use (1a) even when the neither of the two limits above apply. As discussed in Pegler (2016), this is a simplification that works reasonably well and allows 5 at least semi-analytical progress to be made. The simplicity of the model has also led to a large number of authors adopting versions of it. We proceed in that spirit, analysing the model at face value.
We denote the glacier terminus by x c (t); this is the location where ice cover ends. Since x c is a free boundary, two boundary conditions are required. One is needed to close the elliptic problem (1a) and another to determine the evolution of x c . The former is a condition on extensional stress at the ice front (e.g. Schoof, 2007b, appendix B): We take the second condition to be a 'calving law'. The process of calving remains relatively poorly understood, but several calving laws have been developed on theoretical grounds. Our aim is to illustrate how different calving laws can lead to qualitatively, rather than quantitatively, different dynamics in the outlet glacier. We consider two possible calving laws. The first is the 'CD' model from Nick et al. (2010) and the second is a calving law that states that ice breaks off when the glacier reaches its flotation thickness.
The CD model works based on the assumption that water in surface crevasses affects the depth to which those crevasses can penetrate. When they penetrate deeply enough to connect with basal crevasses, calving occurs. When they do not, there is no calving and the ice front simply moves at the velocity of the ice. Algebraic manipulation of the Nick et al. (2010) CD model 5 shows that connections with basal crevasses occur instantly in the model when ice thickness is at (or below) a value h c . In other words, the evolution of the calving front x c satisfies where the dot indicates differentiation with respect to time. Note that the domain lies to the left of x c , soẋ c < u implies that 10 ice is removed by calving. h c itself can be written as a function of the crevasse water depth and of local bedrock depth. This function can be expressed as and ν is The Cryosphere Discuss., doi:10.5194/tc-2017-42, 2017 Manuscript under review for journal The Cryosphere Discussion started: 28 March 2017 c Author(s) 2017. CC-BY 3.0 License.
The function φ is the ratio of calving front thickness to depth to bedrock; the form of φ is illustrated in panel a of Fig. 2.
As an alternative to the CD model, in which the function φ is defined through (1i), we consider a glacier that calves 'at flotation.' This means that calving occurs when h = −(ρ w /ρ i )b at the calving front, but not when h is larger. This is easy to incorporate into the calving framework (1f)-(1h) above: we simply have to replace the definition of φ in (1i) with the simpler This condition is effectively what applies in previous work on marine ice sheet flow without sidewall drag as considered in e.g. Schoof (2007a Once the calving front is afloat, h c no longer depends on bedrock depth: when −d w /b < 1/2, we simply have In other words, for a fixed water depth parameter d w and sufficiently large bedrock depths the Nick et al. (2010) CD model is actually a calving law that simply states that ice breaks off a floating glacier shelf at a critical thickness that is determined purely by the parameter d w (panel b of Fig. 2).

15
To complete the notation for our model, we also define the grounding line position x = x g . For a glacier with a floating shelf, this is the point x = x g (t) at which θ changes discontinuously: For calving at flotation, there is no floating shelf; the grounding line and calving front coincide. The CD calving model goes further, permitting cases where there is no floating ice shelf and the flotation condition h = −(ρ i /ρ w )b is attained nowhere 20 along the glacier, including at the calving front. To keep our terminology as simple as possible, we identify the grounding line in that case with the terminus location, For later convenience, we also define the ice thickness h g at the grounding line and the flotation thickness h f at the grounding line through For glaciers with a floating shelf, we always have h g = h f , and for that reason, existing theories for marine ice sheets generally do not make a distinction between h g and h f (e.g. Schoof, 2007a). The distinction becomes relevant when there is no floating ice shelf, in which case we only have h g ≥ h f . We will use h f frequently below, as it is a simple function of bedrock depth at the grounding line, and therefore determined purely by the geometry of the glacier channel. h g additionally depends on the One of the biggest pitfalls of the calving model in Nick et al. (2010) is that it predicts no calving at all if d w = 0 and surface crevasses are free of water. We contrast the Nick et al. (2010) calving model, which is based on tensile failure, with the shear failure model of Bassis and Walker (2011). The tensile failure model requires d w > 0 and predicts calving for any h below the value given by (1f), instantaneously removing all parts of the glacier shelf that are too thin. By contrast, the shear failure model of Bassis and Walker (2011) predicts that calving will start at a critical calving front thickness and not occur below 5 that thickness, so the inequality in (1g) would need to be reversed. It also predicts that once initiated, the calving front will continue to fracture as it moves into thicker ice inland. This is the basis of the catastrophic calving cliff instability mechanism for marine ice sheet collapse advocated by Pollard et al. (2015), but cannot be captured by an analogue of (1f). It is clear that ice sheets whose calving cliff is larger than the critical thickness for shear failure simply cannot persist: they are guaranteed to disintegrate completely or to stabilize in some shape where the calving front thickness is below the critical thickness for shear 10 failure, and the shear failure model by itself does not provide a timescale for that disintegration. We exclude such shear failure from consideration here and focus purely on the CD calving model. 3 Solution of the model 20

Non-dimensionalisation
In the remainder of this paper, we will consider the problem (1) in dimensionless form. The purpose of doing so is two-fold. We define dimensionless variables as u = Dropping asterisks on the dimensionless variables immediately, we obtain and the boundary conditions at the terminus being and either 10 with φ given by (1i) for the CD calving model, or by (1j) (which states that φ ≡ r −1 ) for calving at flotation.

Direct numerical solution
The system (4) can be solved numerically as posed. In this paper, we focus on solutions of the steady-state version of the problem by a shooting method, which provides a straightforward alternative to a solution by more established time-stepping methods. As our method has not been used previously in this context, we sketch it here for completeness; results are presented 15 at the end of this section and in Fig. 3.
We can write the steady state problem as a four-dimensional, first-order autonomous system of differential equations if, in addition to h, we define the phase space variables q, σ and χ through For technical reasons associated with singular behaviour of the steady state problem near ice divides, we also define a new independent variable η through to obtain a first order system of differential equations from (4): We assume there is an ice divide at x = 0, where u = q = 0. Technically, the ice divide then becomes a fixed point of the The trick is to determine the value of h 0 for which the boundary conditions at the glacier terminus are satisfied 5 by means of a shooting method. Given h 0 , the fixed point has a unique orbit that emerges from it. In other words, if h 0 is known, then the solution to (6) can be computed uniquely. A constraint on h 0 therefore arises from imposing the boundary conditions at the glacier terminus, of the form  as a increases, and (ii) a small ice sheet for which x g increases with a. The larger solution branch also contradicts existing understanding of marine ice sheet dynamics, precisely because an increase in surface mass balance causes the grounding line to retreat into shallower water.
Our aim in what follows is to explain this behaviour using the same boundary layer approach as in Schoof (2007a).
4 Approximation: a small lateral aspect ratio 25

A local force balance version of the model
If we take ε as defined in (2) withB given by (1d), we have In other words, ε is a measure of the lateral aspect ratio [w]/[x]. A narrow channel ensures that ε is small, which is the basis for our approximation scheme. With ε small, we can neglect gradients of the depth-integrated extensional stress in (4a) (that Panel ( is, gradients of 4h|u x | 1/n−1 u x ) everywhere except close to the terminus, and find For the case of m = 1/n (which arises naturally from theories of hard-bed sliding, see e.g. Weertman (1957), Fowler (1981), we get a diffusive model for ice thickness evolution, This is essentially analogous to 'shallow ice' models in ice sheet flow (Fowler and Larson, 1978): we have a local balance of forces, and an ice flux that depends on ice thickness and ice surface slope (see also Kowal et al., 2013). Other choices of m > 0 also imply an ice flux uh that is an increasing function of width w, thickness h and surface slope −(h x + b x ), but that flux cannot be computed in closed form.

The grounding line boundary layer
The model (7) holds everywhere except near the grounding line and in the floating ice shelf. Following Schoof (2007a), we can use the method of matched asymptotic expansions (Holmes, 1995) to capture the behaviour of ice flow in that region. This requires us to rescale the dimensionless model (4) to bring back extensional stress at leading order while maintaining an O(1) ice flux q = uh. The appropriate rescaling turns out to be By contrast with Schoof (2007a), we also have to include a potentially non-zero ice shelf length here, so we put We treat H and U as functions of (X, T ) and X c as a function of T . The rescaling in H implies that ice thickness at the grounding line must be small compared with the interior of the glacier. If there is a floating portion, the glacier must however also reach its flotation thickness at the grounding line. We assume that the glacier is at least near flotation if it has a calving front that remains above flotation. This implies that bed elevation must be small compared with ice thickness in the interior.
Specifically, we rescale and assume that B ∼ O(1); the analogous case of laterally unconfined flow discussed in Schoof (2007a) also requires shallow bed topography. In addition, we assume that thickness b and width w change significantly only over length scales associated with the glacier length as a whole. Over the short length scale associated with the boundary layer coordinate X, we treat B = ε −n 2 /(n+1) 2 b(x g ) and width W = w(x g ) as constant. These additional constraints are again analogous to those made in  Schoof (2007a), and imply that we should treat b and b x as small in the outer problem (7). With B constant at leading order in the boundary layer, we can also define a scaled flotation thickness which we will use throughout the rest of the paper as a proxy for water depth to bedrock.
The result is a boundary layer model at leading order in ε. We do not give the detailed derivation here but merely state its form: The additional boundary condition at the calving front takes the form Note that the boundary layer is in a pseudo-steady state. This is again analogous to Schoof (2007a); the time scale for dynamic adjustment of ice thickness in the boundary layer and of calving front position relative to the grounding line is much shorter than the time scale relevant to the evolution problem (7) (see also Pattyn et al., 2012). We emphasize that there is no assumption here that the glacier as a whole is in steady state.
Despite working at leading order in ε 1, we have retained two terms that contain factors of ε in (9). We now make 15 further assumptions about the physics of the flow near the grounding line. Our fundamental assumption will be that lateral drag −w −1/n−1 h|u| 1/n−1 u plays a leading order role in force balance at the grounding line, but that the floating ice shelf, if it exists, is not so long as to fully buttress the grounding line. By this, we mean that the depth-integrated extensional stress 4h|u x | 1/n−1 u x is comparable in magnitude to h 2 all the way up to the grounding line, as is the case at the terminus x c by dint of the boundary condition (4d), but that sidewall drag cannot be neglected. 20 In that physical regime (termed a 'distinguished limit', in which all physical processes are potentially active), we have to assume that the basal drag coefficient in (9a) and the calving coefficient in (9e) are both of O(1), meaning that the two parameters Γ and Λ defined through are of O(1). We confine our analysis to parameter regimes where this is the case. Note that with m > 0 and n ≥ 1, this implies 25 strictly speaking that γ 1, and basal friction upstream of the boundary layer is formally small in the parameter regime we are considering. Similarly, λ must also be small. Asymptotic matching is the mathematical formalism by which the boundary layer problem and the 'outer' problem (7) for the dynamics of the bulk of the glacier are connected (Holmes, 1995). With the assumptions on Γ in place, this leads to so-called matching conditions between the boundary layer and the outer problem: Here Q is the flux at the boundary of the domain of the outer problem, to be determined through the solution of the boundary layer problem. Physically, the first condition states that the flux near the grounding line in the 'outer' problem is the flux that enters the boundary layer at its upstream end, the second states that near that upstream end of the boundary layer, extensional 5 stress gradients have become insignificant and flux is given by a shallow ice type formula (with U = Q/H, the condition can be re-written as Q ∼ −H|H X | n−1 H X , the appropriate local-force-balance formula in our case), while the third states that velocities in the interior of the boundary layer are large compared with those in the rest of the glacier. Structurally, the boundary layer problem above is very similar to that in Schoof (2007a), with additional physics due to lateral shearing and calving accommodated at the cost of a more complicated formulation.

10
From the perspective of the model (8) for the dynamics of the glacier as a whole, the purpose of the boundary layer is to provide the relevant boundary conditions at x = x g . As (8) is a diffusion model for h, it requires two boundary conditions at any moving boundary. Where a floating portion exists, a condition on h at the grounding line arises straightforwardly from (1k); for a grounded calving front, an equivalent condition is provided by (4e).
As in previous work (Schoof, 2007a), we can show that the second boundary condition takes the form of a flux condition that 15 can be found by solving the boundary layer problem: the problem (9)  We give additional detail on how to compute that relationship between flux, geometry and model parameters in appendix B and in the supplementary material, and the code used to solve the problem is also included in the supplementary material.
Importantly, we are able to show the relationship takes the form The practical use of this form is that it reduces the complexity of the flux formula: for a given set of constants n, m and r, values of H f /Λ, the CD model has a grounded calving front, while the calving front becomes the terminus of a floating ice shelf at larger values of H f /Λ. Note that flux in the CD model is always a decreasing function of H f /Λ for H f /Λ slightly less than the critical value of 2 for changeover from a grounded to a floating terminus. This observation will be key to our interpretation of the physics involved in the anomalous flux-thickness relationship.
Other features of our solutions are also shown in panels a-c of Fig. 4. Each panel isolates one parameter (H f , Λ, W ) on the 30 horizontal axis, but normalizes it as dictated by (10), and plots it against an also normalized flux (again as dictated by (10) There are two aspects of the CD model flux solution that we still need to explain in more detail: (i) why flux decreases with increasing flotation thickness H f in some circumstances and (ii) why flux approaches a constant limit for large H f and so 15 becomes independent of depth to bedrock in the channel, depending instead only on the calving parameter Λ. We turn to these problems next.
Q ≈ 1 4 This formula is essentially a modification of formula (29) in Schoof (2007b), and its derivation is a translation of appendix A of Schoof (2007a) to our modified boundary layer problem. The omission of the extensional stress gradient can also be formalized 5 on the basis that the density difference (1 − r) between ice and water is small, leading to gradients of HΣ being negligible in the balance of forces (see the supplementary material).
For the case of calving at flotation, it is easy to extract an analytical formula for flux as a function of channel width and depth to bedrock from (11). Specifically, we have H g = H f and Σ g = (1 − r)H f /2, so we get 10 which simply generalizes formula (3.51) in Schoof (2007a) with m = 1/n to the case of lateral as well as basal drag. With Γ = 0, we can also immediately recognize a version of formula (10b) with Panels b and c of Fig. 6 shows that (12)

Grounded calving fronts
Here, we are interested in the anomalous relationship between Q and the flotation thickness H f . Recall that H f = −B/r is given by depth to bedrock B, and is therefore prescribed for a given grounding line location. The actual ice thickness H g at the grounding line is equal to H f when the glacier has a floating ice shelf or calves at flotation; for a grounded calving cliff, H f may exceed H g . However, as Fig. 2 shows, H g always increases with H f . Equation (11) further shows that flux Q increases 5 with ice thickness H g and extensional stress Σ g at the grounding line. The anomalous relationship must therefore hinge on Σ g decreasing sufficiently rapidly as flotation thickness H f increases.
Note that the anomalous decrease in flux with increasing flotation thickness is most pronounced around the critical value

calving front thickness as
where the prime on φ denotes an ordinary derivative. For the CD model, (1i) shows that φ(Λ/H f 0 ) = φ(1/2) = r −1 , and In other words, at linear order, a small perturbation in bedrock depth has no effect on ice thickness at the calving front in the CD model. Note that because φ in the CD model is continuously differentiable, this holds regardless of whether H f is positive or negative, that is, regardless of whether the perturbation causes a calving cliff thicker than flotation or a floating ice shelf to form. calving front thickness H g = H c ≈ H f 0 . The stress condition (9d) can also be approximated to first order in H f as As we have assumed that the flotation thickness perturbation H f is negative, we see an The extensional stress perturbation occurs because the calving cliff ice thickness has not changed at first order, but bedrock depth is shallower. The calving cliff now protrudes further above the water line and the depth-averaged normal stress exerted on it by the water is smaller. As a result, the extensional stress in the ice has to increase. This increase in stress is what leads to 10 the increase in flux caused by the decrease in flotation thickness H f . In fact, for small H f < 0, (11) then becomes and Q increases as H f becomes more negative.
This is consistent with the behaviour shown in Fig. 4. For grounded calving, we always find the anomalous relationship between Q and H f , regardless of the basal friction parameter. Panels b and c of Fig. 6 also show that (16) is accurate only for 15 very small H f ; this is presumably a result of the fact that 1 − r = 0.1 is not extremely small, and of the fact that the quadratic term in (14) starts to become large enough to affect results (indeed, panel a of Fig. 2 indicates that a linearization of φ is unlikely to be accurate for grounded calving fronts except very close to H f /Λ = 2.)

Floating calving fronts
We can conversely take the case of H f > 0, which leads to the formation of a floating ice shelf. As the calving front thickness 20 does not change to first order, the extensional stress at the calving front remains equal to Suppose that basal drag is not so large as to render lateral drag insignificant on the grounded portion of the boundary layer.
In that case, even though the ice at the grounding line is slightly thicker than at the calving front (by H f ), the driving stress in the floating ice shelf is small compared with the other forces acting on the shelf. In particular, the surface slope of the 25 ice shelf is small because of the small density difference between ice and water. Most of the reduction in thickness between grounding line and calving front is accounted for by the bottom of the ice shelf sloping upwards (see panel b1 of Fig. 7). The surface slope of the ice shelf (which causes the driving stress) is only (1 − r)/r times the bottom slope. Since the driving stress is weak, the dominant balance of forces on the ice shelf is then between the gradient of depth-integrated extensional stress (HΣ) X = 4(H|U X | 1/n−1 U X ) X and the lateral drag W −(n+1)/n H|U | 1/n−1 U : in other words, (9a) becomes approximately with the driving stress an O(1 − r) correction. It follows that the floating ice shelf acts to reduce extensional stress at the grounding line relative to its value at the calving front; this is the 'buttressing' effect of the ice shelf.
For small H f , we obtain a short ice shelf and can treat H ≈ H f 0 and U ≈ Q/H f 0 in the shelf as constant, so that X c being the length of the floating ice shelf. This effect -the linear reduction in extensional stress in the floating shelf with 5 distance from the calving front -is illustrated in panel b of Fig. 7.
The shelf length is dictated by H f . A bigger flotation thickness requires a longer ice shelf before the calving front thickness H f 0 is reached, potentially leading to more buttressing. The ice shelf thickness gradient is −H X = HU X /U ≈ H 2 f 0 Σ n /(4 n Q), and the shelf length X c is constrained by the fact that the decrease in ice thickness between grounding line and calving front is This allows us to compute X c and hence the stress at the grounding line Substituting this in (11), flux satisfies for small H f > 0, At first glance, it does not seem that (20) is much use -it defines Q implicitly. However, from (19), it is not difficult to show that an increase in H f leads to a decrease in extensional stress Σ g at the grounding line and, therefore, to a decrease in flux 15 Q. The stress decreases because the ice shelf lengthens as H f increases and the total amount of lateral drag on the ice shelf increases.
Again, we have given an ad hoc derivation for (20). We can formalize that derivation as shown in the supplementary material, once more based on the small density difference 1−r. Panel b of Fig. 6 shows that (20) is more qualitatively than quantitatively accurate for the case of no basal friction Γ = 0. Again, this is presumably the result of 1 − r not being extremely small, and of 20 higher order terms in the approximation scheme above becoming important.
However, as panel a of Fig. 6 also shows, the anomalous behaviour disappears entirely for floating ice shelves when the basal friction coefficient Γ becomes large. In such cases, the argument above must become qualitatively incorrect. The change-over from the anomalous behaviour for flow dominated by lateral drag to the 'normal' behaviour obtained with significant basal drag occurs because when the basal friction coefficient is large, ice velocities near the grounding line become small. This has 25 two effects: it (i) reduces the lateral drag term W −(n+1)/n H|U | 1/n−1 U in (9a) and (ii) increases the thickness gradient and, therefore, the driving stress. Specifically, conservation of mass in the floating shelf dictates that H X = −HU X /U = −H|Σ| n−1 Σ/(4 n U ), so that H X becomes large when U is small. The driving stress −(1 − r)HH X can then no longer be ignored in (9a); (18) is no longer applicable, and neither is (20). An increase in flotation thickness can now potentially cause an increase in ice flux, at least when the calving front is afloat. This is shown in panel c of Fig. 7, and is also described in more formal detail in the supplementary material.
It is relatively straightforward to estimate how large Γ needs to be in order for driving stress to appear at leading order 5 in the shelf. Note that thickness, velocity, and stress are continuous across the grounding line. As a result of (21), so is the thickness gradient H X , but not the surface slope. With large Γ on the grounded side, driving stress balances basal shear, −HH X ∼ Γ|U | 1/n−1 U . In order for driving stress to appear at leading order in the shelf, it should be comparable to lateral drag, so −(1 − r)HH X ∼ W −(n+1)/n H|U | 1/n−1 U . It follows that 10 With r = 0.9, this corresponds to ΓW (n+1)/n /Λ ∼ 20, consistent with panel a in Fig. 6, where the green line corresponds to Finally, consider the limiting case of very large basal friction coefficient (meaning, ΓW (n+1)/n /H (1 − r) −1 ) combined with an ice shelf that has limited extent. By an extension of the argument above, this corresponds to a large driving stress and to lateral drag playing an insignificant role in force balance. In this case, we can make our theory agree with previous work for laterally unconfined flow in Schoof (2007a) by simply ignoring lateral drag in (9a), (HΣ) X − (1 − r)HH X = 0.
Integrating and applying the boundary condition (17) to find that extensional stress at the grounding line is simply given as This is the same extensional stress at the grounding line as we would expect in the case of calving at flotation. The flux increases monotonically with floatation thickness when this is substituted into (11) and Γ is assumed large, Since we are assuming that m = 1/n, this is actually nothing more than a scaled version of formula (3.51) in Schoof (2007a).
As panel a of Fig. 6 shows, we do get agreement between the CD calving law results and calving at flotation for large basal friction coefficients, at least while H f /Λ remains close enough to the critical value of 2: the flux curves then agree well with each (as indicated by the arrow labelled 'driving stress dominant'). 20 For larger H f /Λ, this agreement ceases. The shelf gets long enough that, even with large enough basal friction on the grounded portion, lateral drag on the floating shelf cannot be ignored. The next section describes in more detail the mechanics of a very long ice shelf.

The finite flux limit for large flotation thickness
For a fixed value of Λ and large flotation thickness H f , the flux Q appears to approach a finite limit in panel a of Fig. 4. That limit is of the form with C ≈ 8.67 × 10 −3 for n = 3, regardless of the choice of basal drag parameter Γ. The physics behind this is relatively This situation was previously explored by Hindmarsh (2012) and Pegler (2016). These authors find that ice flux through 15 the calving front is determined in a boundary layer around the calving front in which extensional stress is significant. In our notation, the boundary layer takes exactly the same form as (9) for floating ice (θ = 0), but with different matching conditions: where H c is the prescribed calving front thickness, and The analysis of this boundary layer (a formal derivation of which is included in the supplementary material) is much the same 25 as for (9) Consider as the special case of no basal drag on the grounded part of the glacier. We can show how (24) confirms that we expect an anomalous flux-depth-to-bedrock relationship due to buttressing in the ice shelf. Take a grounded calving cliff just at flotation, with thickness H f = H c . The flux is given by (12) with Γ = 0. Compare this with the flux through a long floating 10 ice shelf that terminates in a calving cliff of the same thickness H c . The solution (24) then predicts that the flux through the floating shelf is smaller than through the grounded calving front, even though the two have the same thickness. This is true at least when the density difference 1 − r is small, because the exponent n 2 /(n + 1) on (1 − r) in the formula for flux for the grounded cliff in (13) is smaller than the exponent n for the floating cliff in (24). This provides further evidence for the buttressing action of the ice shelf leading to an anomalous flux-flotation-thickness relationship. In this paper, we have applied the boundary layer analysis of Schoof (2007a) to a model for channelized outlet glacier flow, incorporating a parameterised description of lateral drag (Dupont and Alley, 2005) and a simple calving law due to Nick et al. (2010). The purpose of this work has been to show how calving and lateral drag can potentially combine to produce a very different relationship between ice flux at the grounding line and glacier bed geometry from that for laterally unconfined marine 5 ice sheet flow. For the latter, ice flux is an increasing function of depth to bedrock, while for a channelized outlet glacier, we find that an 'anomalous' relationship in which flux decreases with increasing depth to bedrock is possible (panel a of Fig. 4 above).
Such an anomalous relationship has significant consequences for stable glacier margin positions. Consider the model (8) for the flow of the glacier as a whole. Two boundary conditions apply at the free boundary x = x g . One of these is a thickness condition, while the second is the flux condition (10a), which can be written in the form x is ice flux, and Q g is the flux Q predicted by the boundary layer problem, written in terms of the original dimensionless parameters and the channel geometry at the grounding line.  (2012))

15
If Q g does not necessarily increase with B, then steady grounding lines located on a downward-sloping beds can become unstable. This is illustrated in Fig. 3, where steady grounding line positions on a downward-sloping bed are plotted against the accumulation rate a over the ice sheet, which is assumed to be spatially uniform, as is channel width w. The steady state grounding line position of (8) is defined implicitly by wQ (−b(x g ), w, λ, γ, m, n, r) = awx g .
Treating x g as a function of a and differentiating both sides with respect to a, we have and when the stability condition (26) is satisfied, dx g /da > 0, so that a stable grounding line must advance when a is increased: a stable ice sheet gets larger when it receives more surface snowfall. Figure 3 shows a solution branch for which this is not the case.
Conversely, we may see grounding lines attain stable steady state positions on upward-sloping beds if Q g decreases with depth to bedrock −b. A second mechanism by which this can happen is the dependence of discharge wQ g on width w: a sufficiently narrow bottleneck in the channel could stabilize a grounding line on an upward slope even if Q g did increase with depth −b, because wQ g is an increasing function of w (this argument is due to Jamieson et al., 2012). It is worth noting that simulations of Greenland outlet glaciers using the CD calving law (Nick et al., 2010) have similarly produced steady states located on upward-sloping beds. Our work suggests that this may be due not only to narrowing of the channel but also to the calving law.

5
Our aim has not been to be authoritative in establishing the existence of an anomalous flux-depth relationship: our model contains at least two components that can be improved upon. First, the parameterised description of lateral drag should eventually be dispensed with, replacing our model with one that resolves the cross-channel dimension. The scaling that underlies our boundary layer model should still be applicable in that case, but the actual boundary layer model will consist of a set of coupled partial differential equations (as opposed to ordinary differential equations) and is likely to be much more onerous to 10 solve for a large number of parameter combinations, as we have been able to do here.
Second, the calving law we have employed is relatively poorly constrained by observation and is based on a number of simple assumptions about how cracks form near a calving front. Furthermore, it relies entirely on water depth in surface crevasses as a control parameter that should itself be determined by additional physics governing the drainage of surface melt water.
We have chosen to take the calving model at face value, simply prescribing the crevasse water depth control parameter. This 15 is worth emphasizing as the dependence of calving cliff height on flotation thickness predicted by the calving model turns out to be key to the anomalous flux-depth relationship.
For a floating ice shelf, calving cliff height is simply proportional to crevasse water depth and independent of depth to bedrock. In other words, the CD model can then be thought of as a generic calving model that imposes a fixed thickness at the glacier terminus. Moving the grounding line to a location with greater flotation thickness (or equivalently, depth to bedrock) therefore leads to a longer ice shelf forming before it can reach the prescribed calving cliff height. If the mechanical effect of the ice shelf is primarily to provide lateral drag, then a longer shelf leads to a greater reduction in extensional stress between calving front and grounding line, and therefore to lower ice flux despite a greater depth to bedrock at the grounding line. Whether this occurs or not is a function of basal drag on the grounded part of the glacier: if basal drag upstream of the grounding line is moderate compared with lateral drag, then the surface slope and driving stress of the floating shelf will be small, so the effect of the shelf is mostly to generate lateral drag. By contrast, if basal drag is large upstream of the grounding line, then the floating shelf will be relatively steeply sloped and lateral drag will play a lesser role in force balance there, leading to the possibility that the floating shelf does not cause a reduction in extensional stress and hence flux through the grounding line. The changeover between the two regimes happens when, in the notation of section 3.1, the basal drag coefficient is approximately (see equation For a grounded calving front, the Nick et al. (2010) calving law predicts that calving cliff height decreases relatively slowly when the calving front is moved to a location with shallower depth to bedrock. In turn, this leads to more of the calving cliff 20 26 The Cryosphere Discuss., doi:10.5194/tc-2017-42, 2017 Manuscript under review for journal The Cryosphere Discussion started: 28 March 2017 c Author(s) 2017. CC-BY 3.0 License.
being exposed above the water line, and consequently to larger extensional stresses acting on the calving front, and these larger extensional stresses cause ice flux to increase as depth to bedrock is decreased.
This contrasts with an alternative 'calving at flotation' calving law, in which calving front height is always proportional to depth to bedrock and no floating shelf forms. In that case, extensional stress at the grounding line increases with depth to bedrock, and so does ice flux. 5 We close by noting that our approach can potentially be used to study the effect of other calving laws relatively simply in future, by replacing the function that specifies ice thickness at the calving front. Since an anomalous flux-depth-to-bedrock relationship may be possible at least in principle has significant consequences for stable outlet glacier configurations, and it may be worth testing this before embarking on simulations of actual glaciers using different calving laws.
7 Code availability 10 The MATLAB code used in the computations reported is included in the supplementary material.
The sign of this last eigenvalue is negative when is small enough. In that case, the fixed point has a two-dimensional centre manifold, a stable and an unstable manifold (Wiggins, 2003). The centre manifold has no dynamics (it consists of other fixed points, corresponding to different values of h 0 , or orbits that do not satisfy lim η→−∞ χ = 0), so that the fixed point must be 15 approached in the limit η → −∞ along the unique unstable manifold.
The parameters in this re-scaled version of the model are r, n, m, H f , C and L, while U, H and X are dependent and independent variables, respectively. It then can be shown that the transformed boundary layer problem has a solution if and only if the parameters r, n, m, H f , C and L satisfy some functional relationship with each other. Using this fact, it is easy to show that the simplified flux laws (10) must hold.
and switching to an independent variable ζ defined through This transforms the boundary layer problem into a non-singular dynamical system in which the matching conditions (9g) correspond to a fixed point (Ψ, ξ, Q) = (1, 0, 1) being attained as ζ → −∞, and there is a unique orbit along which this can happen; to prove the uniqueness of that orbit, an additional transformation to ν = ξ (n 2 (m+1)−n)/(2(n+1) 2 ) may be required, 10 but the basic argument remains the same as in Schoof (2011). The boundary conditions (9d)-(9e) then provide two further constraints: one on the location of the point along the orbit that corresponds to the calving front, and another to relate the parameters of the model, meaning a functional relationship between r, n, m, H f , C and L. In practice, this must be solved for numerically by integrating the transformed dynamical system. More complete details on the solution of the boundary layer problem can be found in the supplementary text and the numerical code provided as supplementary material for this paper.