- Email: [email protected]
Foreland basin—a FORTRAN program to model the formation of foreland basins resulting from the flexural deflection of the lithosphere caused by a time-varying distributed load
Computers & Geosciences Vol. 19, No. 9, pp. 1321-1332, 1993 Printed in Great Britain. All rights reserved
0098-3004/93 $6.00 + 0.00 Copyright © 1993 Pergamon Press Ltd
F O R E L A N D B A S I N - - A F O R T R A N PROGRAM TO MODEL THE FORMATION OF F O R E L A N D BASINS RESULTING FROM THE FLEXURAL DEFLECTION OF THE LITHOSPHERE CAUSED BY A TIME-VARYING DISTRIBUTED LOAD JUAN HOMERO HINOJOSAl and KEVIN L. MICKUS2 tDivision of Arts & Sciences, Laredo State University, Laredo, TX 78040-9960 and 2Department of Geoseiences, Southwest Missouri State University, Springfield, MO 65804, U.S.A. (Received 17 November 1992; accepted 17 April 1993)
A~traet--A FORTRAN program is described that calculates the deflection of a thin elastic plate resulting from a time-varying distributed load. The program is useful for modeling the subsidence of a foreland basin where the elastic plate is the lithosphere and the distributed load is a wedge-shaped thrust sheet. The deflection is calculated for a specific load distribution, and hence at a specific time but one can prescribe various loads at different times in the loading history so that the basin subsidence history can be shown. Key Words: Flexure, Elastic plates, Foreland basins, Thrust sheets, Finite differences, Gauss-Seidel iteration, Successive overrelaxation.
INTRODUCTION
Foreland basins have long been recognized as important elements of the tectonic cycle [and geosynclinal cycle (Aubouin, 1965)]. Dickinson (1974) introduced the term foreland basin and proposed two endmembers: (1) retroarc basins (e.g. Late MesozoicCenozoic basins of the Rocky Mountain region) which are situated behind a magrnatic arc and are associated with subduction of oceanic lithosphere, and (2) peripheral basins (e.g. modern Indo-Gangetic Basin south of the Himalayan Mountains) which form beside foreland thrust belts that are adjacent to a suture zone. Whatever their origins, foreland basins are elongated subsiding troughs that form on the cratonic side of thrust belts throughout the world. Foreland basins range greatly in size, shape, and particularly in the rheological properties of the lithosphere, so one must be careful of the terminological problems involved in describing how foreland basins evolve (Allen and Allen, 1990). The driving mechanism for foreland basin formation is tectonic thickening of the crust resulting from the thrusting of sediments onto the craton (Price, 1973; Dickinson, 1974). To better understand the mechanisms of sedimentary basin formation, several authors (e.g. Sleep and Snell, 1976; Beaumont, 1978; Steckler and Watts, 1978) began quantitative modeling of basin formation using elastic and viscoelastic bending theory. Beaumont (1981) and Jordan (1981) modeled foreland basin formation by having thrust sheets act as loads on an elastic litho-
sphere. Beaumont and Jordan's models only predict the overall basin geometry so Fleming and Jordan (1990) incorporated time-dependent erosion of mass and the deposition of this mass into a foreland basin to predict the basin's stratigraphy. There has been much debate about the state of the lithosphere during foreland basin formation (Beaumont, 1978; Jordan, 1981). Beaumont used a lithosphere that behaves viscoelastically and showed that this theory predicted the subsidence of the Alberta Basin in Canada. Using more geological control (stratigraphic and structural), Jordan explained the subsidence of the Rocky Mountain foreland basin in Wyoming using elastic plates. Although both theories may explain foreland basin subsidence in different situations, one may have to use a model that evokes both situations where there is an initial softening of the lithosphere followed by a period of near constant rigidity (Kusznir and Karner, 1985). An additional complexity of foreland basin formation is determining the elastic thickness of the lithosphere at the time of loading which is controlled by its thermal state which in turn is controlled by the tectonic history of the region (Courtney and Beaumont, 1983). Many foreland basins, for example the Val Verde in West Texas (Wuellner, Lehtonen, and James, 1986) and Arkoma in Oklahoma and Arkansas (Houseknecht, 1986), were formed initially on a thinned passive margin with lithospheric thickness increasing with time. Stockman, Beaumont, and Boutilier (1986) examined this problem of differing lithospheric properties in the formation of foreland basins.
1321
1322
J.H. HINOJOSAand K. L. M1CKUS
Although there are many arguments on how foreland basins form, elastic plate theory is adequate for understanding many significant features of their formation. We present a program that uses elastic theory to solve for the deflection of a thin plate resulting from a time-varying distributed load. As the plate is being deformed, the load can be expressed by any number of suitable geologic events (e.g. structural or sediment loading). We use an advancing thrust sheet as a load to explain the formation of foreland basins. The program calculates the iithospheric deflection resulting from a given load at a specific time and the history of basin formation can be shown by calculating the deflection at specific times in the thrusting history. THEORY
The flexural deflection of a thin elastic plate resulting from a distributed, time-varying load has been solved by many authors (Turcotte and Schubert, 1982; Watts, Karner, and Steckler, 1982). We solve how an elastic plate is deflected while a time-varying load in the form of an advancing thrust sheet is emplaced [e.g. Appalachian Mountains (Quinland and Beaumont, 1984) or the Idaho-Wyoming thrust belt (Jordan, 1981)]. The deflection of the lithosphere results in the formation of a foreland basin. To solve for the deflection of a uniformly thick elastic lithosphere we use a fourth-order, linear elastic equation:
D V4W(X, t) "~ (6p)gw(x, t) = P(x, t),
(1)
where
6p = (Pm - P,)
if P(x, t) 4: O,
fip=(pm--ps)
ifP(x,t)=0,
where D is the flexural rigidity of the elastic lithosphere, V4 is the biharmonic operator, w(x, t) is the amount of elastic deflection as a function of spatial coordinates (x) and time (t), Pm is the mantle density, Ps is the sediment density within the foreland basin, Pt is the thrust sheet density, g is the acceleration of gravity, and P(x, t) is the distributed and time-varying load resulting from advancement of the thrust sheets. D is related to the thickness of the lithosphere (T,) by
ET~ D = - 12(l - v~) '
(2)
where E is Young's modulus and v is Poisson's ratio. The geometry of foreland basins generally is twodimensional (Beaumont, 1981; Allen, Homewood, and Williams, 1986), so Equation (l) can be simplified to a function of only one spatial variable (e.g. x). Equation (1) then can be written as 04w D ~x 4 + (6p )gw = P(x, t),
(3)
where w = w(x, t). The load, P(x, t), can be written as
P(x, t) = ptgh(x, t),
(4)
where h(x, t) is the variable topography of the thrust sheet. Equation (3) now becomes t~4w D ~x 4 + (6p)gw = ptgh(x, t).
(5)
To simplify the determination of h(x, t), the advancing thrust sheet is modeled as a triangular shaped wedge advancing at a uniform rate v. The simplification of the thrust sheet geometry is justified because numerous published cross sections of thrust sheets [e.g. the Appalachian Mountains (Hatcher, 1989) or the Ouachita-Marathon orogeny (Ross, 1986; Viele and Thomas, 1989)] can be considered to be roughly triangular in cross section. Most thrust sheets do not advance at a uniform rate but occur in several pulses or episodes in which the advancement rate can be modeled as uniform (Hatcher, 1989; Reed and Strickler, 1990). So, our model can be applied to the specific episodes and the cumulative effects can be summed together. The expression for a triangular shaped wedge advancing at a uniform rate is given by
h(x,t)=(vt-x)tanO,
O~t<<.TM,
O~x<<.vt
(6) where TM is the duration of the thrusting event in millions of years. Note that Equation (6) satisfies t?h/dt + v Oh/Ox = 0. This condition states that the wedge thickness does not change with time at a point following the motion of the thrust sheet. We therefore are not allowing for local deformation of the advancing thrust wedge with time. Equation (5) was solved numerically using finite differences with Gauss-Seidel iteration and successive overrelaxation (SOR). The initial and boundary conditions used to solve Equation (5) are w=0 ~2W Ox 2
att=0, ~3W =0 Ox 3
Ow w=~--~x=0
O~X<.Xm, x at x = 0 , t > 0
atX=Xmax,
t>0
(7) (8) (9)
where Xm,x is the total length of the plate. PROGRAM USAGE The program F O R E L A N D BASIN as given in the Appendix was written in F O R T R A N 77 on a Digital Equipment Corporation (DEC) 5000/200 computer using the ultrix operating system version 4.2A. A sample input and output is given in the Appendix to illustrate program usage. The program consists of a main program and three subroutines (INPUT, OUTPUT, FLEX). A brief description of each subroutine is given in the main program listing.
1323
Modeling the formation of foreland basins Foreland basin model geometry
Wa
L 0 -2
= thrust sheet length = a n g l e o f thrust
h 0 = m a x i m u m load t o p o g r a p h y T T = maximum thrust wedge thickness x 0 = foreland basin width x a = arc position
-3
w 0 = f o r e l a n d b a s i n d e p th w a = arch a m p l i t u d e
-4
-5
0
50
I
I
I
I
i
L
J
k
I
100
150
200
250
300
350
400
450
5O0
Distance from edge of plate (km) Figure 1. Geometry of elastic plate (lithosphere) and distributed load (thrust sheet). The program calculates the deflection of the lithosphere along any given profile along the surface of the Earth (Fig. 1). The input parameters describing the properties of the lithosphere, the amount of the initial load, where to calculate the deflections along the profile, and factors involved in numerically solving Equation (5) are described briefly in the main program. These parameters will be discussed in more detail next. The distance units are in kilometers (km), densities are in kg/m 3, loads are in megaPascals (MPa), Young's modulus is in Pascals (Pa), and accelerations are in m/s 2. Equation (5) is solved numerically using an iterative Gauss-Seidel technique. Most iterative techniques for solving differential equations converge slowly to a solution where the most important problem becomes determining a technique to speed up the convergence (Forsythe, Malcolm, and Moler, 1977). To speed up the convergence, we used the successive overrelaxation (SOR) method. The SOR method requires the use of an acceleration parameter, o9, where o9 has the value 1 < o9 < 2. If ~o is 1, then the SOR method reverts back to the regular Gauss-Seidel technique and the convergence is slow. For any other number in the given range, the convergence is accelerated. The optimum selection of o9 is problem-dependent and is difficult to determine "a priori" (Carnahan, Luther, and Wilkes, 1969) but we have determined that to = 1.5 works best for our problem. The user may want to vary o~ to understand the complete usage of the program. The parameters im~xand dx are used to describe the grid of nodal points used to solve numerically for the
deflection of the elastic plate. The length of the elastic plate must be greater than the region of deformation. This is so that the deflection goes to zero far from the thrust sheet front. /max is the number of nodal points used to describe the elastic plate. For resolution and accuracy purposes, /max should be greater than about 60. Also, /maxmust be less than the dimensions of the program's arrays (presently set at 201). dx is the spatial increment of the nodal points and is determined by dividing the length of the elastic plate by / m a x - 1. In our work in the Ouachita orogenic belt, we use 500km for the length of the plate and /max= 101 which would make dx = 5 km. The lithospheric thickness, To, is the effective elastic thickness of the plate we are deflecting. The determination of this thickness for continental regions is not trivial because the composition of the lower crust and the lithospheric temperatures are not well known (McNutt, Diament, and Kogan, 1988). Subsequently, there are many techniques (e.g. seismic refraction, gravity modeling, admittance studies, flexural studies) that are used to determine To. In forward gravity modeling, the gravity field of the deflected plate is compared with the observed Bouguer gravity anomalies. The lithospheric thickness and density are differed until the observed and predicted gravity anomalies match. In the admittance techniques, the averaged ratios of the Bouguer gravity anomaly to topography in selected wavelength ranges are compared to those predicted by the deflected plate models. The problem with admittance studies in continental areas, unlike oceanic areas, is that the load causing the deflection may not have a
1324
J.H. HINOJOSAand K. L. MICKUS
corresponding topographic expression (Karner and Watts, 1983; Forsyth, 1985) which is a requirement for the admittance technique. Some authors (e.g. Royden and Karner, 1984) have used flexure modeling of foreland basins constrained by seismic reflection and drill core data to determine the lithospheric thickness. Another technique we have used to determine To is to invert magnetic anomaly data for Curie point depths (Hinojosa and Mickus, 1991; Mickus and Hinojosa, 1992). The Curie point depth is defined as the depth to the 550°C isotherm (e.g. Mayhew, 1982, 1985) assuming the magnetic mineral is magnetite. Once the Curie point depths were calculated, Te was obtained by linear interpolation to the 450°C isotherm. This is justified because the geotherm at shallow depths (crustal depths) increases in a nearly linear fashion. Even though there are numerous techniques to estimate Te one cannot be sure that the value being used is correct. The best technique may be to estimate Te by more than one technique and average the values obtained. The maximum load (load 0) used in the program is determined "a posteriori" in a trial and error iterative scheme so as to satisfy the constraint imposed by the maximum thrust wedge thickness ( T T ) obtained from cross sections (cf. Fig. 1). The parameter T T is related to the maximum plate deflection [w (x -- 0)] and the maximum load topography (h0) by T T = w ( x = O) + h o. Note that the maximum load topography (h0) is that part of the thrust sheet that is above the plane z = 0, where z = 0 is the surface of the Earth. The parameter load 0 then can be calculated from
(1989) and Reed and Strickler (1990). The maximum thrust length and thickness we estimated are 50 and 6 km, respectively. Frankel and others (1988) estimated the duration of thrusting to be 40-50 Myr. Reed and Strickler (1990) recognized at least four distinct thrusting events with periods of no deformation between events in the Marathon Uplift. Even so, to simplify our models we will assume that the thrust rate was constant. The densities we used were determined from gravity modeling (Smith, 1986). The densities for the thrust sheet, basin sediments, and the mantle are 2820, 2600, and 3300 kg/m 3, respectively. We used two lithospheric thickness: (1) 12 km which was estimated from Curie point depth calculations (Hinojosa and Mickus, 1991), and (2) 29 km which was estimated from gravity modeling (Smith, 1986). We used 500 km as the length of the elastic plate, because published maps of the Val Verde Basin region (Frankel and others, 1988) show that the basin is narrower than 500 km. All the parameters are shown in Table 1. We calculated the deflection of the lithosphere using elastic thicknesses of both 12 and 29 km at 10Myr intervals for a total duration of 40 Myr (Figs. 2 and 3). Using a lithospheric thickness of Table 1. Flexure model parameters
Symbol
Value
Length of model
Xmax
500 km
Effective elastic thickness
T,
12 km, 29 km
where h0 = T T - w (x = 0). Most of the parameter values used in the program have to be estimated and depend on the basin being modeled. The thrust sheet length (L) and the thickness ( T T ) can be estimated from published geological data and cross sections. To model the basin's formation, the user can estimate the thrust sheet's rate of advancement, determine the thrust length and thickness, and subsequently the loads for specific time periods during the thrusting history.
Maximum thrust
L
50 km
Maximum thrust sheet thickness
T]"
6 km
Sediment density
p,
2600 Kg m'3
Thrust sheet density
p,
2820 Kg m'3
Mantle density
p.
3300 Kg m a
Total duration of thrusting
At
40 My
EXAMPLE
Average rate of advance of thrust sheet
v
1.25 km My a
Acceleration of
g
9.8 m s 2
Young's modulus
E
0.7x10 u Pa
Poisson's ratio
v
0.25
SOR acceleration parameter
oJ
1.5
Number of nodal points
imax
101
Convergence parameter
eps
.003
load 0 = ptgho,
The application of the computer program is illustrated for the Pennsylvanian-Permian Ouachita orogenic event in West Texas. Specifically, we will show some preliminary results for formation of the Val Verde Basin which is the southeastern part of the oil-rich Permian Basin. There have been many geological studies (mainly for petroleum exploration) of the Val Verde Basin but few geophysical, and no flexural studies. The geological studies are summarized by Frankel and others (1988). We calculated maximum thrust lengths and thicknesses from cross sections through the basin by Nicholas and Waddell
Parameter
length
gravity
1325
Modeling the formation of foreland basins Foreland basin evolution
l 0 Myr
..~-"U'~"~'~:Z:'--_-..'~-~ ,~
20 Myr / ,
-I
E
_// / "t
-
/ // ]/
30 l~yr / " /
T e = 12 km
-2 ..' 40 Myr ,' /
-3
/ / /
-4
-5 0
I
J
I
I
I
I
I
I
J
I
50
100
150
200
250
300
350
400
450
500
Distance from edge of plate (km) Figure 2. Deflection of 12-km thick lithosphere in 10 Myr increments resulting from thrust loads emplaced on lithosphere during Ouachita orogenic event in Val Verde Basin region of West Texas. 12km (Fig. 2), the maximum plate deflection at 40 Myr is about 4.7 km. This should not be confused with the basin depth. The basin depth is the plate deflection at the thrust wedge front, which is approximately 1.6 km for the 12-km thick plate at 40 Myr. The corresponding basin width is about 53 km. There
is a forebulge, or arch, positioned at about 140 km from the edge of the plate, or about 40 km from the zero crossover point, with an amplitude of about 0.21 km. Using a thicker lithosphere (Fig. 3), the basin depth at 40 Myr is about 2.4 kin, with a basin width of about 120 km. The arch is located about
Foreland basin evolution
~ - . ' 7 . . v/ / 20 Myr
..- / /
..
f//
..
;.>
j
-1
E
30 Myr -2
_.""
/
/
.,'
T e = 29 km
// 40 Myr
/ -3
-
/
/ / -4
-5
_
I
I
50
I00
150
I
I
I
L
200
250
300
350
i
I
I
400
450
500
Distance from edge of plate ( k i n ) Figure 3. Deflection of 29-km thick lithosphere in 10 M y r increments resulting from thrust loads emplaced on lithosphere during Ouachita orogenic event in Val Verde Basin region of West Texas.
1326
J . H . HINOJOSA and K. L. MICKUS Table 2. Deflection of lithosphere resulting from various thrust loads. Parameters include basin width (x0), arch position (xa), basin depth (w0), arch amplitude (wa), lithospheric thickness (Te), thrust length (L), and time (My). All distance units are in kilometers and loads are in megaPascals Te
Load0
L
Wo
Wo
xo
xo
Time
12
8.75
12.5
-0.27
0.02
69.6
37.9
10
12
17.50
25.0
-0.78
0.07
62.6
39.1
20
12
26.25
37.5
-1.30
0.14
57.6
39.0
30
12
35.00
50.0
-1.64
0.21
52.5
38.7
40
29
13.25
12.5
-0.21
0.01
134.2
58.3
10
29
26.50
25.0
-0.76
0.06
131.5
73.5
20
29
39.75
37.5
-1.55
0.13
126.4
76.1
30
29
53.00
50.0
-2.35
0.22
120.1
74.9
40
75 k m from the zero crossover point, a n d has a n amplitude o f a b o u t 0.22 km. Table 2 shows the loads, basin widths, arch positions, a n d arch amplitudes as a function o f time. The input a n d o u t p u t for the situation o f T e = 12 k m a n d L = 50 k m is shown in the Appendix. CONCLUSIONS The p r o g r a m described is useful for modeling the subsidence of a foreland basin formed by loading of thrust sheets o n an elastic lithosphere. Even t h o u g h the p r o g r a m calculates the deflection at a specific time, the basin's subsidence history can be determined by k n o w i n g the loads at various times t h r o u g h o u t the orogenic event. The use o f the p r o g r a m is a first step in modeling the f o r m a t i o n of foreland basins. I m p r o v e m e n t s to the technique would be to allow variable elastic thickness, t o p o g r a p h y on the thrust sheet, erosion of the thrust sheet a n d subsequent sediment deposition into the basin, the i n c o r p o r a t i o n of viscoelastic bending theory, a n d the possible inclusion of two-dimensional effects. REFERENCES
Allen, P. A., and Allen, J. R., 1990, Basin analysis: principles and applications: Blackwell Scientific Publications, Oxford, 451 p. Allen, P. A., Homewood, P., and Williams, G. D., 1986, Foreland basins: an introduction, in Allen, P. A., and Homewood, P., eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, Blackwell Scientific, Oxford, p. 3-12. Aubouin, J., 1965, Geosynclines, in Developments in geotectonics: Elsevier, Amsterdam, 335 p. Beaumont, C., 1978, The evoluation of sedimentary basins on a viscoelastic lithosphere: theory and examples: Geophys. Jour. Roy. Astr. Soc., v. 55, no. 2, p. 471~497.
Beaumont, C., 1981, Foreland basins: Geophys. Jour. Roy. Astr. Soc., v. 65, no. 2, p. 291-329. Carnahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied numerical methods: John Wiley & Sons, New York, 604 p. Courtney, R. C., and Beaumont, C., 1983, Thermally activated creep and flexure of the oceanic lithosphere: Nature, v. 305, no. 5931, p. 201-204. Dickinson, W. R., 1974, Plate tectonics and sedimentation, in Dickinson, W. R., ed., Tectonics and sedimentation: Soc. Economic Paleontologists and Mineralogists, Spec. Publ. vol. 22, p. 1-27. Fleming, P. B., and Jordan, T. E., 1990, Stratigraphic modeling of foreland basins: interpreting thrust deformation and lithosphere rheology: Geology, v. 18, no. 5, p. 430-434. Frankel, H. N., and others, 1988, The Permain Basin region, in Sloss, L. L., ed., Sedimentary cover-North American Cration: U.S.: The Decade of North American Geology, v. D-2, Geol. Soc. America, p. 261 306. Forsyth, D. W., 1985, Subsurface loading and estimates of the flexural rigidity of continental lithosphere: Jour. Geophys. Res., v. 90, no. BI4, p. 12623-12632. Forsythe, G. E., Malcolm, M. A., and Moler, C. B., 1977, Computer methods for mathematical computations: Prentice-Hall Inc., New York, 259 p. Hatcher, R. D., 1989, Tectonic synthesis of the U.S. Appalachians, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian-Ouachita Orogen in the United States: The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 511-536. Hinojosa, J. H., and Mickus, K. L., 1991, Comparative geophysical study of the Ouachita foreland basins in Texas: EOS: Trans. Am. Geophys. Union, v. 72, no. 44, p. 430. Houseknecht, D. W., 1986, Evolution from passive margin to foreland basin: the Atoka Formation of the Arkoma basin, south-central U.S.A., in Allen, P. A., and Homewood, P.. eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, Blackwell Scientific, Oxford, p. 327 345. Jordan, T. E., 1981, Thrust loads and foreland basin evolution, Cretaceous, western United States: Am. Assoc. Petroleum Geologists Bull. v. 65, no. 12, p. 250(~2520.
Modeling the formation of foreland basins Karner, G. D., and Watts, A. B., 1983, Gravity anomalies and flexure of the lithosphere at mountain ranges: Jour. Geophys. Res. v. 88, no. B1 I, p. 10449-10477. Kusznir, N., and Karner, G. D., 1985, Tectonic subsidence calculated from lithified basin strata (abst.): Geol. Soc. America, Abstracts with Programs, v. 14, no. 5, p. 534. Mayhew, M. A., 1982, Application of satellite magnetic anomaly data to Curie isotherm mapping: Jour. Geophys. Res., v. 87, no. 6, p. 48464854. Mayhew, M. A., 1985, Curie isotherm surfaces inferred from high-altitude magnetic anomaly data: Jour. Geophys. Res., v. 90, no. 3, p. 2647 2654. McNutt, M. K., Diament, M., and Kogan, M. G., 1988, Variations of elastic plate thickness at continental thrust belts: Jour. Geophy. Res., v. 93, no. B8, p. 8825-8838. Mickus, K. L., and Hinojosa, J. H., 1992, Flexure origin of the Kerr basin, Texas (abst.): EOS: Trans. Am. Geophys. Union, v. 73, no. 24, p. 289. Nicholas, R. L., and Waddell, D. E., 1989, The Ouachita system in the subsurface of Texas, Arkansas, and Louisiana, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian~)uachita Orogen in the United States:-The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 661~571. Price, R. A., 1973, Large-scale gravitational flow of supracrustal rocks, southern Canadian Rockies, in DeJong, K. A., and Schotten, R., eds., Gravity and tectonics: John Wiley & Sons, New York, p. 491-502. Quinlan, G. M., and Beaumont, C., 1984, Appalachian thrusting, lithospheric flexure, and the Paleozoic stratigraphy of the Eastern Interior of North America: Can. Jour. Earth Sciences, v. 21, no. 9, p. 973-995. Reed, T. A., and Strickler, D. L., 1990, Structural geology and petroleum exploration of the Marathon thrust belt, West Texas, in Laroche, T. M., and Higgins, L., eds., Marathon thrust belt: structure, stratigraphy, and hydrocarbon potential: West Texas Geol. Soc. Field Seminar, p. 39-82. Ross, C. A., 1986, Paleozoic evolution of southern margin of Permian Basin: Geol. Soc. America Bull., v. 97, no. 5, p. 536-554.
Royden, L. H., and Karner, G. D., 1984, Flexure of the continental lithosphere beneath the Appenines and Carpathians foredeep basins: evidence for an insutficient topographic load: Am. Assoc. Petroleum Geologists Bull., v. 68, no. 6, p. 704-712. Sleep, N. H., and Snell, N. S., 1976, Thermal contraction and flexure of mid-continent and Atlantic marginal basins: Geophys. Jour. Roy. Astr. Soc., v. 45, no. 1, p. 125-154. Smith, K. J., 1986, A gravity and tectonic study of the southwestern portion of the Ouachita system: unpubl. masters thesis, Univ. Texas, E1 Paso, 99 p. Steckler, M. S., and Watts, A. B., 1978, Subsidence of the Atlantic-type continental margin off New York: Earth and Planet. Sci. Lett., v. 41, no. 1, p. 1 13. Stockmal, G. S., Beaumont, C., and Boutilier, R., 1986, Geodynamic models of convergent margin tectonics: transition from rifted margin to overthrust belt and consequences for foreland basin development: Am. Assoc. Petroleum Geologists Bull., v. 70, no. 2, p. 181-190. Turcotte, D. L., and Schubert, G., 1982, Geodynamics: applications of continuum mechanics to geological problems: John Wiley & Sons, New York, 450 p. Viele, G. W., and Thomas, W. A., 1989, Tectonic synthesis of the Ouachita orogenic belt, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian-Ouachita Orogen in the United States: The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 695-728. Watts, A. B., Karner, G. D., and Steckler, M. S., 1982, Lithospheric flexure and the evolution of sedimentary basins: Phil. Trans. Roy. Soc. London, v. 305A, no. 2, p. 249-281. Wuellner, D. E., Lehtonen, L. R., and James, W. C., 1986, Sedimentary-tectonic development of the Marathon and Val Verde basins. West Texas, U.S.A.: a PermianCarboniferous migrating foredeep, in Allen, P. A., and Homewood, P., eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, BlackweU Scientific, Oxford, p. 327 345.
APPENDIX Program Listing ******************************************************************** PROGRAM
FORELAND
BASIN
PROGRAM TO COMPUTE THE FLEXURAL DEFORMATION OF AN ELASTIC PLATE PLATE DUE TO A DISTRIBUTED, TIME VARYING THRUST SHEET THAT WILL FORM A FORELAND BASIN. THE PROBLEM IS SOLVED NUMERICALLY USING A ONE-DIMENSIONAL, FOURTH-ORDER, LINEAR ELASTIC FLEXURE EQUATION SUBJECT TO A DISTRIBUTED LOAD, WITH A UNIFORM EFFECTIVE ELASTIC THICKNESS USING THE GAUSS-SEIDEL ITERATION METHOD WITH SUCCESSIVE OVER-RELAXATION (SOR).
INPUT PARAMETERS : RHOM RHOT RHOS L
= = = =
EPS
=
E G NU OMEGA
= = = =
1327
MANTLE DENSITY (ALL DENSITIES ARE IN KG/M**3) THRUST-WEDGE DENSITY SEDIMENT DENSITY WITHIN THE FORELAND BASIN THRUST DISTANCE IN THE DIRECTION OF THE ADVANCING WEDGE (KILOMETERS) CONVERGENCE PARAMETER FOR THE GAUSS-SEIDEL ITERATION METHOD (FOR BEST RESULTS WE USE 1.E-3 TO I.E-4) Y O U N G ' S M O D U L U S (Pa) GRAVITATIONAL ACCELERATION (M/S**2) POISSON'S RATIO (DIMENSIONLESS) SOR ACCELERATION PARAMETER, MUST BE BETWEEN 1 A N D 2. W E U S U A L L Y U S E 1 . 5 B U T T H E S E L E C T I O N IS TRIAL AND ERROR DEPENDING ON THE PROBLEM
328
J . H . HINOJOSA and K. L. MICKUS = MAXIMUM NUMBER OF NODAL POINTS. FOR RESOLUTION C IMAX c T H E N U M B E R S H O U L D B E > 60. B U T I T M U S T B E L E S S THAN PROGRAM'S ARRAYS DIMENSION SIZE (PRESENTLY SET C C AT 201). IF ONE WANTS MORE THAN 201 THE DIMENSION C OF THE PROGRAM'S ARRAYS MUST BE INCREASED. c = SPATIAL INCREMENT OF NODAL POINTS (KILOMETERS) DX c = EFFECTIVE ELASTIC THICKNESS OF THE LITHOSPHERE TE c (KILOMETERS) C LOAD0 = M A X I M U M T H R U S T W E D G E L O A D (MPa). C c OUTPUT PARAMETERS : c c W(I) = DEFLECTION OF THE ELASTIC PLATE AS A FUNCTION OF c DISTANCE. (KILOMETERS) C DI S T (I ) = DISTANCE FROM COMPUTATIONAL ORIGIN (KILOMETERS) C L O A D (I ) = D I S T R I B U T E D L O A D (MPa) c C OTHER INTERNAL PARAMETERS: C = FLEXURAL RIGIDITY (NEWTON-M) C D = DEFLECTION OF THE ELASTIC PLATE AS A FUNCTION OF C Wl(I) DISTANCE FROM THE PREVIOUS ITERATION C = LENGTH OF THE ELASTIC PLATE, XMAX=DX*(IMAX-1) C XMAX C SUBROUTINES USED: C c INPUT - READS IN THE VARIOUS PARAMETERS C c OUTPUT - WRITES OUT THE VARIOUS PARAMETERS, THE LOAD AS A c FUNCTION OF DISTANCE, AND THE DEFLECTION OF THE c LITHOSPHERE AS A FUNCTION OF DISTANCE c FLEX - CALCULATES THE DEFLECTION OF THE LITHOSPHERE c DUE TO AN ADVANCING THRUST SHEET c c WRITE UNITS: C C 52 - T H E I N P U T P A R A M E T E R S A N D T H E D E F L E C T I O N O F T H E C LITHOSPHERE AS A FUNCTION OF DISTANCE. c 53 - T H E L O A D D U E T O T H E A D V A N C I N G T H R U S T S H E E T A S C A FUNCTION OF DISTANCE. c 6 - QUESTIONS ASKING FOR THE FILE CONTAINING THE C INPUT PARAMETERS, AND OUTPUT FILES THAT CONTAIN c CALCULATED DEFLECTION AND LOAD. ALSO A BANNER c A N N O U N C I N G T H E N A M E O F T H E P R O G R A M IS P R I N T E D . c READ UNITS : C c c 51 - T H E I N P U T P A R A M E T E R S c 5 - T H E N A M E S O F T H E I N P U T (OLD) A N D O U T P U T (NEW) F I L E S c c THIS PROGRAM WAS WRITTEN ON A DEC 5000/200 USING AN ULTRIX c OPERATING SYSTEM VERSION 4.2A c ******************************************************************** REAL W(201),WI(201),LOAD(201),L,DIST(201),LOAD0,OMEGA,NU,M CHARACTER*20 FNAMEI,FNAME2,FNAME3 CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD C WRITE(6, ****************************************************** WRITE(6, *)'* F O R E L A N D B A S I N -- C O M P U T E S T H E F L E X U R A L *' W R I T E ( 6 , *) '* DEFORMATION OF AN ELASTIC PLATE TO FORM A *' WRITE(6, *)'* FORELAND BASIN *' WRITE(6, ****************************************************** WRITE(6, *)'ENTER INPUT DATA FILE NAME:' R E A D ( 5 , ' (A) ' ) F N A M E I O P E N (51, F I L E = F N A M E I , S T A T U S = 'O L D ' ) WRITE(6,*) 'ENTER OUTPUT DATA FILE NAME:' R E A D (5, " (A) ' ) F N A M E 2 O P E N (52, F I L E = F N A M E 2 , S T A T U S = ' N E W ' ) W R I T E ( 6 , * ) 'DO Y O U W A N T T O W R I T E O U T T H E L O A D ( Y / N ) ? ' R E A D ( 5 , ' (A) ' ) A N S IF ( A N S . E Q . ' Y ' ) T H E N W R I T E (6, *) 'E N T E R N A M E O F L O A D F I L E : ' R E A D (5, " (A) ' ) F N A M E 3 O P E N (53, F I L E = F N A M E 3 , S T A T U S = ' N E W ' )
Modeling the formation of foreland basins ELSE CONTINUE ENDIF ******************************************************************** C
READ
IN T H E P A R A M E T E R S
CALL INPUT ******************************************************************** C D E T E R M I N E T H E D E F L E C T I O N OF T H E E L A S T I C P L A T E ******************************************************************** CALL FLEX ******************************************************************** C W R I T E O U T T H E C O M P U T E D D E F L E C T I O N A N D LOAD, C AND THE INPUT PARAMETERS ******************************************************************** CALL OUTPUT C C5OSE(51) CLOSE(52)
CLOSE(53) C STOP END ******************************************************************** ******************************************************************** SUBROUTINE INPUT ******************************************************************** C S U B R O U T I N E T O I N P U T P A R A M E T E R S U S E D TO D E S C R I B E T H E R H E O L O G I C A L C P R O P E R T I E S OF A T H I N E L A S T I C P L A T E A N D T H E D I S T R I B U T E D L O A D C (SEE T H E M A I N P R O G R A M F O R A D E S C R I P T I O N OF T H E S E PARAMETERS) ******************************************************************** R E A L L, LOAD0, OMEGA, NU, M CHARACTER* 1 ANS C O M M O N / P A R M / L O A D 0 , DX, IMAX, D, OMEGA, EPS, M, L, TE, A N S C O M M O N / D E N S / R H O M , RHOT, RHOS, G READ(51,*) READ(51,*) READ(51,*)
E,NU,OMEGA,IMAX,L RHOM,RHOT,RHOS,G DX,TE,LOAD0,EPS
DX=DX*I.E3 LOAD0=LOAD0*I.E6 M = (LOAD0/(L*I.E3)) TE=TE*I.E3 ******************************************************************* C CALCULATE THE FLEXURAL RIGIDITY ******************************************************************* D = E*(TE**3)/(12.*(I.-(NU)**2)) C RETURN END ******************************************************************** ******************************************************************** SUBROUTINE OUTPUT ******************************************************************** C S U B R O U T I N E TO W R I T E O U T T H E I N P U T PARAMETERS, D I S T R I B U T E D LOAD C AND THE DEFLECTION OF A THIN ELASTIC PLATE ******************************************************************** REAL W(201),LOAD(201),DIST(201),LOAD0,L CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD C LOAD0=LOADO/I.E6 TE=TE/I.E3 WRITE(52,100) TE,LOAD0,L,RHOM,RHOT,RHOS IF(ANS.EQ.'Y') T H E N WRITE(53,200) WRITE(53,'(FIO.3,5X,E9.4)') (DIST(I),LOAD(I),I=I,IMAX) ELSE CONTINUE ENDIF WRITE(52,300) W R I T E ( 5 2 , ' ( 2 ( I X , F I O . 3 ) ) ' ) (DIST(I),W(I),I=I,IMAX) C i00 F O R M A T ( ' T H E L I T H O S P H E R I C T H I C K N E S S IS ', F5.1,/, 1"THE M A X I M U M L O A D IN M P A IS ', F8.3,/,
1329
1330
J . H . HINOJOSA and K. L. MICKUS
200 300 C
2 ' T H E M A X I M U M T H R U S T D I S T A N C E IS " , F 5 . 1 , / , 3'THE DENSITIES OF THE MANTLE, THRUST SHEET AND 4 S E D I M E N T S A R E ", 3 ( F 6 . 1 , 1 X ) , / ) FORMAT(3X,'DISTANCE',8X,'LOAD') FORMAT(3X,'DISTANCE',3X,'DEFLECTION')
RETURN END ******************************************************************** ******************************************************************** SUBROUTINE FLEX ******************************************************************** C DETERMINES THE FLEXURAL DEFLECTION OF A UNIFORM, THIN C ELASTIC PLATE DUE TO A TIME VARYING, DISTRIBUTED LOAD C (SEE T H E M A I N P R O G R A M F O R A D E S C R I P T I O N O F T H E P A R A M E T E R S ) ******************************************************************** C REAL W(201),WI(201),LOAD(201),L,DIST(201),LOAD0,OMEGA,M CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD ******************************************************************** C DEFINE THE DISTRIBUTED LOAD ********************************************************************
10
DO i0 I I 0 = I , I M A X DIST(II0)=DX*FLOAT(II0-1)*I.E-3 IF ( D I S T ( I 1 0 ) . L E . L ) T H E N LOAD(I10)=LOAD0-M*(FLOAT(I10-1)*DX) ELSE LOAD(II0)=0. ENDIF CONTINUE
C BETA1 = D/(DX**4.) BETA2 = (RHOM - RHOT)*G B E T A 3 = (RHOM - RHOS)*G ******************************************************************** C START GAUSS-SEIDEL ITERATION TO COMPUTE DEFLECTION W(I): DO 20 I 2 0 = I , I M A X W ( I 2 0 ) = 0. 20 CONTINUE ******************************************************************** C RETURN HERE FOR NEXT ITERATION: ******************************************************************** 30
D O 40 I 4 0 = I , I M A X WI(I40) = W(I40) 40 CONTINUE ******************************************************************** C C
BOUNDARY CONDITIONS: W'' = W ' ' ' = 0 A T X = 0 (AT M A X I M U M
LOAD)
W(1) = ( I . / ( 2 . * B E T A I + B E T A 2 ) ) * ( L O A D ( 1 ) & - 2.*BETAI*W(3) + 4.*BETAI*W(2)) ******************************************************************** C A P P L Y SOR: ******************************************************************** = (i. - OMEGA)*WI(1) + OMEGA*W(1) IF (LOAD(2) .NE. 0) T H E N W(2) = ( I . / ( 5 . * B E T A I + B E T A 2 ) ) * ( L O A D ( 2 ) & - BETAI*W(4) + 4.*BETAI*W(3) & + 2.*BETAI*W(1)) ELSE W(2) = ( I . / ( 5 . * B E T A I + B E T A 3 ) ) * ( L O A D ( 2 ) & - EETAI*W(4) + 4.*EETAI*W(3) & + 2.*BETAI*W(1)) ENDIF ******************************************************************** W(1)
C A P P L Y SOR: ******************************************************************** W(2)
=
(I.
-
OMEGA)*WI(2)
+ OMEGA*W(2)
C DO 50 I 5 0 = 3 , I M A X - 2 IF(LOAD(I50) .ME. 0.) THEN W(I50) = (I./(6.*BETAI+BETA2))*(LOAD(I50)-BETAI*W(I50-2) & + 4.*BETAI*W(I50-1) + 4.eBETAI*W(IS0+I) & BETAI*W(I50+2))
Modeling the formation of foreland basins ELSE W(I50)
=
& &
133
(1./(6.*BETAI+BETA3))*(LOAD(I50)-BETAI*W(I50-2) + 4.*BETAI*W(I50-1) + 4.*BETAI*W(I50+I) BETAI*W(I50+2))
ENDIF ******************************************************************** C
APPLY
SOR:
******************************************************************** W ( I 5 0 ) = (1. - O M E G A ) * W l ( I 5 0 ) + O M E G A * W ( I 5 0 ) 50 CONTINUE ******************************************************************** C C
BOUNDARY CONDITIONS: W = W' = 0 A T X = X M A X
(FAR FROM THRUST
WEDGE
FRONT)
******************************************************************** W(IMAX-I)
=
&
(I./(7.*BETAI+BETA3))*(LOAD(IMAX-I) + 4.*BETAI*W(IMAX-2) - BETAI*W(IMAX-3))
******************************************************************** C
APPLY
SOR:
******************************************************************** W ( I M A X - I ) = (1. - O M E G A ) * W I ( I M A X - 1 ) + O M E G A * W ( I M A X - I ) ******************************************************************** C
COMPARE
WITH
PREVIOUS
ITERATION
TO TEST
FOR CONVERGENCE:
******************************************************************** D O 60 I 6 0 = I , I M A X IF (ABS(W(I60)-WI(I60)).GT.EPS) CONTINUE D O 70 I T 0 = I , I M A X W(I70) = -W(I70)*I.E-3 CONTINUE
60
70
GOTO
30
C RETURN END COMMAND ! ! ! ! ! !
IN.FLEXURE OUT.DEFLECTION ¥ OUT.LOAD
FILE
FOR RUNNING
FORELAND
BASIN
INPUT FILE DESCRIBING THE RHEOLOGICAL PROPERTIES PLATE AND THE DISTRIBUTED LOAD OUTPUT FILE CONTAINING THE PLATE DEFLECTION AS A OF DISTANCE YES OR NO ANSWER IF USER WANTS A FILE CONTAINING OUTPUT FILE CONTAINING THE LOAD AS A FUNCTION OF
OF THE FUNCTION THE LOAD DISTANCE
INPUT .7eli
.25 1 . 5
3300.
2820.
2600.
35.0
1.E-3
5.
12.
i01
50.
9.80
l Y O U N G ' S M O D U L U S (E), P O I S S O N ' S R A T I O (NU), S O R ! ACCELERATION PARAMETER (OMEGA), MAXIMUM NUMBER OF ! N O D A L P O I N T S (IMAX), T H R U S T D I S T A N C E (L) ! D E N S I T I E S O F M A N T L E (RHOM), T H R U S T - W E D G E (RHOT), ! AND FORELAND BASIN SEDIMENTS (RHOS), GRAVITATIONAL ! C O N S T A N T (G) ! S P A T I A L I N C R E M E N T O F N O D A L P O I N T S (DX), E L A S T I C ! T H I C K N E S S (TE), M A X I M U M L O A D ( L O A D 0 ) , C O N V E R G E N C E ! P A R A M E T E R (EPS) OUTPUT
THE THE THE THE
LITHOSPHERIC THICKNESS IS 12.0 M A X I M U M L O A D I N M P A IS 35.000 M A X I M U M T H R U S T D I S T A N C E IS 50.0 DENSITIES OF THE MANTLE, THRUST SHEET AND DISTANCE 0.000 5.000 10.000 15.000 20.000 25.000 30.000 35.000 40.000 45.000 50.000 55.000 60.000 65.000
CAGEO 19,9--1
DEFLECTION -4.736 -4.405 -4.075 -3.746 -3.420 -3.098 -2.784 -2.479 -2.186 -1.906 -1.643 -1.398 -1.172 -0.966
SEDIMENTS
ARE
3300.0
2820.0
2600.0
1332
J . H . HINOJOSAand K. L. MICKUS 70.000 75.000 80.000 85.000 90.000 95.000 100.000 105.000 Ii0.000 115.000 120.000 125.000 130.000 135.000 140.000 145.000 150.000 155.000 160.000 165.000 170.000 175.000 180.000 185.000 190.000 195.000 200.000 205.000 210.000 215.000 220.000 225.000 230.000 235.000 240.000 245.000 250.000 255.000 260.000 265.000 270.000 275.000 280.000 285.000 290.000 295.000 300.000 305.000 310.000 315.000 320.000 325.000 330.000 335.000 340.000 345.000 350.000 355.000 360.000 365.000 370.000 375.000 380.000 385.000 390.000 395.000 400.000 405.000 410.000 415.000 420.000 425.000 430.000 435.000 440.000 445.000 450.000 455.000
-0.780 -0.613 -0.466 -0.337 -0.225 -0.129 -0.049 0.018 0.073 0.116 0. 149 0.174 0.191 0.201 0.206 0.206 0.203 0.196 0.187 0.177 0.165 0.152 0.139 0 125 0 112 0 099 0 087 0 075 0 065 0 055 0 045 0 037 0.029 0.023 0.017 0.012 0.007 0.003 0.000 -0.002 -0.004 -0.006 -0.007 -0.008 -0.008 -0.009 -0.009 -0.009 -0.009 -0.008 -0.008 -0.007 -0.007 -0.006 -0.006 -0 0 0 5 - 0 004 -0 004 -0 003 -0 003 -0 003 -0.002 -0.002 -0.001 -0.001 -0.001 -0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
460.000 465.000 470.000 475.000 480.000 485.000 490.000 495.000 500.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
DISTANCE LOAD 0.000 .3500E+08 5.000 .3150E+08 10.000 .2800E+08 15.000 .2450E+08 20.000 .2100E+08 25.000 .1750E+08 30.000 .1400E+08 35.000 .I050E+08 40.000 .7000E+07 45.000 .3500E+07 50.000 .0000E+00 55.000 .0000E+00 60.000 .0000E+00 65.000 .0000E+00 70.000 .0000E+00 75.000 .0000E+00 80.000 .0000E+00 85.000 .0000E+00 90.000 .0000E+00 95.000 .0000E+00 100.000 .0000E+00 FOR All THE OTHER LOCATIONS
THE LOAD
IS 0 . 0
0098-3004/93 $6.00 + 0.00 Copyright © 1993 Pergamon Press Ltd
F O R E L A N D B A S I N - - A F O R T R A N PROGRAM TO MODEL THE FORMATION OF F O R E L A N D BASINS RESULTING FROM THE FLEXURAL DEFLECTION OF THE LITHOSPHERE CAUSED BY A TIME-VARYING DISTRIBUTED LOAD JUAN HOMERO HINOJOSAl and KEVIN L. MICKUS2 tDivision of Arts & Sciences, Laredo State University, Laredo, TX 78040-9960 and 2Department of Geoseiences, Southwest Missouri State University, Springfield, MO 65804, U.S.A. (Received 17 November 1992; accepted 17 April 1993)
A~traet--A FORTRAN program is described that calculates the deflection of a thin elastic plate resulting from a time-varying distributed load. The program is useful for modeling the subsidence of a foreland basin where the elastic plate is the lithosphere and the distributed load is a wedge-shaped thrust sheet. The deflection is calculated for a specific load distribution, and hence at a specific time but one can prescribe various loads at different times in the loading history so that the basin subsidence history can be shown. Key Words: Flexure, Elastic plates, Foreland basins, Thrust sheets, Finite differences, Gauss-Seidel iteration, Successive overrelaxation.
INTRODUCTION
Foreland basins have long been recognized as important elements of the tectonic cycle [and geosynclinal cycle (Aubouin, 1965)]. Dickinson (1974) introduced the term foreland basin and proposed two endmembers: (1) retroarc basins (e.g. Late MesozoicCenozoic basins of the Rocky Mountain region) which are situated behind a magrnatic arc and are associated with subduction of oceanic lithosphere, and (2) peripheral basins (e.g. modern Indo-Gangetic Basin south of the Himalayan Mountains) which form beside foreland thrust belts that are adjacent to a suture zone. Whatever their origins, foreland basins are elongated subsiding troughs that form on the cratonic side of thrust belts throughout the world. Foreland basins range greatly in size, shape, and particularly in the rheological properties of the lithosphere, so one must be careful of the terminological problems involved in describing how foreland basins evolve (Allen and Allen, 1990). The driving mechanism for foreland basin formation is tectonic thickening of the crust resulting from the thrusting of sediments onto the craton (Price, 1973; Dickinson, 1974). To better understand the mechanisms of sedimentary basin formation, several authors (e.g. Sleep and Snell, 1976; Beaumont, 1978; Steckler and Watts, 1978) began quantitative modeling of basin formation using elastic and viscoelastic bending theory. Beaumont (1981) and Jordan (1981) modeled foreland basin formation by having thrust sheets act as loads on an elastic litho-
sphere. Beaumont and Jordan's models only predict the overall basin geometry so Fleming and Jordan (1990) incorporated time-dependent erosion of mass and the deposition of this mass into a foreland basin to predict the basin's stratigraphy. There has been much debate about the state of the lithosphere during foreland basin formation (Beaumont, 1978; Jordan, 1981). Beaumont used a lithosphere that behaves viscoelastically and showed that this theory predicted the subsidence of the Alberta Basin in Canada. Using more geological control (stratigraphic and structural), Jordan explained the subsidence of the Rocky Mountain foreland basin in Wyoming using elastic plates. Although both theories may explain foreland basin subsidence in different situations, one may have to use a model that evokes both situations where there is an initial softening of the lithosphere followed by a period of near constant rigidity (Kusznir and Karner, 1985). An additional complexity of foreland basin formation is determining the elastic thickness of the lithosphere at the time of loading which is controlled by its thermal state which in turn is controlled by the tectonic history of the region (Courtney and Beaumont, 1983). Many foreland basins, for example the Val Verde in West Texas (Wuellner, Lehtonen, and James, 1986) and Arkoma in Oklahoma and Arkansas (Houseknecht, 1986), were formed initially on a thinned passive margin with lithospheric thickness increasing with time. Stockman, Beaumont, and Boutilier (1986) examined this problem of differing lithospheric properties in the formation of foreland basins.
1321
1322
J.H. HINOJOSAand K. L. M1CKUS
Although there are many arguments on how foreland basins form, elastic plate theory is adequate for understanding many significant features of their formation. We present a program that uses elastic theory to solve for the deflection of a thin plate resulting from a time-varying distributed load. As the plate is being deformed, the load can be expressed by any number of suitable geologic events (e.g. structural or sediment loading). We use an advancing thrust sheet as a load to explain the formation of foreland basins. The program calculates the iithospheric deflection resulting from a given load at a specific time and the history of basin formation can be shown by calculating the deflection at specific times in the thrusting history. THEORY
The flexural deflection of a thin elastic plate resulting from a distributed, time-varying load has been solved by many authors (Turcotte and Schubert, 1982; Watts, Karner, and Steckler, 1982). We solve how an elastic plate is deflected while a time-varying load in the form of an advancing thrust sheet is emplaced [e.g. Appalachian Mountains (Quinland and Beaumont, 1984) or the Idaho-Wyoming thrust belt (Jordan, 1981)]. The deflection of the lithosphere results in the formation of a foreland basin. To solve for the deflection of a uniformly thick elastic lithosphere we use a fourth-order, linear elastic equation:
D V4W(X, t) "~ (6p)gw(x, t) = P(x, t),
(1)
where
6p = (Pm - P,)
if P(x, t) 4: O,
fip=(pm--ps)
ifP(x,t)=0,
where D is the flexural rigidity of the elastic lithosphere, V4 is the biharmonic operator, w(x, t) is the amount of elastic deflection as a function of spatial coordinates (x) and time (t), Pm is the mantle density, Ps is the sediment density within the foreland basin, Pt is the thrust sheet density, g is the acceleration of gravity, and P(x, t) is the distributed and time-varying load resulting from advancement of the thrust sheets. D is related to the thickness of the lithosphere (T,) by
ET~ D = - 12(l - v~) '
(2)
where E is Young's modulus and v is Poisson's ratio. The geometry of foreland basins generally is twodimensional (Beaumont, 1981; Allen, Homewood, and Williams, 1986), so Equation (l) can be simplified to a function of only one spatial variable (e.g. x). Equation (1) then can be written as 04w D ~x 4 + (6p )gw = P(x, t),
(3)
where w = w(x, t). The load, P(x, t), can be written as
P(x, t) = ptgh(x, t),
(4)
where h(x, t) is the variable topography of the thrust sheet. Equation (3) now becomes t~4w D ~x 4 + (6p)gw = ptgh(x, t).
(5)
To simplify the determination of h(x, t), the advancing thrust sheet is modeled as a triangular shaped wedge advancing at a uniform rate v. The simplification of the thrust sheet geometry is justified because numerous published cross sections of thrust sheets [e.g. the Appalachian Mountains (Hatcher, 1989) or the Ouachita-Marathon orogeny (Ross, 1986; Viele and Thomas, 1989)] can be considered to be roughly triangular in cross section. Most thrust sheets do not advance at a uniform rate but occur in several pulses or episodes in which the advancement rate can be modeled as uniform (Hatcher, 1989; Reed and Strickler, 1990). So, our model can be applied to the specific episodes and the cumulative effects can be summed together. The expression for a triangular shaped wedge advancing at a uniform rate is given by
h(x,t)=(vt-x)tanO,
O~t<<.TM,
O~x<<.vt
(6) where TM is the duration of the thrusting event in millions of years. Note that Equation (6) satisfies t?h/dt + v Oh/Ox = 0. This condition states that the wedge thickness does not change with time at a point following the motion of the thrust sheet. We therefore are not allowing for local deformation of the advancing thrust wedge with time. Equation (5) was solved numerically using finite differences with Gauss-Seidel iteration and successive overrelaxation (SOR). The initial and boundary conditions used to solve Equation (5) are w=0 ~2W Ox 2
att=0, ~3W =0 Ox 3
Ow w=~--~x=0
O~X<.Xm, x at x = 0 , t > 0
atX=Xmax,
t>0
(7) (8) (9)
where Xm,x is the total length of the plate. PROGRAM USAGE The program F O R E L A N D BASIN as given in the Appendix was written in F O R T R A N 77 on a Digital Equipment Corporation (DEC) 5000/200 computer using the ultrix operating system version 4.2A. A sample input and output is given in the Appendix to illustrate program usage. The program consists of a main program and three subroutines (INPUT, OUTPUT, FLEX). A brief description of each subroutine is given in the main program listing.
1323
Modeling the formation of foreland basins Foreland basin model geometry
Wa
L 0 -2
= thrust sheet length = a n g l e o f thrust
h 0 = m a x i m u m load t o p o g r a p h y T T = maximum thrust wedge thickness x 0 = foreland basin width x a = arc position
-3
w 0 = f o r e l a n d b a s i n d e p th w a = arch a m p l i t u d e
-4
-5
0
50
I
I
I
I
i
L
J
k
I
100
150
200
250
300
350
400
450
5O0
Distance from edge of plate (km) Figure 1. Geometry of elastic plate (lithosphere) and distributed load (thrust sheet). The program calculates the deflection of the lithosphere along any given profile along the surface of the Earth (Fig. 1). The input parameters describing the properties of the lithosphere, the amount of the initial load, where to calculate the deflections along the profile, and factors involved in numerically solving Equation (5) are described briefly in the main program. These parameters will be discussed in more detail next. The distance units are in kilometers (km), densities are in kg/m 3, loads are in megaPascals (MPa), Young's modulus is in Pascals (Pa), and accelerations are in m/s 2. Equation (5) is solved numerically using an iterative Gauss-Seidel technique. Most iterative techniques for solving differential equations converge slowly to a solution where the most important problem becomes determining a technique to speed up the convergence (Forsythe, Malcolm, and Moler, 1977). To speed up the convergence, we used the successive overrelaxation (SOR) method. The SOR method requires the use of an acceleration parameter, o9, where o9 has the value 1 < o9 < 2. If ~o is 1, then the SOR method reverts back to the regular Gauss-Seidel technique and the convergence is slow. For any other number in the given range, the convergence is accelerated. The optimum selection of o9 is problem-dependent and is difficult to determine "a priori" (Carnahan, Luther, and Wilkes, 1969) but we have determined that to = 1.5 works best for our problem. The user may want to vary o~ to understand the complete usage of the program. The parameters im~xand dx are used to describe the grid of nodal points used to solve numerically for the
deflection of the elastic plate. The length of the elastic plate must be greater than the region of deformation. This is so that the deflection goes to zero far from the thrust sheet front. /max is the number of nodal points used to describe the elastic plate. For resolution and accuracy purposes, /max should be greater than about 60. Also, /maxmust be less than the dimensions of the program's arrays (presently set at 201). dx is the spatial increment of the nodal points and is determined by dividing the length of the elastic plate by / m a x - 1. In our work in the Ouachita orogenic belt, we use 500km for the length of the plate and /max= 101 which would make dx = 5 km. The lithospheric thickness, To, is the effective elastic thickness of the plate we are deflecting. The determination of this thickness for continental regions is not trivial because the composition of the lower crust and the lithospheric temperatures are not well known (McNutt, Diament, and Kogan, 1988). Subsequently, there are many techniques (e.g. seismic refraction, gravity modeling, admittance studies, flexural studies) that are used to determine To. In forward gravity modeling, the gravity field of the deflected plate is compared with the observed Bouguer gravity anomalies. The lithospheric thickness and density are differed until the observed and predicted gravity anomalies match. In the admittance techniques, the averaged ratios of the Bouguer gravity anomaly to topography in selected wavelength ranges are compared to those predicted by the deflected plate models. The problem with admittance studies in continental areas, unlike oceanic areas, is that the load causing the deflection may not have a
1324
J.H. HINOJOSAand K. L. MICKUS
corresponding topographic expression (Karner and Watts, 1983; Forsyth, 1985) which is a requirement for the admittance technique. Some authors (e.g. Royden and Karner, 1984) have used flexure modeling of foreland basins constrained by seismic reflection and drill core data to determine the lithospheric thickness. Another technique we have used to determine To is to invert magnetic anomaly data for Curie point depths (Hinojosa and Mickus, 1991; Mickus and Hinojosa, 1992). The Curie point depth is defined as the depth to the 550°C isotherm (e.g. Mayhew, 1982, 1985) assuming the magnetic mineral is magnetite. Once the Curie point depths were calculated, Te was obtained by linear interpolation to the 450°C isotherm. This is justified because the geotherm at shallow depths (crustal depths) increases in a nearly linear fashion. Even though there are numerous techniques to estimate Te one cannot be sure that the value being used is correct. The best technique may be to estimate Te by more than one technique and average the values obtained. The maximum load (load 0) used in the program is determined "a posteriori" in a trial and error iterative scheme so as to satisfy the constraint imposed by the maximum thrust wedge thickness ( T T ) obtained from cross sections (cf. Fig. 1). The parameter T T is related to the maximum plate deflection [w (x -- 0)] and the maximum load topography (h0) by T T = w ( x = O) + h o. Note that the maximum load topography (h0) is that part of the thrust sheet that is above the plane z = 0, where z = 0 is the surface of the Earth. The parameter load 0 then can be calculated from
(1989) and Reed and Strickler (1990). The maximum thrust length and thickness we estimated are 50 and 6 km, respectively. Frankel and others (1988) estimated the duration of thrusting to be 40-50 Myr. Reed and Strickler (1990) recognized at least four distinct thrusting events with periods of no deformation between events in the Marathon Uplift. Even so, to simplify our models we will assume that the thrust rate was constant. The densities we used were determined from gravity modeling (Smith, 1986). The densities for the thrust sheet, basin sediments, and the mantle are 2820, 2600, and 3300 kg/m 3, respectively. We used two lithospheric thickness: (1) 12 km which was estimated from Curie point depth calculations (Hinojosa and Mickus, 1991), and (2) 29 km which was estimated from gravity modeling (Smith, 1986). We used 500 km as the length of the elastic plate, because published maps of the Val Verde Basin region (Frankel and others, 1988) show that the basin is narrower than 500 km. All the parameters are shown in Table 1. We calculated the deflection of the lithosphere using elastic thicknesses of both 12 and 29 km at 10Myr intervals for a total duration of 40 Myr (Figs. 2 and 3). Using a lithospheric thickness of Table 1. Flexure model parameters
Symbol
Value
Length of model
Xmax
500 km
Effective elastic thickness
T,
12 km, 29 km
where h0 = T T - w (x = 0). Most of the parameter values used in the program have to be estimated and depend on the basin being modeled. The thrust sheet length (L) and the thickness ( T T ) can be estimated from published geological data and cross sections. To model the basin's formation, the user can estimate the thrust sheet's rate of advancement, determine the thrust length and thickness, and subsequently the loads for specific time periods during the thrusting history.
Maximum thrust
L
50 km
Maximum thrust sheet thickness
T]"
6 km
Sediment density
p,
2600 Kg m'3
Thrust sheet density
p,
2820 Kg m'3
Mantle density
p.
3300 Kg m a
Total duration of thrusting
At
40 My
EXAMPLE
Average rate of advance of thrust sheet
v
1.25 km My a
Acceleration of
g
9.8 m s 2
Young's modulus
E
0.7x10 u Pa
Poisson's ratio
v
0.25
SOR acceleration parameter
oJ
1.5
Number of nodal points
imax
101
Convergence parameter
eps
.003
load 0 = ptgho,
The application of the computer program is illustrated for the Pennsylvanian-Permian Ouachita orogenic event in West Texas. Specifically, we will show some preliminary results for formation of the Val Verde Basin which is the southeastern part of the oil-rich Permian Basin. There have been many geological studies (mainly for petroleum exploration) of the Val Verde Basin but few geophysical, and no flexural studies. The geological studies are summarized by Frankel and others (1988). We calculated maximum thrust lengths and thicknesses from cross sections through the basin by Nicholas and Waddell
Parameter
length
gravity
1325
Modeling the formation of foreland basins Foreland basin evolution
l 0 Myr
..~-"U'~"~'~:Z:'--_-..'~-~ ,~
20 Myr / ,
-I
E
_// / "t
-
/ // ]/
30 l~yr / " /
T e = 12 km
-2 ..' 40 Myr ,' /
-3
/ / /
-4
-5 0
I
J
I
I
I
I
I
I
J
I
50
100
150
200
250
300
350
400
450
500
Distance from edge of plate (km) Figure 2. Deflection of 12-km thick lithosphere in 10 Myr increments resulting from thrust loads emplaced on lithosphere during Ouachita orogenic event in Val Verde Basin region of West Texas. 12km (Fig. 2), the maximum plate deflection at 40 Myr is about 4.7 km. This should not be confused with the basin depth. The basin depth is the plate deflection at the thrust wedge front, which is approximately 1.6 km for the 12-km thick plate at 40 Myr. The corresponding basin width is about 53 km. There
is a forebulge, or arch, positioned at about 140 km from the edge of the plate, or about 40 km from the zero crossover point, with an amplitude of about 0.21 km. Using a thicker lithosphere (Fig. 3), the basin depth at 40 Myr is about 2.4 kin, with a basin width of about 120 km. The arch is located about
Foreland basin evolution
~ - . ' 7 . . v/ / 20 Myr
..- / /
..
f//
..
;.>
j
-1
E
30 Myr -2
_.""
/
/
.,'
T e = 29 km
// 40 Myr
/ -3
-
/
/ / -4
-5
_
I
I
50
I00
150
I
I
I
L
200
250
300
350
i
I
I
400
450
500
Distance from edge of plate ( k i n ) Figure 3. Deflection of 29-km thick lithosphere in 10 M y r increments resulting from thrust loads emplaced on lithosphere during Ouachita orogenic event in Val Verde Basin region of West Texas.
1326
J . H . HINOJOSA and K. L. MICKUS Table 2. Deflection of lithosphere resulting from various thrust loads. Parameters include basin width (x0), arch position (xa), basin depth (w0), arch amplitude (wa), lithospheric thickness (Te), thrust length (L), and time (My). All distance units are in kilometers and loads are in megaPascals Te
Load0
L
Wo
Wo
xo
xo
Time
12
8.75
12.5
-0.27
0.02
69.6
37.9
10
12
17.50
25.0
-0.78
0.07
62.6
39.1
20
12
26.25
37.5
-1.30
0.14
57.6
39.0
30
12
35.00
50.0
-1.64
0.21
52.5
38.7
40
29
13.25
12.5
-0.21
0.01
134.2
58.3
10
29
26.50
25.0
-0.76
0.06
131.5
73.5
20
29
39.75
37.5
-1.55
0.13
126.4
76.1
30
29
53.00
50.0
-2.35
0.22
120.1
74.9
40
75 k m from the zero crossover point, a n d has a n amplitude o f a b o u t 0.22 km. Table 2 shows the loads, basin widths, arch positions, a n d arch amplitudes as a function o f time. The input a n d o u t p u t for the situation o f T e = 12 k m a n d L = 50 k m is shown in the Appendix. CONCLUSIONS The p r o g r a m described is useful for modeling the subsidence of a foreland basin formed by loading of thrust sheets o n an elastic lithosphere. Even t h o u g h the p r o g r a m calculates the deflection at a specific time, the basin's subsidence history can be determined by k n o w i n g the loads at various times t h r o u g h o u t the orogenic event. The use o f the p r o g r a m is a first step in modeling the f o r m a t i o n of foreland basins. I m p r o v e m e n t s to the technique would be to allow variable elastic thickness, t o p o g r a p h y on the thrust sheet, erosion of the thrust sheet a n d subsequent sediment deposition into the basin, the i n c o r p o r a t i o n of viscoelastic bending theory, a n d the possible inclusion of two-dimensional effects. REFERENCES
Allen, P. A., and Allen, J. R., 1990, Basin analysis: principles and applications: Blackwell Scientific Publications, Oxford, 451 p. Allen, P. A., Homewood, P., and Williams, G. D., 1986, Foreland basins: an introduction, in Allen, P. A., and Homewood, P., eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, Blackwell Scientific, Oxford, p. 3-12. Aubouin, J., 1965, Geosynclines, in Developments in geotectonics: Elsevier, Amsterdam, 335 p. Beaumont, C., 1978, The evoluation of sedimentary basins on a viscoelastic lithosphere: theory and examples: Geophys. Jour. Roy. Astr. Soc., v. 55, no. 2, p. 471~497.
Beaumont, C., 1981, Foreland basins: Geophys. Jour. Roy. Astr. Soc., v. 65, no. 2, p. 291-329. Carnahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied numerical methods: John Wiley & Sons, New York, 604 p. Courtney, R. C., and Beaumont, C., 1983, Thermally activated creep and flexure of the oceanic lithosphere: Nature, v. 305, no. 5931, p. 201-204. Dickinson, W. R., 1974, Plate tectonics and sedimentation, in Dickinson, W. R., ed., Tectonics and sedimentation: Soc. Economic Paleontologists and Mineralogists, Spec. Publ. vol. 22, p. 1-27. Fleming, P. B., and Jordan, T. E., 1990, Stratigraphic modeling of foreland basins: interpreting thrust deformation and lithosphere rheology: Geology, v. 18, no. 5, p. 430-434. Frankel, H. N., and others, 1988, The Permain Basin region, in Sloss, L. L., ed., Sedimentary cover-North American Cration: U.S.: The Decade of North American Geology, v. D-2, Geol. Soc. America, p. 261 306. Forsyth, D. W., 1985, Subsurface loading and estimates of the flexural rigidity of continental lithosphere: Jour. Geophys. Res., v. 90, no. BI4, p. 12623-12632. Forsythe, G. E., Malcolm, M. A., and Moler, C. B., 1977, Computer methods for mathematical computations: Prentice-Hall Inc., New York, 259 p. Hatcher, R. D., 1989, Tectonic synthesis of the U.S. Appalachians, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian-Ouachita Orogen in the United States: The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 511-536. Hinojosa, J. H., and Mickus, K. L., 1991, Comparative geophysical study of the Ouachita foreland basins in Texas: EOS: Trans. Am. Geophys. Union, v. 72, no. 44, p. 430. Houseknecht, D. W., 1986, Evolution from passive margin to foreland basin: the Atoka Formation of the Arkoma basin, south-central U.S.A., in Allen, P. A., and Homewood, P.. eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, Blackwell Scientific, Oxford, p. 327 345. Jordan, T. E., 1981, Thrust loads and foreland basin evolution, Cretaceous, western United States: Am. Assoc. Petroleum Geologists Bull. v. 65, no. 12, p. 250(~2520.
Modeling the formation of foreland basins Karner, G. D., and Watts, A. B., 1983, Gravity anomalies and flexure of the lithosphere at mountain ranges: Jour. Geophys. Res. v. 88, no. B1 I, p. 10449-10477. Kusznir, N., and Karner, G. D., 1985, Tectonic subsidence calculated from lithified basin strata (abst.): Geol. Soc. America, Abstracts with Programs, v. 14, no. 5, p. 534. Mayhew, M. A., 1982, Application of satellite magnetic anomaly data to Curie isotherm mapping: Jour. Geophys. Res., v. 87, no. 6, p. 48464854. Mayhew, M. A., 1985, Curie isotherm surfaces inferred from high-altitude magnetic anomaly data: Jour. Geophys. Res., v. 90, no. 3, p. 2647 2654. McNutt, M. K., Diament, M., and Kogan, M. G., 1988, Variations of elastic plate thickness at continental thrust belts: Jour. Geophy. Res., v. 93, no. B8, p. 8825-8838. Mickus, K. L., and Hinojosa, J. H., 1992, Flexure origin of the Kerr basin, Texas (abst.): EOS: Trans. Am. Geophys. Union, v. 73, no. 24, p. 289. Nicholas, R. L., and Waddell, D. E., 1989, The Ouachita system in the subsurface of Texas, Arkansas, and Louisiana, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian~)uachita Orogen in the United States:-The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 661~571. Price, R. A., 1973, Large-scale gravitational flow of supracrustal rocks, southern Canadian Rockies, in DeJong, K. A., and Schotten, R., eds., Gravity and tectonics: John Wiley & Sons, New York, p. 491-502. Quinlan, G. M., and Beaumont, C., 1984, Appalachian thrusting, lithospheric flexure, and the Paleozoic stratigraphy of the Eastern Interior of North America: Can. Jour. Earth Sciences, v. 21, no. 9, p. 973-995. Reed, T. A., and Strickler, D. L., 1990, Structural geology and petroleum exploration of the Marathon thrust belt, West Texas, in Laroche, T. M., and Higgins, L., eds., Marathon thrust belt: structure, stratigraphy, and hydrocarbon potential: West Texas Geol. Soc. Field Seminar, p. 39-82. Ross, C. A., 1986, Paleozoic evolution of southern margin of Permian Basin: Geol. Soc. America Bull., v. 97, no. 5, p. 536-554.
Royden, L. H., and Karner, G. D., 1984, Flexure of the continental lithosphere beneath the Appenines and Carpathians foredeep basins: evidence for an insutficient topographic load: Am. Assoc. Petroleum Geologists Bull., v. 68, no. 6, p. 704-712. Sleep, N. H., and Snell, N. S., 1976, Thermal contraction and flexure of mid-continent and Atlantic marginal basins: Geophys. Jour. Roy. Astr. Soc., v. 45, no. 1, p. 125-154. Smith, K. J., 1986, A gravity and tectonic study of the southwestern portion of the Ouachita system: unpubl. masters thesis, Univ. Texas, E1 Paso, 99 p. Steckler, M. S., and Watts, A. B., 1978, Subsidence of the Atlantic-type continental margin off New York: Earth and Planet. Sci. Lett., v. 41, no. 1, p. 1 13. Stockmal, G. S., Beaumont, C., and Boutilier, R., 1986, Geodynamic models of convergent margin tectonics: transition from rifted margin to overthrust belt and consequences for foreland basin development: Am. Assoc. Petroleum Geologists Bull., v. 70, no. 2, p. 181-190. Turcotte, D. L., and Schubert, G., 1982, Geodynamics: applications of continuum mechanics to geological problems: John Wiley & Sons, New York, 450 p. Viele, G. W., and Thomas, W. A., 1989, Tectonic synthesis of the Ouachita orogenic belt, in Hatcher, R. D., Thomas, W. A., and Viele, G. W., eds., The Appalachian-Ouachita Orogen in the United States: The Decade of North American Geology, v. F-2, Geol. Soc. America, p. 695-728. Watts, A. B., Karner, G. D., and Steckler, M. S., 1982, Lithospheric flexure and the evolution of sedimentary basins: Phil. Trans. Roy. Soc. London, v. 305A, no. 2, p. 249-281. Wuellner, D. E., Lehtonen, L. R., and James, W. C., 1986, Sedimentary-tectonic development of the Marathon and Val Verde basins. West Texas, U.S.A.: a PermianCarboniferous migrating foredeep, in Allen, P. A., and Homewood, P., eds., Foreland basins: Intern. Assoc. Sedimentologists, Spec. Publ., v. 8, BlackweU Scientific, Oxford, p. 327 345.
APPENDIX Program Listing ******************************************************************** PROGRAM
FORELAND
BASIN
PROGRAM TO COMPUTE THE FLEXURAL DEFORMATION OF AN ELASTIC PLATE PLATE DUE TO A DISTRIBUTED, TIME VARYING THRUST SHEET THAT WILL FORM A FORELAND BASIN. THE PROBLEM IS SOLVED NUMERICALLY USING A ONE-DIMENSIONAL, FOURTH-ORDER, LINEAR ELASTIC FLEXURE EQUATION SUBJECT TO A DISTRIBUTED LOAD, WITH A UNIFORM EFFECTIVE ELASTIC THICKNESS USING THE GAUSS-SEIDEL ITERATION METHOD WITH SUCCESSIVE OVER-RELAXATION (SOR).
INPUT PARAMETERS : RHOM RHOT RHOS L
= = = =
EPS
=
E G NU OMEGA
= = = =
1327
MANTLE DENSITY (ALL DENSITIES ARE IN KG/M**3) THRUST-WEDGE DENSITY SEDIMENT DENSITY WITHIN THE FORELAND BASIN THRUST DISTANCE IN THE DIRECTION OF THE ADVANCING WEDGE (KILOMETERS) CONVERGENCE PARAMETER FOR THE GAUSS-SEIDEL ITERATION METHOD (FOR BEST RESULTS WE USE 1.E-3 TO I.E-4) Y O U N G ' S M O D U L U S (Pa) GRAVITATIONAL ACCELERATION (M/S**2) POISSON'S RATIO (DIMENSIONLESS) SOR ACCELERATION PARAMETER, MUST BE BETWEEN 1 A N D 2. W E U S U A L L Y U S E 1 . 5 B U T T H E S E L E C T I O N IS TRIAL AND ERROR DEPENDING ON THE PROBLEM
328
J . H . HINOJOSA and K. L. MICKUS = MAXIMUM NUMBER OF NODAL POINTS. FOR RESOLUTION C IMAX c T H E N U M B E R S H O U L D B E > 60. B U T I T M U S T B E L E S S THAN PROGRAM'S ARRAYS DIMENSION SIZE (PRESENTLY SET C C AT 201). IF ONE WANTS MORE THAN 201 THE DIMENSION C OF THE PROGRAM'S ARRAYS MUST BE INCREASED. c = SPATIAL INCREMENT OF NODAL POINTS (KILOMETERS) DX c = EFFECTIVE ELASTIC THICKNESS OF THE LITHOSPHERE TE c (KILOMETERS) C LOAD0 = M A X I M U M T H R U S T W E D G E L O A D (MPa). C c OUTPUT PARAMETERS : c c W(I) = DEFLECTION OF THE ELASTIC PLATE AS A FUNCTION OF c DISTANCE. (KILOMETERS) C DI S T (I ) = DISTANCE FROM COMPUTATIONAL ORIGIN (KILOMETERS) C L O A D (I ) = D I S T R I B U T E D L O A D (MPa) c C OTHER INTERNAL PARAMETERS: C = FLEXURAL RIGIDITY (NEWTON-M) C D = DEFLECTION OF THE ELASTIC PLATE AS A FUNCTION OF C Wl(I) DISTANCE FROM THE PREVIOUS ITERATION C = LENGTH OF THE ELASTIC PLATE, XMAX=DX*(IMAX-1) C XMAX C SUBROUTINES USED: C c INPUT - READS IN THE VARIOUS PARAMETERS C c OUTPUT - WRITES OUT THE VARIOUS PARAMETERS, THE LOAD AS A c FUNCTION OF DISTANCE, AND THE DEFLECTION OF THE c LITHOSPHERE AS A FUNCTION OF DISTANCE c FLEX - CALCULATES THE DEFLECTION OF THE LITHOSPHERE c DUE TO AN ADVANCING THRUST SHEET c c WRITE UNITS: C C 52 - T H E I N P U T P A R A M E T E R S A N D T H E D E F L E C T I O N O F T H E C LITHOSPHERE AS A FUNCTION OF DISTANCE. c 53 - T H E L O A D D U E T O T H E A D V A N C I N G T H R U S T S H E E T A S C A FUNCTION OF DISTANCE. c 6 - QUESTIONS ASKING FOR THE FILE CONTAINING THE C INPUT PARAMETERS, AND OUTPUT FILES THAT CONTAIN c CALCULATED DEFLECTION AND LOAD. ALSO A BANNER c A N N O U N C I N G T H E N A M E O F T H E P R O G R A M IS P R I N T E D . c READ UNITS : C c c 51 - T H E I N P U T P A R A M E T E R S c 5 - T H E N A M E S O F T H E I N P U T (OLD) A N D O U T P U T (NEW) F I L E S c c THIS PROGRAM WAS WRITTEN ON A DEC 5000/200 USING AN ULTRIX c OPERATING SYSTEM VERSION 4.2A c ******************************************************************** REAL W(201),WI(201),LOAD(201),L,DIST(201),LOAD0,OMEGA,NU,M CHARACTER*20 FNAMEI,FNAME2,FNAME3 CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD C WRITE(6, ****************************************************** WRITE(6, *)'* F O R E L A N D B A S I N -- C O M P U T E S T H E F L E X U R A L *' W R I T E ( 6 , *) '* DEFORMATION OF AN ELASTIC PLATE TO FORM A *' WRITE(6, *)'* FORELAND BASIN *' WRITE(6, ****************************************************** WRITE(6, *)'ENTER INPUT DATA FILE NAME:' R E A D ( 5 , ' (A) ' ) F N A M E I O P E N (51, F I L E = F N A M E I , S T A T U S = 'O L D ' ) WRITE(6,*) 'ENTER OUTPUT DATA FILE NAME:' R E A D (5, " (A) ' ) F N A M E 2 O P E N (52, F I L E = F N A M E 2 , S T A T U S = ' N E W ' ) W R I T E ( 6 , * ) 'DO Y O U W A N T T O W R I T E O U T T H E L O A D ( Y / N ) ? ' R E A D ( 5 , ' (A) ' ) A N S IF ( A N S . E Q . ' Y ' ) T H E N W R I T E (6, *) 'E N T E R N A M E O F L O A D F I L E : ' R E A D (5, " (A) ' ) F N A M E 3 O P E N (53, F I L E = F N A M E 3 , S T A T U S = ' N E W ' )
Modeling the formation of foreland basins ELSE CONTINUE ENDIF ******************************************************************** C
READ
IN T H E P A R A M E T E R S
CALL INPUT ******************************************************************** C D E T E R M I N E T H E D E F L E C T I O N OF T H E E L A S T I C P L A T E ******************************************************************** CALL FLEX ******************************************************************** C W R I T E O U T T H E C O M P U T E D D E F L E C T I O N A N D LOAD, C AND THE INPUT PARAMETERS ******************************************************************** CALL OUTPUT C C5OSE(51) CLOSE(52)
CLOSE(53) C STOP END ******************************************************************** ******************************************************************** SUBROUTINE INPUT ******************************************************************** C S U B R O U T I N E T O I N P U T P A R A M E T E R S U S E D TO D E S C R I B E T H E R H E O L O G I C A L C P R O P E R T I E S OF A T H I N E L A S T I C P L A T E A N D T H E D I S T R I B U T E D L O A D C (SEE T H E M A I N P R O G R A M F O R A D E S C R I P T I O N OF T H E S E PARAMETERS) ******************************************************************** R E A L L, LOAD0, OMEGA, NU, M CHARACTER* 1 ANS C O M M O N / P A R M / L O A D 0 , DX, IMAX, D, OMEGA, EPS, M, L, TE, A N S C O M M O N / D E N S / R H O M , RHOT, RHOS, G READ(51,*) READ(51,*) READ(51,*)
E,NU,OMEGA,IMAX,L RHOM,RHOT,RHOS,G DX,TE,LOAD0,EPS
DX=DX*I.E3 LOAD0=LOAD0*I.E6 M = (LOAD0/(L*I.E3)) TE=TE*I.E3 ******************************************************************* C CALCULATE THE FLEXURAL RIGIDITY ******************************************************************* D = E*(TE**3)/(12.*(I.-(NU)**2)) C RETURN END ******************************************************************** ******************************************************************** SUBROUTINE OUTPUT ******************************************************************** C S U B R O U T I N E TO W R I T E O U T T H E I N P U T PARAMETERS, D I S T R I B U T E D LOAD C AND THE DEFLECTION OF A THIN ELASTIC PLATE ******************************************************************** REAL W(201),LOAD(201),DIST(201),LOAD0,L CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD C LOAD0=LOADO/I.E6 TE=TE/I.E3 WRITE(52,100) TE,LOAD0,L,RHOM,RHOT,RHOS IF(ANS.EQ.'Y') T H E N WRITE(53,200) WRITE(53,'(FIO.3,5X,E9.4)') (DIST(I),LOAD(I),I=I,IMAX) ELSE CONTINUE ENDIF WRITE(52,300) W R I T E ( 5 2 , ' ( 2 ( I X , F I O . 3 ) ) ' ) (DIST(I),W(I),I=I,IMAX) C i00 F O R M A T ( ' T H E L I T H O S P H E R I C T H I C K N E S S IS ', F5.1,/, 1"THE M A X I M U M L O A D IN M P A IS ', F8.3,/,
1329
1330
J . H . HINOJOSA and K. L. MICKUS
200 300 C
2 ' T H E M A X I M U M T H R U S T D I S T A N C E IS " , F 5 . 1 , / , 3'THE DENSITIES OF THE MANTLE, THRUST SHEET AND 4 S E D I M E N T S A R E ", 3 ( F 6 . 1 , 1 X ) , / ) FORMAT(3X,'DISTANCE',8X,'LOAD') FORMAT(3X,'DISTANCE',3X,'DEFLECTION')
RETURN END ******************************************************************** ******************************************************************** SUBROUTINE FLEX ******************************************************************** C DETERMINES THE FLEXURAL DEFLECTION OF A UNIFORM, THIN C ELASTIC PLATE DUE TO A TIME VARYING, DISTRIBUTED LOAD C (SEE T H E M A I N P R O G R A M F O R A D E S C R I P T I O N O F T H E P A R A M E T E R S ) ******************************************************************** C REAL W(201),WI(201),LOAD(201),L,DIST(201),LOAD0,OMEGA,M CHARACTER*I ANS COMMON/PARM/LOAD0,DX,IMAX,D,OMEGA,EPS,M,L,TE,ANS COMMON/DENS/RHOM,RHOT,RHOS,G COMMON/DEFL/W,DIST,LOAD ******************************************************************** C DEFINE THE DISTRIBUTED LOAD ********************************************************************
10
DO i0 I I 0 = I , I M A X DIST(II0)=DX*FLOAT(II0-1)*I.E-3 IF ( D I S T ( I 1 0 ) . L E . L ) T H E N LOAD(I10)=LOAD0-M*(FLOAT(I10-1)*DX) ELSE LOAD(II0)=0. ENDIF CONTINUE
C BETA1 = D/(DX**4.) BETA2 = (RHOM - RHOT)*G B E T A 3 = (RHOM - RHOS)*G ******************************************************************** C START GAUSS-SEIDEL ITERATION TO COMPUTE DEFLECTION W(I): DO 20 I 2 0 = I , I M A X W ( I 2 0 ) = 0. 20 CONTINUE ******************************************************************** C RETURN HERE FOR NEXT ITERATION: ******************************************************************** 30
D O 40 I 4 0 = I , I M A X WI(I40) = W(I40) 40 CONTINUE ******************************************************************** C C
BOUNDARY CONDITIONS: W'' = W ' ' ' = 0 A T X = 0 (AT M A X I M U M
LOAD)
W(1) = ( I . / ( 2 . * B E T A I + B E T A 2 ) ) * ( L O A D ( 1 ) & - 2.*BETAI*W(3) + 4.*BETAI*W(2)) ******************************************************************** C A P P L Y SOR: ******************************************************************** = (i. - OMEGA)*WI(1) + OMEGA*W(1) IF (LOAD(2) .NE. 0) T H E N W(2) = ( I . / ( 5 . * B E T A I + B E T A 2 ) ) * ( L O A D ( 2 ) & - BETAI*W(4) + 4.*BETAI*W(3) & + 2.*BETAI*W(1)) ELSE W(2) = ( I . / ( 5 . * B E T A I + B E T A 3 ) ) * ( L O A D ( 2 ) & - EETAI*W(4) + 4.*EETAI*W(3) & + 2.*BETAI*W(1)) ENDIF ******************************************************************** W(1)
C A P P L Y SOR: ******************************************************************** W(2)
=
(I.
-
OMEGA)*WI(2)
+ OMEGA*W(2)
C DO 50 I 5 0 = 3 , I M A X - 2 IF(LOAD(I50) .ME. 0.) THEN W(I50) = (I./(6.*BETAI+BETA2))*(LOAD(I50)-BETAI*W(I50-2) & + 4.*BETAI*W(I50-1) + 4.eBETAI*W(IS0+I) & BETAI*W(I50+2))
Modeling the formation of foreland basins ELSE W(I50)
=
& &
133
(1./(6.*BETAI+BETA3))*(LOAD(I50)-BETAI*W(I50-2) + 4.*BETAI*W(I50-1) + 4.*BETAI*W(I50+I) BETAI*W(I50+2))
ENDIF ******************************************************************** C
APPLY
SOR:
******************************************************************** W ( I 5 0 ) = (1. - O M E G A ) * W l ( I 5 0 ) + O M E G A * W ( I 5 0 ) 50 CONTINUE ******************************************************************** C C
BOUNDARY CONDITIONS: W = W' = 0 A T X = X M A X
(FAR FROM THRUST
WEDGE
FRONT)
******************************************************************** W(IMAX-I)
=
&
(I./(7.*BETAI+BETA3))*(LOAD(IMAX-I) + 4.*BETAI*W(IMAX-2) - BETAI*W(IMAX-3))
******************************************************************** C
APPLY
SOR:
******************************************************************** W ( I M A X - I ) = (1. - O M E G A ) * W I ( I M A X - 1 ) + O M E G A * W ( I M A X - I ) ******************************************************************** C
COMPARE
WITH
PREVIOUS
ITERATION
TO TEST
FOR CONVERGENCE:
******************************************************************** D O 60 I 6 0 = I , I M A X IF (ABS(W(I60)-WI(I60)).GT.EPS) CONTINUE D O 70 I T 0 = I , I M A X W(I70) = -W(I70)*I.E-3 CONTINUE
60
70
GOTO
30
C RETURN END COMMAND ! ! ! ! ! !
IN.FLEXURE OUT.DEFLECTION ¥ OUT.LOAD
FILE
FOR RUNNING
FORELAND
BASIN
INPUT FILE DESCRIBING THE RHEOLOGICAL PROPERTIES PLATE AND THE DISTRIBUTED LOAD OUTPUT FILE CONTAINING THE PLATE DEFLECTION AS A OF DISTANCE YES OR NO ANSWER IF USER WANTS A FILE CONTAINING OUTPUT FILE CONTAINING THE LOAD AS A FUNCTION OF
OF THE FUNCTION THE LOAD DISTANCE
INPUT .7eli
.25 1 . 5
3300.
2820.
2600.
35.0
1.E-3
5.
12.
i01
50.
9.80
l Y O U N G ' S M O D U L U S (E), P O I S S O N ' S R A T I O (NU), S O R ! ACCELERATION PARAMETER (OMEGA), MAXIMUM NUMBER OF ! N O D A L P O I N T S (IMAX), T H R U S T D I S T A N C E (L) ! D E N S I T I E S O F M A N T L E (RHOM), T H R U S T - W E D G E (RHOT), ! AND FORELAND BASIN SEDIMENTS (RHOS), GRAVITATIONAL ! C O N S T A N T (G) ! S P A T I A L I N C R E M E N T O F N O D A L P O I N T S (DX), E L A S T I C ! T H I C K N E S S (TE), M A X I M U M L O A D ( L O A D 0 ) , C O N V E R G E N C E ! P A R A M E T E R (EPS) OUTPUT
THE THE THE THE
LITHOSPHERIC THICKNESS IS 12.0 M A X I M U M L O A D I N M P A IS 35.000 M A X I M U M T H R U S T D I S T A N C E IS 50.0 DENSITIES OF THE MANTLE, THRUST SHEET AND DISTANCE 0.000 5.000 10.000 15.000 20.000 25.000 30.000 35.000 40.000 45.000 50.000 55.000 60.000 65.000
CAGEO 19,9--1
DEFLECTION -4.736 -4.405 -4.075 -3.746 -3.420 -3.098 -2.784 -2.479 -2.186 -1.906 -1.643 -1.398 -1.172 -0.966
SEDIMENTS
ARE
3300.0
2820.0
2600.0
1332
J . H . HINOJOSAand K. L. MICKUS 70.000 75.000 80.000 85.000 90.000 95.000 100.000 105.000 Ii0.000 115.000 120.000 125.000 130.000 135.000 140.000 145.000 150.000 155.000 160.000 165.000 170.000 175.000 180.000 185.000 190.000 195.000 200.000 205.000 210.000 215.000 220.000 225.000 230.000 235.000 240.000 245.000 250.000 255.000 260.000 265.000 270.000 275.000 280.000 285.000 290.000 295.000 300.000 305.000 310.000 315.000 320.000 325.000 330.000 335.000 340.000 345.000 350.000 355.000 360.000 365.000 370.000 375.000 380.000 385.000 390.000 395.000 400.000 405.000 410.000 415.000 420.000 425.000 430.000 435.000 440.000 445.000 450.000 455.000
-0.780 -0.613 -0.466 -0.337 -0.225 -0.129 -0.049 0.018 0.073 0.116 0. 149 0.174 0.191 0.201 0.206 0.206 0.203 0.196 0.187 0.177 0.165 0.152 0.139 0 125 0 112 0 099 0 087 0 075 0 065 0 055 0 045 0 037 0.029 0.023 0.017 0.012 0.007 0.003 0.000 -0.002 -0.004 -0.006 -0.007 -0.008 -0.008 -0.009 -0.009 -0.009 -0.009 -0.008 -0.008 -0.007 -0.007 -0.006 -0.006 -0 0 0 5 - 0 004 -0 004 -0 003 -0 003 -0 003 -0.002 -0.002 -0.001 -0.001 -0.001 -0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
460.000 465.000 470.000 475.000 480.000 485.000 490.000 495.000 500.000
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
DISTANCE LOAD 0.000 .3500E+08 5.000 .3150E+08 10.000 .2800E+08 15.000 .2450E+08 20.000 .2100E+08 25.000 .1750E+08 30.000 .1400E+08 35.000 .I050E+08 40.000 .7000E+07 45.000 .3500E+07 50.000 .0000E+00 55.000 .0000E+00 60.000 .0000E+00 65.000 .0000E+00 70.000 .0000E+00 75.000 .0000E+00 80.000 .0000E+00 85.000 .0000E+00 90.000 .0000E+00 95.000 .0000E+00 100.000 .0000E+00 FOR All THE OTHER LOCATIONS
THE LOAD
IS 0 . 0