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 bx(x) = 1. Of special interest is the volume perturbation
which determines the volume time–scale TV (Equation (17a)) and therefore the adjustment time TM of the model glacier (Equation (21)). The principal question concerns the asymptotic behaviour of V1
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 C0 = x(l – x) and diffusion coefficient D0 = Ex2(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 b1(X) = 1) leads to the differential equation
for the steady–state change in thickness h1(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
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 h1 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
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 h1(X) from Equation (A2) can hardly be distinguished by eye from plots of At(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 V1
deduced from Equation (A2) is
The factor f = VJUnh1(I0)) in Equations (14) is evidently / = 1/(1 + )/E) + O(S) if E = 0(1).
The shape of A1 is essentially unaffected as δ → 0 and A1 is inversely proportional to 5 over almost the entire length of the glacier. V1
Tv 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 A1 and V1
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 D0
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 O0, or more importantly to the form of the ratio D0/C0, 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/∂α
. The partial derivative (∂h/∂α
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
D0ZC0 = –(QhZQa)cl =
(rZs)h0Za0. Since A0 decreases towards the terminus and O0
probably increases, D0
must decrease faster than C0
as the terminus zone is approached. This may be crudely quantified for a glacier on a flat bed by assuming that A0 decreases as some power of the distance from the terminus U — X (a square–root shape is commonly assumed to be realistic). Then D0ZC0 = (r/.s)h0/ (–dA0/dx) = O(l0 – x). We see that in this case the ratio D0ZC0
will decrease proportional to I0 – x. The distributions of D0
and C0 for the ε model presented in this paper (Fig. 8) clearly show this behavior. So do published calculations of D0
and C0 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 C0 and D0
is realistically modeled outside the terminus zone where we assume a flux relationship of this form to be valid.
Both D0 and C0 in Equation (Ala) decrease linearly as x → l0
. Based on the above discussion, this is not realistic and we replace D0 = Ex2(1 − δ − x) in Equation (Ala) with D0 = Ex(1 − δ − x)2
in order to investigate the effect of the shape of D0 near the terminus on the solution. In addition to replacing (1 − δ − x) with (1 − δ − x)2 to produce the required ratio of D0ZC0
near the terminus, we have also replaced x2
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 V1
is therefore O(l).
Away from the boundary layer, the asymptotic behavior of A1O) 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 Zi1 outside the boundary layer can be found by integrating A1O–) in the interior of the glacier (Equation (A6)). This is the most important contribution to V1
and we can derive the following order–of–magnitude estimate
The factor ƒ = V1/(l0h1(l0)) 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 h1(Z0) = 1/8. For E in the range 0 E 1, Zi1 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 V1
TV and TM go to infinity as δ → 0 as was the case for the £,8 model. The dependency of It1
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 D0 and C0, especially near the terminus. Therefore, it must be crucial that the distributions of D0 and C0 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 V1
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 D0 and C0 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. D0 and C0) 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 V1
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 D0 and C0
leads to the results discussed above. Both D0 and C0 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)). h0 in the ablation area is then given by
when A0 is uniform A0(x) = –1 in the ablation area (Paterson, 1980). The corresponding D0 and C0 (Equations (6)) are
where we have added the term 6 in the expression for C0in order to provide a non–zero C0 at the terminus χ = Z0–
We see that D0 decreases faster than C0 as χ – /0 as expected.
Adding 6 to C0 is admittedly an arbitrary way of specifying the terminus dynamics. As discussed above, it is the consistency of A0, D0, and C0 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 A0 is different from Equation (A8) and consequently the corresponding expressions for D0 and C0 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 A1(Jt) is confined to the ablation area, since then the thickness change A1(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 b1(X) = 1, this model is expressed by the following equation
with D0 and C0 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 = l0 − 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 A1. The correct scaling is . Using Equations (A9), Equation (AlO) transforms to
and the equation determining the lowest–order term A1’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 V1
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
which is finite in the limit 6 = 0. the factor f =
V1/(l0h1(l0) in Equations (14) is O(δ) as δ → 0.
Equation (A13) shows that, away from the terminus, A1 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 V1
and therefore rv and
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 A0, D0, and
C0. We expect other models with carefully specified
and Cn (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 C0 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 A0, O0, and C0, based on a specific flux relationship and datum mass balance, is essential for successful analytical models.