Solving the Vialov equation of glaciology in terms of elementary functions

TL;DR

This paper derives many new exact solutions to the Vialov equation in glaciology using elementary functions, expanding the limited set of known solutions for glacier profile modeling.

Contribution

It introduces a novel application of Chebyshev's theorem to find elementary solutions for the Vialov equation, broadening analytical tools in glaciology.

Findings

Many new exact solutions derived

Solutions expressed in elementary functions

Enhanced understanding of glacier profile equations

Abstract

Very few exact solutions are known for the non-linear Vialov ordinary differential equation describing the longitudinal profiles of alpine glaciers and ice caps under the assumption that the ice deforms according to Glen's constitutive relationship. Using a simple, yet wide, class of models for the accumulation rate of ice and Chebysev's theorem on the integration of binomial differentials, many new exact solutions of the Vialov equations are obtained in terms of elementary functions.

Equations

$${\dot{ϵ}}_{ij}=\mathcal{A}{\sigma }_{\text{eff}}^{n-1}{s}_{ij}$$

$${\sigma }_{\text{eff}}=\sqrt{\frac{1}{2}\text{mbox}Tr\left({\hat{s}}^2\right)}$$

$$xc(x)=\frac{2\mathcal{A}}{n+2}{\left(\rho gh\left|\frac{dh}{dx}\right|\right)}^n{h}^2,$$

$$h(x)={\left\{\mp \frac{2\left(n+1\right)}{n\rho g}{\left(\frac{n+2}{2\mathcal{A}}\right)}^{1/n}\int dx{\left[xc(x)\right]}^{1/n}\right\}}^{\frac{n}{2(n+1)}}\equiv A{\left[V(x)\right]}^{\frac{n}{2(n+1)}},$$

$$A\equiv {\left[\frac{2\left(n+1\right)}{n\rho g}{\left(\frac{n+2}{2\mathcal{A}}\right)}^{1/n}\right]}^{\frac{n}{2(n+1)}},$$

$$V(x)\equiv \int dx{\left[xc(x)\right]}^{1/n}.$$

$$h(x)$$

$$H$$

$$h(x)=h_0{x}^{\frac{n+m+1}{2(n+1)}}$$

$$c(x)=a+b{x}^r,$$

$$I\left(x;p,q,r\right)=\int dx{x}^p{\left(a+b{x}^r\right)}^q,p,q,r\in \mathbb{Q},r\ne 0$$

$$\int dx{x}^{1/3}{\left(a+b{x}^r\right)}^{1/3}=\frac{3x^{4/3}}{4\left(4+r\right){\left(a+b{x}^r\right)}^{2/3}}$$

$$\cdot \left[ar{\left(\frac{b{x}^r}{a}+1\right)}^{2/3}{}_2{F}_1\left(\frac{2}{3},\frac{4}{3r};1+\frac{4}{3r};\frac{-b{x}^r}{a}\right)+4\left(a+b{x}^r\right)\right],$$

$$\frac{p+1}{r},q,\frac{p+1}{r}+q$$

$$\frac{p+1}{r}$$

$$\frac{p+1}{r}+q$$

$$\frac{4}{3}\simeq 1.33,\frac{2}{3}\simeq 0.667,\frac{4}{9}\simeq 0.444,\frac{1}{3}\simeq 0.333,\frac{4}{15}\simeq 0.267,\frac{2}{9}\simeq 0.222,$$

$$\frac{4}{21}\simeq 0.190,\frac{1}{6}\simeq 0.167,\frac{4}{27}\simeq 0.148,\frac{2}{15}\simeq 0.133,...,0.$$

$$2,\frac{4}{5}=0.8,\frac{1}{2}=0.5,\frac{4}{11}\simeq 0.364,\frac{2}{7}\simeq 0.286,\frac{4}{17}\simeq 0.235,\frac{1}{5}=0.2,$$

$$\frac{4}{23}\simeq 0.174,\frac{4}{23}\simeq 0.174,\frac{2}{13}\simeq 0.154,...,0.$$

$$V(x)=\frac{n{a}^{1/n}}{n+1}{x}^{1+1/n}+D,$$

$$h(x)=A{\left[V(x)\right]}^{\frac{n}{2(n+1)}}=A{\left[\frac{n{a}^{1/n}}{n+1}{x}^{1+1/n}+D\right]}^{\frac{n}{2(n+1)}},$$

$$V(x)=\frac{3b^{1/3}}{4+r}{x}^{(4+r)/3}+D.$$

$$h(x)=A{\left[\frac{3b^{1/3}}{4+r}{x}^{(4+r)/3}+D\right]}^{3/8}.$$

$$r$$

$$V(x)$$

$$h(x)$$

$$h(x)=h_0{\left(a+b{x}^{4/3}\right)}^{3/8}$$

$$r$$

$$V(x)$$

$$h(x)$$

$$r$$

$$V(x)$$

$$h(x)$$

$$r$$

$$V(x)$$

Full text

Solving the Vialov equation of glaciology in terms of

elementary functions

Valerio Faraoni1

1 Department of Physics and Astronomy, Bishop’s University

2600 College Street, Sherbrooke, Québec Canada J1M 1Z7

vfaraoni@ubishops.ca

Very few exact solutions are known for the non-linear Vialov ordinary differential equation describing the longitudinal profiles of alpine glaciers and ice caps under the assumption that the ice deforms according to Glen’s constitutive relationship. Using a simple, yet wide, class of models for the accumulation rate of ice and Chebysev’s theorem on the integration of binomial differentials, many new exact solutions of the Vialov equations are obtained in terms of elementary functions.

PACS:92.40.vk, 92.40.Cy, 92.40.vv

1 Introduction

A part of the mathematical modelling of alpine glaciers and polar ice sheets and ice caps is the description of their longitudinal profiles, which is based on non-linear differential equations. The microphysics and the rheology of ice play a crucial role in determining the shape of glaciers. A good model for the response of glacier ice to stress is Glen’s law relating the strain rate tensor $\dot{\epsilon}_{ij}$ to stresses in the ice (Glen 1955)

[TABLE]

where $s_{ij}$ is the deviatoric stress tensor,

[TABLE]

is the effective stress, and ${\cal A}$ is a (temperature-dependent) constant (Paterson 1994; Cuffey and Paterson 2010; Hooke 2005; Greve and Blatter 2009). The value $n=3$ is adopted for glacier flow in most theoretical and modelling work.

Let $x$ be a coordinate along the glacier bed in the direction of the ice flow. Assuming incompressible and isotropic ice, steady state, a flat bed (this means that the bed is a plane which, in general, has non-zero slope), and Glen’s law, the longitudinal glacier profile (or local ice thickness) $h(x)$ obeys the Vialov ordinary differential equation (Vialov 1958; Paterson 1994; Cuffey and Paterson 2010; Hooke 2005; Greve and Blatter 2009)

[TABLE]

where $c(x)$ is the accumulation rate of ice, that is, the flux density of ice volume in the $z$-direction perpendicular to $x$, with the dimensions of a velocity. The absolute value in Eq. (1.3) is introduced when one looks for solutions in the finite interval $x\in\left[0,L\right]$. If $x=0$ and $x=L$ denote the glacier summit and terminus, respectively, then the local surface slope $dh/dx$ is negative and its absolute value must be taken. If instead $x=0$ denotes the glacier terminus while $x=L$ is the summit, it is $dh/dx>0$. For ice caps and ice sheets, once a solution for the longitudinal profile of half of a glacier is found in $\left[0,L\right]$, it is extended to the interval $\left[-L,L\right]$ (or to $\left[0,2L\right]$, respectively) by reflection about the vertical line $x=0$ (or $x=L$, respectively) passing through the summit. A consequence of this procedure is that the surface profile $h(x)$ of an ice cap or ice sheet is not differentiable at the summit, where the left and right derivatives of $h$ are finite and opposite and, usually, also at the terminus where the slope $dh/dx$ and the basal stress $\tau_{b}=-\rho gh\,dh/dx$ diverge (here $\rho$ is the ice density and $g$ is the acceleration of gravity). This is, however, common procedure in the literature (Paterson 1994; Cuffey and Paterson 2010; Hooke 2005; Greve and Blatter 2009). The non-linearity of the Vialov equation (1.3) is a direct consequence of the non-linearity of Glen’s law (1.1).

The formal solution of the Vialov equation (1.3) can be expressed as the integral

[TABLE]

where the upper sign applies if the summit is at $x=0$ and $dh/dx<0$, while the lower sign applies if $x=L$ is the summit,

[TABLE]

and the integral

[TABLE]

is determined up to an arbitrary integration constant $D$. A function $c(x)$ modelling the accumulation rate of ice must be prescribed. Even for simple choices of $c(x)$, the integral (1.6) can rarely be computed in terms of elementary functions, which has led to stagnation in the literature on this subject, but a few analytic solutions of the Vialov equation are known111Other analytic profiles (Nye 1951a; Nye 1951b; Faraoni and Vokey 2015) follow from the rather unrealistic assumption of perfecly plastic ice used in the early days of theoretical modelling and when the deformation of the ice is irrelevant. (Böðvardsson 1955; Vialov 1958; Weertman 1961; Paterson 1972; Bueler 2003; Bueler et al. 2005). Solutions in $\left[0,L\right]$ with $x=0$ and $x=L$ denoting the position of the glacier summit and terminus, respectively, include:

• •

$c=$ const., which yields the Vialov profile ((Vialov 1958), see also (Paterson 1994; Cuffey and Paterson 2010; Hooke 2005; Greve and Blatter 2009))

[TABLE]

• •

If $c(x)$ is chosen as a step function, the Weertman-Paterson profile is obtained by matching two Vialov solutions (Weertman 1961; Paterson 1972).

Assuming instead $x=0$ at the glacier terminus and $x=L$ at the summit, one obtains the following solutions.

• •

The model $c(x)=c_{m}x^{m}$ is used in the literature, with the value $m=0$ believed to be appropriate for ice caps and $m=2$ for alpine glaciers. In scaling theory, according to the Buckingham Pi theorem (Buckingham 1914), the mass balance rate is supposed to scale as $l^{m}$, while the characteristic thickness $h$ of a glacier or ice cap is assumed to scale with its characteristic length $l$ as $h\sim l^{s}$. The exponents $s=\frac{m+n+1}{2(n+1)}$ and $s=\frac{m+1}{n+2}$ are predicted by scaling theory for ice caps and for alpine glaciers, respectively (Bahr et al. 2015). The power law

[TABLE]

(with $h_{0}$ a constant) solves the Vialov equation (1.3) with $c(x)=c_{m}x^{m}$. The exponent $\frac{n+m+1}{2(n+1)}$ was deduced in scaling theory (Bahr et al. 2015; Faraoni 2016). This solution includes the case $c=$ const. and also the profile $h(x)=h_{0}\,\sqrt{x}$ which reproduces the parabolic profile first obtained by Nye under the simplifying assumption of perfectly plastic ice (Nye 1951a; Nye 1951b), which is very different from the more realistic Glen law (1.1) but is obtained as the limit $n\rightarrow+\infty$ of Eq.  (1.9). As shown in Sec. 2.1, the profile $h(x)=h_{0}\,\sqrt{x}$ is not restricted to the unrealistic assumption of perfectly plastic ice but is also a solution of the Vialov equation following from the realistic Glen law. This fact is significant because this profile is currently used in a number of applications (e.g., (Benn and Hulton 2010; Ng et al. 2010)) and is appropriate when the internal deformation of the ice is irrelevant.

In Sec. 2 a simple, yet broad, model of the function $c(x)$ describing the ice accumulation rate is postulated and the Chebysev theorem on the integration of differential binomials is applied to the search of exact solutions of the Vialov equation in terms of elementary functions, in the form (1.4). Infinitely many new solutions in terms of elementary functions can be obtained, some of which are reported in appendix A, while known solutions are re-derived. Sec. 3 contains a discussion of these solutions and of the method employed.

2 Chebysev theorem and Vialov equation

Let us return to Eq. (1.3) and let us search for solutions of the form (1.4) when the integral (1.6) can be expressed in terms of elementary functions. A wide class of reasonable models for the accumulation rate function is the choice

[TABLE]

where $a,b$, and $r$ are constants and where $r$ is chosen to be rational for reasons explained below. Special cases include:

$a=0$, $b>0$, $r>0$. In this case $x=0$ is the location of the glacier terminus corresponding to zero accumulation rate, while the glacier summit is at $x=L$, where the accumulation rate of ice assumes its largest value $c_{\text{max}}$. Then it follows that $b=c_{\text{max}}/L^{r}$. The choice $r=2$ is appropriate to describe alpine glaciers (Bahr et al. 2015). Although it is not done in the literature, a better model would assume $a<0$ to describe ablation at the glacier terminus. 2. 2.

$c(x)=a-|b|x^{r}$ with $a>0,b<0$, and $r>0$. In this case it is appropriate to locate the summit at $x=0$ and the terminus at $x=L$, with $c(x)$ a decreasing function of $x$ in $\left[0,L\right]$ vanishing at $x=L$ and with the constants assuming the values $a=c_{\text{max}}$, $b=-c_{\text{max}}/L^{r}$. An alternative choice consists of having $c(L)<0$ in order to describe ablation at the terminus.

With the choice $r\in\mathbb{Q}$, the integral $V(x)$ falls into the category

[TABLE]

(if $r=0$ the integral is trivial). In practice, for glacier flow it is $p=q=1/n=1/3\in\mathbb{Q}$. The integral (2.11) can be expressed in terms of an hypergeometric function,

[TABLE]

but this representation is of little use for practical purposes, for example when, in statistics, one needs a simple model of longitudinal glacier profile $h(x)$ to fit a large number of glaciers. For numerical studies of a single glacier, it is convenient to integrate numerically Eq. (1.3) but for other problems a simple analytic formula for $h(x)$ is required. A necessary and sufficient condition for the integral (2.11) to be expressed in terms of elementary functions is the

Chebysev theorem (Chebysev 1853; Marchisotto and Zakeri 1994):

the integral (2.11) admits a representation in terms of elementary functions if and only if at least one of

[TABLE]

is an integer.

Since $n=3$, it is $p=q=1/3\in\mathbb{Q}$, and $r$ in Eq. (2.10) is chosen to be rational (in the glaciological literature $r$ is usually the integer 0 or 2). Atmospheric models which could provide hints to fix the function $c(x)$ are not currently coded to have the ability to discriminate between a real number $r$ and a rational approximation of it. One then has

[TABLE]

Given the freedom in the choice of the parameters $a,b$, and $r$ of the model (2.10), one requires that $r\in\mathbb{Q}$ and searches for values of $r$ such that $(p+1)/r$ or $q+(p+1)/r$ are integers.

• •

By imposing that $\frac{4}{3r}\equiv m_{0}\in\mathbb{Z}$, one obtains $r\equiv\frac{4}{3m_{0}}$,  $m_{0}=1,2,3,\,...\,,+\infty$. This choice produces the sequence of values of $r$

[TABLE]

• •

Imposing $\frac{4+r}{3r}\equiv m_{0}\in\mathbb{Z}$ gives $r=\frac{4}{3m_{0}-1}$,  $m_{0}=1,2,3,\,...\,,+\infty$ and the sequence of values of $r$

[TABLE]

Not all these values of $r$ are appropriate from the glaciological point of view to describe the ice accumulation rate (2.10). However, the values 0 and 2 universally used in the literature, and many values of potential interest lying between these two extremes, are reproduced. The value $r=0$ is usually suggested for ice caps and ice sheets while the value $r=2$ is suggested for alpine glaciers (Böðvardsson 1955; Vialov 1958; Weertman 1961; Paterson 1972; Bueler 2003; Bueler et al. 2005). The representation of the integral $V(x)$ in terms of elementary functions falling into the range covered by the Chebysev theorem include the following special cases.

2.1 Choice $c(x)=$ constant

The choice $c(x)=$ constant can be obtained by setting $b=0$ (in which case $r$ drops out of the discussion) or when $b\neq 0$ with $r=0$ (in which case the Chebysev theorem as stated does not apply). In both cases the integration is trivial and, in the first case, one obtains

[TABLE]

where $D$ is an integration constant, and the longitudinal glacier profile

[TABLE]

which reproduces the Vialov profile (1.7), (1.8) always associated with the choice $c=$ const. in the glaciological literature (Paterson 1994; Cuffey and Paterson 2010; Hooke 2005; Greve and Blatter 2009). Setting $D=0$ yields the parabolic profile $h(x)=h_{0}\,\sqrt{x}$ irrespective of the value of $n$.

2.2 Choice $c(x)=b\,x^{r}$

In this case, with $a=0$ and $b>0$, and without choosing a specific value of $r$, one obtains the integral222Strictly speaking, in this degenerate case there is no need to assume that $r\in\mathbb{Q}$ and use the Chebysev theorem. In fact, these ingredients were not assumed in the recent work (Bahr et al. 2015) deriving this power law solution.

[TABLE]

The corresponding longitudinal glacier profile is

[TABLE]

Setting $r=0,D\neq 0$ reproduces $c=$ const. and gives the Vialov profile (1.7) and (1.8). Setting instead the integration constant $D$ to zero yields $h(x)=h_{0}\,x^{(4+r)/8}$. As already noted, the value $r=2$ is appropriate to describe alpine glaciers (Bahr et al. 2015; Faraoni 2016) and gives $h(x)\propto x^{3/4}$. Setting instead $r=0$, which is appropriate for ice caps, yields the well known profile $h(x)\propto\sqrt{x}$ (Paterson 1994; Bahr et al. 2015).

The choice $c(x)=b\,x^{r}$, usually written as $c(x)=c_{m}x^{m}$, reproduces the power law solution $h(x)\propto x^{\frac{n+m+1}{2(n+1)}}$ of (Bahr et al. 2015; Faraoni 2016). In fact, setting $r=m$ and $n=3$ yields $h\sim x^{(4+r)/8}$.

As $x$ becomes large the highest order term is dominant and, in all of these solutions, the profile then approaches $h(x)\sim\sqrt{x}$.

Other examples of plausible models of the accumulation rate $c(x)=a+b\,x^{r}$ leading to representations of the integral (1.6) in terms of elementary functions and to relatively simple exact profiles are reported in appendix A.

3 Discussion

Analytic expressions describing longitudinal glacier profiles are needed in several problems of glaciology (e.g., (Thorp 1991; Ng et al. 2010; Benn and Hulton 2010)). However, under the realistic and well tested assumption that glacier ice deforms according to Glen’s constitutive relationship (1.1), the Vialov ordinary differential equation (1.3) ruling these longitudinal glacier profiles is non-linear and obtaining analytic solutions in closed form in terms of elementary functions is difficult. Only a few exact solutions are known in the literature (Böðvardsson 1955; Vialov 1958; Weertman 1961; Paterson 1972; Bueler 2003; Bueler et al. 2005). By assuming a simple, yet general, model for the accumulation rate of ice appearing in Eq. (1.3), the Chebysev theorem provides a necessary and sufficient condition for the integral (1.6) expressing a formal solution of the Vialov equation to be represented in terms of elementary functions. The solutions provided by the Chebysev theorem include the known solutions, with the exception of the Böðvardsson, Vialov, and Bueler profiles (Böðvardsson 1955; Vialov 1958; Bueler 2003; Bueler et al. 2005).

The initial condition $h=0$ of the Vialov equation (1.3) is imposed at the glacier terminus ($x=0$ or $x=L$, depending on the geometry adopted, which determines also the sign of $dh/dx$), which is a singular point of the equation corresponding to divergent surface slope $dh/dx$. In this situation, the usual uniqueness theorems for ordinary differential equations (e.g., (Brauer and Noel 1986)) do not hold and this is the reason why one can find multiple solutions of the Vialov equation, and why the solutions obtained by using the Chebysev theorem do not always generate the well known Vialov (1958) profile, and do not reproduce other profiles (Böðvardsson 1955; Bueler 2003; Bueler et al. 2005).

An infinite number of solutions in terms of elementary functions is guaranteed by the Chebysev theorem, corresponding to rational values of the constant $r$, and they can be found easily with computer algebra. The current models for the ice accumulation rate $c(x)$ are very unsophisticated ($c=$ const. being perhaps the most popular choice) and the 3-parameter choice $c(x)=a+b\,x^{r}$, $r\in\mathbb{Q}$ allows freedom to extend these models. Of course, other functional choices may be appropriate to model the accumulation rate $c(x)$ and, at the same time, provide analytic profiles $h(x)$. However, exact solutions of the Vialov equation (1.3) in simple form have been hard to find and sometimes they correspond to unintuitive choices of $c(x)$ which make the corresponding analytic profile $h(x)$ more of a toy model achieving one desired physical property than a realistic description of the shape of alpine glaciers and ice caps. This is the case of the Bueler profile (Bueler 2003; Bueler et al. 2005; Greve and Blatter 2009), which exhibits a finite basal stress $\tau_{b}=-\rho gh\,dh/dx$ at the glacier terminus, contrary to the Vialov and other profiles. The old Chebysev (1853) theorem extends the scope of existing analyses. The values of the parameters $a,b$, and $r$ in Eq. (2.10) appropriate to particular geographic locations have to be determined by data-fitting and are expected to be different for different situations (alpine glaciers, polar ice caps, cirque glaciers, etc.).

Acknowledgments

This work is supported by Bishop’s University.

Appendices

A Exact solutions of the Vialov equation for

some rational values of $r$

In the range of parameters $\left(p,q,r\right)$ in which the Chebysev theorem is satisfied, computer algebra easily provides the integral (1.6) for the choice (2.10) of $c(x)$. Here some of these integrals and the corresponding longitudinal glacier profiles are reported for various values of the parameter $r$ listed in Eq. (2.10), which correspond to values of $m_{0}$ reported in Eq. (2.15).

[TABLE]

where $D$ is, as usual, an integration constant. For $D=0$ one obtains

[TABLE]

and, if also $a=0$, one obtains again $h(x)=h_{0}\sqrt{x}$. Another possibility is

[TABLE]

where $D_{0}$ is another constant. Other possibilities are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]