Maximizing friction in the erosion of glacial valleys

Abstract

The cross-sectional profile of a glacial valley can be obtained with a variational principle in which the friction against the valley walls and the glacier bed is extremized, subject to a Lagrangian constraint. We show that the actual valley profile maximizes the friction, thus settling an old debate.

Introduction

Since the early days of geomorphology, it has been known that valleys carved by glaciers are U-shaped, as opposed to valleys carved by rivers which are V-shaped (Campbell, 1865; McGee, 1894). Focusing on the former, it is intuitive that following step by step the detailed erosion process of the valley walls and glacier bed through time requires many assumptions and details. This process can be modelled numerically (e.g. Harbor, 1995; Seddik and others, 2009; Yang and Shi, 2015). However, it is simpler to describe the final product of the glacier action, which can be done with an analytical approach first proposed by Hirano and Aniya (1988). From the phenomenological point of view, crude models of transverse glacial valley profiles used parabolas (Svensson, 1959) or power-law profiles, but there was no theory behind this so-called ‘power-flaw’, which amounted to mere data-fitting (Pattyn and Van Huele, 1998).

The Hirano and Aniya (1988) approach (refined by Morgan (2005)) uses a variational principle in which the friction of the ice against the valley walls, which is a functional of the transverse profile of the glacial valley, is extremized subject to an appropriate constraint. To fix the geometry, denote with x ∈ ( − x ₀, x ₀) a coordinate transverse to the glacier flow, while the margins of the ice have coordinates ±x ₀. Let the ice thickness at x be H(x). Correspondingly, the profile of the glacial valley is given by the bed elevation B(x) = H _s − H(x), where the constant H _s denotes the maximum ice thickness reached at x = 0 (see Fig. 1).

Fig. 1. The glacial valley is contained between −x ₀ and x ₀, the ice thickness (blue, dashed curve) is H(x) and is maximum at x = 0, where H = H _s. The bed elevation (red, solid curve) over the deepest point in the ice (valley profile) is B(x) = H _s − H(x).

As is clear from Figure 1, the cross-sectional area of the valley is

(1)$$A = 2x_0H_{\rm s}-\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} B\lpar x\rpar = \int_{{-}x_0}^{ + x_0} H {\rm d}x\comma \;$$

and the friction force is modelled by Coulomb's law,

(2)$$f = \mu N = \mu PA_{{\rm contact}} = \mu \rho _{{\rm ice}}gHA_{{\rm contact}}\comma \;$$

where N = ρ _icegHA _contact is the normal force caused by the cryostatic pressure, μ is the friction coefficient, ρ _ice is the ice density, g is the acceleration of gravity, and A _contact is the area of contact between the ice and the wall.

A word of caution about the assumptions of the Hirano–Aniya–Morgan theory is mandatory. The friction model may not be entirely appropriate, since friction could depend on the velocity and the Coulomb law adopted in this theory could be an oversimplification. Indeed, friction is taken to depend on the velocity in the numerical simulations of the dynamical development of glacial valley profiles of Seddik and others (2009).

Moreover, because normal stress should be used in the calculation of friction, the assumption of hydrostatic pressure may fail outside of the shallow ice approximation. It is possible that the model stretches the validity of its assumptions, but since the shallow ice approximation holds well for so many valley glaciers and is required extensively in glacier modelling anyway, perhaps this is not a big issue.

A second potential problem is that glacier hydrology is neglected in the Hirano–Aniya–Morgan scenario. Hirano and Aniya (1988, 1990, 2005) mention that the effective pressure may not be entirely cryostatic but proceeds by assuming that it is, citing older work in support. The issue was discussed in a debate (Harbor, 1990; Hirano and Aniya, 1988, 1990). Arguments supporting the model point out that water pressure is not significant along the entire length of the glacier, quoting similar results for a landslide block model, but the issue was not settled. Here we work in the context of the Hirano–Aniya–Morgan model, but the basic assumptions should be questioned further and the potential problems studied by dedicated work.

The a priori expectation that friction should be maximized invoking the second law of thermodynamics is too naive. In fact, the second law of thermodynamics says that the entropy of a closed system (or of the universe consisting of an open system plus its surroundings) never decreases. Apart from the fact that entropy and friction do not coincide, the frictional system considered here is not closed: it is clearly driven by climate forcing from the atmosphere, and energy is exchanged as heat within the glacier and the valley walls, through subglacial hydrology, etc. Driven systems have an energy source that allows the entropy to decrease. Moreover, although only the final state is considered here through a variational principle, the dissipative process itself involves non-equilibrium thermodynamics.

For an element of contact length ds between ice and rock and unit length downstream, the elementary friction force is

(3)$${\rm d}f = \mu \rho _{{\rm ice}}gH\;{\rm d}s = \mu \rho _{{\rm ice}}gH\sqrt {1 + {\lpar {{H}^{\prime}} \rpar }^2} \;{\rm d}x\;\comma \;$$

where a prime denotes differentiation with respect to x. The total friction is

(4)$$\mu \rho _{{\rm ice}}g\int_{{-}x_0}^{ + x_0} H \sqrt {1 + {\lpar {H}^{\prime}\rpar }^2} {\mkern 1mu} {\rm d}x{\mkern 1mu} \comma \;$$

a functional of the ice thickness H(x) (or, alternatively, of the glacial valley profile B(x)). A Lagrangian constraint must be imposed to keep the area of the cross-section constant (Morgan, 2005); considering a unit length downstream, this means that the volume of the ice is kept constant.¹ This requirement amounts to imposing $\int _{{-}x_0}^{ + x_0} {\rm d}x\;H\lpar x \rpar = {\rm const} .$ Extremizing the total friction subject to this constraint then yields

(5)$$\delta J\equiv \delta \int_{{-}x_0}^{ + x_0} {\rm d} x\left[{H\sqrt {1 + {\lpar {H}^{\prime}\rpar }^2} -\lambda H} \right]\equiv \delta \int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} L\lpar {H\comma \;{H}^{\prime}} \rpar = 0{\mkern 1mu} \comma \;$$

where λ is a Lagrange multiplier, J[H(x)] is a functional of the ice thickness H(x), and the constants μ, ρ _ice and g have been subsumed into λ.

The Lagrangian L does not depend explicitly on x and the corresponding Hamiltonian is conserved:

(6)$${\cal H} = \displaystyle{{-H} \over {\sqrt {1 + {\lpar {H}^{\prime}\rpar }^2} }} + \lambda H = C{\mkern 1mu} \comma \;$$

where C is an integration constant. We will refer to Eqn (6) as the Morgan equation (Morgan, 2005). The existence of smooth solutions H(x) in ( − x ₀, x ₀) with the properties

(7)$$H\lpar { \pm x_0} \rpar = 0\;\comma \;\; \; \; \; \; H\lpar x \rpar = H\lpar {-x} \rpar \;\comma \;\; \; \; \; \; {H}^{\prime}\lpar 0 \rpar = 0\;\comma \;$$

requires H ^′′(x) < 0, λ > 0 and C > 0 (Morgan, 2005). Analytical solutions of the Morgan equation (6) have been found in Harbor (1990); Morgan (2005); further formal solutions can be found in Chen and others (2015a) using the fact that the Morgan equation is analogous to the Friedmann equation of cosmology (Faraoni and Cardini, 2017). Moreover, all solutions of (6) are roulettes (Chen and others, 2015b) (a roulette is the trajectory described by a given point lying on a curve that rolls without slipping along another given curve).

In their original paper, Hirano and Aniya (1988) argued that friction should be minimized during the erosion process. Their method and conclusions were criticized in Harbor (1990) (see also Hirano and Aniya (1990); Morgan (2005); Hirano and Aniya (2005)). In particular, Harbor (1990) argued that friction should be maximized, not minimized. Although (after imposing the correct Lagrangian constraint) this change of perspective does not affect the first-order variational principle δJ[H(x)] = 0, it would seem intuitive that friction should be maximized instead of minimized (perhaps entropy production rate is maximized, but the previous authors do not provide sound arguments to support a maximum or a minimum). However, in problems involving geophysical flows in the earth sciences, sometimes the opposite point of view is supported, based on the argument that minimizing friction optimizes the flow. This is the case, for example, of equilibrium beach profiles in oceanography (see Jenkins and Inman, 2006; Faraoni, 2019; Maldonado and Uchasara, 2019; Maldonado, 2020). The issue of whether the final profile for a glacial valley corresponds to a minimum or a maximum of the friction (subject to the constraint above), which could help understanding better the physics behind the model, has never been settled in the literature on glacial valley transverse profiles. Here we show that the final configuration indeed corresponds to a maximum of the friction by considering the second variation of the functional J[H(x)] of the valley profile and establishing its sign.

Maximizing friction

As seen in the previous section, the action functional is (Morgan, 2005)

(8)$$J = \int_{{-}x_0}^{ + x_0} {\rm d} x\left({H\sqrt {1 + {{H}^{\prime}}^2} -\lambda H} \right)\equiv \int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} L\lpar {H\comma \;{H}^{\prime}} \rpar {\mkern 1mu} .$$

Let us consider variations η(x) around the path that extremizes the action J, which vanish at the endpoints, η( − x ₀) = η(x ₀) = 0 and are parametrized by a parameter α. More precisely, the varied paths are given by (Weber and Arfken, 2004)

(9)$$H\lpar {x\comma \;\alpha } \rpar = H\lpar {x\comma \;0} \rpar + \alpha \;\eta \lpar x \rpar \;\comma \;$$

where H(x, 0) is the path extremizing J[H(x, α)]. Using

(10)$${H}^{\prime}\lpar {x\comma \;\alpha } \rpar = {H}^{\prime}\lpar {x\comma \;0} \rpar + \alpha \;\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}\;\comma \;$$

(11)$$\displaystyle{{\partial H} \over {\partial \alpha }} = \eta \;\comma \;\; \; \; \displaystyle{{\partial {H}^{\prime}} \over {\partial \alpha }} = \displaystyle{{{\rm d}\eta } \over {{\rm d}x}}\;\comma \;\; \; \; \displaystyle{{\partial ^2{H}^{\prime}} \over {\partial \alpha ^2}} = \displaystyle{{\partial ^2H} \over {\partial \alpha ^2}} = 0\;\comma \;$$

one obtains (Weber and Arfken, 2004)

(12)$$\displaystyle{{\partial ^2J} \over {\partial \alpha ^2}} = \int_{{-}x_0}^{ + x_0} {\rm d} x\left[{\displaystyle{{\partial^2L} \over {\partial {{H}^{\prime}}^2}}{\left({\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}} \right)}^2 + 2{\mkern 1mu} \displaystyle{{\partial^2L} \over {\partial H\partial {H}^{\prime}}}{\mkern 1mu} \eta {\mkern 1mu} \displaystyle{{{\rm d}\eta } \over {{\rm d}x}} + \displaystyle{{\partial^2L} \over {\partial H^2}}{\mkern 1mu} \eta^2} \right]{\mkern 1mu} .$$

Substituting the derivatives

(13)$$\displaystyle{{\partial L} \over {\partial H}} = \sqrt {1 + H{\prime}^2 } -\lambda \;\comma \;$$

(14)$$\displaystyle{{\partial ^2L} \over {\partial H^2}} = 0\;\comma \;$$

(15)$$\displaystyle{{\partial ^2L} \over {\partial H\partial {H}{\prime}}} = \displaystyle{{{H}{\prime}} \over {\sqrt {1 + H{\prime}^2 } }}\;\comma \;$$

(16)$$\displaystyle{{\partial L} \over {\partial {H}{\prime}}} = \displaystyle{{H{H}{\prime}} \over {\sqrt {1 + H{\prime}^2 } }}\;\comma \;$$

(17)$$\displaystyle{{\partial ^2L} \over {\partial H{\prime}^2 }} = \displaystyle{{H\lsqb {{H}^{\prime \prime}\lpar {1 + H{\prime}^2 } \rpar -H{\prime}^2 } \rsqb } \over {{\lpar {1 + H{\prime}^2 } \rpar }^{3{\rm /}2}}}\;\comma \;$$

one obtains

$$\displaystyle{{\partial ^2J} \over {\partial \alpha ^2}} = \int_{{-}x_0}^{ + x_0} {\rm d} x\left\{{\displaystyle{{H\lsqb {{H}^{\prime \prime}\lpar {\lpar 1 + {{H}^{\prime}}^2} \rpar -{{H}^{\prime}}^2} \rsqb } \over {{\lpar {1 + {{H}^{\prime}}^2} \rpar }^{3/2}}}{\mkern 1mu} {\left({\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}} \right)}^2\!\!\! + \displaystyle{{{H}^{\prime}} \over {\sqrt {1 + {{H}^{\prime}}^2} }}{\mkern 1mu} \displaystyle{{{\rm d}\lpar \eta^2\rpar } \over {{\rm d}x}}} \right\}{\mkern 1mu} \comma \;$$

where the integral is now computed along the trajectories extremizing J, i.e. those that realize the condition ∂J/∂α = 0. We know that $H\lpar x \rpar \ge 0\; \forall x\in \lsqb {-x_0\comma \;x_0} \rsqb$ and H = 0 only at the boundaries x = ±x ₀, and that ${H}^{\prime \prime}\lpar x \rpar \lt 0\; \forall x$ (see Fig. 1), therefore HH ^′′(1 + H ^′2) < 0 for x ∈ ( − x ₀, x ₀). Using also ${-}H\lpar x \rpar H{\prime}^2 \lpar x \rpar \le 0\; \forall x$, one obtains

(19)$$\displaystyle{{H\lsqb {{H}^{\prime \prime}\lpar {\lpar 1 + {{H}^{\prime}}^2} \rpar -{{H}^{\prime}}^2} \rsqb } \over {{\lpar {1 + {{H}^{\prime}}^2} \rpar }^{3/2}}}{\mkern 1mu} \left({\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}} \right)^2 \le 0\;\;\;\;\;\;\forall x\in \lpar {-x_0\comma \;x_0} \rpar .$$

Therefore, the first integral in Eqn (18) is non-positive,

(20)$$\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} \displaystyle{{H\lsqb {{H}^{\prime \prime}\lpar {\lpar 1 + {{H}^{\prime}}^2} \rpar -{{H}^{\prime}}^2} \rsqb } \over {{\lpar {1 + {{H}^{\prime}}^2} \rpar }^{3/2}}}{\mkern 1mu} \left({\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}} \right)^2 \le 0{\mkern 1mu} .$$

Let us evaluate now the sign of the second integral

(21)$$\int _{{-}x_0}^{ + x_0} {\rm d}x\;\displaystyle{{{H}^{\prime}} \over {\sqrt {1 + H{\prime}^2 } }}\;\;\displaystyle{{{\rm d}\lpar {\eta^2} \rpar } \over {{\rm d}x}} = \displaystyle{{\eta ^2\;{H}^{\prime}} \over {\sqrt {1 + H{\prime}^2 } }}\vert {} _{{-}x_0}^{ + x_0} -\int _{{-}x_0}^{ + x_0} {\rm d}x\;\eta ^2\;\displaystyle{{\rm d} \over {{\rm d}x}}\left({\displaystyle{{{H}^{\prime}} \over {\sqrt {1 + H{\prime}^2 } }}} \right)\;\comma \;$$

where we have integrated by parts. The first term on the right-hand side vanishes because the variation η(x) vanishes at the endpoints x = ±x ₀. To proceed, use the equation of motion (6) satisfied by the trajectories H(x) that extremize J, where the constant C is positive (Morgan, 2005; Faraoni and Cardini, 2017). Substituting into the integral (21) yields

(22)$$\displaylines{I\equiv{-}\int_{{-}x_0}^{ + x_0} {\rm d} x\eta ^2\displaystyle{{\rm d} \over {{\rm d}x}}\left({\displaystyle{{{H}^{\prime}} \over {\sqrt {1 + {{H}^{\prime}}^2} }}} \right)\cr = {-}\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} \eta ^2\displaystyle{{\rm d} \over {{\rm d}x}}\left[{\left({\lambda -\displaystyle{C \over H}} \right){H}^{\prime}} \right]\cr = {-}\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} \eta ^2\displaystyle{{C{{H}^{\prime}}^2} \over {H^2}} + \int_{{-}x_0}^{ + x_0} {\rm d} x\left({\displaystyle{C \over H}-\lambda } \right){H}^{\prime \prime}\equiv I_1 + I_2.} $$

Since the integrand appearing in the first integral is manifestly non-negative, it is I ₁ ≤ 0. Then,

(23)$$I_2 = \int_{{-}x_0}^0 {\rm d} x\left({\displaystyle{C \over H}-\lambda } \right){H}^{\prime \prime} + \int_0^{ + x_0} {\rm d} x\left({\displaystyle{C \over H}-\lambda } \right){H}^{\prime \prime}{\mkern 1mu} \semicolon \;$$

since C/H → +∞ as the valley boundaries ± x ₀ are approached and H(x) → 0⁺, both integrals on the right-hand side of Eqn (22) (and, therefore, their signs) are dominated by the regions around ±x ₀, where ${C \over H}-\lambda \simeq {C \over H}$ and λ is bounded. Since ${H}^{\prime \prime}\lpar x \rpar \lt 0\; \forall x\in \lpar {-x_0\comma \;x_0} \rpar$, both integrals are negative and it is I ₂ < 0.

Putting everything together, we finally have

(24)$$\displaylines{\displaystyle{{\partial ^2J} \over {\partial \alpha ^2}} = \underbrace{{\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} \displaystyle{{H\lsqb {{H}^{\prime \prime}\lpar {1 + {{H}^{\prime}}^2} \rpar -{{H}^{\prime}}^2} \rsqb } \over {{\lpar {1 + {{H}^{\prime}}^2} \rpar }^{3/2}}}{\left({\displaystyle{{{\rm d}\eta } \over {{\rm d}x}}} \right)}^2}}_{{\;\lt 0\;}} \cr \underbrace{{-\int_{{-}x_0}^{ + x_0} {\rm d} x{\mkern 1mu} \displaystyle{{C{{H}^{\prime}}^2\eta ^2} \over {H^2}}}}_{{\;\lt 0\;}} \cr \underbrace{{ + \int_{{-}x_0}^{ + x_0} {\rm d} x\left({\displaystyle{C \over H}-\lambda } \right){H}^{\prime \prime}}}_{{\;\lt 0\;}} \cr \hskip-6pc \lt 0{\mkern 1mu} .} $$

This is our main result: ∂²J/∂α ² < 0 and the path extremizing J corresponds to a maximum of the friction functionalJ.

Conclusions

The final profile of a glacial valley due to the erosion by a glacier over periods of time of the order of 10 000 years can be obtained through a variational principle that extremizes the friction of the ice against the valley walls and bed, subject to the constraint of fixed ice volume (Hirano and Aniya, 1988). The correct variational principle leading to the Morgan equation (6) and its analytical solutions were obtained only after a few attempts (Harbor, 1990; Hirano and Aniya, 1990, 2005; Morgan, 2005; Chen and others, 2015a; Faraoni and Cardini, 2017). A problem debated in Harbor (1990); Hirano and Aniya (1990); Morgan (2005); Hirano and Aniya (2005) and still open after all these years is whether the friction is maximized or minimized. This issue, which appears in similar situations involving geophysical flows, generates conflicting intuitive answers and can only be settled once and for all by calculation. By studying the second variation of the relevant action functional (8), we have established that the glacial valley profile indeed corresponds to a maximum of the friction.