<< . .

. 29
( : 51)



. . >>





ice thickness
1 bar

300“400 m
200“300
100“200
0“100

boundary conditions are not required is that the margins. We have assumed a uniform threshold
algorithm refers only to downslope grid points. shear stress of 1 bar for this reconstruction.
This is useful for applications in ice-marginal ar- In the alpine environment, the bed-slope com-
eas because only a portion of the ice sheet need ponent of the basal shear stress is dominant. As
be considered. such, variations in ice thickness are small and are
One of the bene¬ts of the threshold-sliding associated with areas of bed curvature: thicken-
model is that it can be broadly applied to the re- ing occurs in areas of bed concavity and thinning
construction of any ice sheet or glacier in the con- occurs in areas of bed convexity. Rapid changes
tinuum from large ice sheets to alpine glaciers. in glacier thickness are, therefore, restricted to
As a robust and broadly applicable model, it can slope breaks and glacier snouts.
provide preliminary reconstructions for any area Figure 6.19b presents the results of our re-
and requires the constraint of just a single input construction based on the input data in Figure
variable, the basal shear stress. 6.19a. Ice is up to several hundred meters thick
As an example of the model application in glacial valleys and is generally thicker where
in an alpine environment, we will recon- bed slopes are gentler, except for areas close to
struct the Wisconsin glaciation of the central glacier snouts where thicknesses taper. Modify-
Brooks Range of Alaska using input data com- ing the assumed basal shear stress from 1 bar
piled by the Alaska PaleoGlacier Group (http:// to a greater or lesser value leads to thicker or
instaar.colorado.edu/QGISL/data intro.html). This thinner glaciers but does not signi¬cantly change
data includes a 100-m resolution DEM of bed to- the relative thicknesses illustrated in Figure
pography and a vector dataset of Wisconsin ice 6.19b.
6.5 THRUST SHEET MECHANICS 147


In this reconstruction, steep areas in the high- ground surface
(a)
lands of the Brooks Range are predicted to have
thin ice, as indicated by the dark tones in the
L
grayscale map in Figure 6.19b. However, it is likely
that below some minimum thickness, ice cannot
form and locally the area will be ice free. decollement
For example, in steep terrain the snow depth
H
may be insuf¬cient to insulate deeper layers from
summer melting or it could be too thin to initi-
ate compaction. To account for this, the alpine re- (b)
s
constructions could be post-processed to remove H
ice from areas that have insuf¬cient thickness q
t=0
for ice formation. In such cases, the reconstruc-
tion should then be recomputed with the revised
mask to provide an improved reconstruction.



t>0
6.5 Thrust sheet mechanics
Fig 6.20 Schematic diagram of the 2D thrust sheet model.
Deformation within Earth™s crust takes place by
brittle fracture and slip along faults and joints. the shape of the fold depend on the yield stress
As such, it is not appropriate, in detail, to model and fault geometry? In this example we prescribe
the crust as a ¬‚uid. At large scales, however, plas- the tectonic geometry and slip rate, but it is also
tic ¬‚ow models do a surprisingly good job at re- possible to prescribe the stress state along the
producing the observed geometries of structural fault and allow the slip rate to respond to both
folds. The reason why plastic ¬‚ow models are suc- the tectonic stress and the overburden stress, al-
cessful is that they represent the collective ef- lowing the fault to slow down as overburden ac-
fects of slip along a complex system of faults and cumulates.
joints. As a rock is stressed, it deforms very little The thickness of rock above the fault surface
until the yield stress of the rock is exceeded and is governed by Eq. (6.46) with a uniform yield
a new fault or joint is formed. Subsequent defor- stress:
mation is accommodated by slip along faults and
„0
h|∇h + ∇b| =
joints when the threshold frictional force is ex- (6.49)
ρg
ceeded. Rock fracture and slip along a network
of faults are both threshold-controlled motions. Equation (6.49) can be solved with the sandpile
As such, a perfectly plastic or threshold-sliding algorithm using a very similar approach to the
model is often a good model for the large-scale one used for modeling 2D ice sheets over com-
deformation of the crust. plex topography. The advance of the fault tip is
In this section we consider a simple exam- modeled by adding mass to the system directly
ple of thrust-sheet displacement using a per- above the fault tip at a rate equal to the vertical
fectly plastic model. We assume the existence of component of slip. Appendix 5 provides a C code
a crustal ramp of angle θ and width L separating that implements this approach.
the ground surface from a decollement surface at Figure 6.21 illustrates the fold evolution for
a case with „0 = 20 bars, ρ = 2700 kg/m3 , H =
depth H (Figure 6.20). Motion along the thrust
3 km, L = 300 km, θ = 2—¦ , and s = 1 cm/yr at
fault takes place with slip rate s . The key ques-
three time increments t = 10 Myr, 20 Myr, and
tions we will ask are: (1) what is the 2D shape of
the fold above the thrust fault, and (2) how does 30 Myr following the initiation of thrusting. The
148 NON-NEWTONIAN FLOW EQUATIONS



(a)
(a) 3
4
100 km 100 km
3
h2
h
(km) 2 (km) 1
1
linear
0
forelimb
0
edge
(curvature
increases
through
time)

t = 10 Myr t0 = 5 bar
(b)
4 3
(b)
3
h2
h
(km) 2 (km) 1
1
0
0




t = 20 Myr L = 150 km
3
4
(c) (c)
3
h2
h
(km) 2 (km) 1
1
0
0




t = 30 Myr q = 1° 100 200
0 300
200 400
0
x (km) x (km)

Fig 6.21 Shaded-relief images and cross-sectional plots of Fig 6.22 Sensitivity of the thrust sheet system to variations
in the (a) yield stress („0 = 5 bars compared to 20 bars for
the fold topography developed over a thrust fault for a
the reference case), (b) fault width (L = 150 km instead of
perfectly-plastic model of crustal deformation, assuming
„0 = 20 bar, ρ = 2700 kg/m3 , H = 3 km, L = 300 km, 300 km), and (c) fault angle (θ = 1—¦ ).
θ = 2—¦ , and s = 1 cm/yr at three time increments t = 10 Myr,
20 Myr, and 30 Myr following the initiation of thrusting.
total shortening equals the width of the thrust
sheet does the fault achieve a true steady-state
maximum structural relief in the model has also taper.
been restricted to be 2H (i.e. the fault becomes The sensitivity of fold geometry to model pa-
horizontal, forming a large syncline, once the rameters is shown in Figure 6.22. Each panel
fault tip has achieved an elevation equal to H shows the shaded relief and cross section of the
fold at t = 10 Myr for a case with (a) a low yield
above the ground surface). The shaded-relief im-
stress („0 = 5 bars), (b) a low yield stress („0 =
ages show the topography of the fold (given by
5 bars) and a narrow fault (L = 150 km), and (c)
the sum of the fault-surface elevation b and the
a low yield stress („0 = 5 bars) and thrust angle
rock thickness h). Early in the simulation, fault
(θ = 1—¦ ). Lower yield stresses lead to thinner rock,
displacement creates an asymmetric fold charac-
terized by a steeper backlimb than frontlimb. The as expected. Out in front of the fault, this thin-
forelimb progrades out in front of the fault to a ning has the effect of creating a broader zone of
maximum length controlled by the structural re- deformation and hence a higher degree of oro-
lief and the yield stress. Early in the simulation, genic curvature for faults with otherwise equal
the forelimb has a linear planform shape. Over geometries and shortening amounts.
time, the curvature migrates from the edge of the To date, most ¬‚uvial landform evolution
orogen into the center. Only when the amount of models have assumed relatively simple tectonic
6.6 GLACIAL EROSION BENEATH ICE SHEETS 149


forcings. When more realistic tectonic forcings are calculated based on the reconstructed ice-
have been considered (e.g. Willett et al., 2000), surface topography. Our approach is to model the
the effects on mountain-belt geometry have been geometry and ¬‚ow in ice sheets and the resultant
shown to be very important. Although simpli¬ed, subglacial erosion iteratively through time. Fig-
the perfectly plastic model of topographic evolu- ure 6.23 illustrates the basic steps of the model
tion above a crustal-scale fault provides a start- and the input and output data at each step us-
ing point to quantify rock uplift rates above pre- ing the New York State Finger Lakes as the ex-
scribed structural elements. Coupling this model ample application. In this example, we use the
with bedrock erosion processes could help to fur- modern Finger Lakes topography as input to de-
ther unravel the relationships between tectonics termine the bedrock-erosion pattern that would
and erosion in thrust systems. occur if an ice sheet were to advance over the
area today. Input data for step 1 of the model
includes (1) the bedrock topography, (2) the po-
sition of the ice margin, and (3) the thresh-
6.6 Glacial erosion beneath ice
old basal shear stress for motion. These inputs
sheets are used to reconstruct the geometry of the ice
sheet using the threshold-sliding model and the
Bedrock erosion beneath ice sheets is the domi- sandpile method. Following the ice-surface recon-
nant agent of Quaternary erosion and landform struction, the lithospheric de¬‚ection under the
evolution in most of the Northern Hemisphere. ice load must be considered if the study area
Despite the dominance of subglacial erosion in is of large extent. If the study area is less than
≈ 100 km in width, the lithosphere is assumed
many parts of the world, our understanding of
subglacial erosional landforms and the feedback to be perfectly rigid. If the area is signi¬cantly
processes between ice ¬‚ow, bedrock erosion, and larger than 100 km the de¬‚ection of the litho-
bedrock topography beneath ice sheets is limited. sphere should be determined using the 2D ¬‚ex-
For example, what processes and feedbacks con- ure equation and estimates for the local elastic
trol the morphology of overdeepenings beneath thickness of the crust. If de¬‚ection is computed,
glaciers and ice sheets? How do these processes the model returns to step 1 where the ice-surface
differ on the margins of ice sheets from their in- topography is recomputed taking into account
teriors? What processes control the size and dis- the lithospheric de¬‚ection in the bed topogra-
tribution of regularly spaced glacial lakes such phy. Although 100 km is given as a rule-of-thumb,
as the New York State Finger Lakes? Although the precise threshold length scale where ¬‚exu-
considerable work has focused on modeling land- ral compensation becomes signi¬cant depends
form evolution in the alpine-glacial environment on the local elastic thickness.
(e.g. Harbor, 1992; Braun et al., 1999; MacGregor In step 2, basal-¬‚ow velocities are determined
et al., 2000), previous work has just touched on with the balance-velocity method of Budd and
landform evolution beneath ice sheets. In ad- Warner (1996). The balance-velocity method re-
dition to providing a better understanding of quires that the pattern of accumulation and abla-
glacial-landform evolution, determining the con- tion be speci¬ed. For many paleo ice sheets, this
trols on glacial-lake formation could help to con- pattern may not be well constrained. For these
strain paleo ice-sheet thickness, a critical param- cases, the accumulation/ablation pattern may be
eter that is dif¬cult to determine using exist- modeled using an equilibrium line altitude (ELA)
ing methods such as relative sea-level inversion and an accumulation/ablation gradient. For ex-
(Peltier, 1994). ample, assuming an ELA of 1 km and an accumu-
In this example application we use the lation gradient of 0.1 m/yr per km, the ablation
rate at an elevation of 500 m is ’0.5 m/yr. In step
threshold-sliding model as the ¬rst step in a
glacial-landform evolution model for erosion ice 3, the basal velocities are used to determine the
sheets. Following the ice-surface reconstruction, bedrock-erosion rate using Hallet™s model. Hal-
basal-sliding velocities and bedrock erosion rates let™s model estimates erosion rates as a power
150 NON-NEWTONIAN FLOW EQUATIONS



(a) (c)




(b)



1 km
is area of study
> 100 km?




30 km


(d)


0.5 bars




0.5 bars



Fig 6.23 A ¬‚ow chart illustrating the model inputs and outputs with erosion of the modern New York State Finger Lakes as the
example application. (a) Inputs to the model include the initial bedrock topography (from US Geological Survey DEMs), the
ice-margin position (given by the position of the Valley Heads Moraine (New York State Geological Survey, 1999)), and the
threshold basal shear stress (assumed to be 0.5 bars). (b) In step 1 these inputs are used to reconstruct the ice-surface
topography. Following this step, the de¬‚ection of the lithosphere under the ice load should be determined if the study area is
equal to or larger than 100 km in extent. For this example we have neglected ¬‚exural subsidence. (c) In step 2 the ice-surface
topography and pattern of accumulation and ablation are used with the balance-velocity method of Budd and Warner (1996) to
determine the basal-sliding velocities. (d) In step 3, Hallet™s model, which relates the bedrock erosion rate to a power function of
the sliding velocity, is used to determine the spatial pattern of erosion. The model then returns to step 1 where the erosion
pattern is used to determine a new bed topography for the next iteration.




function of the basal-sliding velocity. In step 4, Any study of glacial erosion must contend
the pattern of erosion rate is used, together with with two complicating factors. First, the initial,
the model time step, to modify the bedrock to- preglacial surface is generally unknown. This
pography according to b = E t, where b is the complication applies to individual glacial ad-
bedrock topography, E is the erosion rate, and t vances as well as to erosion during the full dura-
is the time step. The new bedrock topography, re- tion of the Quaternary. As a speci¬c example, it
¬‚ecting the bedrock scour during the time step is commonly assumed that the ice ¬‚ow respon-
t, is then used as input for the next model it- sible for the carving of the New York State Fin-
eration, which begins back at step 1. ger Lakes took advantage of preexisting ¬‚uvial
6.6 GLACIAL EROSION BENEATH ICE SHEETS 151


valleys of the Allegheny Plateau formed by north- the abrading particles, embedded in the ice, is
¬‚owing rivers (von Engeln, 1965). It is dif¬cult the most important controlling factor for ero-
to test this hypothesis, however, since we have sion. Hallet™s model can be summarized as
no constraints on the depth of pre-glacial ¬‚u-
‚b
= ’c V n (6.50)
vial valleys. If valleys were deep, it may still be
‚t
possible that ice ¬‚ow and bedrock erosion were
where ‚b/‚t is the erosion rate, V is the basal
controlled principally by the ice-surface topogra-
sliding velocity, and n and c are empirical co-
phy of the Ontario Lobe rather than by the un-
ef¬cients. n may taken to be equal to 2 (e.g.
derlying bed topography. Numerical modeling of
MacGregor et al., 2000), although this value is
glacial erosion is ideal for coping with this un-
not well-constrained empirically. The value of
certainty because multiple working hypotheses
c may be constrained using sediment-¬‚ux esti-
can be tested and evaluated for their impact on
mates from modern glaciers (e.g. Hallet et al.,
model results. For example, a suite of initial to-
1996). In this chapter we use Hallet™s model, re-
pographies can be used to determine a suite of
¬‚ecting the consensus that sliding velocity is the
corresponding model results to assess the sensi-
dominant control on bedrock erosion rates. Nu-
tivity of the model results to the initial topogra-
merical modeling in glacial geomorphology is
phy. Results that are insensitive to variations in
still in its infancy. Three of the most important
initial topography can be considered to be robust.
contributions have been made by Harbor (1992),
A second complicating factor in studying ero-
MacGregor et al. (2000), and Braun and his col-
sion beneath ice sheets is the dif¬culty of separat-
leagues (Braun et al., 1999; Tomkin and Braun,
ing variations in the intensity of glacial erosion
2002). All of these contributions focused on the
from variations in the erodibility of bedrock. In
alpine glacial environment, but they represent
some cases, structural control is easy to recog-
the previous work most similar to the work in
nize. For example, the topography of the Finger
this chapter. Harbor™s model illustrated how U-
Lakes Region is dominated by a weathering scarp
shaped glacial valleys may form from pre-existing
separating the ¬ne-grained shale of the Ontario
¬‚uvial valleys. However, because he assumed a
Lowlands from the coarser-grained sandstone and
simple form for the glacier (i.e. a constant-area
shale of the Allegheny Plateau. In other areas,
cross-sectional ¬ll with a level surface), his model
however, minor differences in deformation may
cannot readily be generalized to more complex
have in¬‚uenced erosion more greatly than varia-
2D cases in which a level surface cannot be as-
tions in ice thickness or basal sliding. This is not
sumed. MacGregor et al. (2000) developed a 1D
a great concern for distinctive landforms such as
model to study the formation of overdeepenings
¬nger lakes. However, if bedrock geologic maps
using a model based on mass-balance, ice-creep,
and structural data are available for a region,
and stress-controlled basal sliding. They showed
structural control should be evaluated through
that several classic alpine-glacial landforms, in-
correlations between topography and structure.
cluding overdeepenings and hanging tributaries,
Determining the controls on glacial bedrock-
could be reproduced by their model. Their results
erosion rates is an active area of research. Boul-
were limited to 1D, however.
ton (1974) emphasized the role of ice thickness
Braun and colleagues (Braun et al., 1999;
in determining erosion rates, while Hallet (1979,
Tomkin and Braun, 2002) constructed the ¬rst
1996) emphasized the importance of sliding ve-
2D glacial-landform evolution model. These au-
locity. In Boulton™s model, the hydrostatic pres-
thors used a continuity equation to model local
sure induced by the ice overburden pressure in-
changes in ice thickness:
creases the normal stress acting on subglacial de-
bris, in turn increasing their ability to abrade the
‚h
=∇·F +M
bed. In Hallet™s model, ice thickness is not a fac- (6.51)
‚t
tor because hydrostatic pressures are assumed to
where the ¬‚ux F was dependent on creep and
act equally on the tops and bottoms of rounded
sliding velocity and M is the mass balance or rate
subglacial particles. In this case, the velocity of
152 NON-NEWTONIAN FLOW EQUATIONS


of accumulation and ablation. Sliding velocities the basal shear stress can also be varied to simu-
in their model were given by late ˜˜sticky™™ spots on the bed. Both assumptions
yield nearly identical results.
2A s βc (ρgh)n
us = |∇(h + b)|n’1 ∇(h + b) (6.52) The ice-surface topography and basal-sliding
N+P
velocities above a ¬‚at surface with white noise
where h is the ice thickness, b is the bedrock variations (shown in Figure 6.24c) are illustrated
elevation, and A s , βc , N , n, and P are empiri- in Figures 6.24a and 6.24b for a basal shear stress
of 1 bar. The grid for this example is 50 km —

<< . .

. 29
( : 51)



. . >>