Annals of Glaciology 37 2003
# International Glaciological Society

Numerical investigation of the effects of temporal variations
in basal lubrication on englacial strain-rate distribution
Shin SUGIYAMA,1 * G. Hilmar GUDMUNDSSON,2 Jakob HELBING3
1

Institute of LowTemperature Science, Hokkaido University, Sapporo 060-0819, Japan
E-mail: sugishin@pop.lowtem.hokudai.ac.jp
2
British Antarctic Survey, Natural Environment Research Council, Madingley Road, Cambridge CB3 0ET, England
3
Versuchsanstalt fu« r Wasserbau, Hydrologie und Glaziologie, Eidgeno« ssischeTechnische Hochschule, ETH-Zentrum,
CH-8092 Zu«rich, Switzerland

ABSTRACT. The effects of spatial and temporal variations in basal lubrication on

the englacial strain rate and surface velocity distribution are investigated with anumerical ice-flow model. General aspects of the solutions are compared to measurements made on Lauteraargletscher, Switzerland, in 2001, that showed diurnal fluctuations in both surface velocity and englacial vertical strain. We find that spatial
gradients in basal lubrication can set up variations in the deviatoric stress field that
increases with distance to the bed and has a maximum value near the glacier surface.
This stress field produces a significant strain rate near the surface. The temporal evolution of a slippery zone is identified as a possible cause of the observed diurnal variations in the vertical strain rate. Although general aspects of the measurements can be
explained in this way, the calculated vertical strain rates are too small, suggesting thatthe modeled effective viscosity values using Glen’s flow law are too large near thesurface.

INTRODUCTIONMETHODS

Temperate glaciers are known to change their flow velocity
over a time-scale of a week, a day and in some cases hours(Willis, 1995). Such relatively short-term fluctuations areunderstood to be controlled by changes in subglacial conditions, particularly by subglacial water-pressure fluctuations(Iken and Bindschadler, 1986; Jansson, 1995). The transmission of the basal variability in flow speed to the surface iscomplex and at first sight counterintuitive (Balise and Raymond, 1985), and the influence of temporal changes in basal
lubrication on surface velocities and glacial deformation difficult to predict in detail.
In this study, we use a numerical finite-element model
to determine how a non-linear viscous medium (ice) moving over a bed with spatially non-uniform basal lubrication causes flow perturbations. The aim of this work is
not to comprehensively study the general aspect of thetransmission of basal slipperiness perturbations to thesurface of glaciers, but rather to find basal conditions thatcan reproduce the particular diurnal fluctuations in surface flow velocity and vertical strain observed during theablation season in 2001 on Lauteraargletscher, Bernese
Alps, Switzerland.

Field equations*

Present address: Versuchsanstalt fu«r Wasserbau, Hydrologie und Glaziologie, Eidgeno«ssische Technische
Hochschule, ETH-Zentrum, CH-8092 Zu«rich, Switzerland.

The problem to be solved is the two-dimensional flow ofincompressible viscous material governed by Stokes equations:r¢ u ˆ 0;

rp ¡ ·¢u ˆ »g ;

…1†
…2†

where u ˆ …u; w† is the velocity vector, p is the hydrostaticpressure, · is the viscosity, » is the density of ice, andg ˆ …gx ; gz † is the gravity vector. We employ parallel-sided
slab geometry and use a Cartesian coordinate system withthe x axis along the bed pointing down-glacier, and the zaxis normal to the bed pointing upward. The constitutive
relationship used is
1 n…3†"_ ij ˆ
½ ˆ A½en¡1 ½ij ;2· ij
where "_ ij and ½ij are the components of the strain-rate anddeviatoric-stress tensors, respectively, and ½e is the effective
stress. The rate factor A and the flow-law exponent n arematerial parameters. We use the commonly accepted valueof n ˆ 3 for the flow-law exponent, and use a rate factor ofA ˆ 10 MPa ^3 a ^1, which gives modeled surface velocities
similar to those observed. Taking into account the reduction
in shear stress due to the valley shape by using a shape factorof 0.5 (Nye, 1965), this value of the rate factor compares
favorably with 75 MPa ^3 a ^1 which was used by
Gudmundsson (1999) in a three-dimensional numericalmodeling of Unteraargletscher, Bernese Alps, Switzerland.
49

Sugiyama and others: Basal lubrication and englacial strain-rate distribution
and midpoint, and the maximum resolution of the velocity
field is 250 m longitudinally and 25 m vertically.
The discretization error was estimated by comparing
numerical results with the analytical solution for plane-slab
flow for no slip. The largest error of 1.2% occurred in thehorizontal surface velocity. By comparing the numericaland the analytical vertical velocity profiles, we found thatthe main source of error was near the bed, where the largest
spatial gradients exist.
Boundary conditions
Fig. 1. Finite-element mesh used in this study. (a) All elements
and coordinates used for the longitudinal cross-section modelof a conceptional glacier with its head at x ˆ 0 m. (b)Enlargement of the mesh near the glacier bed. Shaded region
indicates the 1m thick subglacial thin layer used to introduce
basal sliding.
Numerical scheme and model description
Equations (1^3) were solved with a finite-element methodthat was coded for this study. Figure 1a shows the finiteelement mesh used in the computation. The geometry of the
6 km central part of the finite-element grid is a simplified
longitudinal cross-section of Lauteraargletscher. We assumed
an ice thickness of 400 m and a surface slope of 4³ (Funk andothers,1995; Bauder, 2001). To minimize boundary effects, thenumerical grid was extended an additional 6 km upstream
and 6 km downstream from this central section. We used
Galerkin’s method with quadratic shape functions (Zienkiewicz,1977) to compute the velocity and pressure fields. Starting from the solution of a linearly viscous flow, a new effective
viscosity distribution was calculated from the previously
determined strain-rate field, and the calculation repeateduntil the horizontal velocity field converged within 10^5 m a^1.
The computation gives the velocity components at each node

Basal motion is simulated by introducing a thin sub-basal
layer within the central part of the finite-element model.
The viscosity of this layer can be regarded as the effective
viscosity of subglacial till, or, which in the case of Lauteraargletscher is more appropriate, as being related to the form
drag generated by sub-grid bed roughness. Spatially varying basal slipperiness is generated by giving different viscosity values to some of the elements of the sub-basal layer. A
no-slip boundary condition is applied for the bottom nodesof the sub-basal layer (Fig. 1b).
Ignoring the effects of longitudinal stresses within the
sub-basal layer, the deformational velocity of the layer isrelated to the viscosity of the sub-basal layer, its thickness
and the shear stress at the top of the layer. This gives rise toa sliding law of the form
00

ub ˆ 2dA0 ½bn ˆ c…x†½bn ;where ub is the basal sliding velocity, d is the layer thickness,½b is the basal shear traction, and A 0 and n0 are the parameters of the layer. We assume n0 ˆ 1. This method of introducing basal motion through a sub-basal layer has beenused before (e.g.Vieli and others, 2000). Because longitudinal
stress gradients will cause the thickness of the layer tochange with time, this method is not suitable for transient
calculations. Here we are interested in obtaining snapshotsof the velocity field for a given glacier geometry, so we do

Fig. 2. Computed fields of (a) horizontal velocity, (b) vertical velocity, (c) longitudinal deviatoric stress component ½xx, and(d) vertical strain rate. Only the central domain of the modeled field (5000 < x < 11 000 m) is shown. Contour lines were
established using the method of triangulation with linear interpolation.
50…4†

Sugiyama and others: Basal lubrication and englacial strain-rate distribution
not allow the surface or other boundaries within the finiteelement mesh to evolve with time.
Stress-free conditions are imposed at the surface, andatmospheric pressure is neglected (½ij nj ˆ 0, where nj is asurface-normal vector component). The lower end of themodeled domain (x ˆ16 000 m) was connected to the upperend (x ˆ 0 m) so that the glacier is assumed to be infinitely
long (periodical boundary condition).
RESULTS AND DISCUSSION
Transmission of basal sliding to the surface
Figure 2 shows model results for which the rate factors of theelements in the sub-basal layer are set to 25 MPa ^1 a ^1 for aslippery zone within 7000 < x < 9000 m, and 5610^4 MPa^1a^1 elsewhere. Introducing this slippery zone in the modelincreased the surface horizontal velocities directly above
the center of the zone by about 30% (Fig. 2a).With increasing distance up- and down-glacier, the amplitude of thevelocity perturbation decreases. For example, five ice thicknesses away from the boundaries of the slippery zone (x ˆ5000 and 11000 m), the increase in surface velocities is down
to 8%. The calculated vertical velocity field shows
extending and compressive flows above the boundaries(Fig. 2b). Furthermore, tensile and compressive deviatoricstresses occur throughout the entire glacier but are strongest
near the surface (Fig. 2c). Because the deviatoric stress isrelated to the vertical strain rate by Equations (1) and (3),this stress field generates vertical strain rate as shown in Figure 2d. The computed vertical strain rate near the surface
helps to explain the observed vertical strain fluctuations discussed in the next subsection. Analytical solutions for smallamplitude slipperiness perturbations show a similar increase in deviatoric stress and strain with distance awayfrom the source of the disturbance (i.e. the slippery zone)(Vonmoos,1999).
Figure 3a shows the longitudinal distribution of basal
and surface horizontal velocities calculated from the samebasal condition as used in Figure 2 (A0 ˆ 25 MPa^1 a^1 at
7000 < x < 9000 m).We define the transmission rate of horizontal speed as
¢us
…5†
fu ˆ
;
ub
where ¢us is the difference in horizontal surface velocity
between non-slip and slip conditions and ub is the basal
velocity at the same position (Fig. 3a). Calculations were performed for various sets of viscosity of the sub-basal layer anddifferent lengths of the slippery zone (L). Results are summarized in Figure 3b, which shows the horizontal surface flow
speed us and transmission rate of horizontal velocity fu at
x ˆ 8000 m as functions of L and the basal velocity ub. As
expected, the surface velocity increases with L.The transmission rate of horizontal speed fu also increases with L, andbecomes larger than 0.75 for ub > 5 m a^1 when L ˆ 6000 m.
In general, fu varies with L similar to an analysis for a linear
medium (Balise and Raymond,1985), although our calculated
values of fu are slightly smaller. This discrepancy might becaused by an effective material softening close to the glacier
bed that introduces non-linearity into the flow law. Blatterand others (1998) used the boundary condition of zero basal
shear traction (free slip) in their numerical model, which
resulted in much higher transmission rates than ours.

Fig. 3. (a) Longitudinal distribution of the horizontal velocity
obtained for the sliding condition described in the text. Surface(solid circle) and basal velocities (open circle) were computed
at nodes and midpoints of the finite-element mesh. Solid line
indicates the surface velocity for the non-slip condition. (b)Horizontal surface velocities at x ˆ 8000 plotted against thebasal velocity ub at x ˆ 8000 and the length of the slipperyzone L. The surface velocity is indicated for each experimentwith a unit of m a^1 and contoured with solid lines.The shaded
contours show the transmission rate of horizontal velocity fu.
In order to examine the transfer of a velocity gradientfrom the bed to the surface, we assumed the rate factor of
the sub-basal layer to vary linearly with x within a regionof length L in further numerical experiments. Figure 4ashows the resulting surface and basal horizontal velocity
distribution for A0 given by
A0 …x† ˆ ¡0:01x ‡ 105 …MPa¡1 a¡1 †5500 < x < 10 500 …m† :…6†Both the longitudinal gradients of basal and surface
velocities were calculated at x ˆ 8000 m, as well as thetransmission rate of horizontal velocity gradient defined as
fgrad ˆdus dub
=
:
dx dx

…7†

Figure 4b shows a contour plot of the surface velocity gradient dus =dx and fgrad for a range of lengths L and basal
velocity gradients dub =dx. The surface velocity gradientincreases as the basal gradient and L increase, but the transmission rate is rather small. In this case, the maximum valueof fgrad is about 0.5, and almost no transmission occurs for L
5 3000 m.
51

Sugiyama and others: Basal lubrication and englacial strain-rate distribution
Furthermore:(a) The surface velocity and vertical strain measurements
were made within the ablation area with the snowline
about 1km upstream.
(b) Subglacial water pressure at the study site fluctuateddiurnally throughout the summer, with an averagepressure decreasing from July to August.
(c) Surface velocity varied diurnally at least as far as§1.5 km up- and downstream from the study site, andthe amplitude of the diurnal variation increased in the
upstream direction.

Fig. 4. (a) Longitudinal distribution of the surface (solid
circle) and basal velocities (open circle) obtained for the sliding
condition given by Equation (5). The solid line is the surface
velocity for the non-slip condition. (b) Longitudinal gradientsof the horizontal surface velocity at x ˆ 8000 plotted against thebasal velocity gradient dub=dx at x ˆ 8000 and the length of
the slippery zone L.The surface velocity gradient is denoted foreach experiment with a unit of 10^3 a^1 and contoured with solidlines.The shaded contours are the transmission rates of horizontal velocity gradient fgrad.
Diurnal fluctuation in the surface flow speed andvertical strain
We now compare our modeling predictions to measurements of diurnal velocity and strain fluctuations made on
Lauteraargletscher from 22 to 27 August 2001. Lauteraargletscher, a tributary of Unteraargletscher in the Bernese
Alps of Switzerland, is a temperate valley glacier that isabout 5 km long and 1km wide. Details of the measurements
are reported by Sugiyama and Gudmundsson (2003), but
the main findings can be summarized as follows.
(A) Surface flow velocity at the center of the glacier varied
diurnally, with an approximate minimum speed of 30
m a^1 in the morning and a maximum of 60 m a^1 in theafternoon.
(B) From 24 to 26 August, diurnal fluctuations in verticalstrain occurred in a 174 m deep borehole where thetotal ice thickness was about 400 m. The depth of theborehole increased by about 30 mm from morning toevening, and then decreased again by a similar amounttowards the next morning. This vertical stretching/compression corresponds to a strain of 1.7610^4.
52

If we assume no transverse component in the flow
velocity, observation (B) suggests a strongly compressive flowat the study site in the daytime. This compressive flow could
be due to the comparatively larger increase in flow velocities
in the upstream direction in the early morning (c). The factthat the surface flow velocity was well correlated with thesubglacial water pressure during the study period suggeststhat the surface velocity variation can be attributed to thediurnal water-pressure variations. The decrease in the average water pressure in the ablation season (b) seems to berelated to the development of the subglacial drainage systembeneath the study site. Because the drainage system starts todevelop in the lower reach of the glacier and then progresses
upward during the ablation season, water pressure fluctuateswith larger amplitude and causes greater basal sliding in theupper reach. Therefore, we ascribe the observed verticalstrain to the longitudinal gradient of basal velocity.
From Figure 4b, one can determine the basal velocity
gradient necessary to produce the observed vertical strain.
However, the observed strain rate was about 0.1a^1, which is
too large to be generated by the velocity gradients in Figure
4b. This suggests that the perturbation of slipperiness waseven larger than the one used to produce Figure 4b.
The observed diurnal variation in flow speed and its connection to basal pressure variations suggests that the glacierbecame partly decoupled from the bed during daytime. Thismay be due to increased basal cavity formation or due to thefailure of a thin till layer. If such a decoupled zone begins todevelop upstream from the measurement site and thenincreases in spatial extent and sliding magnitude with moving down-glacier in the course of the day, large deviatoricstresses and vertical strain rates near the surface can resultthat are consistent with the results in Figure 2c and d. Sucha basal condition can be simulated by prescribing a slipperyzone as described in the previous section. The sliding magnitude, spatial size and position of the slippery zone were varied
to fit the observed velocity and strain rate at the study site aswell as possible. Because the observed diurnal changes inthickness are only a few centimeters at most, changes in thickness and surface slopes can be ignored in the modeling calculations, and the diurnal variations can be analyzed withoutchanging the model geometry.
We introduced a slippery zone upstream of the observation site by changing the rate factor of the sub-basal layer
as explained above. By subsequently changing the magnitude, size and position of the slippery zone, the observed
spatial and temporal pattern in vertical strain rate could be
reproduced qualitatively. Figure 5a and b show surface andbasal velocity distributions calculated for the determinedslippery zones. The conditions of the slippery zone at eachtime were determined by adjusting the distribution of A0 to

Sugiyama and others: Basal lubrication and englacial strain-rate distribution

Fig. 5. Diurnal evolutions of (a) surface and (b) basal velocity
distributions computed with the basal condition scenario
described in the text.
fit the observed surface velocity at the study site. The temporal evolution of the slippery zone is explained as follows:
The slippery zone starts to develop upstream of the observation site at about 0800 h when the surface melting starts.
Then it increases in extent and magnitude, and forms a sliding boundary beneath the study site at 1500 h. In the evening, when the surface-melting rate decreases, the basal
velocity starts to decrease and the zone moves down-glacierduring the period 1800^0000 h. Eventually, this downglacier propagation of the slippery zone replaces the compressive flow at the study site with a weaker extending flow
(0300^0600h).

Fig. 6. (a) Diurnal variations of (a) horizontal surface
velocity and (b) vertical strain rate. Open circles with thedashed line are model results, and the solid lines are measurement data taken on Unteraargletscher on 25 August 2001.Vertical strain rate is computed as the mean for the upper 200 m ofthe glacier, and measurements were taken from a 174 m deep
borehole and plotted with the scale magnified 10 times.
Figure 6 compares the field measurements with the computed temporal variations of horizontal surface velocity andvertical strain rate in the upper 200 m of the ice at the studysite (x ˆ 8000 m). The model reproduces the observed temporal variation in velocity and vertical strain on 25 August
2001. Even though the computed strain rates are still considerably smaller than the measured ones, they are an orderof magnitude larger than those obtained in Figure 4b.
The fact that the computed strain rate accounts for only
about 10% of the measured borehole-depth changes showsthat some important aspects of the flow are not described

Fig. 7. Computed fields of (a) horizontal velocity, (b) vertical velocity, (c) longitudinal deviatoric stress component ½xx, and (d)vertical strain rate. Modified flow law (Equation (7)) was used with the basal conditions identical to those in Figure 2.
53

Sugiyama and others: Basal lubrication and englacial strain-rate distribution
adequately by the model. To obtain larger strain rates, themodel could assume higher flow speeds in the upper region,
but this assumption would contradict observed surface
velocities. Another possibility is that the flow law used isincorrect. In the present model, stress concentration near
the surface in Figure 2c does not generate clear strain-rate
concentration in Figure 2d. This is because Glen’s flow law(Equation (3)) gives much higher effective viscosity as theeffective stress decreases near the glacier surface. To assess
the influence of effective viscosity near the surface, anotherexperiment was carried out with modified Glen’s flow law,"_ ij ˆ A…½en¡1 ‡ ½0n¡1 †½ij :…8†Here, as elsewhere, a constant value ½0 has been introducedin the flow law to avoid mathematical singularity where theeffective viscosity is infinite when the effective stress vanishes(Hutter, 1983). We took ½0 ˆ 0.1MPa so that ice becomes lessviscous under a low-stress condition, and applied the basal
conditions identical to those used for Figure 2. Computed
flow velocities, stress and strain-rate fields are shown in Figure 7 and compared with Figure 2. It is interesting to see thevertical strain increases near the surface while the deviatoricstress decreases (Fig. 7c and d). The results suggest that themodeled vertical strain will be improved by Equation (7)but will still be insufficient due to the reduction in the stress.
Further study on the flow-law exponent or flow law itself isrequired to interpret the field data adequately. The model
can also be improved by including the average longitudinalstrain and transverse stresses over the area in question.
The studied velocity variation accompanied by a strainrate anomaly is similar in nature to the ``mini-surges’’ orpropagation of enhanced motion waves, observed in Variegated Glacier, Alaska, U.S.A. (Raymond and Malone, 1986;
Kamb and Engelhardt, 1987). Balise and Raymond (1985)
applied their analysis of a linear viscous flow to interpret themeasurements on the surface by a propagation of basal
velocity anomaly that was much sharper than the one we used
for Figure 5b. They also suggested that the vertical strainingof ice might explain uplift events, which are increases in surface elevation during a fast-flow period. In order to apply ourmodel to uplift events, quantitative discussion described aboveis crucial because the uplift rates measured atVariegated Glacier (Kamb and Engelhardt,1987), or at Unteraargletscher by
Iken and others (1983), are in the same order of ¹10^4 day^1
that is higher than our modeling result in Figure 6b.
SUMMARY AND CONCLUSION
We used a finite-element model with Glen’s flow law to study
the effects of a local perturbation in basal slipperiness on thesurface velocity and the englacial strain-rate distribution.
The model results were compared to recent measurementsof diurnal strain-rate variations on Lauteraargletscher.
Our model of down-glacier propagation of a slippery zone
was able to qualitatively explain the observed temporal variation in flow velocities and vertical strain rates. However,quantitative comparison between model predictions and fieldmeasurements showed that calculated vertical strain rates
were about an order of magnitude too small. The reason for

54

this discrepancy remains unclear. A possible explanation isthat the modeled ice near the surface is too stiff, which may
arise from ignoring the average longitudinal strain or
because the flow-law exponent of ice varies with stress.

ACKNOWLEDGEMENTS
We are grateful to R. Naruse of the Institute of Low Temperature Science, Hokkaido University, for his support forthis study. All the members of Section of Glaciology, VAWETHZ, Zu«rich, assisted us in the fieldwork at the glacier.We
also thank L. Tarasov, H. Blatter, R. C. A. Hindmarch, D.
Vaughan and an anonymous reviewer for their beneficial
comments on the manuscript. The Swiss National Science
Foundation through grant No. 2100-063770.00 and the Inoue
Scientific Field Study Foundation funded this research.

REFERENCES
Balise, M. J. and C. F. Raymond.1985.Transfer of basal sliding variations tothe surface of a linearly viscous glacier. J. Glaciol., 31(109), 308^318.
Bauder, A. 2001. Bestimmung der Massenbilanz von Gletschern mit Fernerkundungsmethoden und Fliessmodellierungen: eine Sensitivita«tsstudie
auf dem Unteraargletscher. Eidg. Tech. Hochschule, Zu«rich. Versuchsanst.
Wasserbau, Hydrol. Glaziol. Mitt. 169.
Blatter, H., G. K. C. Clarke andJ. Colinge.1998. Stress and velocity fields inglaciers: Part II. Sliding and basal stress distribution. J. Glaciol., 44(148),457^466.
Funk, M., G. H. Gudmundsson and F. Hermann.1995. Geometry of the glacier
bed of the Unteraarglacier, Bernese Alps, Switzerland. Z. Gletscherkd.
Glazialgeol., 30,1994,187^194.
Gudmundsson, G. H. 1999. A three-dimensional numerical model of theconfluence area of Unteraargletscher, Bernese Alps, Switzerland. J. Glaciol.,45(150), 219^230.
Hutter, K. 1983.Theoretical glaciology; material science of ice and the mechanics ofglaciers and ice sheets. Dordrecht, etc., D. Reidel Publishing Co.; Tokyo,
Terra Scientific Publishing Co.
Iken, A. and R. A. Bindschadler. 1986. Combined measurements of subglacial water pressure and surface velocity of Findelengletscher, Switzerland: conclusions about drainage system and sliding mechanism. J.
Glaciol., 32(110),101^119.
Iken, A., H. Ro«thlisberger, A. Flotron and W. Haeberli. 1983. The uplift of
Unteraargletscher at the beginning of the melt season ö a consequenceof water storage at the bed? J. Glaciol., 29(101), 28^47.
Jansson, P. 1995. Water pressure and basal sliding on Storglacia«ren, northern
Sweden. J. Glaciol., 41(138), 232^240.
Kamb, B. and H. Engelhardt.1987.Waves of accelerated motion in a glacierapproaching surge: the mini-surges of Variegated Glacier, Alaska,
U.S.A. J. Glaciol., 33(113), 27^46.
Nye, J. F. 1965. The flow of a glacier in a channel of rectangular, elliptic orparabolic cross-section. J. Glaciol., 5(41), 661^690.
Raymond, C. F. and S. Malone.1986. Propagating strain anomalies duringmini-surges of Variegated Glacier, Alaska, U.S.A. J. Glaciol., 32(111),
178^191.
Sugiyama, S. and G. H. Gudmundsson. 2003. Diurnal variations in verticalstrain observed in a temperate valley glacier. Geophys. Res. Lett., 30(2),1090. (10.1029/2002GL016160.)
Vieli, A., M. Funk and H. Blatter. 2000. Tidewater glaciers: frontal flow
acceleration and basal sliding. Ann. Glaciol., 31, 217^221.
Vonmoos, M. 1999. Auswirkungen basaler Sto«rungen auf das Geschwindigkeitsfeld und die Oberfla«che eines Gletschers. (Diplomarbeit, Versuchsanstalt fu« rWasserbau, Hydrologie und Glaziologie, ETH, Zu«rich.)
Willis, I. C. 1995. Intra-annual variations in glacier motion: a review. Prog.
Phys. Geogr., 19(1), 61^106.
Zienkiewicz, O. C. 1977. The finite element method.Third edition. London, etc.,
McGraw-Hill Book Co.