Analogues of glacial valley profiles in particle mechanics and in cosmology

It has long been acknowledged in glaciology (since its inception) (Campbell 1865; McGee 1894), and it is common knowledge in elementary geography, that valleys carved by glaciers are U-shaped, whereas valleys carved by rivers are V-shaped. Here we focus on the former. The detailed and continued process of reshaping a valley by a glacier via erosion of the valley walls and floor over time is complex and is best modelled with numerical techniques (Harbor 1995; Seddik et al. 2009; Yang and Shi 2015). If one is interested only in the final result of the glacier action, simpler analytical approaches can be used. Given the scarcity of analytical models in the literature, theoretical approaches to this problem are valuable. Hirano and Aniya (1988) proposed a variational principle extremizing the friction of ice against the valley walls, subject to an appropriate constraint. Let the cross-sectional profile of a glacial valley be described by a function y(x), where x is a coordinate transverse to the glacier flow. Hirano and Aniya (1998) argued that friction (a functional of the cross-profile y(x)) should be minimal at the end of the erosion process, subject to the constraint that the contact length of the cross-profile of the ice is constant. This contact length between two endpoints x ₁ and x ₂ of the transverse profile iswhere a prime denotes differentiation with respect to x. The friction force is modelled by Coulomb’s law as f = μN, where μ is the friction coefficient and the normal force is N = ρghA ₁. Here ρ is the ice density, g is the acceleration of gravity, h is the ice thickness, and A ₁ is the area of contact between the ice and the bed.¹ By considering a unit width of ice in the longitudinal direction, the friction force due to an element of contact length ds is df = μPds. Further, P = η(y_s  − y), where y_s denotes the ice surface and η is a constant. Extremizing the frictionsubject to the constraint (eq. 1.1) leads towhere λ is a Lagrange multiplier and L is the Lagrangian. The Euler–Lagrange equationyields the ordinary differential equation (ODE) (Hirano and Aniya 1988)where c ₁ is a constant. This equation is also obtained by solving the classic catenary problem of mechanics (Goldstein 1980), and therefore, it is not surprising that Hirano and Aniya (1988) obtained catenary solutions of eq. (1.5). Their method and conclusions have been criticized by Harbor (1990; see also Hirano and Aniya 1990; Hirano and Aniya 2005; Morgan 2005). First, the assumptions in the model of Hirano and Aniya (1988) are inconsistent with common assumptions of glaciology (Harbor 1990). Second, friction should be maximized, not minimized (Harbor 1990; although this change does not affect the first-order variational principle, which only requires the friction integral to be extremized). In a reply, Hirano and Aniya (1990) agreed on this point but stood by the validity of application of the variational principle and of their previous result. Further critique (Morgan 2005) pointed out that there is no physical basis for the isoperimetric constraint (eq. 1.1), which should be replaced by the requirement that the area of the cross-section of the glacial valley is fixed. The rationale is that by considering a unit width of ice in the direction of longitudinal flow, the ice volume is thus kept constant (Harbor 1990).

The new Lagrangian constraint in the variational principle leads (Morgan 2005) to the ODEwhere y(x) is now the ice thickness at transverse coordinate² x, λ is again a Lagrange multiplier, and C is a constant, with λ > 1 and C > 0 required to have a smooth symmetric solution y(x) on the interval [−x ₀, x ₀] with y′(0) = 0 (Morgan 2005). Equation (1.6) is adopted in the rest of this work to describe the transverse profiles of glacial valleys. An exact solution of eq. (1.6) was provided by Morgan (2005):where

Another formal solution of eq. (1.6) for λ ² < 1 has been given recently by Chen et al. (2015a):where D is an integration constant and |w| > 1. For C = 0 and |λ| < 1, the solutions are linear.

The constraint of fixed cross-sectional area of the valley may perhaps seem questionable but no better constraint has been proposed in the literature thus far. The requirement of fixed cross-sectional area is also used in the numerical modelling of the erosion process leading to U-shaped valleys (Seddik et al. 2009). In any case, some constraint must be imposed, because if the friction integral is maximized without constraints, the first-order variationwhere the Lagrangian is nowproduces an equation that only admits unphysical solutions. In fact, as the Lagrangian (eq. 1.12) does not depend explicitly on x, the corresponding Hamiltonian $\mathcal{\mathscr{H}}$ ₀= p _y y′ − L ₀ is conserved, whereis the momentum canonically conjugated to y. The Euler–Lagrange equation (eq. 1.4) for L ₀ has the first integralwhere C is a constant. Using the variable ξ ≡ y _s  − y > 0, eq. (1.14) is equivalent to $C\sqrt{1+{({\xi }^{\prime })}^2}+\xi =0$ , which requires that C < 0; hence we set C ≡ −C ₂, where C ₂ > 0 has the dimensions of a length. One obtainswhich requires ξ ≥ C ₂ > 0. All the solutions ξ(x) of eq. (1.15) are unbounded from above and are given by(with x ₀ another integration constant) and require |x| ≥ C ₂, which does not describe a valley geometry. Therefore, some Lagrangian constraints must be imposed when extremizing the friction integral f [y(x)].

The fact that eq. (1.6) has an analogue in point particle mechanics was missed in the glaciology literature, whereas the fact that it has an analogue in the Friedmann equation of cosmology was recently noted in passing by Chen et al. (2015a, 2015b). Equation (1.6) is a special case of Friedmann-type equations, which are of fundamental importance in cosmology. A mathematical peculiarity of this type of equations demonstrated by Chen et al. (2015b) is that the graphs of all the solutions (in our case, of the transverse valley profiles y(x)) are roulettes. A roulette is the locus of a point that lies on, or inside, a curve that rolls without slipping along a straight line.³

Our goal is to explore the analogues of the ODE (eq. 1.6) in point particle mechanics and in cosmology, obtaining insight into the properties of this equation and establishing the basics of the formal analogies for future reference. We uncover a type of cosmological singularity, which was studied recently (Bouhmadi-López et al. 2008) for spatially flat universes in the now abundant literature on cosmological singularities (Wald 1984; Barrow et al. 1986; Caldwell 2002; Shtanov and Sahni 2002; Barrow 2004; Nojiri et al. 2005; Dabrowski et al. 2007; Fernandez-Jambrina 2007; Dabrowski and Denkiewicz 2009; Frampton et al. 2011; Bouhmadi-López et al. 2015; Beltran Jimenez et al. 2016).

Following Morgan (2005), we assume C > 0 (but see the discussion below). The ODE (eq. 1.6) can be rewritten aswhere (y′)²/2 is the kinetic energy of a particle of unit mass in one-dimensional motion if x and y(x) are the analogues of time and position, respectively, whereasis an effective potential energy and E = −1/2 is the total mechanical energy of the particle. Newton’s second law y″ = −dV/dy then rules the motion of the particle and eq. (2.1) is a first integral corresponding to conservation of energy (with the value E = −1/2). The possible motions y(x) are candidate analogues of the cross-profiles of glacial valleys. A qualitative understanding of the motion can be obtained from the potential V(y) and its intersections with the horizontal line E = −1/2 (Weierstrass approach; Bochicchio and Laserra 2007; Destrade et al. 2007; Bochicchio et al. 2011).

It is V(y) < 0 ∀y ≠ 0, V(0) = 0, V(y) → −1/(2λ ²) as y → ±∞, there is a vertical asymptote y = C/λ > 0, and V′(y) = Cy/(λy − C)³, therefore the function V(y) increases for y < 0 and for y > C/λ, decreases for 0 < y < C/λ, and is maximum at y = 0 (Fig. 1). We look for regions of bounded motions y(x) > 0 corresponding to finite ice thickness. We restrict to the situation C > 0, λ > 0 for which the vertical asymptote y = C/λ of V(y) lies in the y > 0 region. However, eq. (1.6) and the potential V(y) are invariant under the exchange (C, λ) → (−C, −λ).

Relativistic cosmology (Wald 1984; Liddle 2003; Carroll 2004) is obtained by assuming that the four-dimensional spacetime of general relativity is spatially homogeneous and isotropic about every point of three-space, which leads uniquely to the Friedmann–Lemaître–Robertson–Walker (FLRW) line element (Wald 1984; Liddle 2003; Carroll 2004)in comoving polar coordinates (t, r, θ, φ), where the curvature index is normalized to k = 0,±1, in units in which the speed of light c is unity. The matter content of the universe is assumed to be in the form of a single perfect fluid with energy density ρ and pressure P. The functions a(t), ρ(t), and P(t) satisfy the Einstein–Friedmann ODEs (Wald 1984; Liddle 2003; Carroll 2004)where an overdot denotes differentiation with respect to t and $H(t)\equiv \dot{a}/a$ is the Hubble function, and G is Newton’s constant. A third convenient (but not independent) equation expressing covariant conservation of energy is (Wald 1984; Liddle 2003; Carroll 2004)

The Friedmann equation (eq. 3.2) is the analogue in cosmology of eq. (1.6) describing the ice thickness in glacial valley transverse profiles. The analogy was recently noted, but not pursued, by Chen et al. (2015b). The most straightforward identification between eqs. (1.6) and (3.2) is achieved by setting k = 1, which corresponds to a closed universe. Then the cosmological analogue of eq. (1.6)can be rewritten aswhereare positive constants. Equation (3.2) gives the energy density of the analogue cosmic fluid asWe stress that we study the analogy with possible theoretical models of the future universe dominated by an exotic dark energy, which abound in the literature, and not the actual universe at the present time, which is indeed dominated by dark energy but also contains a non-negligible proportion of dark matter. Other details of the present universe cannot be reproduced by a single equation in the analogy with glacial valley profiles.

Two properties of the density ρ(t) are relevant. First, the energy density is always positive, which is expected of “reasonable” forms of matter but is not at all guaranteed in any formal analogy. Second, the density diverges if a → a ₀, corresponding to a spacetime singularity.

“Reasonable” matter in general relativity is supposed to satisfy energy conditions that prohibit negative energy densities and energy flows faster than light (Wald 1984; Liddle 2003; Carroll 2004). When applied to a perfect fluid, the energy conditions are the weak energy condition (ρ ≥ 0 and ρ + P ≥ 0), the strong energy condition (ρ + P ≥ 0 and ρ + 3P ≥ 0), the null energy condition (ρ + P ≥ 0) and the dominant energy condition (ρ ≥ 0 and ρ ≥ |P|; Wald 1984; Liddle 2003; Carroll 2004).

The effective pressure P of the analogue cosmic fluid is deduced by imposing eq. (3.4), which yieldsand can be rewritten aswhere the upper sign applies when a > a ₀ and the lower sign when a < a ₀. Using eqs. (3.8)–(3.10), one concludes easily that for a > a ₀, the weak, strong, and null energy conditions are always satisfied and the dominant energy condition is satisfied only when a ≥ 3a ₀/2. When a < a ₀, all the energy conditions are violated. Equations of state of the cosmic fluid corresponding to the lower sign in eq. (3.10) and violating the weak energy condition have been discussed by Barrow (2004). Equations of state of the cosmic fluid of the form $P={\sum }_{k=1}^m{c}_k{\rho }_{(k)}^k$ have been studied by Chen et al. (2015a, 2015b), Stefancic (2005), Frampton et al. (2011) and Bouhmadi-López et al. (2015). Quadratic equations of state, in particular, have been the subject of further attention (Nojiri and Odintsov 2004, 2005; Ananda and Bruni 2006a, 2006b; Capozziello et al. 2006; Silva e Costa 2009). Pressures depending on fractional powers of the density have been studied by Nojiri et al. (2005). As discussed below, they give rise to a peculiar type of singularity. Although traditional cosmology and relativity textbooks report only linear barotropic equations of state P = P ₀ + A ρ, following the discovery of the acceleration of the universe in 1998, the literature abounds with exotic non-linear equations of state for the dark energy fluid postulated to explain this acceleration. Let us consider now the acceleration equation (eq. 3.3), which, using eq. (3.9), reduces toClearly, the universe is accelerating if a < a ₀ (corresponding to ρ + 3P < 0) and decelerating if a > a ₀ (corresponding to ρ + 3P > 0). The value a ₀ of the scale factor corresponds to a spacetime singularity. In fact, the Einstein equation R _ab  − (1/2)g _ab R = 8π GT _ab (where R _ab is the Ricci tensor) gives R ≡ g ^ab R _ab  = −8π G T, where T is the trace of the fluid energy momentum tensor (Wald 1984; Liddle 2003; Carroll 2004)with four-velocity u ^a . The Ricci scalar is, therefore,which diverges in the limit a → a ₀ signalling a spacetime singularity. It is not yet clear whether the value a ₀ of the scale factor can actually be approached during the dynamics. To answer this question, note that eq. (3.2) can be written asexcluding the orbits of the solutions of the dynamical system (eqs. 3.2 and 3.3) from a certain volume of the (a, $\dot{a}$ ) phase space.

If a > a ₀, the dynamical constraint (eq. 3.14) can be written as a ≤ λ a ₀/(λ−1). The coefficient of a ₀ is λ/(λ−1) = 1 + 1/(λ−1) > 1, and therefore, in this regime we have a ₀ < a < λ a ₀/(λ − 1); the scale factor a(t) is bounded from above but it can get arbitrarily close to the value a ₀ corresponding to the singularity.

If instead a < a ₀, then the constraint (eq. 3.14) becomes a ≥ λ a ₀/(λ + 1) and one has 0 < λ a ₀/(λ + 1) ≤ a < a ₀. In both cases, the scale factor is bounded from below by a positive constant (hence one cannot have a Big Bang- or Big Crunch-type singularity corresponding to a → 0; Wald 1984; Liddle 2003; Carroll 2004) but it can reach the singularity⁵ a ₀. It corresponds to the fact that y = 0 (equivalently a = 0) lies in the region of the phase space forbidden by the constraint (eq. 3.14).

The boundary values a = λ a ₀/(λ∓1) are formal solutions of eq. (3.5) obtained by setting a = constant, but they do not satisfy eq. (3.3) (they would be meaningless as analogues of glacial valley profiles).

Consider again the acceleration equation (eq. 3.3): if a > a ₀, then $\ddot{a}$  < 0 and the curve representing the scale factor has concavity facing downwards. As this curve is continuous, it always decreases and eventually crosses the horizontal line a = a ₀.

Using eq. (3.7), eq. (3.5) is written asas a → a ₀. This asymptotic equation is easily integrated towhere the integration constant t ₀ has the meaning of time at which the singularity occurs and the positive sign of the square root must be chosen because a > a ₀. This situation constitutes a physically meaningful analogue of glacial valleys because y is the thickness (maximum at x = 0 and minimum at the valley boundaries; Morgan 2005) and it is interesting in cosmology because it provides an example of a finite-time singularity even when ρ > 0 and ρ + P > 0 (i.e., without violating the weak energy condition), an interesting situation discussed by Barrow (2004). Vice versa, if a < a ₀ then $\ddot{a}$  > 0 and a(t) always increases, eventually crossing the horizontal line a = a ₀. The asymptotic equation (eq. 3.15) now integrates tochoosing the negative sign in front of the square root because now a < a ₀. This situation is not a meaningful analogue of a glacial valley cross-profile. The slope of this function becomes infinite where a(t) → a ₀, corresponding to a spacetime singularity in the cosmological analogue. In this case, the universe has a minimum size a _min = λ a ₀/(λ + 1) and it bounces upon reaching it. The bounce occurs when the weak energy condition is violated, which is exactly what is happening here because ρ + P = 2ρ ₀ a/3(a − a ₀)³ < 0. The violation of the weak energy condition signals a very exotic phantom fluid, which causes the universe to accelerate with Hubble function increasing according to $\dot{H}=-4\pi G(\rho +P)+k/a^2>0$ if k ≥ 0. This superacceleration causes the universe to expand superexponentially and reach a singularity at a finite time. In the standard cosmological literature, the scale factor of a phantom-dominated universe diverges at a finite time in the future at a Big Rip singularity (Caldwell 2002), but here the situation is different because the scale factor stays finite while the Hubble function H, energy density ρ, pressure P, and Ricci scalar R all diverge as a → a ₀. This situation corresponds instead to a type of singularity studied recently in spatially flat universes and called a Big Freeze singularity (Bouhmadi-López et al. 2008) or Type III singularity in the classifications of Nojiri et al. (2005) and Bamba et al. (2012). A Big Freeze singularity also appears in cosmology in the context of Palatini f(R) gravity (Borowiec et al. 2012, 2016; Szydlowski et al. 2016; Stachowski et al. 2016). A Big Freeze singularity was not previously reported for positively curved universes. The situation is essentially the same as for spatially flat universes because, as a → a ₀ in eq. (3.6), the divergent term proportional to (a − a ₀)⁻² in the energy density dominates over the curvature term −1/a ², which stays finite. Finite-time singularities, including Big Rip (Caldwell 2002) and sudden future singularities (Barrow et al. 1986; Barrow 2004), have been the subject of a significant amount of work in cosmology (Shtanov and Sahni 2002; Kofinas et al. 2003; Calcagni 2004; Gorini et al. 2004; Nojiri et al. 2005; Stefancic 2005; Dabrowski et al. 2007; Fernandez-Jambrina 2007; Bouhmadi-López et al. 2008; Dabrowski and Denkiewicz 2009; Frampton et al. 2011; Bamba et al. 2012; Bouhmadi-López et al. 2015; Beltrán Jiménez 2016).

Finally, although not interesting for the original glacial valley problem, the cosmic analogue of eq. (1.6) for C = 0 (noted also by Chen et al. (2015b)) corresponds to a ₀ = 0 andand requires |λ| ≤ 1. Then $\dot{a}=\mathrm{constant}$ and the solutions are linear in time and include a static universe with a = constant as a special case.

In glacial morphology studies, researchers content themselves with fitting data of glacial valley ice thickness with parabolas y(x) = a x ² + b x + c (following an early practice initiated by Svensson (1959), which is not free of critique (Pattyn and Van Huele 1998)). Other fitting curves used include power-law profiles y(x) = a x ^b , possibly with different powers b for each half-profile going from the bottom at x = 0 to each valley side. There is a deep disconnect between theory and practice here. As a parabola is just the second-order Taylor expansion of any sufficiently regular function with a minimum, which could solve infinitely many ODEs, data fitting with parabolas is of no help when one attempts to test models and to discriminate among theoretical approaches to the problem of the cross-sectional profiles of glacial valleys, which predict different ODEs for the profile y(x). In this sense, fitting parabolas or power-laws and determining best-fit parameters is deeply unsatisfactory from the theoretical point of view and does not contribute to understanding the mechanism that generated these profiles.

Here we focused on the variational principle approach to the problem of glacial valley erosion, in the form given by Morgan (2005). The resulting equation (eq. 1.6) for the ice thickness y(x) has analogues in point particle mechanics and in cosmology, which we explored in detail. These analogies contribute to a better understanding of this ODE, its solutions, and the conditions (on the parameters λ and C) for their physical viability. The previous analysis of Morgan (2005) is augmented by the graphical study of the effective potential V(y) in the mechanical analogy. The problem of glacial valley profiles provides also an interesting example of a finite-time singularity in current theoretical cosmological models without violating the weak energy condition. The finite-time singularity is caused by the peculiar effective equation of state (eq. 3.10), which falls into the broader category P = w ρ − α ρ ^m (with w, α, and m constants), which has been the subject of wide interest in cosmology⁶ (Shtanov and Sahni 2002; Kofinas et al. 2003; Barrow 2004; Calcagni 2004; Gorini et al. 2004; Nojiri et al. 2005; Stefancic 2005; Frampton et al. 2011; Bouhmadi-López et al. 2015).

The present work also highlights two open problems in glaciology. First, the friction of glacier ice against the valley walls and floor is unlikely to be described purely by Coulomb’s law, but it should include viscous friction that depends on the velocity (a friction model quadratic in the velocity is used, for example, in the numerical work of Seddik et al. (2009)). Second, the numerical analyses of the formation of glacial valleys ignore the variational approach (but not its Lagrangian constraint) and should be compared with it. These issues will be revisited in future work.

This work is supported by the Natural Sciences and Engineering Research Council of Canada.