### Three Linearized Glacier Models; Asymptotic Behaviour In The Limit Of Small Sliding At The Terminus

In the following we will use perturbation methods to analyze three different linear steady–state models of the form of Equation (7), with the time derivative omitted and with a uniform mass–balance perturbation *b*_{x}(x) = 1. Of special interest is the volume perturbation

which determines the volume time–scale T_{V} (Equation (17a)) and therefore the adjustment time T_{M} of the model glacier (Equation (21)). The principal question concerns the asymptotic behaviour of *V*_{1}
as sliding at the terminus goes to zero. Sliding at the terminus had little effect on the response of the numerical glacier models in section V. Equation (la), on the other hand, seems to predict that sliding at the terminus plays a key role in determining the response. Our goal is to explain in analytical terms the perhaps unexpected result of the numerical glacier models in section V.

Perturbation models (Cole, 1968; Kevorkian and Cole, 1981) are often used to find approximate solutions to problems involving a small parameter 6. The solution of the problem is expanded as an asymptotic expansion in 6 (the “outer expansion”). Often this expansion fails to satisfy a boundary or a regularity condition. In that case, a different asymptotic expansion (the “inner expansion”) has to be used in a so–called boundary layer. In the following we will make frequent use of the “big *O* notation”, i.e. *f(x) = 0(g(x))* as *O(g(x))* as x → 0 means that the limit of *(f/g)* is finite as *x → 0*.

1. **Nye’s (1963a) £,δ model**

The model is described by the (non–dimensional) kinematic wave velocity C_{0} = x(l – x) and diffusion coefficient *D*_{0} = Ex^{2}(1 − δ − x)

107] with 0 > δ >> 1 and 0 ≼ *x* ≼ *l*
_{0} = (1-δ) (see section VII for discussion). Integration of Equation (7) with respect to *x* (omitting the time derivative and putting b_{1}(X) = 1) leads to the differential equation

for the steady–state change in thickness h_{1}(X). By putting *y = 1 − δ − x*, this equation is more conveniently expressed as

There is no boundary condition, but the solution is required to be regular at the singular point y = 0 of the differential equation (Nye, 1963a, p. 437–40). Obviously, *h*
_{1}(y = 0) = 1/δ. This suggests the expansion *h*
_{1} = . Substitution into Equation (Alb) and collection of terms with equal powers of & leads to differential equations for each of . The equation for is

Then

where *k* is a constant of integration. This solution is regular for all *y*, 0 ≼ *x* ≼ *l*
_{0} = (1 − δ).

This equation has a regular solution at *у =* 0 if *K* = 1.

Consequently, the constant of integration in *h[°>* can be determined from the equation for . Similarly, a constant of integration in can be determined by requiring that be regular and in this way the expansion for h_{1} can be built up to any order without the introduction of a boundary layer. The lowest order, however, is sufficient for our purposes. It is given by

Thus

This approximation is very good for small 8. For 6 as small as 0.006, as used by Nye (1963a), and *E* in the range 0.2 ≼ *E* ≼ 2.0, plots of h_{1}(X) from Equation (A2) can hardly be distinguished by eye from plots of A_{t}(x) based on exact solutions of Equation (Ala). For very small £, the approximation fails since it is implicitly assumed in the preceding analysis that *E* = 0(1). For *E =* Ο(δ), a different scaling of the variables is necessary to produce a consistent approximation.

The volume perturbation *V*_{1}
deduced from Equation (A2) is

The factor *f* = *VJUnh*_{1}(I_{0})) in Equations (14) is evidently / = 1/(1 + *)/E) +* O(S) if *E* = 0(1).

The shape of A_{1} is essentially unaffected as δ → 0 and A_{1} is inversely proportional to 5 over almost the entire length of the glacier. *V*_{1}
and therefore

T_{v} and т_{м} go to infinity as S – 0. This means that the amount of sliding at the terminus determines the response of the whole glacier. It is difficult to imagine a physical process that can bring about this result. Therefore, we will now analyze a slightly different model in order to find out whether the 1/6 dependency of A_{1} and *V*_{1}
only applies to this particular model or whether it is likely to be true in general for models of this kind.

2 A slightly modified £,6 model

The diffusion coefficient *D*_{0}
must go to zero at the terminus (see section III for discussion). Thus, the steady–state form of Equation (7) will in general have a singular point at the terminus. This is indeed the case for the £,6 model which has a singular point at *χ =* I – 6. We may expect the steady–state solution to be sensitive to the form of O_{0}, or more importantly to the form of the ratio *D*_{0}/C_{0}, as the singular point at the terminus is approached. Assuming *q = q(h,a)* and using the cyclic rule of partial derivatives, we can write −(∂*h*/∂*α*
_{
q
}. The partial derivative (∂*h*/∂*α*
_{
q
} expresses how small changes in A and α at a fixed location must be related to each other for the flux to remain constant at that location.

In the interior of the glacier, outside a zone near the terminus where special considerations must apply, we expect a flux relationship of the form of Equation (4b) to be valid. A flux relationship of this form predicts that

*D*_{0}ZC_{0} = –(QhZQa)_{cl} =

(rZs)h_{0}Za_{0}. Since A_{0} decreases towards the terminus and *O*_{0}
probably increases, *D*_{0}
must decrease faster than *C*_{0}
as the terminus zone is approached. This may be crudely quantified for a glacier on a flat bed by assuming that A_{0} decreases as some power of the distance from the terminus *U — X* (a square–root shape is commonly assumed to be realistic). Then *D*_{0}ZC_{0} = (r/.s)h_{0}/ (–dA_{0}/dx) = *O(l*_{0} – x). We see that in this case the ratio *D*_{0}ZC_{0}
will decrease proportional to *I*_{0} – x. The distributions of *D*_{0}
and C_{0} for the ε model presented in this paper (Fig. 8) clearly show this behavior. So do published calculations of *D*_{0}
and C_{0} based on real glacier profiles (e.g. Nye, 1965b, fig. 5; Reeh, 1983, p. 38–53; however, figure 7 in Nye (1963b) does not show this as well).

The above argument does not apply to the terminus zone, where a flux relationship of the form of Equation (4b) breaks down. However, it turns out to be crucial that the form of C_{0} and *D*_{0}
is realistically modeled outside the terminus zone where we assume a flux relationship of this form to be valid.

Both D_{0} and C_{0} in Equation (Ala) decrease linearly as *x → l*_{0}
. Based on the above discussion, this is not realistic and we replace *D*_{0} = Ex^{2}(1 − δ − x) in Equation (Ala) with *D*_{0} = Ex(1 − δ − x)^{2}
in order to investigate the effect of the shape of D_{0} near the terminus on the solution. In addition to replacing (1 − δ − x) with (1 − δ − x)^{2} to produce the required ratio of *D*_{0}ZC_{0}
near the terminus, we have also replaced *x*^{2}
with *Y* This is done in order to simplify the analysis and does not change the nature of the problem. With this modification, Equation (Ala) becomes

or defining as before y = 1 − δ − x

]
One quickly finds that the method used for Equation (Alb) does not work here. It is not possible to derive an outer expansion which is regular at *у* = 0. This indicates the occurrence of a boundary layer at the terminus.

A differential equation for the boundary layer is deduced by re–scaling both the dependent and the independent variables. In this case, the correct scaling is given by *y** = *y*/δ, , which leads to

We see that the small parameter 6 has been scaled out of the problem. This is not unusual for perturbation problems and means that at least in the boundary layer the full equation has to be solved. The unique regular solution of Equation (A5) is readily found to be

From this solution, both the nature of the boundary layer and the asymptotic behavior of the solution in the interior of the glacier can be derived.

The boundary layer is of width 0(6) and within it the solution is

0(1/8). Its contribution to *V*_{1}
is therefore O(l).

Away from the boundary layer, the asymptotic behavior of A_{1}O) as 6 – O depends on the behavior of
*y*/δ → ∞ for a fixed value of *y.* Depending on the numerical value of E, we find that as *y*
^{* → ∞}

where Γ is the gamma function, .

Then is given by

as δ → 0 for a fixed *y,* in the interior of the glacier outside the boundary layer at the terminus.

The contribution to *V–* from Zi_{1} outside the boundary layer can be found by integrating A_{1}O–) in the interior of the glacier (Equation (A6)). This is the most important contribution to *V*_{1}
and we can derive the following order–of–magnitude estimate

The factor ƒ = *V*_{1}/(l_{0}h_{1}(l_{0})) in Equations (14) goes to zero as 6 –* 0. This arises because the thickness change is only partially spread up–glacier from the terminus where h_{1}(Z_{0}) = 1/8. For *E* in the range 0 *E* 1, Zi_{1} is even *independent* of 8 outside the boundary layer at the terminus. This is similar to what we found for the Ε model in section V (Fig. 9).

Equation (A7) predicts that *V*_{1}
and therefore

T_{V} and T_{M} go to infinity as δ → 0 as was the case for the £,8 model. The dependency of *It*_{1}
and *V*_{1}
on 8 as 8 – 0 is, however, quite different from the £,8 model and this indicates that analytical models of this kind are highly sensitive to the details of the distributions of D_{0} and C_{0}, especially near the terminus. Therefore, it must be crucial that the distributions of D_{0} and C_{0} are realistic in models of this kind. Otherwise, almost any result may be obtained. In particular, the prediction of the £,5 model that the volume change *V*_{1}
is inversely proportional to sliding at the terminus, cannot be regarded as a valid theoretical prediction for real glaciers, since other seemingly more realistic distributions of D_{0} and C_{0} yield quite different predictions.

3. **A model derived from a steady–state profile**

In numerical models of glacier flow, the terminus has traditionally been treated in a very crude way (e.g. Budd and Jenssen, 1975; Bindschadler, 1982; Waddington, unpublished). The grid size of the models is usually so coarse that detailed modeling of the terminus is out of the question. Instead, it is assumed that the terminus plays no active role in the dynamics of the glacier as a whole. The position of the terminus is determined from simple continuity requirements, i.e. the shape of the terminus is assumed constant and the rate of change of its position is determined by the ice flow past the last grid node up–stream from the terminus. These models cannot produce any effect of the terminus on the glacier dynamics, since the expected passive behavior of the terminus is built into the models as an assumption.

The assumption of a passive terminus is a physically appealing one. One would expect realistic analytical models to have a similar property, even if it is not included at the outset as an assumption. In this view, the crucial part of the model is the correct specification of the model parameters (i.e. D_{0} and C_{0}) in the interior of the glacier, where the proposed flux relationship is assumed valid. The region very close to the terminus, where the flux relationship breaks down, must be included in the model in some way, but from physical grounds one would expect this region to be isolated from the rest of the model. Mathematically, this could for example arise if an isolated boundary layer is located at the terminus. Needless to say, the behavior of *V*_{1}
in the limit of zero sliding at the terminus must be non–singular. Otherwise, the terminus region determines the response of the whole glacier.

We will now derive a simple glacier model based on the flux relationship expressed by Equation (4b) and investigate whether consistent specification of D_{0} and *C*_{0}
leads to the results discussed above. Both D_{0} and C_{0} are defined as partial derivatives of the same quantity (Equations (6)). Therefore, they cannot be specified independently, but must be based on a specific flux relationship, datum mass balance, and the corresponding steady–state profile. We will investigate a two–dimensional ice sheet on a flat bed and take *K* – 1, г = 3, and *s* = 5 in the flux relationship (Equation (4b)). h_{0} in the ablation area is then given by

when A_{0} is uniform A_{0}(x) = –1 in the ablation area (Paterson, 1980). The corresponding D_{0} and C_{0} (Equations (6)) are

where we have added the term 6 in the expression for C_{0}in order to provide a non–zero C_{0} at the terminus χ = Z_{0}–

We see that D_{0} decreases faster than C_{0} as *χ* – /_{0} as expected.

Adding 6 to C_{0} is admittedly an arbitrary way of specifying the terminus dynamics. As discussed above, it is the consistency of A_{0}, D_{0}, and C_{0} in the interior of the glacier which is of primary concern. Given consistency, the details of the terminus dynamics should not be important for the dynamics of the glacier as a whole. Therefore, this *ad hoc* way of specifying the terminus dynamics is adequate for our purpose.

In the accumulation area, the expression for A_{0} is different from Equation (A8) and consequently the corresponding expressions for D_{0} and C_{0} are different from Equations (A9). Arithmetic simplicity makes it worthwhile to analyze the ablation area without considering the accumulation area. This can be done if it is assumed that that mass–balance perturbation A_{1}(Jt) is confined to the ablation area, since then the thickness change A_{1}(Ac) in the ablation area will be independent of the accumulation area. Assuming an ablation area extending from x = 0 to x = 1 with a uniform mass–balance perturbation b_{1}(X) = 1, this model is expressed by the following equation

with D_{0} and C_{0} expressed by Equations (A9). The accumulation area is located to the left of χ = 0 and is not considered.

After a change of variables to y = l_{0} − x = 1 − x for convenience, the outer expansion to lowest order satisfies the equation

The solution of this equation is

This solution is singular at *у =* 0 for any value of the integration constant *K.* Therefore, a boundary layer at *у* = 0 is needed. The integration constant *K* must be determined by matching to the boundary layer (Kevorkian and Cole, 1981, p. 24).

The boundary–layer equation is deduced by re–scaling *у* and A_{1}. The correct scaling is . Using Equations (A9), Equation (AlO) transforms to

and the equation determining the lowest–order term A_{1}’^{0}’, in the inner expansion is

This equation has the unique non–singular solution.

In very simple terms, matching of Equations (All) and (A 12) is done by requiring that, when expressed in the same variables, the equations agree to highest order “outsidejust outside” the boundary layer. When *y*
^{* → ∞} (by letting 8–0 for a fixed *y),* the lowest–order term in the boundary layer (Equation (A 12)) becomes

Expressed in outer variables, this expression for the boundary layer becomes

as δ → 0 for a fixed *у.* This means that we must have *K* = O if the boundary layer (Equation (A 12)) is to merge smoothly with the outer solution (Equation (All)). Thus,

The boundary layer (Equation (A12)) is of width *O(δ*^{2}) and within it the solution is *O(1/δ)*. Its contribution to *V*_{1}
is therefore 0(6) and goes to zero as 6–0. The outer solution (Equation (A 13)) can be integrated without difficulty to find *V*
_{1}

which is finite in the limit 6 = 0. the factor f =

V^{1}/(l^{0}h^{1}(l^{0}) in Equations (14) is *O(δ)* as δ → 0.

Equation (A13) shows that, away from the terminus, A_{1} is determined by the dynamics in the interior of the glacier and independent of the amount of sliding at the terminus. Equation (A 14) shows that *V*_{1}
and therefore r_{v} and

*T*_{M}
are also determined by the dynamics in the interior of the glacier, and finite in the limit 8–0. This is what we expected to find for an analytical model with consistent specification of A_{0}, *D*_{0}, and

C_{0}. We expect other models with carefully specified

*D*_{0}
and C_{n} (for example, based on a flux relationship expressed by Equation (4b) with other values of *r* and *s* than used here, or models with different formulation of C_{0} near the terminus), to exhibit the same properties as the model presented here. In section VI, we found that geometrical similarity of steady–state glacier profiles leads to the conclusion that *f* is proportional to the sliding rate at the terminus. This is the case here since *f* = 0(8). Therefore, the ratio *ƒ/u(l) = ƒ/δ* in Equation (la) is approximately constant for small 8, and Equation (la) is consistent with Equation (lb). This model agrees with the numerical results of the paper, which are derived from a more complete formulation, which does not permit an analytical analysis.

We conclude that consistent specification of A_{0}, O_{0}, and C_{0}, based on a specific flux relationship and datum mass balance, is essential for successful analytical models.