1. Introduction
Lewis Glacier on Mount Kenya (lat. 0°09′S., long. 37°19′E.) is the subject of a climate-monitoring program in which both recent climatic change in East Africa and the interaction between climate and glacier are being studied (Reference HastenrathHastenrath, 1975; Reference Caukwell and HastenrathCaukwell and Hastenrath, 1977; Reference Hastenrath and CaukwellHastenrath and Caukwell, 1979; Reference Hastenrath and KrussHastenrath and Kruss, [1981], Reference Hastenrath and Kruss1982; Reference Bhatt, Bhatt, Hastenrath and KrussBhatt and others, [1982]; Reference KrussKruss, in press). For the most important earlier work on Lewis Glacier the reader should refer to Reference Troll and WienTroll and Wien (1949) and Reference CharnleyCharnley (1959). Of considerable interest in the context of the present project is the response of Lewis Glacier to given climatic forcings, an effect which is dependent on the time scale and magnitude of the climatic event as well as on the physical characteristics of the glacier itself. In view of the link between climate, net balance, and the glacier, general characteristics can be investigated by examining net balance variation. A calculation is made of the response of the extent of the Lewis Glacier terminus to repeated sinusoidal netbalance oscillations of various periods and amplitudes.
The technique used involves first the modeling of a steady-state glacier of a certain length. This is accomplished by translating the vertical netbalance profile by the amount necessary to produce the required length. Sinusoidal perturbations are then applied about the displaced net balance and the terminus response is monitored.
For Lewis Glacier, two separate steady states of the correct length are computed; one by translating the curve of net balance versus elevation with respect to the balance axis; the other by moving this curve along the orthogonal elevation axis. The former translation is essentially a precipitation effect with the glacier receiving a snow-fall increase over its entire area. The latter is most directly described as a change in equilibrium-line altitude (Reference Allison and KrussAllison and Kruss, 1977). The repeated sinusoidal fluctuations are then applied along the same axis as the initial translation, and the pseudo-sinusoidal glacier reaction is quantitatively analyzed for terminus-respone amplitude and terminus time lag. Terminus-response amplitude is defined as half the difference between the maximum and minimum glacier lengths measured along a modeled central line. Time lag is the time elapsing between a net-balance event and the corresponding glacier reaction and is measured from applied net-balance extreme to terminus-response extreme.
2. Glacier Dynamics
The numerical glacier model developed for this study is an extension of the deformational model of Reference Budd and JenssenBudd and Jenssen ([1975]). It is essentially a twodimensional, short-time-step model, but a threedimensional parameterization of the continuity equation is included. The basic glacier variables which must be computed are ice velocity and depth and glacier surface elevation and extent. Velocity, depth, and surface elevation are determined at grid points spaced evenly along a longitudinal modeled line, whilst terminus position is followed between grid paints. All values are re-calculated over repeated short time intervals (0.1 year for Lewis Glacier). Earlier versions of this basic model have been employed for studies of various non-surging and surging glaciers (Reference BuddBudd, 1975; Reference Allison and KrussAllison and Kruss, 1977; Reference Budd, McInnes, Pittock, Pittock, Frakes, Jenssen, Peterson and ZillmanBudd and Mclnnes, [c1978]). Close matching has been found between calculated and observed values of ice depth, surface velocity, and glacier extent.
For the Lewis Glacier calculations, a simplified flow scheme derived from the general model is used. Reference Budd and JenssenBudd and Jenssen ([1975]) suggest that for ordinary glaciers the large-scale average-velocity distribution may be computed assuming deformation only with no basal slip. Reproduction of the velocity with this assumption is sufficient for this study. For nonsurging, temperate glaciers in middle and high latitudes, Reference Budd and JenssenBudd and Jenssen ([1975]) and Reference SmithSmith (unpublished) employed computation schemes which includes only a deformational ice-velocity formulation. Reference Allison and KrussAllison and Kruss (1977) modeled the retreat of the (tropical) Carstensz Glacier also with a deformationai model. The matching between computed results and observed velocity values found for a range of valley glaciers using a non-sliding velocity scheme suggests the possibility of employing such an approach for this study of Lewis Glacier, which is a slow-moving temperate glacier (maximum velocity less than 5 m a−1 in 1978).
The present model, excluding the basal-sliding segment, is outlined in Table I: Table I defines one complete computation cycle in the sequence programmed, each pass through this sequence constituting one time step. The equation numbers in this table refer to equations in the text. The total surface velocity U
s (see total velocity block of Table I) is the sum of the surface deformationai velocity V
s (Equation (2) in Table I) and the sliding velocity V
b. The deformational-velocity scheme follows from Reference Budd and JenssenBudd and Jenssen ([1975]). The change in ice depth ∆Z over a time step, and hence the ice depth Z itself, is determined from considerations of glacier net balance and internal ice re-distribution, whilst ice-surface elevation E is the bedrock elevation b plus the ice depth (see ice-depth block of Table I). The exact glacier lengthy X is also determined from continuity (Equations (16) to (18)). Reference BuddBudd and Jenssen (1975) during their modeling found the second right-hand term of Equation (4) in Table I to be very small. This also proved to be true for Lewis Glacier during some initial computations. Hence, this longitudinal stress “correction” term is not included in the computations reported here. The specific flow model employed in this study is described by the Table I equations with 
and V
b both set to zero.
Table I. The Calculation Routine of the Lewis Glacier Flow Model
The importance of 
in the calculation of velocity can also be estimated using an iterative approach to the solution of the set of equations (modified by a stress shape factor) involving longitudinal stress outlined in Paterson (1981, p. 89–91) (personal communication from I.M. Whillans). Solutions were obtained for the regions of mean positive and negative down-glacier velocity gradients. Setting 
to zero was found to be quite an acceptable approximation within the limits of ice modeling.
Regarding the surface velocity, this model was employed by Reference Hastenrath and KrussHastenrath and Kruss ([1981], Reference Hastenrath and Kruss1982) in computing the velocity of Lewis Glacier since late in the nineteenth century; it was found that surface velocities could be well reproduced. Reference Bhatt, Bhatt, Hastenrath and KrussBhatt and others ([1982]) also used this model in estimating the bedrock topography beneath Lewis Glacier in 1978 from known values of the surface velocity and details of the topography. The numerical calculations gave ice depths comfortably within the error ranges of results from seismic and gravimetric analyses. Reference Bhatt, Bhatt, Hastenrath and KrussBhatt and others ([1982] and Reference Hastenrath and KrussHastenrath and Kruss ([1981], Reference Hastenrath and Kruss1982) both discuss the particular suitability of this basic model to the specific situation of Lewis Glacier.
Calculation of the deformational ice velocity is of primary importance in this study. Many researchers have utilized flow laws of ice to relate velocity and shear stress. Reference BuddBudd (1969) discusses several flow laws applicable to differing stress ranges which involve the deformational velocity V i and the basal shear stress τb. In this model, a straightforward power flow law is used as only a comparatively small range need be covered, i.e.
where V i is the vertically-averaged mean deforma- tional velocity, Z is ice depth, and k and n are empirically-determined constants. Reference GlenGlen (1955) and Reference MellorMellor (1959) found n to be between 3 and 4 for stresses greater than one bar. For lower stresses, up to approximately 0.5 bar, Reference MellorMellon and Smith (1966) suggest that n is close to one. Present-day basal stresses of Lewis Glacier are of the order of one bar and less, and hence an intermediate value for n is used. This value is used by Reference BuddBudd (1975) and is near to that of Reference Budd and JenssenBudd and Jenssen ([1975]), However, the value of n is not strictly defined for a valley glacier; if anything, n may be larger than 2. The effect on the computed retreat of Lewis Glacier of changing n from 2 to 3 is explored by Reference KrussKruss (in press). The values of k input for Lewis Glacier are 0.16 and 0.14 bar−n a−1 for n equals 2 and 3, respectively, from a formulation for temperate ice by Reference Budd and JenssenBudd and Jenssen ([1975]).
The surface deformational velocity V s is obtained from V i using
which follows for a power law definition of the velocity vz at depth z from the surface (Reference BuddBudd, 1969) in the expression
Thus, the calculation of surface velocity is dependent on the derivation of a basal-stress solution (see Equations (1) and (2)).
The form of the τb solution is determined by the wavelength of the surface features to be modeled. For wavelengths greater than of order ten times the ice depth, Reference BuddBudd (1971) has suggested that
where 
is the vertical average of the longitudinal stress deviator and τC is the centerline down-slope stress. However, this equation may rather be representational than quantitative (cf. Reference NyeNye, 1979; Reference Hutter, Hutter, Legerer and SpringHutter and others, 1981). Reference NyeNye (1965[a]) found that the frictional effect of valley walls must be included when treating valley glaciers. This effect is embodied in the stress shape factor s of the centerline stress equation
where ρ is ice density and g is gravitational acceleration. The ice surface slope α is defined by
where E is the ice-surface elevation.
The stress deviator 
is found via a flow law applicable to a wider stress range than the power approximation (Reference Butkovich and LandauerButkovich and Landauer, 1960)
where τ1 and ε1 are constants of value 0.3 bar and 0.685 a−1, respectively, for temperate ice (Reference Budd and JenssenBudd and Jenssen, [1975]). V is the vertical average of the total velocity, i.e.
where V b is the basal sliding velocity. The surface deformational velocity may be determined from Equations (1), (2), and (4) to (8). The computation sequence programmed is included in Table I.
The estimation of both ice depth and glacier length is dependent on mass continuity and must include mass addition at the ice surface and internal mass re-distribution due to glacier flow. If we take a volume element from bedrock to surface across the glacier
where Δz is the change in ice depth of the element over the time interval Δt, A is the net balance, and W s and δx are mean width and length of the volume element, respectively. The net ice inflow ΔI is a function of the change, over δx, in crosssection mass flux, i.e.
where V is cross-section mean velocity and Ω is cross-sectional area. In a similar way to Reference NyeNye (1965[a]), V is taken proportional to the centerline velocity, i.e.
where U s is the total surface velocity and C v, the cross-sectional velocity ratio, is dependent on valley shape. If the shape of the valley cross-section is parameterized by a simple power fit of order m, i.e.
where m is the valley power, then
Combining Equations (9) through (11) and (13) for δx small leads to
This equation is employed in the glacier-flow model for the computation of ice depth change over a time step. Thus, the ice depth Z and the surface elevation E are defined at ary time (see ice depth block of Table 1).
The parameters C v and m are constant for each grid point. However, the surface width is recalculated after each time step from input values of reference ice depth and width, Zref, and Wref, respectively. From (Equation 12)
The glacier extent is determined more accurately than the grid spacing by calculating the volume of ice past the last grid point in use, and then computing the length necessary to contain this volume within a specified longitudinal snout shape (Equations (16) to (18) of Table I). The volume of ice at the current time t (M a (t)) in this terminus region is the volume at the preceding time step (Ma (t–1) plus the change over a time step resulting from net balance and ice inflow past the final grid point, i.e. in a similar way to Equation (14),
where r is the distance from the last grid point in use to the glacier terminus. The net balance is a mean over r whilst the remaining values are at the final grid point. A simple power fit of order J to this snout volume and the ice depth at the last grid point gives a snout length
Finally, the total glacier length X is the distance along the modeled line to the last grid point plus r, i.e.
where N g is the number of grid points in use and δX is the grid spacing.
Much of the data input to this model must be defined at grid points spaced 50 m apart for Lewis Glacier along a central, longitudinal modeled line. An important parameter is the bedrock elevation along this line. The bedrock topography is known below the 1978 terminus (4600 m) and at the head of the glacier (4980 m) as the Lewis Glacier commences on a rock slope. For grid points covering the present ice, the bedrock elevation has been determined using three distinct techniques (Reference Bhatt, Bhatt, Hastenrath and KrussBhatt and others, [1982]), a maximum ice depth in 1978 of less than 50 m being found.
The flow cross-sectional area and the glacier surface width at each grid point are approximated by a mathematical power-law representation of the crosssection shape, i.e. the valley power m and reference width W ref and depth Z ref values. As cross-section profiles at grid points within the present ice were not well defined, a valley power m of 2, which was constant for the entire glacier, was used; this value was based on the value appropriate to the exposed valley walls. Definition of m at a value other than unity allows for a modeling of cross-sectional area and surface-width variation with time. The reference width W ref and depth Z ref values used correspond approximately to the surface widths and ice depths at the late nineteenth-century maximum glacier.
Neither the cross-section velocity ratio C v nor the stress shape factor s are well known although they are discussed by Reference NyeNye (1965[a]). Hence, a model tuning process was followed for the final definition of these variables. Employing the Reference NyeNye (1965[a]) information as a guide, C v and a values were initially defined for each grid point. Model calculations were then carried out for a range of C v and s values until the best fit to 1978 ice depths and surface velocities was found. For Lewis Glacier, C v values were about 0.7 and s values were derived in the range 0.8–0.9.
Net balance is defined by altitude bands for Lewis Glacier rather than at each grid point, thus allowing a feedback response between net balance and changes in surface elevation. The net balance employed is that constructed for the 1978/79 balance year (March to March); it exhibits maximum ablation of about 4 m a−1 water equivalent near the glacier terminus and increases to a maximum accumulation rate of about 1 m a−1 for the 4850–4950 m band before decreasing slowly once more at higher elevations.
3. Results
The response of Lewis Glacier is modeled about a steady state of length 1.32 km which is the mean of the late-nineteenth century maximum of 1.60 km and the 1978 minimum of 1.04 km (Reference KrussKruss, in press). Netbalance sinusoidal oscillation periods from 20 to 1000 years are covered, as are applied amplitudes of 0.1 to 0.5 m a−1 of ice (with respect to the balance axis) and 10 to 50 m (with respect to the elevation axis). Figures 1 and 2 summarize the time lag and amplitude of the terminus response for these sinusoidal fluctuations. The plotted time lags are means over the applied net-balance amplitude range and also between the lags at the computed maximum and minimum extents, the lag behavior of the glacier at these two extremes being somewhat different.
Fig. 1. Time lag between net-balance variation and Lewis Glacier terminus response. The mean time lag (ordinate) is plotted against the period of the applied sinusoidal balance oscillation (abscissa) and is expressed in years (squares, left scale) and as a percentage of the applied period (dots, right scale). Dashed and solid lines denote translation with respect to the balance and elevation axes, respectively.
Fig. 2. Amplitude of the Lewis Glacier terminus response. The terminus response amplitude (ordinate) is plotted against the period of the applied sinusoidal balance oscillation (abscissa). A dashed line indicates a response to variation in the balance axis of amplitude 0.3 m a−1 of ice. A solid line indicates a translation with respect to the elevation axis of amplitude 30 m, which is most directly described as a 30 m change in equilibrium-line altitude.
The time lag between sinusoidal net balance event and glacier response (Fig. 1) exhibits a very similar pattern for translation along the balance and elevation axes, with the latter on average about 20% lower in magnitude. The time lag expressed in years increases from a minimum of about 10 years for an applied period of 20 years to a maximum at a 1000 year period of about 30 years. The time lag at 1000 years for elevation-axis translation is not given in Figure 1 because there are instability problems associated with the very small glacier minimum lengths obtained at this long period. Also included in Figure 1 is the time lag as a percentage of applied period. These curves show a rapid decrease from a maximum of more than 50% for a 20 year period to a minimum at 1000 years of only a few per cent. In contrast to this high dependence on applied period, time lag is virtually independent of net balance amplitude. At a given period, the average variation about the mean values given in Figure 1 (see above) is ±5%.
The terminus-response amplitude (Fig. 2) increases rapidly with period for the higher frequencies but approaches a maximum value at periods longer than about 200 years. For the results included in Figure 2, which correspond to translations along the balance and elevation axes of 0.3 m a–1 of ice and 30 m respectively, the terminus-response amplitude varies from about 10 m at a 20 year period to a maximum of more than 150 m. Further, in contrast to time lag, the response amplitude was found to be linearly dependent on applied amplitude (coefficient of determination R 2 equals 1.0). However, a doubling of the applied net balance amplitude does not result in a two-fold increase in terminus amplitude, rather the change is by a mean factor of about 1.9 for Lewis Glacier.
It is relevant here to compare the characteristic results obtained for the present balance-axis study with similar results for other alpine glaciers. The terminus reaction of Hintereisferner, Austria, to sinusoidal net-balance oscillations in the balance axis has been computed for mean lengths of 9.7 and 7.7 km (Reference KrussKruss, unpublished). Reference NyeNye (1965[b]) estimated the frequency response of Storglaciären, Sweden, and South Cascade Glacier, U.S.A.; both glaciers have similar lengths (about 3.5 km) but. Storglaciären moves at only half the speed. The frequency response calculated at 8.0 km for Berendon Glacier, Canada, is included by Reference Untersteiner and NyeUntersteiner and Nye (1968). In the Hintereisferner study, changes in the terminus extent were modeled in the same way as in the present work. For the remaining glaciers, however, the glacier length was held constant and the change in the ice depth near the glacier terminus in response to harmonic net balance fluctuations was calculated. Hence, whilst it is appropriate to compare many of the results of these studies, the Hintereisferner and Lewis Glacier response amplitudes are not compatible In terms of absolute magnitude with the amplitude results for Berendon Glacier, Storglaciären, and South Cascade Glacier.
The five glaciers all exhibit time lags of order 10 years at a 20 year period, while at the longer periods the time lags approach a maximum. At a 1000 year period, the lags for Hintereisferner, Berendon Glacier, Storglaciären, South Cascade Glacier, and Lewis Glacier are about 110, 55, 55, 40, and 30 years, respectively.
The curve of time lag expressed as a percentage of applied period versus this time period exhibits a peak for each glacier except Lewis Glacier. This maximum is at about 55% and 90 years for Hintereisferner, and about 45%, 50 years; 45%, 55 years; and 45%, 25 years for Berendon Glacier, Storglaciären, and South Cascade Glacier, respectively. The Lewis Glacier maximum is evidently at a period less than or equal to 20 years (see Fig. 1) and is at least 5%. At periods below this peak, the time lag asymptotically approaches a value apparently about 25% of applied period (equivalent to 90° phase lag), whilst the curve tends towards zero at longer periods.
The time lag was found to be largely independent of applied net-balance amplitude over the range 0.1 to 0.5 m a–1 of ice for both Lewis Glacier and Hintereisferner. This is an inherent feature of Nye′s analysis (Reference NyeNye, 1965[b]).
The response amplitude approaches a maximum asymptotically at the longer periods. For a period of 1000 years, the terminus amplitudes for Hintereisferner and Lewis Glacier are about 1700 and 160 m, respectively, for an applied amplitude of 0.3 m a–1 of ice. At the shorter periods the terminus response is very small, being about one metre for 0.3 m a–1 of ice applied amplitude and 20 year period. A basic feature of Nye′s approach is a linear dependence between the amplitude of the response near the terminus and the applied net-balance amplitude. Such a dependence was also found for the Hintereisferner and Lewis Glacier terminus responses.
It is relevant to discuss the importance in these. Lewis Glacier calculations of the parameters (C V, s, and J (the cross-section velocity factor, stress shape factor, and the power of the longitudinal snout shape) which are dependent on the glacier dimensions. The terminus response of the Lewis is quite stable to changes in C v and s; there exists a stabilizing feedback between the effect on the integrated mass flux caused by varying these parameters and the opposing influences of re-adjustment in the glacier dimensions and velocity. Consider, for example, the extreme case of a steady-state glacier with net balance and surface width fixed in time. Changing C v or s will not bring about any change in terminus position in this case; rather, re-adjustment is in the velocity and depth at each point within the glacier such that the integrated net balance remains balanced by cross-section flux. A similar though less extreme argument may be applied to this response environment, where only the integrated response of the whole glacier at the terminus is being studied.
Numerical experiments support this assertion. Varying s by 0.2 (25%) results in less than 10% changes in both computed time lag and amplitude responses. Similarly, changing C v by 0.2 (25%) produces variations in lag and length response of at most 5%.
The terminus power J is somewhat different as it must, in part, be set to ensure smooth growth and retreat of the ice past grid points. The parameterization of the near-terminus region embodied by Equations (16) and (17) compensates for the necessarily finite and limited number of grid points. Significant changes of J from its best value will not affect the major features of the response but will result in computationally unacceptable discontinuities in terminus movement through grid points.
4. Conclusions
A number of characteristic features concerning the reaction at the glacier terminus to sinusoidal net-balance fluctuations of various frequencies and amplitudes are apparent. Time lag is dependent on applied net-balance period but essentially independent of applied amplitude. At periods below about 10 years this lag approaches a minimum of 25% of the applied period, whilst at the longer applied periods an approach to a maximum in absolute terms is found. For Lewis Glacier this maximum is about 30 years for variation in the balance axis and about 5 years less for change in equilibrium line altitude. This difference is representative of variations in time lag to be expected depending on the particular climatic forcing involved. It occurs here because the Lewis Glacier net-balance elevation profile is very steep in the ablation zone and relatively quite flat in the accumulation zone; for such curves translation parallel to the elevation axis produces net-balance change concentrated in the lower glacier.
The response amplitude is dependent on both applied period and amplitude, the dependence in the latter case being linear. For a given applied amplitude, the response amplitude is comparatively very small at the shorter periods, i.e. periods on the order of decades and less for Lewis Glacier, and approaches a maximum asymptotically for the longer periods, greater than about 200 years for Lewis Glacier.
Acknowledgements
The guidance of S. Hastenrath in regard to this particular project and W.F. Budd, U. Radok, D. Jenssen, B.J. McInnes, and Ian Whillans concerning model development is gratefully acknowledged. This work was supported by U.S. National Science Foundation grants EAR76−18881, EAR77−13I30, and EAR79−23897 and by University of Wisconsin, Wisconsin Alumni Research Foundation, and Vilas Fellowships.
