A theory of glacier dynamics and instabilities Part 2: Flatbed ice streams

Abstract

In Part 1, we have considered the dynamics of topographically confined glaciers, which may undergo surge cycles when the bed becomes temperate. In this Part 2, we consider the ice discharge over a flatbed, which would self-organize into alternating stream/ridge pairs of wet/frozen beds. The meltwater drainage, no longer curbed by the bed trough, would counter the conductive cooling to render a minimum bed strength at some intermediate width, toward which the stream would evolve over centennial timescale. At this stationary state, the stream width is roughly twice the geometric mean of the stream height and length, which is commensurate with its observed width. Over a flatbed, streams invariably interact, and we deduce that the neighboring ones would exhibit compensating cycles of maximum velocity and stagnation over the centennial timescale. This deduction is consistent with observed time variation of Ross ice streams B and C (ISB/C), which is thus a manifestation of the natural cycle. Moreover, the model uncovers an overlooked mechanism of the ISC stagnation: as ISB widens following its reactivation, it narrows ISC to augment the loss of the meltwater, leading to its stagnation. This stagnation is preceded by ice thickening hence opposite to the thinning-induced surge termination.

1. Introduction

In Part 1 (Ou, 2021), we have considered the dynamics of a glacier confined in a topographic trough, which may exhibit cyclic surges when the bed becomes temperate. We establish that a prognostic model of the phenomenon must entail global momentum balance and subglacial hydrology in its closure. The requirement is met by the undrained plastic bed (UPB) formulated by Tulaczyk and others (2000) and numerical calculations employing UPB have produced realistic surge cycles (Bougamont and others, 2011; Robel and others, 2013). Since these models have fixed glacier margins, they cannot address the free-boundary problem of the ice streams, as they have intended. Bougamont and others (2011), for example, liken the simulated surge cycle to the time variability of the Ross ice streams, including the recent stagnation of Ice Stream C (ISC, also named after Kamb). Although observed streams do exhibit time variation, it differs qualitatively from surge cycles characterized by short surge followed by long period of lull (Meier and Post, 1969; Shabtaie and others, 1988). Robel and others (2013) note in addition that the observed stream width lies near the stability boundary between cyclic and steady regimes, which however serves no prognostic purpose unless the stream width can be internally constrained.

In contrast to topography-confined glaciers, there is no apparent spatial linkage of Ross ice streams to the bed topography (Retzlaff and others, 1993; Bennett, 2003), so as a more relevant approach, numerical calculations have been carried out over a flatbed, which show that the ice discharge would self-organize into alternating fast/slow bands, but they remain inadequate in addressing the ice stream problem, as noted next. Hindmarsh (2009), for example, shows that thermal instability due to the strain heating would speed the local flow and draw in the adjacent ice, and the accompanying cold advection would curb the thermal instability to result in alternating warm/cold bands. Since his fast flow remains a viscous creep, it only slightly perturbs the global momentum balance, whose inclusion has nonetheless lessened the resolution dependence of the solution. His model, however, is deficient in simulating Ross ice streams characterized by fast sliding motion (Joughin and others, 2002), which would strongly perturb the global momentum balance, so the inclusion of the latter in a model is critical. Indeed, without such balance, the sliding streams are only one to two gridpoints wide (Payne and Dongelmans, 1997; Hulton and Mineter, 2000), a narrowness that is thus an artifact of the model resolution, and the singular stress at the stream edge obviously cannot be accommodated by the global momentum balance.

Remedying the above shortfall, Brinkerhoff and Johnson (2015) have included the global momentum balance in their model, which shows that the streams attain finite widths, thus removing the singular stress at their margins. In addition, the streams exhibit surge cycle, much like that observed for a topography-confined glacier when it is sufficiently wide (Part 1). The modeled variability, thus, is due to the thermal switch considered in Part 1 (Ou, 2021), and the stream width is essentially that divides the cyclic-surge and steady-sliding regimes. In other words, the additional degree of freedom associated with the unknown stream width can be countered by the Hopf bifurcation when the unstable stationary state vaults into a stable limit cycle (Robel and others, 2013). Being a thermal-induced surge cycle, however, it cannot explain the observed variability, as alluded to earlier; then their neglect of the meltwater dispersal cannot be justified over a flatbed. Since the stream is predicated on a wet bed, one expects the efficiency of the meltwater drainage to play a key role in regulating its width.

Kyrke-Smith and others (2015) have included the meltwater dispersal in their model, and to counter the unknown stream width, they posit a constraint on the viscous closure of the water film by the strain heating. The constraint, however, involves an arbitrary parameter (their Λ_∞) limiting the weakness of the bed, and a stationary state is attainable only when this constraint is less severe. Then, when the streams do settle into this stationary state, there is no longer time variability, so the model cannot address the observed temporal variation, including the stagnation of ISC.

With the above, we discern the needed expansion of the theoretical framework presented in Part 1 (Ou, 2021) to address the ice stream problem, as expounded next. First of all, without the topographic trough confining the glacial mass flux, we must allow inflowing ice flux from the adjacent ridge (termed ‘entrainment’), which indeed is comparable to the upper-glacier catchment in observations (Shabtaie and others, 1988). The accompanying cold advection is the facilitator of the self-organized stream/ridge pair, which thus must be considered a single physical entity. Second, without the topographic constriction, an undrained substrate is no longer tenable, and the efficiency of the meltwater drainage necessarily regulates the stream width defined by a wet bed. Finally, the closure must allow interaction of the streams over a flatbed, which should not be stymied by extraneous constraints. These expansions notwithstanding, we are seeking a minimal model to address the genesis and interaction of the flatbed ice streams. The striving for simplicity uncluttered by unwarranted sophistication should aid the understanding of the governing physics.

Since we are concerned with the physical closure of the whole stream/ridge system, it differs in scope from recent studies addressing the migrating stream boundary due to local concentration of the shear stress hence strain heating (Jacobson and Raymond, 1998; Schoof, 2012). As such concentration is caused by the weak stream bed (Echelmeyer and others, 1994; Harrison and others, 1998), which in turn is a consequence of the substream hydrological balance, the latter must be considered in assessing the stream behavior. Inferring the boundary migration from the englacial heat balance alone by fixed or diagnosed internal variables offers no prognosis of the boundary motion or stability (Jacobson and Raymond, 1998; Haseloff and others, 2018). In fact, the margin need not move if the heightened shear heating is countered by the advective cooling, the latter sometimes neglected in the aforementioned models. As we shall see later (Section 2.2), the relative warmth near the stream margin has only minor effect on the stream-wide hydrological balance, so it is neglected in our minimal model.

This Part 2 is organized as follows: In Section 2, we consider an isolated stream/ridge pair and derive its stationary state. In Section 3, we consider the interaction between neighboring pairs and examine their equilibrium configuration, which will be compared with observations. We provide additional discussion in Section 4 and summarize the paper in Section 5.

2. Isolated stream/ridge

We first consider an isolated stream/ridge pair shown in Figure 1, which depicts the middle section between the ice divide and terminus where the driving stress is broadly peaked to endow the primary stream behavior (Bindschadler and others, 2001). The state variables pertain to this middle section, and the longitudinal dimension enters only through the stream half-length that defines the upper-glacier catchment distance and the longitudinal surface slope. All symbols and their prescribed values (referred to as ‘standard’ for convenience) are listed in Appendix. Both the stream and bed conditions are shown with the external variables boxed while internal ones unboxed, and they are denoted as subscripts 0 and 1, respectively, for the ridge and stream.

Fig. 1. Isolated stream/ridge pair, with boxed-in variables prescribed and others, prognostic. l is the longitudinal distance scale, h ₀ and h ₁ are the ridge and stream height, w ₀ and w ₁ are the catchment width and the stream half-width, u ₀, v _i and u ₁ are the longitudinal creep, the inward creep and the mean stream velocity, τ ₀ and τ _b are the driving and the basal stress, their difference is half the effective pressure N, which is the difference between the overburden p _i and water pressures p _w, $\dot{a}$ and $\dot{g}$ are the accumulation and the geothermal flux, respectively. The stream boundary is defined by the basal stress equaling the driving stress, and the water pressure is reduced over a narrow transition under the ridge to zero, which divides the frozen and wet bed.

The external parameters are the constant accumulation $\dot{a}$ and the geothermal flux $\dot{g}$, the latter specifies the ice height h ₀ when the bed becomes temperate to initiate the streams (Part 1). This height combined with the stream half-length l then specifies the driving stress τ ₀ and the longitudinal creep u ₀ over the ridge outside its interacting zone with the stream. The internal variables are the height h ₁, the half-width w ₁, the mean basal stress τ _b and sliding velocity u ₁ of the stream, as well as the inflow v _i supplied over the catchment width w ₀ of the ridge.

As drawn on the till surface, the stream boundary is defined by a basal stress equaling the driving stress, which thus separates the free slip from the non-slip beds. The substream basal stress is approximately half the effective pressure N (Tulaczyk and others, 2000), a difference between the overburden (p _i) and the water pressure (p _w). Because of the narrowness of the shear zone, we have equated the mean and the centerline basal stresses, as denoted by τ _b. Since the finite water pressure and the outward meltwater flux at the stream boundary can only be depleted by freezing over a finite distance, the stream boundary contains no discontinuity in the temperature (at the pressure-melting point), the basal stress (at the driving stress), the longitudinal velocity (~0) or bed condition (wet).

The external parameters are justifiable as the ice-sheet properties unaffected by the streams. The internal stream variables, other than its width, are the same as those of Part 1 (Ou, 2021), which represent minimal descriptors of the stream. The addition of the inflow velocity and the catchment width as internal variables are necessitated by the mass coupling of the stream/ridge, which however can be closed by the viscous flow law and the mass balance over the ridge. Therefore, compared with Part 1 (Ou, 2021), Part 2 indeed has just one additional degree of freedom associated with the stream width, whose determination is the original aim of this study.

Figure 1 illustrates the essence of the self-organization: as the ice sheet thickens, the increasing sequestration of the geothermal heat would warm the bed to the pressure-melting point to trigger the sliding motion (Part 1); the ensuing surface lowering would induce an inward flow from the adjacent ridge and the accompanying cold advection would freeze the ridge bed, resulting in self-organized stream/ridge pairs of wet/frozen beds. This contrast in the bed condition is well-observed over the Siple coast (Bentley and others, 1998; Kamb, 2001), which thus can be the outcome of the self-organization rather than a precondition that differentiates the disparate flow, as suggested by some authors (Blankenship and others, 2001).

2.1 Mass coupling

The ice entrainment from the ridge into the stream would couple the mass balances of the two, which are discussed successively below. For simplicity, we assume the surface displacement to be sufficiently slow that steady-state mass balance holds. Over the ridge, the longitudinal driving stress is

(1)$$\tau _0 = \rho _{\rm i}gh_0^2 /l, \;$$

and, applying a linear flow law (see Part 1: Ou, 2021), the vertical-averaged longitudinal creep is (van der Veen, 2013, Eqn (5.30))

(2)$$u_0 = h_0\tau _0/( {3\nu } ) , \;$$

where ν is the ice viscosity (Doake and Wolff, 1985). The lowering of the stream surface would induce an inward stress τ _i of

(3)$$\tau _{\rm i} = \rho _{\rm i}gh_0^2 \eta /w_0, \;$$

where w ₀ is the catchment width over the ridge and

(4)$$\eta \equiv ( h_0-h_1) /h_0\ll 1$$

is the fractional surface depression of the stream, assumed small at the outset for simplicity. Since this depression is explicitly calculated in the model, this approximation will be checked a posteriori for consistency. The inward stress would drive a vertical-averaged inward flow:

(5)$$v_{\rm i} = h_0\tau _{\rm i}/( {3\nu } ) , \;$$

accompanied by a mass flux

(6)$$q_{\rm_ i} = h_0v_{\rm i}.$$

Because of the assumed steady-state mass balance, this influx is supplied by the accumulation over the catchment width of the ridge minus the longitudinal flux divergence, or

(7)$$q_{\rm _i} = w_0( {\dot{a}-h_0u_0/l} ) \sim w_0\dot{a}.$$

We have neglected the longitudinal flux divergence as it is more than an order smaller than the accumulation for the standard case (Appendix). Combining Eqns (3)–(7), we link the two internal variables η and w ₀:

(8)$$\eta = \displaystyle{{\dot{a}} \over {u_0}}\displaystyle{{w_0^2 } \over {h_0l}}.$$

Physically, a greater depression implies a greater inflow, which needs to be supplied by the accumulation over a wider catchment zone. Based on the functional form of this equation, we define the width scale

(9)$$[ w ] \equiv ( {h_0l} ) ^{1/2}, \;$$

so this equation becomes

(10)$$\eta = \omega w_0 ^{ { \prime ^{\vskip 2pt 2}}} , \;$$

where primed variables are non-dimensionalized and

(11)$$\omega = \dot{a}/u_0$$

is a dimensionless parameter measuring the strength of the accumulation hence referred to as the accumulation parameter.

We next consider the mass balance of the stream, which states

(12)$$\dot{a}w_1 + \dot{a}w_0 = w_1q_1/l, \;$$

or the accumulation and the influx from the ridge (the second term, see Eqn (7)) are balanced by the longitudinal divergence of the stream flux. This equation can also be interpreted as the mass balance of the combined stream/ridge pair: that is, owing to the negligible longitudinal flux of the ridge, the total accumulation over the pair is discharged through the fast stream. As discussed in Part 1 (Ou, 2021), applying the global momentum balance and the linear flow law to the lateral strain, the stream velocity, averaged across the stream, is (Raymond, 1996)

(13)$$u_1 = \displaystyle{{\tau _0-\tau _{\rm b}} \over {3\nu h_0}}w_1^2 , \;$$

where τ _b is the (unknown) mean basal stress. Note that because of the small surface depression (4), the same driving stress τ ₀ applies over the stream/ridge pair. Defining additional scales

(14)$$[ \tau ] = \tau _0, \;$$

(15)$$[ u ] = u_0\displaystyle{{{[ w ] }^2} \over {h_0^2 }}, \;$$

and

(16)$$[ q ] = h_0[ u ] , \;$$

the velocity and stream flux are non-dimensionalized to

(17)$$u_1 ^ {\hskip1pt{\prime}} = q_1^ {\hskip1pt{\prime}} = ( {1-\tau_{\rm b} ^{\hskip2pt {\prime}} } ) w_1 ^{\hskip2pt {\prime} ^{\hskip1pt2} } , \;$$

and Eqn (12) becomes

(18)$$w_0^ {\hskip1pt{\prime}} = w_1^ {\hskip1pt{\prime}} ( {q_1^ {\hskip1pt{\prime}} -\omega } ) /\omega , \;$$

which links the catchment width to the stream width.

Eliminating the catchment width $w_0^ {\hskip1pt{\prime}}$ from Eqns (10) and (18), we arrive at

(19)$$\eta = w_1^ {\hskip1pt{\prime}\hskip1pt ^2}[ {( {1-\tau_{\rm b}^ {\hskip1pt{\prime}} } ) w{_1 ^ {\hskip1pt{\prime}} }^2-\omega } ] ^2/\omega , \;$$

which contains only stream variables. As alluded to earlier, the additional degrees of freedom introduced by the entrainment have been closed by the viscous flow law and the mass balance over the ridge.

2.2. Meltwater balance

For a temperate bed, the basal heat balance is the same as the meltwater balance, and its difference from Part 1 is that the lateral dispersal of the meltwater, no longer curbed by the topographic trough, must be included. Since the stream is defined by an effective pressure lower than half the driving stress (Tulaczyk and others, 2000, see Fig. 1), one expects the efficiency of the meltwater dispersal to play a key role in constraining the stream width. To maintain the high water-pressure necessary for the sliding motion, the drainage system must be of the distributed type, which may take the form of linked cavities, a thin water film or their hybrid (Fountain and Walder, 1998), but for our purpose, the distinction is immaterial since it can be absorbed into the effective hydraulic conductivity. For simplicity and indeed with stronger observational evidence (Engelhardt and Kamb, 1997), we take the water film of uniform depth d as the proxy medium for the water transport.

Since the hydrological timescale is short compared with the timescale of the ice variability (Fricker and others, 2007; see Section 2.3), we assume a steady-state heat balance of the form (Kyrke-Smith and others, 2014)

(20)$$u_1\tau _{\rm b} = \dot{g}\left({\displaystyle{{h_0} \over {h_1}}-1} \right) + \displaystyle{{\rho _{\rm i}Ld^3} \over {12\mu }}\displaystyle{{2( {\tau_0-\tau_{\rm b}} ) } \over {w_1^2 }}, \;$$

which states that the frictional heating (the first term) is balanced by the conductive cooling in excess of the geothermal flux (the second term) and the drainage cooling associated with the outward meltwater flux (the last term). For the conductive cooling, we have assumed an englacial conductive equilibrium (Part 1) and noted that the marginal warmth only slightly perturbs the stream-wide conductive flux because of its narrowness (Harrison and others, 1998) hence neglected. For the meltwater drainage, we have applied a scaled Darcy's law neglecting the surface slope because the stream resembles a plug flow (Joughin and others, 2002, Fig. 6) and noted that the effective pressure has end values as indicated in Figure 1. It is noted that the outward meltwater flux is wholly determined by the stream properties, and it can be accommodated by the ridge because of the slight set-back of the frozen bed (Fig. 1). Since both freezing and melting are allowed in the hydrological balance, there is no need to distinguish the expanding versus the shrinking stream (Haseloff and others (2018).

When non-dimensionalized, the heat balance (20) becomes

(21)$$0 = \underbrace{{\alpha u_1^ {\hskip1pt{\prime}} \tau _{\rm b}^ {\hskip1pt{\prime}} }}_{{\dot{\,f}}}-\underbrace{\eta }_{{\dot{c}}}-\underbrace{{\beta ( {1-\tau_{\rm b}^ {\hskip1pt{\prime}} } ) /w{_1^ {\hskip1pt{\prime}} }^2}}_{{\dot{d}}}, \;$$

where

(22)$$\alpha \equiv [ u ] [ \tau ] /\dot{g}$$

is the ‘heating’ parameter measuring the strength of the frictional heating and

(23)$$\beta = \displaystyle{{\rho _{\rm i}Ld^3} \over {12\dot{g}\mu }}\displaystyle{{2\;[ \tau ] } \over {{[ w ] }^2}}$$

is the ‘drainage’ parameter measuring the efficiency of the meltwater dispersal. Since this drainage parameter differs from the hydraulic conductivity only by a numerical factor, they are referred to interchangeably. We have used the symbols $\dot{f}$, $\dot{c}$ and $\dot{d}$ for the ‘frictional’ heating, ‘conductive’ cooling and ‘drainage’ cooling, respectively. Substituting from Eqns (17) and (19), we combine the mass and heat balances into a single equation governing the two prognostic variables: the half-width $w_1^ {\hskip1pt{\prime}}$ and the basal stress $\tau _{\rm b}^ {\hskip1pt{\prime}}$ of the stream:

(24)$$ 0 = \underbrace{{\alpha w{_1 {^{\hskip-1pt\prime}} }^2( {1-\tau_{\rm b}{^{\hskip-1pt\prime}} } ) \tau _{\rm b}{^{\hskip-1pt\prime}} }}_{{\dot{\,f}}}-\underbrace{{w{_1 {^{\hskip-1pt\prime}} }^2{[ {( {1-\tau_{\rm b} {^{\hskip-1pt\prime}} } ) w{_1 {^{\hskip-1pt\prime}} }^2 - \omega } ] }^2/\omega }}_{{\dot{c}}} - \underbrace{{\beta ( {1-\tau_{\rm b} {^{\hskip-1pt\prime}} } ) /w{_1 {^{\hskip-1pt\prime}} }^2}}_{{\dot{d}}}, \;$$

which contains three dimensionless parameters: ω, α and β.

To tidy up Eqn (24), we use the shorthands

(25)$$\tau \equiv \tau _{\rm b}^ {\hskip1pt{\prime}} , \;$$

(26)$$\zeta \equiv w_1 ^ {\hskip1pt{\prime} \hskip1pt^2} , \;$$

to yield

(27)$$0 = \underbrace{{\alpha \zeta ( {1-\tau } ) \tau }}_{{\dot{\,f}}}-\underbrace{{\zeta {[ {( {1-\tau } ) \zeta -\omega } ] }^2/\omega }}_{{\dot{c}}}-\underbrace{{\beta ( {1-\tau } ) /\zeta }}_{{\dot{d}}}.$$

Since this equation is quadratic in τ, it can be easily solved for τ(ζ), which is the solid curve drawn in Figure 2 for the standard case of (ω, α, β = 0.21,  2.3,  0.24). This curve bisects the net cooling and heating regions, as indicated by the encircled signs, and it has a minimum at some intermediate width marked by box C. The dashed line corresponds to zero surface depression η, which intersects the solid curve at box B. Since the surface depression (10) cannot be negative, it is set to zero to the left of the dashed line wherein the heat balance (27) with $\dot{c} = 0$ would extend the solid curve from boxes B to A. As we shall see later, box A marks both the onset and termination of the stream, whose half-width is then given by

(28)$$w_2 ^ {\hskip1pt{\prime}} = ( {\beta /\alpha } ) ^{1/4}.$$

Fig. 2. Steady-state heat balance (solid curve) on (ζ,τ) phase space for the standard case of (ω, α, β = 0.21, 2.3, 0.24). For a wide stream (of greater ζ), the balance is between the frictional heating $( \dot{f})$ and conductive cooling $( \dot{c})$ whereas for a narrow stream between the frictional heating $( \dot{f})$ and the drainage cooling $( \dot{d})$. To the left of the dashed line, the surface depression (η) is set to zero. Solid arrow indicates the stream evolution from its onset at box A to box B when the conductive cooling commences and to box C when the basal stress is at its minimum and the stream attains stationarity. The dash-dotted line is the entropy production (σ) whose maximum roughly aligns with the minimum basal stress.

It is observed that the solid curve representing Eqn (27) is U-shaped with a minimum in the basal stress, which can be understood by noting that: (a) because the three terms have different power dependence on the stream width, the drainage cooling would dominate at small width, the conductive cooling at large width, and the frictional heating possibly at intermediate width; and (b) because the frictional heating contains an extra appearance of the basal stress, it is overtaken by cooling when the bed is weak. The minimum basal stress thus occurs at some intermediate width when the drainage cooling and conductive cooling are of comparable strength. We shall argue next that following its onset, the stream would evolve toward this minimum when it attains stationarity.

2.3 Stationary state

A stream is initiated from a creep state when the bed first becomes temperate and the meltwater production has weakened the bed to below the driving stress, that is, along the top boundary in Figure 2 with a small width. A narrow stream is subjected to effective drainage of the meltwater, which would widen the stream toward box A when it attains a steady-state meltwater balance. The process is rapid since it is governed by the short hydrological timescale. For this reason, we shall identify box A as the onset of the stream, as alluded to earlier. Beyond box A, we discern possible positive feedback between the stream widening and the bed softening: a slight widening hence faster sliding motion would increase the frictional heating to soften the bed, which would augment the ice flux hence the entrainment to further widen the stream. This positive feedback would propel the stream toward box C when the bed is at its softest to thwart the positive feedback and the stream attains stationarity. Physically, this is when the increasing conductive cooling eventually overcomes the frictional heating to stabilize the bed strength.

It should be noted that although the solid curve represents a steady-state balance, it is the small deviation from this balance that propels the stream evolution. As this evolution involves the slower mass adjustment than the hydrological processes, its timescale should be that of the former, which stems ultimately from the sliding-induced thinning. It is thus of the form

(29)$$[ t ] = \displaystyle{{\eta h_0} \over {[ q ] /l}}$$

(30)$$ = \eta l/[ u ] , \;$$

or it is simply the longitudinal transit time multiplied by the fractional surface depression. For the standard case (Appendix), this transit time is ~2 ka and with η ≈ 0.1 based on the later solution, we have an evolution timescale that is of order 200 years hence centennial, which is consistent with that characterizes the ice stream variation over the Siple coast (Catania and others, 2012). Having identified this timescale, we can infer the temporal behavior of the stream/ridge pair along the solution continuum without explicit time integration. This is the same approach as that taken in Part 1 (Ou, 2021) wherein its utility in depicting the surge cycle has been demonstrated.

In addition to the positive feedback discerned above, we offer next a complementary rationale based on the non-equilibrium thermodynamics (NT). We note that because of the high Reynolds number, the substream water film is likely turbulent (Kamb, 1987), so the effective hydraulic conductivity would be subjected to microscopic fluctuations, whose fractional amplitude can be taken to be the ratio of the turbulent to the mean water velocity. Although Ou (2018) shows that such fluctuations would propel the system toward maximized entropy production (MEP) in a different context, the MEP has been applied over a wide range of phenomena as a veritable generalization of the second law to the NT system (Kleidon and Lorenz, 2005). In our application, the entropy production σ is the product of the thermodynamic flux $( \dot{d})$ and the thermodynamic force (1 − τ) that drives this flux, or

(31)$$\sigma = \dot{d}( {1-\tau } ) $$

(32)$$ = \beta ( {1-\tau } ) ^2/\zeta .$$

We plot the entropy production along the solid curve in Figure 2 (the dash-dotted line) from which we observe that its maximum roughly aligns with the minimum basal stress. Because of the appearance of the stream width in Eqn (32), there is a slight misalignment of the two, which is of little practical consequence for our purpose. The application of the NT renders a different expression of the evolution timescale: for a fractional fluctuation amplitude ɛ, this timescale should be ɛ ⁻² of the hydrological timescale, based on Ou (2018). As a plausible example, we take the turbulent velocity of the meltwater to be an order smaller than its mean velocity, and a hydrological timescale of order 1 year (Fricker and others, 2007), then the evolution timescale could be centennial, similar to that governing the ice mass balance.

Based on above reasonings, we postulate that the stream would evolve on centennial timescale toward the minimum basal stress when it attains stationarity. We shall next examine this stationary state and compare it with observations.

2.4 Application

With the stationary state being of minimum basal stress, we differentiate Eqn (27) against ζ and set dτ/dζ to zero to derive

(33)$$0 = \alpha \zeta ( {1-\tau } ) \tau -\zeta [ {( {1-\tau } ) \zeta -\omega } ] [ {3( {1-\tau } ) \zeta -\omega } ] /\omega + \beta ( {1-\tau } ) /\zeta $$

which can be combined with Eqn (27) to yield

(34)$$\tau = 1-\omega ( {1 + \beta /\zeta^3} ) /\zeta .$$

The equation, being non-transcendental, allows the following procedure for calculating the stationary state. For a guessed ζ, we calculate τ from this equation, both are then substituted into Eqn (27) to check the equality, and the initial ζ is then adjusted until this equality is satisfied. For the standard case shown in Figure 2, the calculated variables have the values of $w_1^ {\hskip1pt{\prime}} = 0.80, \;\,\,\tau _{\rm b} ^ {\hskip1pt{\prime}} = 0.37$ and hence $\eta = 0.11, \;\,\,w_0 ^ {\hskip1pt{\prime}} = 0.73$ and $u_1^ {\hskip1pt{\prime}} = 0.40$. Since the surface depression is much smaller than unity, it supports the approximation (4), as alluded to earlier. Although, as we shall see later, all these values undergo cyclic changes when streams interact with one another, to the degree that the stationary state represents an extremum of the cycle, it has the observational counterpart hence can be compared. Dimensionally, for the standard case, the driving stress thus is 0.23 b (shorthand for bar), the basal stress is 0.09 b, the stream width is 32 km, the catchment width is 15 km, the stream velocity is 77 m a⁻¹ and the stream surface is 110 m lower than the ridge, all seems quite reasonable except for the stream velocity seems low (Joughin and others, 2002).

We next examine the dependence of the stationary state on the dimensionless parameters, as shown in Figure 3 on the (ω, β) space. The axes are in the binary logarithm spanning half to double the standard values (solid squares); the solid lines are for the standard α and the dashed lines for the doubled α. To understand the dependence, we show in Figure 4 how the stationary state defined by the minimum basal stress (box C) would be displaced (solid squares) when the dimensionless parameters are doubled from their standard values (the altered basal heat balance in the dashed lines). A greater ω, say, due to greater accumulation, would raise the stream surface to reduce the conductive cooling, which, as can be inferred from the double-headed arrow on the right, would lower box C and push it to the right. The stationary state, thus, is characterized by a smaller depression, wider stream and weaker bed whereas at the same time Eqn (10) implies a narrower catchment; the wider stream and weaker bed reinforce each other to augment the stream velocity (17); all are as observed in Figure 3.

Fig. 3. State variables plotted against the accumulation (ω) and drainage (β) parameters: (a) stream half-width; (b) basal stress; (c) surface depression; (d) catchment width and (e) stream velocity. Dashed contours are for the doubled heating parameter (α).

Fig. 4. As shown in Figure 2 but showing the modification of the heat-balance curve (solid to dashed lines) and displacement of the stationary state (box C to solid squares) by the doubling of the dimensionless parameters (single-headed arrows).

A greater β implies stronger drainage cooling, which, as can be inferred from the double-headed arrows on the left, would raise the left limb of the solid curve to push box C upward and to the right. The stationary state, thus, is characterized by deeper depression, wider stream and stronger bed, and from Eqn (10) a wider catchment; the stronger bed and wider stream however counter each other to render little net effect on the stream velocity (17); all are as observed in Figure 3.

A greater α implies stronger frictional heating, which would lower box C to weaken the bed, but with little effect on the stream width (only slightly narrower). The weakening of the bed, however, is observed from Eqn (19) to increase the depression and with that the catchment width (10); the weakening bed also speeds up the stream (17); all are again as observed in Figure 3.

Among the external parameters, the most uncertain are perhaps the ice viscosity ν and the water-film depth d. From scale definitions, the dimensionless parameters have the following power dependence on these two parameters: ω ~ ν, α ~ ν ⁻¹, β ~ d ³, so the ice viscosity has no effect on the drainage parameter whereas the water-film depth only affects the drainage parameter. Quantitatively, the doubling and halving of the ice viscosity falls within the parameter domain of Figure 3, which, however, can only accommodate 30% change in the water-film depth. Although one may appraise the sensitivity of the glacial properties from Figure 3, one should be mindful of the limit of such inferences given the crudeness of the model. Perhaps more significantly, Figure 3 suggests relative robustness of the standard solution over considerable ranges of the external conditions; the stream properties, thus, can be summarized crudely as follows: its width is about twice the geometric mean of its height and length (9), its basal stress is less than half of the driving stress, its surface depression is an order smaller than its height and it moves about two orders faster than the creep; all these are consistent with the observed Ross ice streams (Shabtaie and others, 1988; Joughin and others, 2002).

Because of the mass and thermal coupling, a stream/ridge pair must be considered as a single physical entity. Although the pair could attain stationarity when isolated from other pairs, we shall next see that such isolation is untenable and they invariably interact (Catania and others, 2012).

3. Interacting streams

Although steady-state mass balance holds for a stationary stream/ridge pair, it would be perturbed for the time-dependent problem of interacting pairs; but since we are only deducing their qualitative behavior in the following, this mass balance should remain adequate. For a more refined quantitative model, temporal change of the surface displacement obviously needs to be included in the mass balance.

3.1 Centerline separation

Let us imagine the initiation of the streams when the bed first becomes temperate as the ice sheet thickens. Without the topographic constriction, the streams would be randomly distributed and coalesce to form bigger ones as they grow (Kyrke-Smith and others, 2014). Through entrainment of the adjacent ice, they would self-organize into stream/ridge pair differentiated by wet/frozen beds, and the pair would grow to the stationary state discussed above if it is isolated from the neighboring pair. The situation can be visualized via a transect of the surface height across neighboring pairs as shown in Figure 5a. The left and right pairs are distinguished by solid and dashed lines, respectively, with the flat surfaces representing the streams and the curved surfaces, their catchment zones; the shaded strip is the ridge outside the catchment zones that enables the isolation of the two pairs. Such a configuration, however, is untenable for the following reason: the shaded strip, being outside the catchment zone, is not subjected to cold advection hence remains warm-based; it is thus ripe for initiating an embryonic stream to interact with the existing ones, contradicting their isolation.

Fig. 5. Transverse sections of the surface height of the interacting stream/ridge pairs. The left and right pairs are drawn in solid and dashed lines, with the flat and curve surfaces indicating the stream and the ridge, respectively. The centerline distance w _c (outermost vertical lines) decreases from (a) to (c). (a) The pairs are isolated from each other by an ice strip (shaded). (b) The streams vary between boxes D (slow) and C (fast) of Figure 6. (c) The streams vary between boxes A (stagnant) and C (fast) of Figure 6.

Fig. 6. Interaction between an active growing stream (solid arrow) and a passive narrowing one (dashed arrow), resulting in a cycle of: (a) smaller amplitude between fast/slow streams and (b) maximum amplitude between fast/stagnant streams.

If we denote the ‘centerline’ separation by $w_{\rm c} ^ {\hskip1pt{\prime}}$, it then has an upper bound of

(35)$$w_{\rm c} ^ {\hskip1pt{\prime}} = 2( w_1 ^ {\hskip1pt{\prime}} + w_0^ {\hskip1pt{\prime}} ) , \;$$

the right-hand side pertaining to the stationary stream/ridge pair derived in the previous section. Because of this upper bound, there is inherent interaction of the pairs with the widening one squeezing the other one. Because of this compensating behavior, the two stream beds obviously may not be weakening at the same time, so how does this reconcile with the minimization of the basal stress? For this, we posit the following hierarchy: the minimizing tendency applies only for the self-propelled stream growth; once it attains the maximum width, the tendency ceases and the stream can shrink passively if forced by the neighboring stream.

We visualize the interaction as shown in Figure 5b, which depicts the two extrema of the cycle within the stationary centerlines. As the left pair widens (the solid arrow), propelled by the weakening stream bed, it would narrow the right pair (the dashed arrow); when the left stream attains its maximum width, the right stream would be at its narrowest, which then begins to grow, reversing its prior passive role. One may also visualize the interaction on the phase space (Fig. 6a) whereby boxes D and C represent the two extrema of the cycle, and the solid and dashed arrows correspond to that shown in Figure 5b, namely, the active growth of the left stream from boxes D to C causes the passive shrinkage of the right stream from boxes C to D.

For a further decrease of the centerline separation, the narrowing stream may reach box A (Fig. 6b) of zero flow, which is characterized by half-width $w_2^ {\hskip1pt{\prime}}$ (28) and zero depression hence catchment width, as depicted in Figure 5c. The centerline separation thus is

(36)$$w_{\rm c}^ {\hskip1pt{\prime}} = w_1 ^ {\hskip1pt{\prime}} + w_0^ {\hskip1pt{\prime}} + w_2^ {\hskip1pt{\prime}} .$$

We shall next argue that this centerline separation, in fact, constitutes its lower bound. If this was not the case, the narrower stream would have stagnated before the wider stream reaches its maximum width, that is, the wider stream is not yet passive hence must be continually growing; and since the stagnation of the narrower stream would initiate an embryonic stream, both streams would be growing at the same time to push their centerlines apart, which cannot be stationary.

With the above, we have deduced upper and lower bounds of a stationary centerline separation given by Eqns (35) and (36), respectively; but for an initial value problem of the stream generation, once the lower bound is reached, the centerline separation has attained stationarity hence there is no tendency for its further expansion, or the system has reached an equilibrium.

3.2 Equilibrium state

With the above rationale, we thus define the equilibrium state of the interacting streams as defined by the stationary centerlines of minimum separation. Within the centerlines, the neighboring stream/ridge pairs exhibit compensating cycles of maximum amplitude, varying between the fastest (and widest) and stagnated (and narrowest) streams. To aid the visualization, we draw the time progression of the interacting streams in Figure 7, which will be compared later with that of the Ross ice streams B (left, also named after Whillans) and C (right).

Fig. 7. Time progression of the interacting streams, with their surfaces marked by stripes and maximum velocity by arrows. Following its reactivation, the left stream widens, deepens and speeds up. Squeezed by the left stream/ridge pair, the right stream narrows, shoals and slows to stagnation. Dashed line tracks the migration of the ridge, which is of constant height. The figure can be representative of the Ross ice streams B (left) and C (right) from 450 to 150 years ago.

The reactivation of the left stream causes it to widen and weaken its bed (the solid arrow from boxes A to C in Fig. 6b), both acting to speed up the flow (17); the increasing ice flux lowers its surface to entrain greater inward flux from the ridge (19) accompanied by growing catchment zone (10) until the stream bed reaches the minimum strength. The widening of the left stream/ridge pair would narrow the right stream and strengthen its bed (the dashed arrow from boxes C to A in Fig. 6b), being constricted to its right by the stationary centerline. This narrowing slows the stream and reduces the catchment flux until its surface rises to the ridge (box B). With the surface pinned to the ridge height, the narrowing continues to strengthen the bed until it reaches the driving stress and the flow comes to a halt (box A). Since the surface is already flush with the ridge and the bed remains temperate, it is preconditioned for reactivation to repeat the cycle.

To contrast the time progression of the ice streams with the surge cycles of the topography-confined glaciers, we juxtapose them in Figure 8. The widths (striped upper surfaces), the basal stress (striped lower surfaces), the surface height (solid lines) and the ice velocity (dashed lines) are drawn; the shaded columns signify the abrupt activation and termination of the surge.

Fig. 8. Time progression of the ice stream (left) and cyclic surge of the topographic glacier (right). Surface height in solid lines, velocity in dashed lines, striped top and bottom are width and basal stress, respectively. Shaded columns represent onset and termination of the surge characterized by short hydrological timescale. The figure shows the sharp contrast between the two cyclic variations.

At the surge onset, without the meltwater drainage, the basal stress quickly decreases to zero, causing the surge. The ensuing thinning augments the conductive cooling, which strengthens the bed and slows the flow. When the basal stress reaches about half the driving stress, the bed freezes over rapidly to enter the creep phase. There are, thus, three timescales characterizing the surge cycles: the onset and termination of the surge governed by the short hydrological timescale, the surge and creep durations characterized by the surface displacement divided by thinning rate and accumulation, respectively, with the former typically much shorter. For the stream, on the contrary, the meltwater drainage somewhat smooths the basal stress, so there is only one centennial timescale governing its time variation.

The differing causality of the surge from that of the stream manifests in sharply different relation among its state variables. During the surge, the causality goes from thinning to stronger bed to slower flow, resulting in a positive correlation between the surface height and the ice velocity. For the stream, on the contrary, the causality goes from widening to weaker bed to faster flow to lower surface, which results in anticorrelation between the surface height and the ice velocity. Because of these opposite correlations, they should be discernible from observations, as discussed next.

3.3 Observation

Although the Northeast Greenland Ice Stream has the moniker of an ice stream, it is steered by the bed trough (Joughin and others, 2001) to fall in the purview of Part 1 (Ou, 2021) wherein its deduced properties are discussed. Although the Siple coast does not have a flat bed, the relatively low relief does not strongly constrain the Ross ice streams, which represent unique examples of the ‘pure’ ice streams (Bentley, 1987; Bennett, 2003) for which our model may be applicable. Since our model is predicated on coupled stream/ridge pairs, its applications are further limited to ISB and ISC for the reason that ISA is bounded on one side by the Transantarctic Mountains, and that ISD and ISE have poorly developed interstream ridge before they merge into one stream. In addition, only ISB and ISC have their time evolution over the last millennium sufficiently documented (Hulbe and Fahnestock, 2007; Catania and others, 2012) to test our model.

Based on Figure 5c, we have drawn in Figure 7 the projected time evolution of ISB (left) and ISC (right), and the time proceeds from ~450 years to 150 years ago. It is observed that, following the reactivation of ISB, it would widen and thin while ISC would narrow, thicken and eventually stagnate, a prognosis that is indeed consistent with the observed scenario (Catania and others, 2012).

There is extensive discussion of the ISC stagnation, which is mostly viewed as a happenstance with myriad precedents (Anandakrishnan and others, 2001; Bennett, 2003). Among them, the ice and water piracies by ISB are inferred from the current surface relief (Alley and others, 1994), which, however, is likely a result rather than the cause of the ISC stagnation (Price and others, 2001). Additionally, we note that the meltwater cannot be transported across the frozen ridge between streams, and that the upper-glacier supply of the meltwater is not required for maintaining the stream, which after all is initiated by the local production of the meltwater; then the ice piracy is at odds with the observed thickening of ISC before its stagnation. As regards the drainage switch to the channelized system (Retzlaff and Bentley, 1993), there is no observational evidence of conduits under ISC and in fact the bed conditions of ISB and ISC are quite similar (Kamb, 2001). The loss of the meltwater would expose the sticky spots, as inferred from seismic data (Anandakrishnan and others, 2001), which however is a mediating agent rather than a cause of the stagnation. Although the grounding line has retreated in the last hundreds of years, it is likely a response to the stagnation (Hulbe and Fahnestock, 2007) since ice flow is relatively unaffected by the back-stress from the ice shelf (Anandakrishnan and Alley, 1997). However, others have likened the ISC stagnation to the surge termination by thinning-induced freeze-on (Hulbe and Fahnestock, 2004; Bougamont and others, 2011), which is at odds with the observation that the ISC bed remains wet (Bentley and others, 1998) and its stagnation is preceded by thickening – not thinning (Catania and others, 2012).

Instead, based on positive comparison of Figure 7 with observations, we offer the following explanation of the ISC stagnation. First, it is not a happenstance but a recurring event of the natural cycles. Second, the cycles are the outcome of invariable interaction among streams with the neighboring ones exhibiting compensating cycles. As such, the widening of ISB following its reactivation ~450 years ago would narrow ISC to augment the meltwater drainage, which then strengthens the bed to cause the stagnation. Although the loss of lubricant has been mentioned in the literature (Catania and others, 2012), it is largely overlooked, as attested by the above-referenced reviews, and our model provides a concrete affirmation of the mechanism.

There are additional points worth noting about ISC. Although its stagnation can be caused by narrowing-induced loss of the meltwater, as indicated by the dashed arrow in Figure 6b, the bed strengthening might not fully catch up with the narrowing to maintain the heat balance, so the stream might veer slightly into the cooling zone to effectuate partial freezing. This may possibly explain the accretional ice layers recently discovered under ISC (Vogel and others, 2005). To the degree that the meltwater drainage is operative when ISC stagnates, the bed must remain temperate and wet. This differs fundamentally from the surge termination, which, once commences, the positive feedback between the strengthening bed and reduced heating lead unceasingly to a thorough freeze-on (see Part 1); there is no allowance for a temperate and wet bed, contrary to the observed ISC.

4. Discussion

One key difference of flatbed ice streams from topography-confined glaciers considered in Part 1 (Ou, 2021) is the allowance of the meltwater dispersal, which has bended the solution continuum of Figure 2 to render a minimum bed strength. If the meltwater dispersal is artificially suppressed as is the case of Brinkerhoff and Johnson (2015), the stream would simply grow unabated until it bifurcates into limit cycles, which correspond to the cyclic-surge regime considered in Part 1 when the bed trough is sufficiently wide. A narrow trough, on the contrary, would arrest this growth to render a stationary state corresponding to the steady-sliding regime of Part 1. With inclusion of the meltwater dispersal, Kyrke-Smith and others (2014) have availed to them the downward branch of the continuum, and their closure amounts to fixing the bed strength in selecting the stream width, which, thus, is sensitive to the additional parameters related to the water film closure. Our stationary state, on the contrary, is set at the minimum bed strength, which being independent of the water-film mechanics renders a more robust width. Perhaps more significantly, we have postulated self-propelling of the streams toward this stationary state, which allows their interaction to exhibit compensating cycles, as seems the observed case. Since the stream boundary varies not because of instability but the outcome of deterministic cycle, there is no need to invoke drainage switch in stabilizing this boundary (Perol and others, 2015; Meyer and others, 2018) whose observational support remains tentative.

There is much discussion on the mass balances of Ross ice streams because of their effect on the global sea level, but if the interacting streams exhibit compensating cycles, their effect on the net mass balance would be reduced, as is the observed case (Shabtaie and Bentley, 1987). The equilibrium state of stationary centerlines is predicated on the width-compensation of neighboring stream/ridge pairs, which does not imply exact compensation in the ice transport because of its non-linear dependence on the stream width; a quantitative addressing of the problem, however, stretches the capability of our crude model.

Although our model has determined the width and spacing of the ice streams, we have not addressed how the whole pattern is being anchored. Myriad external controls of the stream location have been examined (Winsborrow and others, 2010), but our model suggests that some of the bed differences can be an outcome of the self-organization rather than progenitors of the disparate flows. On the contrary, if there are topographic troughs steering the upper-glacial tributaries, as seemingly the observed case (Joughin and others, 1999), they can provide an anchor for the whole stream/ridge pattern and there need not be external control of individual stream locales.

5. Summary

In Part 1 (Ou, 2021), we have considered the dynamics of a topography-confined glacier, which may undergo surge cycles due to thermal switches. In this Part 2, we have expanded the dynamical framework to address the ice streams over a flatbed. The expansion is twofold: first, without the topographic constriction, the stream would entrain the ambient ice, and the accompanying cold advection would freeze the ambient bed, resulting in self-organized stream/ridge pair; second, without the topographic trough curbing the meltwater dispersal, it would counter the conductive cooling to yield a minimum in the basal stress at some intermediate stream width. Invoking a positive feedback loop and the NT, we posit that the stream would evolve toward this minimum when it attains stationarity.

Although this stationary state exhibits certain sensitivity, it remains relatively robust over considerable ranges of the external condition. Crudely, the stream width is about twice the geometric mean of its height and length, the basal stress is less than half of the driving stress, the surface depression is about an order smaller than the ridge height, and its velocity is about two orders greater than the viscous creep; all are commensurate with the observed ones.

We next consider the invariable interaction of the stream/ridge pairs on a flatbed. We posit that they would attain an equilibrium configuration when the centerlines are the closest allowed whereby the neighboring streams exhibit compensating cycles of maximum flow and stagnation. The model deduction is consistent with the observed interaction of ISB and ISC over the last millennium, including the ISC stagnation ~150 years ago. Based on this positive comparison, the model suggests that the ISC stagnation is not a happenstance but a recurring natural cycle, and it uncovers a heretofore overlooked mechanism: namely, the widening ISB following its reactivation would narrow ISC to augment the loss of the meltwater, leading to its stagnation. Consistent with observations, this stagnation is preceded by thickening – opposite the thinning-induced surge termination. Although there could be incidental freezing as evidenced in the accretionary ice layer, the bed remains temperate and wet, as currently observed under ISC.

Although the model is extremely crude, it has nonetheless abridged more glaring empiricisms of the past models, and to the degree that the model deductions are consistent with observations, we suggest that the model has captured the governing dynamics of the flatbed ice streams. Together with Part 1 (Ou, 2021), we have presented a theoretical framework that integrates the dynamics of surging glaciers and self-organized ice streams while underscoring their generic difference due to the lateral dispersal of the meltwater.