Solution procedure and mesh geometry

The numerical calculations were done with the general-purpose finite-element (FE) program MARC (MARC Analysis Research Corp., 1995). The model uses an eight-node isoparametric, 3-D brick element with trilinear interpolation. An effective viscosity distribution is determined iteratively from a previously calculated strain-rate distribution until the maximum relative change in velocity between two successive iterations over the whole mesh is <0.1%. A detailed discussion of the solution procedure and a comparison of numerical calculations with theoretical solutions of analytically solvable problems can be found in Reference GudmundssonGudmundsson (1994b).

The outlines of the FE mesh are shown in Figure 1. The bedrock geometries of Unteraargletscher and most parts of Finsteraar- and Lauteraargletscher are known from radio-echo soundings (Reference Funk, Gudmundsson. and Hermann.Funk and others, 1995). The surface geometry of the FE mesh is based on a digital terrain model from the year 1991.

Fig. 1. A perspective plot of the FE mesh of Unteraargletscher. Only the outlines of the elements are shown.

The mesh has 7920 nodes and 6424 elements. Typically the edges of the elements have lengths around 60 m in the horizontal plane, and about 40 m in the vertical plane. To estimate the discretization error, a coarser mesh having 3732 nodes and 2800 elements was used for comparison. The coordinates of the FE nodes correspond to the official Swiss coordinate system. The x axis points east, the y axis north and the z axis upwards.

Field equations and constitutive law

The field equations expressing the conservation of mass for an incompressible material, and the conservation of angular and linear momentum can be written as

(1)

and

. (2)

respectively. In the above equations, v _i are the velocity components, ρ is the ice density, f _i, are the components of the body force and σ _ij are the components of the stress tensor. The notation is summarized in Table 1.

Table 1. List of symbols

Glen’s flow law relates the deviatoric stresses to the strain rates and may be written as

. (3)

where are the strain rates, σ ^′_ij are the deviatoric stresses, σ ^′_Π is the second invariant of the deviatoric stress tensor

. (4)

and . Glen’s flow law has traditionally been used for modeling purposes in glaciology, although there are reasons to believe that it may be based on some oversimplification of the actual rheological behavior of glacier ice (e.g. Reference HookeHooke, 1998). Moisture and debris content, ice fabric and ice composition, for example, affect the rheological properties of ice. Since these factors vary from one glacier to another, the values for A and n that are in some sense appropriate for one particular glacier cannot be used with confidence for other glaciers. Furthermore, the flow-law parameters n and A are possibly only approximately independent of the stress regime (Reference HutterHutter, 1983; Reference Pimienta and Duval.Pimienta and Duval, 1987). In fact, the search for a single flow law for ice may be illusory (Reference Lliboutry and Duval.Lliboutry and Duval, 1985). Although Glen’s flow law may not be physically the most appropriate one, it seems to be flexible enough to give a reasonable fit to a number of field observations (Reference PatersonPaterson, 1994). There is also a practical argument for using Glen’s flow law since its use allows an easy comparison with results from other areas. The flow-law parameters A and n are here considered to be purely phenomenological parameters that can be varied freely in an attempt to obtain the best possible agreement between calculated and measured quantities.

Boundary conditions

The upper boundary is a free surface with no forces being exerted on it. The purpose of the modeling calculations is to obtain a snapshot of the velocity and stress fields resulting from the momentary geometry and therefore a mass-balance distribution does not have to be prescribed.

It is known that surface velocities on Unteraargletscher are subjected to temporal variations on various time-scales, ranging from a few hours to a few months (Reference FlotronFlotron, 1973; Reference Iken, Röthlisberger, Flotron. and Haeberli.Iken and others, 1983; Reference GudmundssonGudmundsson, 1996), indicating the occurrence of basal sliding. The velocity fluctuations seem to be limited to the melt season, as measurements on the confluence area showed winter velocities, whose magnitudes are about 75% of summer values, to be stable throughout the winters of 1991–92 and 1992–93 (Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997). The stability of the winter velocities indicates that basal boundary conditions do not change throughout the winter. Basal sliding is then definitely reduced with respect to summer, and possibly even absent. Assuming that basal sliding is insignificant during the winter period, it is ignored and the focus put on the winter velocity regime. A no-slip condition is thus imposed along the lower boundary. This is, admittedly, a somewhat problematic assumption as the possibility that some basal sliding takes place during the winter months cannot be entirely excluded and this can evidently have some effect on the results. On the other hand, the model is in this way kept as simple as possible, having only two adjustable parameters A and n. Including basal sliding in the model by using a sliding law would be possible in principle, but it is not clear how the sliding law could be verified. Although recent results from continuous borehole inclinometry indicate that almost half of the forward motion of Unteraargletscher at a location about 2 km down-glacier from the confluence is due to basal sliding (Reference Gudmundsson, Bauder., Lüthi, Fischer. and Funk.Gudmundsson and others, in press), no similar measurements are available for the confluence area.

Across all three cross-sections of the model (Fig. 1) either stress or velocity distributions must be defined as boundary conditions. It was decided to define the normal and the shear stresses across the transect boundaries rather than to define the velocity distribution, because in that way the applied boundary conditions are independent of the value of the rheological parameter A used for any particular model run. With this approach, only the value of n and not the value of A needs to be adjusted during a model calibration. The value of A that would have led to a best overall agreement between measurements and calculations for some value of n can then be found afterwards through a simple scaling procedure.

An iterative procedure was employed when fixing the applied stress distribution for the cross-sections of the numerical model. In a first step all cross-sections were approximated by ellipses, and the analytical solution for flow in a slightly elliptical channel used to obtain an estimate of the stress tensor across the transect boundaries (Reference NyeNye, 1965). The applied stress distribution was then perturbated in a trial-and-error procedure in order to obtain the smoothest possible longitudinal variation of the calculated stress and velocity fields in the vicinity of the cross-sections. The degree of agreement between calculated and measured velocity gradients was then used as a measure of the appropriateness of the boundary conditions.

Statistical description

A number of available surface velocity measurements were used for the model calibration. The velocity measurements are discussed in Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997 and listed in Reference GudmundssonGudmundsson (1994a). As mentioned above, calculated velocities depend linearly on A. For a model run, with some particular pair of values for the rheological parameters n and A = A _c, the optimal estimate of A, referred to as A _b, is found by minimizing

(5)

with respect to f _A. Here v^m(r _i) are measured velocities at the locations r_i, v^c(r_i) are the corresponding calculated values, and N is the total number of points for which measured velocities are available. A best overall value of A, in a root-mean-square sense, is then given by A _b = f _A A _c, where A _c is the value of A used in the model run. From Equation (5) it can be seen that f _A is given by

. (6)

The normalized root-mean-square error, S _v, is defined as

. (7)

and the normalized difference between the measured and the (scaled) calculated velocity vector at the position of marker i is given by

(8)

where R_i is a dimensionless error vector. The subscript “i” refers to the ith marker located at r_i. Note that is the magnitude of the vector and not the ith component of that vector.

A rough measure of systematic errors is the mean vector sum of all error vectors,

. (9)

If all error vectors were randomly distributed, their sum would represent a random walk. The length of the vector sum, R, would then increase as the square root of N, so that R would tend to zero as N → ∞.

Optimal estimates of flow-law parameters A and n

The optimal values of A, found by minimizing Expression (5) for different values of n, are listed in Table 2 together with the normalized root-mean-square error S _v and the length R of the sum of the error vectors. Only winter measurements obtained from repeated surveying of markers with time intervals larger than 1 month were used in the calculations leading to the values in Table 2. A number of velocity estimates based on repeated surveying with a time interval of 3–4 days were not used, since the errors of the resulting velocity gradients were considered too large. For markers surveyed more than twice, averaged velocities were used.

Table 2. Values of A _b, S _v and R for different values of n

In Figure 2 the value of R is shown as a function of n. As can be seen, R is a convex (concave upward) function of n with a minimum at n = 3, suggesting that n = 3 is the optimal estimate of n. If N random vectors having length l are added together, their sum will, on average, have length . The medial length of R_i was about 0.18 and N = 39, so that a pure statistical variation is expected to give R ≈ 0.18/ = 0.016, a number almost twice as large as the value of R obtained for n = 3. This is only a rough estimate of the contribution of random errors to the length of R, since not all the vectors R_i had the same lengths. A simple estimate of the errors involved can be obtained by using different numbers of marker velocities to calculate the variation of R with respect to n. Several calculations made by using only 80% of the measured velocity vectors revealed that although the shape of the R(n) curve changes somewhat as different subsets of available measurements are used, the minimum of the curve is always found for n = 3.

Fig. 2. R as a function of the parameter n in Glen’s flow law. R is the magnitude of the sum over all error vectors R_i, and is defined in the text. The minimum of R is reached for n = 3, indicating that the best overall agreement between measured and calculated surface velocities is obtained for that value of n.

For n = 3 the value of A that gives the best overall fit between calculated and measured winter velocities is A = 2.37 × 10⁻¹⁵ s⁻¹ kPa⁻³. This value is almost three times smaller than the value recommended by Reference PatersonPaterson, 1994, which is A = 6.8 × 10⁻¹⁵ s⁻¹ kPa⁻³. Including basal sliding in the numerical model would lead to an even smaller optimal value for A and thus make this discrepancy larger. One does not have to resort to a fully 3-D numerical model to realize that the recommended value of A = 6.8 × 10⁻¹⁵ s⁻¹ kPa⁻³ is much too large for Unteraargletscher. Down-glacier of the confluence area the geometry of Unteraargletscher is relatively simple and can be approximately described as a parabolic channel with an aspect ratio of about 2, a slope close to 4° and a thickness of about 300 m. Using the form factors listed in Reference GudmundssonGudmundsson (1997a) indicates that the value of A must be reduced by at least a factor of two from Patersons’s recommended value if calculated medial flow velocities are to be similar to the observed winter-surface velocities.

Spatial error distribution

The variation of the sum of the error vectors with n shown in Figure 2 gives no indication of the spatial variation of the error vectors across the confluence area. Only if the lengths and the orientations of the error vectors are randomly distributed can it be concluded that the numerical model is free of systematic errors and that Glen’s flow law represents an appropriate rheological description of the ice of the confluence.

The error vectors R_i for n = 1, n = 3 and n = 5 can be seen in Figure 3. In addition, the normalized projection (P _i) of the error vectors R_i along the measured flow direction defined as

(10)

is shown as a diamond symbol.

Fig. 3. Relative discrepancy between measured and calculated surface velocities for n = 1 (a), n = 3 (b) and n = 5 (c). Vectors represent a vector quantity R_i defined by , where f_A is defined in the text. Subscript “i” refers to the ith marker. are measured velocities if markers, and are calculated velocities. The size of the diamond symbols is proportional to the length of P_i, defined as the projection of R_i on (Equation (10)). Coordinates are in meters and correspond to those of the official Swiss coordinate system. North points upward.

Markers along two profiles across Lauteraar- and Finsteraargletscher above the confluence, and a profile across the Unteraargletscher below the confluence area, were measured with rather short time intervals of approximately 4 days. The associated measurement errors in the velocities are correspondingly large. Since the surveying errors are estimated to be about 3 cm, the errors of the velocities are about 4 m a⁻¹. On the other hand, markers along the center line down-glacier from the junction point and along two lines from the marginal zones of the confluence towards its center, were measured repeatedly with time intervals of approximately 40 days, leading to much smaller measurement errors. In what follows, the discussion of the spatial variation of the error vectors will be limited to those markers for which long-interval measurements were available.

For n = 1 (Fig. 3a) the largest error vectors are found at the northern margin of the confluence. The error vectors point approximately in the flow direction. The (scaled) calculated velocities are therefore, at that location, smaller than the measured ones. For the marker closest to the junction point the situation is reversed.

For n = 3 (Fig. 3b) the error vectors at the northern and southern margins also point more or less in the flow direction but their magnitudes are overall somewhat smaller than in Figure 3a. For n = 5 (Fig. 3c) the magnitudes of the error vectors close to the northern margins decrease even further with respect to Figure 3a, but towards the southern margin the error vectors in Figure 3c are larger than in Figures 3a and b. The error vectors lying closest to the junction point in a direction almost exactly opposite to the direction of flow for all three values of n. The magnitudes of the error vectors in this area increase with increasing n.

For all values of n there is a systematic variation in both the direction and the magnitude of the error vectors across the confluence area. This is small enough, however, that it can only be seen at those locations where accurate estimates of the surface velocities, based on repeated surveying over a period of several months, are available. For n = 3 the errors are distributed almost symmetrically about the medial line. Calculating the error vectors for further values of n revealed that for n < 3 the errors tend to be larger along the northern than the southern margin; for n > 3 the situation is reversed. For n = 5 the only significant errors are close to the junction point and at the southern margin, but there is still a clear tendency for calculated marginal velocities to be smaller than measured ones.

Despite the non-randomness of the error vectors, their magnitudes are usually less than about 15%, an accuracy that is probably sufficient for most modeling purposes. Since the general flow characteristics of the confluence area form the main focus of this study, an error of this magnitude is considered acceptable.

Velocities

The calculated horizontal velocity distribution along the surface for n = 3 and A = 2.37 × 10¹⁵ a¹kPa³ is shown in Figure 4a. The corresponding vertical velocity distribution can be seen in Figure 4b. The 3-D distribution of the vertical velocities is depicted in Figures 5 and 6 as isosurfaces of constant velocity for n = 1 and n = 3, respectively.

Fig. 4. (a) Calculated horizontal surface velocities (vectors) and calculated horizontal speeds (contour lines). Comparing this velocity distribution with measured winter velocities (see Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 5) reveals a good qualitative agreement. (b) Vertical velocities. The five flow features discussed in the text are indicated by Nos. 1–5. The values of the flow-law parameters are n = 3 and A = 2.37 × 10⁻¹⁵ s⁻¹ kPa⁻³. Units are m a⁻¹.

Fig. 5. Isosurfaces if vertical velocities v_z for A= 1.90 × 10¹¹ s⁻¹ kPa⁻¹ and n = 1. Vertical surface velocities are generally positive. A narrow zone of negative vertical velocities stretches from the junction point (J) down the whole Unteraargletscher. The downward movement if ice is particularly large in the vicinity of J. North is approximately to the right.

Fig. 6. Isosurfaces of vertical velocities v_z, for A = 2.37 × 10¹⁵ s⁻¹ kPa⁻³ and n = 3.

Figure 4a shows a good qualitative agreement with measured surface-winter velocities (Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 5). The calculated velocity distribution has a local maximum close to the center of the confluence, as do the measured winter velocities. Measured summer velocities do not have a maximum at this location (Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 7). The calculated surface velocities of Finsteraargletscher are larger than those of Lauteraargletscher, decrease somewhat as both tributaries enter the confluence, and then increase again. All these qualitative aspects of the calculated velocity pattern are found in the measured velocity fields (Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 5).

As seen in Figures 4b, 5 and 6, the numerical model predicts a complicated spatial pattern of vertical velocity distribution across the surface and with depth for both n = 1 and n = 3. Almost no direct measurements of vertical velocities during the winter period are available, as such measurements have proven logistically difficult to make. In addition, information on the vertical velocity distribution has been considered to be of secondary importance for model calibration. This opinion may be subject to change. As discussed in detail below, it appears that the effects both of local topography (having spatial scales of a few hundred meters) and of the convergence of the two tributaries on the flow field can be seen most clearly in the spatial variation in vertical velocities across the surface.

Vertical flow feature 1

On the basis of the 2-D map-plane model developed and analyzed in Reference GudmundssonGudmundsson, 1997b, it is predicted that the abrupt change in boundary conditions at the junction point (J) leads to the formation of a surface depression located slightly down-glacier of J. The spatial size of the surface depression depends primarily on the value of n, decreasing in size as n increases.

The semicircular zone of highly negative vertical velocities (Figs 4b, 5 and 6) close to J does coincide with a strongly localized surface depression found in nature. This surface depression exists in the numerical model, but is not accurately described by it due to the coarseness of the FE mesh. It does not, as might be expected, give rise to positive vertical velocities along the surface as ice from the sides flows towards the depression to fill it. Instead the surface velocities are downward-oriented, suggesting that the surface depression results from the flow anomaly. A flow anomaly of this type has been predicted to be an essential feature of ice flow at all glacier confluence areas.

The “pressure drop” observed in the 2-D map-plane models (Reference GudmundssonGudmundsson, 1997b) is associated with the sharp velocity increase along the medial line at the junction point, which in turn is caused by the abrupt change in boundary conditions at that location. This pressure drop, which would lead to a surface depression of the type observed on Unteraargletscher, is related to the overall changes in flow regime as ice enters a confluence area. In this respect, the 3-D model confirms the conclusions drawn from the 2-D map model. Evidently, the ice flows faster away from the junction point in the down-glacier direction than it can flow laterally towards it, unless transverse gradients are set up, which then drive the ice towards the junction point.

The 2-D map-plane model predicts not only a localized surface depression close to the junction point, but also two broad surface bulges along the sides of the confluence. As discussed below, flow features 2 presumably represent the upward flow needed to maintain these bulges.

Calculated horizontal and vertical velocities along a profile from the junction point in the down-glacier direction are shown in Figure 7. This may be compared with Figure 4a in Reference GudmundssonGudmundsson, 1997b which is based on the 2-D map-plane model. The increase in velocities at the junction point in Figure 7 is sharp and becomes sharper with increasing n. Consequently, the strain rates and the deviatoric stresses are large. The FE mesh is, however, too coarse to give are exact estimate of the vertical velocities, as can be seen from the fact that the maximum of the vertical velocities is obtained at the FE node lying closest to J. A further mesh refinement would presumably move the maximum even closer to J and give larger negative values. Calculations with a coarser mesh (2800 elements instead of 6424) indeed showed that the position and the magnitude of this maximum depend on the level of mesh refinement. On the other hand, the slope of the horizontal velocity profile (Fig. 7) did not depend significantly on the number of elements used. There is also a tendency for the zone to become more localized with increasing n (Fig. 7), but the coarseness of the FE mesh makes it difficult to make a definitive statement in this respect.

Fig. 7. Calculated horizontal (v_x) and vertical (v_x) velocity profiles along the surface from the junction point (J) towards the center of the confluence. The profile is approximately along a flowline. As discussed in the text, the velocity drop in v_z seen on the left is related to the overall dynamics of the confluence, and leads to the formation of an observed local surface depression. Negative velocities further away from J are, on the other hand, associated with disintegration of the medial moraine caused by its higher than average elevation with respect to the surrounding ice. Symbols represent values at FE nodes.

Around the junction point, the maximum of v _z, as a function of depth, is found at the surface. Several hundred meters further down-glacier (close to the location of a number of boreholes drilled in autumn 1991) the maximum is not at the surface but at a depth corresponding approximately to one-half of the glacier thickness. Even further down-glacier the amplitude of the vertical velocity feature becomes progressively smaller and eventually disappears. Below the confluence center the glacier has adjusted to the new bed geometry, and the divergence of the flowlines over the depth, and the associated vertical velocities, are thus strongly reduced in amplitude.

Strain rates

The isosurfaces of constant vertical strain rates are displayed in Figures 8 and 9 for n = 1 and n = 3, respectively, and for A from Table 2. Negative strain rates correspond to vertical compression, and positive values to vertical extension.

Fig. 8. Vertical strain-rate isosurfaces for n = 1. Vertical strain rates are generally positive (extension). The purple-colored “tube” represents ice being compressed vertically as a result of the weight of the overlying medial moraine. Further south, a similarly shaped yellow-and-orange-colored tube can be seen, which, as explained in the text, is related to the existence if an “inverted” medial moraine directly above it. North is approximately to the right.

Fig. 9. Vertical strain-rate isosurfaces for n = 3. The strain-rate distribution is qualitatively similar to the one seen in Figure 8 calculated for n = 1, but quantitatively large differences can be observed.

The vertical compression resulting from the ice-cored medial moraine can be seen clearly in Figure 8 as an isosurface in the form of a longitudinally aligned cylinder, referred to as “the tube”. This tube is also evident in Figure 9, although less so. Velocity feature 3 is found along the surface area above this tube of vertical compression. Underneath the medial moraine (Figs 8 and 9), the isosurface marking the boundary between positive and negative vertical strain rates (shown in red) is displaced towards the surface with respect to the surrounding ice. Vertical extension is particularly large directly south and north of the medial moraine below velocity feature 4.

Outside of the region directly affected by the medial moraine, vertical strain rates are predominantly positive, as expected for an ablation area. The most noticeable exception to this general rule is found within the confluence area close to the junction point (Fig. 10). As seen in Figures 8 and 9, the tube dips downwards in the up-glacier direction, and reappears at the surface directly down-glacier of the junction point. The complexity of the vertical strain-rate variation with depth along the center line of the confluence is evident in Figure 11. Vertical compression of the basal ice is confined to the area where ice thicknesses increase in the flow direction (Figs 8 and 11). Although the average strain over the depth must be positive where ice thickness increases in flow direction, the no-slip boundary condition forces the vertical velocity component to become zero along the glacier bed, with the result that the zone of vertical compression forms directly above the bed. This at first sight somewhat counter-intuitive point is explained in more detail in Reference GudmundssonGudmundsson, 1997b

Fig. 10. Vertical strain-rate isosurfaces of the confluence area for n = 3, as seen from below the confluence area. The overall flow direction of Unteraargletscher is into the picture plane. Gray arrows indicate flow directions of Finsteraar- and Lauteraargletscher. Vertical strain rates are predominantly negative (compression). A zone of particularly strong vertical compression (–0.15 to –0.05 a⁻¹) is seen extending downward (upward in the figure) from the Junction point (J). Field measurements (Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997) have demonstrated the existence of this zone.

Fig. 11. Vertical strain rates for A = 2.37 × 10⁻¹⁵ s⁻¹ kPa⁻³ and n = 3, across a longitudinal section from the junction point (J) towards the center if the confluence. B indicates the position of a borehole (shown as a thick vertical line) in which vertical strain rates were measured. The section shown is approximately along a flowline. Extension is positive and compression negative. Below point B, strain rates are positive over about the upper two-thirds of the glacier thickness, with maximum strain rates of approximately 0.03 a¹. With increasing depth the bed strain rates become negative, reaching a maximum value if about –0.05 a⁻¹ at bed. Comparison with the measured vertical strain-rate profile (see Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 4) shows a good qualitative agreement. It is difficult to give an accurate estimate if the quantitative agreement between modeled and measured values, because of the small number if available data points.

The vertical strain-rate variation with depth has been measured in a borehole at one location within the confluence (see Reference Gudmundsson., Iken. and Funk.Gudmundsson and others, 1997, fig. 4). The location of the borehole is labeled B in Figure 11. According to measurements, is positive at the surface, decreases with increasing depth, becomes negative at a depth corresponding to about three-fifths of the glacier thickness and becomes increasingly negative towards the bottom. The extension towards the surface, and the compression towards the bottom are nicely reproduced by the numerical model (Fig. 11). For n = 3 the numerical model gives values which are also in reasonably good agreement with measurements. For example, in the vicinity of the bedrock the model gives , whereas measurements give . Since the measured vertical strain-rate variation is based on only four velocity measurements, the qualitative agreement is within observational errors. Vertical compression of the basal ice has been measured only at the confluence area. Measurements made in spring 1997 in boreholes about 2 km down-glacier from the confluence area indicate vertical extension throughout the thickness (Reference GudmundssonGudmundsson, submitted).