Introduction

Columbia Glacier is a large (1100 km²), grounded, iceberg-calving glacier near Valdez, Alaska (Fig. 1). From its first recorded sighting in 1794 (Reference Vancouver, Robinson, Paternoster and Pall MallVancouver, 1798) until 1978, its terminus pushed against a moraine shoal submerged about 20 m below sea-level. It was conjectured (Reference PostPost, 1975) that a tide-water glacier’s retreat would become rapid and irreversible should it disconnect from its shoal. The U.S. Geological Survey began an intensive research program in 1977 to determine whether this would be the case with Columbia Glacier. On the basis of results of three numerical models, such retreat was predicted (Reference MeierMeier and others, 1980). By late 1986, the terminus had retreated 3 km behind the shoal and the ice velocity was twice as great as in 1978 (Reference KrimmelKrimmel, 1987).

Although the initial concern of the research program was the threat that icebergs might pose in the shipping lanes to the Valdez oil port, the data collected were recognized as valuable to glaciological research in general. The chief data source is vertical aerial photography, which has been obtained about five times a year since 1976 over the lowest 18 km of the glacier. The extreme surface roughness there, making even helicopter-assisted surface work difficult and dangerous, is highly amenable to photo-grammetric analysis (Reference Meier, Rasmussen, Krimmel, Olsen and FrankMeier and others, 1985) without the aid of specially placed markers. Airborne radio echo-sounding was used over the lower half of this region in 1978 (Reference Brown, Rasmussen and MeierBrown and others, 1986) to determine the bed topography, the first time ever for a temperate glacier. At several stakes placed in the ice, mass balance was measured on the approximate dates of three of the flights (Reference Mayo, Trabant, March and HaeberliMayo and others, 1979).

The wealth of data collected will benefit research into the poorly understood mechanisms of calving (Reference Meier and PostBrown and others, 1982) and of glacier flow (Reference HookeHooke, 1981), including sliding motion (Reference Alean, Braun, Iken, Schram and ZwostaAlean and others, 1986). Much uncertainty afflicts these endeavors, in two forms: (1) the diversity of proposed functional forms, and (2) the range of numerical values considered for the coefficients appearing in any particular form. Toward these ends, data reduction was undertaken to interpolate values to the nodes of a square grid to ensure that the interpolated values were internally consistent with respect to the principle of conservation of mass. If a glacier-flow model drew its initial conditions from a data set that was not internally consistent, the model would begin by artificially redistributing the glacier mass, not as a realistic projection of the flow but simply to reconcile the inconsistencies.

Although most flow models have heretofore introduced partly theoretical, partly empirical shape factors to parameterize the transverse variations of geometry and flow — with only variations along a single longitudinal axis being treated explicitly — the two-dimensionality of the Columbia Glacier data was deliberately retained in creating the data set. This was done so that it could be used by models, such as ice-sheet models, that explicitly treat the transverse variations of flow. As little as possible was assumed about the flow law itself to enhance the data set’s value in investigating it. The work was motivated by the desire to contribute to glaciology a data set against which flow models and other rheological speculations could be tested, much as the data set of Reference Lettau and DavidsonLettau and Davidson (1957) is used in atmospheric boundary-layer research.

The bed topography had been sounded for only half the area covered by the aerial photography; for a variety of reasons, the resulting accuracy was low. The balance stakes were well distributed spatially but the frequency of their measurement was much lower than that of the aerial photography. Because of these deficiencies, and because of the abundance of photogrammetric data, a method was devised for creating the entire data set solely from the photogrammetric data. That is, aerial photography yielded not only the surface topography and velocity fields but also estimates of the bed topography and mass-balance distributions, making in all an internally consistent data set for each of 21 intervals between successive flights. As a result, observed changes in surface altitude have been separated into their two components without the aid of measurements made from the glacier surface: (1) the contribution from mass-balance processes, and (2) the dynamic thinning or thickening due to redistribution of mass by the ice flow.

Fig. 1. Index map and surface topography. The open arrows show the direction of flow, the main ice stream being indicated by the longer ones and by small dots at 2 km intervals from the head of the glacier. The surface topography was mapped from aerial photography taken 1 September 1981; the large dots are the 120 nodes of the square data grid. Standard matrix notation is followed for the node indices. Both the Universal Transverse Mercator (UTM) and local coordinate systems are shown.

The results, which obey a centered finite-difference approximation of the continuity equation, are tabulated on the nodes of a 762.5 m square grid. The surface altitudes for the dates of each of 22 aerial photography flights have been given in Reference Rasmussen and MeierRasmussen and Meier (1985, appendix B). The velocity components for the 21 intervals between the flights are given in Rasmussen (in press, appendix A); also included in that report are the bed altitudes, linear balance-versus-altitude functions for each interval, and all the tedious detail of creating the data set. This paper describes the method used to estimate the bed topography and mass balance from the aerial photography.

Aerial Photography and Photogrammetry

The coordinates of naturally occurring surface features, in most cases the intersections of crevasses, can be determined within about 3 m from the aerial photography. Because a feature was usually recognizable for several months, despite winter snow cover, its displacement vector could also be measured over the interval from one photography flight to the next. Photogrammetric analysis yielded the coordinates of about 150 surface points in the lower reach on the dates of each of 22 flights during the 4.25 year period from June 1977 to September 1981. There were slightly fewer displacement vectors than surface-altitude points over the intervals between flights, because an individual point could not always be located on the photographs from both flights spanning an interval.

Each displacement vector was converted to the average velocity over its time interval by dividing by the interval’s length, and the velocity was assigned to the mid-point of the trajectory. The intervals were so short and the strain-rates in the region were so low that the error in this approximation is negligible (Reference RasmussenRasmussen, 1983). The spatial distribution of the velocity is shown in component form in Figure 2 as averages over the entire 4.25 year period. The distributions for individual intervals are roughly similar, except that the contour values would all be scaled by a factor ranging, over the intervals, from about 0.7 to 1.4. The direction of the flow at any node was nearly constant over time. The temporal distribution of velocity is shown in Figure 3 for three locations in the region; winter velocities are about twice those of summer.

Fig. 2. Average velocity from 2 June 1977 until 1 September 1981. Both fields are hand contoured, at an interval of 100 m/year, from adjusted values at the grid nodes. Their row and column indices represent them as a 120 node section of the 71 by 63 grid covering the entire glacier. The u-component is in the x-direction, and the v-component in the y-direction.

The advantages of high spatial and temporal density of the photogrammetric data are offset, in part, by the irregular positioning of the natural features used. To permit readier use of the data, the surface altitudes for the date of each flight were interpolated to the nodes of a square grid (Reference Rasmussen and MeierRasmussen and Meier, 1985) by using the method of optimum interpolation (Reference GandinGandin, 1965). The surface topography at the end of the 4.25 year period is shown in Figure 1, and its variations are shown in Figure 4, in which the altitudes are represented as values averaged across the glacier width. The amplitude of the seasonal oscillation is about the same as the annual amount of the long-term lowering over the period, which ranges from about 10 m at the top of the region to about 20 m at the bottom, resulting in a slight steepening of surface slope. Waves of altitude change are not apparent, indicating synchronous change over the entire region.

The velocities derived directly from the displacement vectors were contoured by hand, component by component, at a scale of 1 : 50 000 and a contour interval of 50 m/a. The u(v)-component is measured along the x(y)-axis, which is positive to the east (north). Visual interpolation between the contours gave preliminary estimates of each component, on each of the 120 grid nodes (Fig. 1), for each of the 21 intervals.

The bed topography and mass-balance distributions are inferred from observed changes in the surface topography, which are assumed to be known exactly, and from the initial estimates of the velocity. Because of data error in the velocity estimates, no combination of bed topography and simple mass-balance distribution exists that will be consistent with those initial estimates. Therefore, the initial estimates must be adjusted so that mass-balance ditributions and a bed consistent with the adjusted velocities can be found. To maximize the influence of the original velocity data, the minimum adjustment is sought. The magnitude of the adjustment is not taken to be simply the sum of the unweighted vector differences — between the original and adjusted vectors — nor of the differences scaled by the original vector. Instead, the differences are scaled by independently estimated errors in the original velocity estimates.

For any velocity field, at any node, the standard errors ε^u and ε^v in the preliminary estimates were assumed to be proportional to the gradients |∇u| and |∇v| averaged over about 1 km² area centered on the node. This was not done on glaciological grounds but solely in recognition of the likely human error, both in creating the contour field and in interpolating between them to get the velocity estimates at the grid nodes. The constant of proportionality was an increasing function of the distances from the node to the three nearest displacement vectors. The 4 m standard error in determining a point’s position photogrammetrically was converted to velocity units and was arbitrarily assigned the same direction as the much larger error ascribed to creating and interpolating between the contours, to which it was added. As an illustration, the ε^v at the nodes in the lowest part of the region are shown in Figure 5, along with the v-contours and the locations of the displacement vectors from which the contouring was done, for a mid-1980 interval.

The root-mean-square average over all nodes, all flight intervals, was 27 m/a for ε^u and 67 m/a for ε^v. In accord with the velocity gradients (Fig. 2), the ε^u were roughly uniformly distributed spatially, and the ε^v were about half as big in the interior as on the margins. The combined root-mean-square of both, averaged over all 120 nodes, ranged over the 21 intervals from 36 to 60 m/a, with an average of 49 m/a.

Constrained Interpolation

The condition of internal consistency is expressed by the continuity equation in the form

in which ∇•Q is the horizontal divergence (in the x,y plane) of the flux vector; b is the mass-balance rate, and h is the rate of change of thickness, both measured vertically in ice- equivalent units. The flux vector is the vertical integral of the horizontal velocity vector V:

in which Z and B are the surface and bed altitudes, and the thickness h is their difference. Because the velocity profile is not observed, the factor γ is introduced to relate the flux to the surface velocity, and the product γh is considered for convenience as the characteristic thickness h. Although γ is a feature of the velocity profile, it is more convenient mathematically to associate it with h than with V.

Fig. 3. Down-glacier velocity component as a function of time at three indicated nodes. Velocity is determined as an average value (horizontal line segment) over each of the intervals between successive flights (node 68,23 solid; 60,23 dashed; 49,22 dotted). The hand-drawn smooth curves have these same average values, interval by interval. Behavior of the curves at the end of the period was dictated by reference to the average values (not shown) over the interval between flights 30 and 31.

Fig. 4. Width-averaged surface altitude as a function of time for different locations. See Figure 2 for positions of grid rows. Sea-level (m a.s.l.) refers to the National Geodetic Vertical Datum of 1929.

As shown in Reference RasmussenRasmussen (1985), the factor γ depends on the fraction of the total motion V that is due to the sliding component V _b and on the exponent n in the flow law assumed to govern the deformation component, whose value at the surface is denoted V _d.

Assumed for Columbia Glacier, along with the commonly used value n = 3, was that V _d was equal to half the minimum value of V observed over the 4.25 year period. There is no basis for using one-half other than that it is intermediate between the extremes of assuming either that there is no deformation or that the minimum consists entirely of deformation; as mentioned below, the results are only weakly sensitive to the fraction used. Consequently, as the maximum value attained by V over this period was about twice its minimum, γ varied from node to node and from interval to interval in the range of 0.90–0.95, with an average value of 0.936.

Because the motion is mainly due to sliding, the results, as analyzed by Reference RasmussenRasmussen (in press), are not very sensitive to the value chosen for the flow-law exponent. No attempt is made here to interpret the velocity fields in terms of numerical values for the flow-law parameters or their possible spatial and temporal variation. The only inference made here about the flow regime is that V _d varied little over the period, owing to the slight changes in glacier geometry, and that the variation in V must have been due almost entirely to variation in V _b.

A centered finite-difference approximation of Equation (1) is

substituting Equation (2) for Q in terms of the components u and ν of the velocity vector V; the nodal subscript scheme is indicated in Figure 1. This equation is applied at the 77 interior nodes in the lower reach and involves 205 velocity components: both u and ν at 85 nodes, only u at 29 nodes along the lateral margins, and only ν at six nodes at the top and bottom of the area.

Fig. 5. Estimated errors in the initial estimates of the average v-component over the interval from 22 July to 2 September 1980 for the lowest nine grid rows. The gradient effect is reflected by the large values along the high-gradient margins and the small values in the low-gradient interior. The density effect is reflected by the values at the ends of the third row from the bottom; the effect of the gradient is stronger near the 42 than near the 87, but it is overcome by the differences in the data densities in the vicinities of the two nodes. The dashed contours are intermediate in value between the solid contours.

If values are assumed for the h̄_ijL and F _¡jL for some flight interval L, it is possible to find values of the (u,v) _¡jL that satisfy Equation (4) and that minimize

in which the (u₀,v₀) _ijL are the initial estimates of the velocities. That is, of all the vector fields that satisfy the continuity equation, the (u,v) _ijL are the ones nearest, as measured by D _L; to the initial estimates (u₀,v₀) _ijL. A linear, non-iterative algorithm (Reference RasmussenRasmussen, 1985) introduces a stream function to find the minimizing (u,v) _ijL as the result of four independent, interlacing solutions covering the entire region.

The algorithm can be used as a component of a broader procedure for finding the bed altitudes B _ij and linear balance functions b _L(z) that minimize the 21 interval sum

The h _ijL in Equation (4), as applied to interval L, may be replaced by the observed Z _ijL and the B _ij (from the relation h = Z – В) along with the γ _ijL calculated (Equation (3)) from the initial estimates of velocity. The B _ij are assumed to be constant over the 4.25 year period, and the simple averages of the surface altitudes at the beginning and end of the interval are used for the Z _ijL. The F _ijL are replaced by the observed and . For convenience, the balance function is specified by its values and at z = 100 and 500 m.

What results is the problem of minimizing D, which is non-linear because of the h u and h v terms in Equation (4). The problem is determined by the 4305 observed velocity components, 205 of them for each of the 21 flight intervals. There are M = 162 model parameters, the B _ij at the 120 nodes and the two balance-function parameters for each interval. The observed Z _ijL and Ż _ijL, which are as numerous as the velocity components, are taken to be exact (Rasmussen, in press). The errors in Z _ij on the dates of successive flights are believed to be strongly correlated positively (Rasmussen, in press), so that error in the altitude change tends to be the difference of those errors, which is small; this inference is supported by the regular spatial variation of the changes (Reference Rasmussen and MeierRasmussen and Meier, 1985, appendix B).

Systematic methods of solving the problem were rejected in favor of a more economical parameter-by-parameter method. As an example of the former, Newton’s method may be written as

Here, ΔΡ is the column vector of the adjustments to be made to initial estimates of the M parameters p _I …, p_M, and D' is the column vector of first derivatives of D with respect to the parameters. The elements of the symmetric M×M Hessian matrix D” are the mixed second derivatives of D with respect to pairs of the parameters. That is, in row I of Equation (7): the element of ΔP is ΔP _I, the element of D' is ∂D/∂p _I, and the element in column J of D” is ∂ ² D/∂p _I ∂p _J. After the equation is solved, parameter P _I is revised by adding ΔP _I to it. Both D' and D” are formed from initial estimates of the parameters. The use of Equation (7) to estimate the bed topography and the mass-balance distribution is outlined in Figure 6.

Because M = 162 for the Columbia Glacier data, D” would have 13 203 distinct elements. Of these, only the 840 with respect to two balance parameters from different flight intervals would be identically zero. Numerical finite-difference approximations would have to be used for the elements of D' and D”. Because four values of D are needed to approximate each of the second derivatives, the order of 50 000 evaluations of D would be needed to form the coefficients in Equation (7); each one involves using the linear algorithm to find the minimizing (u,v) _ijL and D_L for each of the 21 flight intervals. Moreover, one solution of Equation (7) does not produce the values of the 162 parameters that minimize D; instead, it must be solved repeatedly, beginning with the initial estimates of the parameters, and revising the parameters each time it is solved, until they converge to their minimizing values.

There is no guarantee that any systematic method will converge or that, if it does, that it will converge to a global minimum. Nor is there a guarantee that the parameter-by-parameter method has these properties. In the case of the Columbia Glacier data, however, the latter method did converge to an apparent global minimum.

Fig. 6. Algorithm for inferring bed altitudes and mass-balance distribution. Denoted by (i,j) are the grid nodes and by L are the intervals between successive dates of photography. The objective is to find the vector P that minimizes D, the resultant of the adjustments (Equation (6)) in applying the linear algorithm (Reference RasmussenRasmussen. ¡985) to the (U₀V₀) for each L. That algorithmn finds the velocity field that is consistent with P according to Equation (4) and that is the closest, as measured by D_L (Equation (5)), to the (u₀,v₀).

Initial estimates of the parameters used with the parameter-by-parameter method could have been used by a systematic method as well. The bed inferred from the airborne radio recho-sounding was used where it existed and was continued into the upper part of the area at about the same depth and shape; along the lateral margins the trend of the exposed rock topography was followed. A zero-flux divergence F was assumed over the entire region, which is equivalent to assuming that the observed changes in surface altitude were due entirely to balance processes.

In contrast to the systematic method, the parameter-by-parameter method required only about 300 evaluations of D in each of the six or so iteration cycles for which it was used. In the first cycle, the bed altitude B_ij at a node was revised by considering what the linear algorithm did there for each of the 21 intervals, as described in the last section of Reference RasmussenRasmussen (1985). If the linear algorithm generally increased the initial estimate of the velocity, B_ij was lowered; if it decreased the initial estimate, then B_ij was raised. The rationale was that increasing the velocity, for example, indicated the need for a greater flux, which could have been attained by using a greater thickness with the existing velocity, so B_ij was lowered. In subsequent cycles, B_ij at a node was revised by taking it to be the B_ij giving the minimum D on a parabola through the three points D(B_ij – 10 m), D(B_ij), and D(B_ij + 10 m), which were obtained by holding the other model parameters constant and varying only that B_ij .

As it turned out, the adjustment of any particular B had little influence on the apparent minimizing value of B at nodes more than two rows or columns away. In addition, adjusting B at a low-velocity node had almost no effect on the minimizing B at any high-velocity node. This behavior, along with the mutual independence of balance parameters from different intervals, indicates that the matrix D would be diagonally dominant.

After all the B_ij had been revised in one iteration cycle, the minimizing balance parameters were estimated as the result of a separate two-parameter search for each of the intervals. In every balance-parameter optimization — for each interval, for each iteration cycle — the required adjustment D _L appeared to be a convex function of the two parameters. A typical example is shown in Figure 7, which illustrates the optimization for interval 26 during the last iteration cycle, in which the finally adopted values of the B_ij were used. After all the balance parameters had been revised, they were used as the basis for further revision of the B_ij in the next cycle.

Fig. 7. Sensitivity of the required velocity adjustment to the balance parameters, for interval 26 (2 September–30 October 1980). The contours show the relative increase in D₂₆ when the parameters are varied experimentally from their adopted values.

The procedure was continued until D was reduced to 0.427, which is probably very near its absolute minimum. That is, the adjustment of the velocity vectors required to make them satisfy continuity was, on the average, less than half the estimated error in the intial estimates of those vectors. No interpretation of the value 0.427 is made here except to suggest that it simply may indicate the amount by which the standard errors ε _u and ε _v were overestimated. The sensitivity of the results to the B _ij is shown in Figure 8 for the case where the finally adopted values are used for all the bed and balance parameters. It does not exhibit the same rough uniformity of response that does exist, from interval to interval, with the balance parameters (Rasmussen, in press). Along the low-velocity margins, Dis almost independent of the bed altitude used. In the higher-velocity interior, however, the adopted altitudes are very near the optimizing values.

Fig. 8. Sensitivity of the velocity adjustment to the bed-topography parameters. The graph at each node shows the effect of changing the bed altitude at only that node, while the adopted bed altitudes are used at all other nodes, and while the adopted values of the mass-balance parameters are used. The curve is a parabola through D(B _ij – 10 m). D(B_ij), and D(B _ij + 10 m). and is scaled by dividing by D(B _ij ). The adopted value of the bed altitude is B _ij and D is the magnitude (Equations (5) and (6)) of the velocity adjustment averaged over all 120 nodes, over all 21 flight intervals. The height of the vertical line through each node corresponds to an increase of D by a ratio of 1.001; for purposes of vertical exaggeration, it occurs in the diagram in three different lengths, in the proportion 1:2: 4.

Results

Although it is economical of computer time, the parameter-by-parameter method demands intimate, often tedious, participation by the person using it. For this reason, and for others, it was not continued to the point of finding the absolute minimum of D. An internally consistent data set for Columbia Glacier has been created that closely resembled the raw data. It had been established that the bed topography and mass-balance distributions could be estimated solely from photogrammetric data. Evidence of incompleteness of the solution exists in Figure 8: the curve for node (62,21), for example, suggests that the actual bed altitude there is 5 m higher than calculated, and the one at (61,21) that it is 3 m lower. The balance-parameter searches (Fig. 7) did not achieve absolute minimization either.

Sensitivity of the results to the assumptions about the value of the flow-law exponent and the fraction of the motion that was due to deformation were not investigated experimentally, but analysis of the equations (Reference RasmussenRasmussen, in press) showed that the inferred mass-balance parameters were nearly independent of those values. It was concluded that a +1 m change in the inferred bed altitudes in the interior would result if either the flow-law exponent was increased to about 3.15 or if the fraction of the minimum velocity assumed to be due to deformation was decreased to about 0.48.

Because the final interpolation of the velocities was close to the intial estimates of them, the results are susceptible to bias that may exist in the initial estimates. Although hand contouring of the photogrammetric displacements was done carefully, that process is partly subjective. If, for example, the velocity at some node was consistently overestimated, the bed altitude there would also be overestimated in order to preserve the product of thickness and velocity, and thus (Equation (2)) the flux. Near the terminus and at the top of the solution region the contouring was hard to do, and thus these places are likely sites for such distortion. Bias may also exist where the contouring was easy to do.

The inferred bed topography is shown in Figure 9. The contouring is not continued to the margins because of the unreliability of the inferred altitudes there. The initial estimates there — made by continuing the trend of the rock topography — were not revised. The contours, therefore, offer only weak evidence that Terentiev Lake (Fig. 1) will not become an arm of the fjord after the glacier terminus retreats past it. The contours generally agree (Rasmussen, in press, fig. 15) with the bed inferred from airborne radio echo-sounding (Reference Brown, Rasmussen and MeierBrown and others, 1986) where it had good coverage, which was in the interior of the lower part of the region. It is not possible to make an absolute comparison because of the large error (30 m) believed to exist in the radio-echo altitudes, and because the two beds were referred to different grids.

The inferred mass-balance distribution is shown in Figure 10, in the form of

for both z = 100 and 500 m. Because b(z) is taken to be linear over each interval, [z,b _z (t)] is collinear with [100, b ₁₀₀(t)] and [500, b ₅₀₀(t)] for any t, and any z. The curves are physically reasonable, overall, but contain several minor irregularities. They exhibit the usual seasonal progression, rising in winter and falling in summer at rates typical of glaciers in this maritime region. The balance gradient db/dz is positive in all but six intervals, three of which occur in winter; when it is negative, it is only mildly so. The most positive gradients occur in spring (intervals 14, 20, 23, and 28), the season when high-altitude ice is still covered by snow, and low-altitude ice is exposed. Average rates are −7.5 m/a at 100m altitude and 0 at 500 m; because the equilibrium line is above 500 m, the latter is slightly too positive.

Two error sources affect the balance parameters far more than they do the bed parameters, as analyzed in Rasmussen (in press). Because the photogrammetry measures only the geometric surface, the accumulation and then ablation of the lower-density snow-pack will induce positive errors in the balances inferred for winter and then negative errors in those inferred for spring, for the balances are defined to be in ice-equivalent units. Therefore, the winter maxima of the b _z (t) curves should be subdued to account for the lower density.

Fig. 9. Inferred bed topography, in m a.s.l.. hand contoured (50 m interval) from bed altitudes at the grid nodes. The contours on the exposed rock (100 m interval) were directly determined photogrammetrically.

If systematic error exists in the photogrammetry, the inferred bed altitudes will be in error by its 21 flight average, and its flight-to-flight variations will nearly all be absorbed in the balance parameters. The first effect is innocuous, for it simply modifies the datum plane to which both bed and surface altitudes are referred. The second, however, can seriously alter the inferred balance change over the interval between two flights with different systematic errors. An example is flight 14, whose systematic error is suspected of being 2 m more negative than the error for flights 13 and 15. The more positive slopes during interval 14–15 in Figure 4 are not supported by the meteorological records at nearby National Weather Service stations Valdez and Cordova.

Ablation was measured at ten stakes ranging in altitude from 130 to 540 m on about the dates of flights 11 and 17 (Reference Mayo, Trabant, March and HaeberliMayo and others, 1979). When the inferred b(z) were integrated over that August-to-August period, they were about 1 m more positive than the stake measurements but showed about the same db/dz gradient. The comparisons for the winter and summer parts of that period (Rasmussen, in press, fig. 20) was poorer. The major difference between the inferred values and the measurements is that the latter did not closely follow any simple dependence on altitude, let alone a linear one, indicating perhaps that the stakes were subject to highly localized conditions of accumulation and ablation.

Discussion

It may be possible to estimate the bed topography and mass-balance distributions of other ice masses by using aerial photography or other remote-sensing instrumentation. If that is undertaken, further automation of the minimization process that produces those estimates may be appropriate. An objective method of forming the initial velocity estimates would eliminate the need for hand contouring, but such a method would have to take into account similarity of the fields, from interval to interval, for the result to be better than that of careful contouring. Using a systematic method of finding the model parameters would eliminate much tedious, time-consuming work but at the cost of substantially more computer time than the parameter-by-parameter method uses. Modifying the linear algorithm to treat zero-flux nodes explicitly would probably not be worth the effort, as they can now be accommodated by making a zero-velocity initial estimate and assigning a small error.

Estimating mass-balance distributions alone is computationally much easier than estimating them simultaneously with bed topography. If the bed had already been determined by independent means, such as echo sounding or by using the method of this paper during the early stages of an aerial photography program, subsequent photography could be used only for estimating the mass balance. The functional dependence of could then be generalized to ten terms, say, so that horizontal gradients and non-linear effects could be estimated.

Fig. 10. Inferred mass balance shown as cumulative change since the beginning of the period, for two indicated altitudes. At any time, the cumulative balance change is linear with altitude.

Estimating both bed topogaphy and mass-balance distributions from aerial photography might have been possible for the lower part of Columbia Glacier because of properties it possesses that some other ice mass may not. Of course, photogrammetry must be feasible for determining both surface altitudes and displacements, which requires a sufficient surface roughness. The immunity of surface features from obliteration by seasonal snow cover is essential to getting the seasonal variations of velocity and, through the minimization method, of the mass-balance distribution. The temporal variation of the results is fixed at the frequency of the photography. Quality of the photogrammetry is also significant: inferred mass-balance parameters are very sensitive to systematic error; the inferred bed parameters, as well as the adjusted velocity fields, are sensitive to the spatial density of the photogrammetric displacements.

The method may not work well on an ice mass with low velocities, because the thickness, and hence the bed altitude, become undefined as the velocity goes to zero. The seaonal variation of velocity at Columbia Glacier may not be important nor may the strength of the divergence flux there be.

The fluctuation of surface altitude, however was essential to getting results for Columbia Glacier. If , the continuity equation would reduce to

and the results would have been good only within a constant factor. If h(x,y) and (z) satisfy the reduced equation, then so do ch(x,y) and (z) for any constant c. The non-zero injected absolute scale into the results. An ice mass in long-term steady state may have sufficient seasonal variation of Z to do this, but photography more frequent than annual would be required.

If long-term Z changes did not exist for some ice mass, and if seasonal Z changes either did not exist or were not measured, then independently determined thicknesses could be used to inject absolute scale. Such thicknesses should be specified at nodes with high velocity, so that they exert a strong influence on the flux fields. Because the linear algorithm decomposes the solution area into four sub-sets, as described in Reference RasmussenRasmussen (1985), thickness must be specified both at a node whose row and column indices add to an even number, and at a node whose indices add to an odd one.

Using the data set to study the ice-flow mechanism will, alas, require much more than just relating the velocity to the geometry, as conventional flow models do (Reference HookeHooke, 1981). In these, the velocity follows a power-law function of the product of thickness and surface slope; n = 3 is usually used for the exponent, as was done here in treating the vertical profile of velocity (Equation (3)). As mentioned earlier, the inferred bed and mass-balance distributions are relatively insensitive to the value chosen. No such law, it would seem, can by itself fit the Columbia Glacier data. The slight changes in thickness and slope (Fig. 4) account for only about 2% (Reference Meier, Rasmussen, Krimmel, Olsen and FrankMeier and others, 1985, p. 34) of the velocity changes (Fig. 3).

Much recent attention has been given to how subglacial water modifies the dependence of velocity on geometry (Reference Alean, Braun, Iken, Schram and ZwostaAlean and others, 1986). Concurrent measurements of the water are, unfortunately, lacking for Columbia Glacier. This makes analysis much more difficult but perhaps not impossible; the water is not a random, independent regime but must obey known physical laws. If two postulates are combined — one relating velocity to the geometry and the water regime, another relating the water regime to the geometry and velocity — it might be possible to eliminate the water regime mathematically and computationally between the two, leaving a relation between velocity and geometry that will probably be much more complicated than the usual power law. It seems formidable.