Step 13: Compute the structure of the lithosphere
Program Shells will require at least 4 nodal data values at each node: elevation, heat-flow, crustal thickness, and mantle lithosphere thickness. (If you look at a new .feg file produced by OrbWin, these values are initially all zero.) Elevation and heat-flow can be interpolated from the .grd files produced in the previous step. My recommendation is start simply, by computing the layer thicknesses through use of 4 assumptions:
The
base of the lithosphere is defined as an isothermal surface.
Heat
conduction through the lithosphere is approximately in steady-state.
All
lithosphere should be isostatic with respect to mid-ocean ridges.
Crust
and mantle lithosphere can each be approximated as laterally homogeneous.
If you agree with these assumptions/approximations, you may compute the necessary layer thicknesses using my older program OrbData, which will interpolate elevation and heat-flow from your two .grd files, and then iteratively solve for crust and mantle lithosphere thickness to achieve isostasy and steady-state heat-flow. This program takes its physical parameters (e.g., density, thermal expansion, thermal conductivity, radioactivity) from an ASCII text file called the input parameter (.in) file, such as the one below:
i2000-01.in = 97001 + Shells 1998 + PB1999 (Earth2P.feg) 0.03 FAULT FRICTION COEFFICIENT 0.85 CONTINUUM FRICTION COEFFICIENT 1.00 BIOT COEFFICIENT (EFFICACY OF PORE PRESSURE) 0.00 BYERLY (0.-.99); FRACTIONAL STRENGTH REDUCTION OF MASTER FAULT 2.3E9, 9.5E4 ACREEP (SHEAR STRESS COEFFICIENT OF CREEP LAW) 4000., 18314. BCREEP (ACTIVATION ENERGY/N/GAS-CONSTANT) (IN KELVIN) 0., 0.0171 CCREEP (DERIVITIVE OF BCREEP WITH RESPECT TO DEPTH; CRUST/MANTLE) 5.E8, 5.E8 DCREEP (MAXIMUM SHEAR STRESS; CRUST/MANTLE) 0.333333 ECREEP (EXPONENT ON STRAIN-RATE IN CREEP-STRESS LAWS)=(1/N) 1412. 6.1E-4 INTERCEPT AND SLOPE OF UPPER MANTLE ADIABAT (K, K/M) 400.E3 ZBASTH: DEPTH OF BASE OF ASTHENOSPHERE AF PLATE WHICH DEFINES VELOCITY REFERENCE FRAME (AF, EU, NA, ...) 4 1.10 ICONVE:0=NONE;1=HAGER/O'CONNELL;2=BAUMGARDNER;3=NV1;4=CONTINENTAL NV1 1.E8 TRHMAX (LIMIT ON BASAL TRACTION) 2.5E+12 TAUMAX (DOWN-DIP INTEGRAL OF SUBDUCTION ZONE TRACTION) 1032. RHOH2O (DENSITY OF WATER, AT P=0 AND T=SURFACE TEMPERATURE) 2816., 3332. RHOBAR (MEAN DENSITY AT P=0 AND T=0; CRUST/MANTLE) 3125. RHOAST (DENSITY OF ASTHENOSPHERE, AT P=0 AND AMBIENT T) 9.8 GRAVITATIONAL ACCELERATION 1000. ONE KILOMETER, EXPRESSED IN CURRENT LENGTH UNITS 6371000. RADIUS OF THE PLANET 2.4E-5, 3.94E-5 VOLUMETRIC THERMAL EXPANSION COEFFICIENT; CRUST/MANTLE 2.7, 3.20 THERMAL CONDUCTIVITY; CRUST/MANTLE 7.27E-7, 3.2E-8 RADIOACTIVE HEAT PRODUCTION, ON VOLUME (NOT MASS) BASIS 273. SURFACE TEMPERATURE, IN KELVIN 1223., 1673. UPPER TEMPERATURE LIMITS, IN KELVIN; CRUST/MANTLE-LITHOSPHERE |
50 MAXIMUM NUMBER OF ITERATIONS 0.0001 ACCEPTABLE CONVERGENCE LEVEL (FRACTIONAL VELOCITY CHANGE) 50.E6 REFERENCE LEVEL OF SHEAR STRESS 1.00E-10 ACCEPTABLE LEVEL OF VELOCITY ERRORS (1 MM/A = 3.17E-11 M/S) F OUTPUT NODE VELOCITIES EVERY ITERATION? (FOR CONVERGENCE STUDIES) |
===== POST-PROCESSING PLOT CONTROL PARAMETERS (NOT USED BY SHELLS) ===== 999 PLOTFILE # ON INPUT DATASET (OR 999 FOR LAST, OR 0 FOR ALL) F --------------------- PLOT GRID OF ELEMENTS F 1000. -2000.+1 PLOT ELEVATIONS F 10.E-3 80.E-3+1 PLOT HEAT-FLOW F 5.E+3 +1 PLOT CRUSTAL THICKNESS F 20.E+3 80.E+3-1 PLOT TOTAL LITHOSPHERE THICKNESS F +1 PLOT TEMPERATURE OF MOHO F +1 PLOT TEMPERATURE OF BASE OF PLATE F +1 PLOT NONLITHOSTATIC PRESSURE ANOMALY AT BASE F +1 PLOT VELOCITY VECTORS BELOW PLATE F +1 PLOT SHEAR TRACTION ON BASE OF PLATE T +1 PLOT SURFACE VELOCITY VECTORS F +1 PLOT VELOCITY CHANGES SINCE LAST TIME F +1 PLOT GREATEST PRIN. STRAIN RATE F --------------------- PLOT HORIZONTAL-VELOCITY DISCONTINUITY F --------------------- PLOT SLIP-RATE OF FAULTS F 3.169E-10 1.0E-15+1 PLOT CRUSTAL THICKENING RATE F +1 PLOT VERTICAL INTEGRAL OF STRESS ANOMALY F --------------------- PLOT MOST COMPRESSIVE STRESS AXES F --------------------- PLOT NET EXTERNAL FORCE ON NODES F 1.00 2 PLOT Log10[Viscosity Integral]&LIMITS(1 or 2) 10 APPROXIMATE NUMBER OF CONTOURS IN EACH PLOT (FOR DEFAULT) 1 COASTLINES? (0=NONE, 1=PRESENT, 2=USER DATA) 12. RMS LENGTH OF VECTORS AND SYMBOLS, IN DEGREES (12. OR 4.) 360. WIDTH OF MAP, IN DEGREES (.LE.360.), E.G., 360. OR 60. 180., 0. (LON,LAT) OF PLOT CENTER 1 PEN WEIGHT FOR LIGHT LINES 2 PEN WEIGHT FOR MEDIUM LINES 3 PEN WEIGHT FOR HEAVY LINES T COLOR? (F GIVES BLACK-AND-WHITE) |
File iEarth5-049.in is another example of such a file, which you may use for starting values. Notice that:
the
first line gives a mnemonic label for the set of parameters;
all
values in the file are in SI units;
where
there are two values on a line, the first is for the crust, and the second is
for the mantle lithosphere;
thermal
expansion is volumetric, not linear;
radioactive
heat production is per unit volume, not per unit mass;
comments
on each line following the values are not read by my programs;
the
values in bold font are the ones particularly important for OrbData;
values
in the blue section will be used by Shells,
but not by OrbData;
values
in the green section will be used by graphics program OrbMapAI, but not by OrbData
or Shells or FiniteMap.
Once you have succeeded with this simple calculation, you might want to consider 2 options for increased realism:
(1) While running OrbData, answer "True" when asked if you want to use a .grd file of seafloor ages to determine model heat-flow and thickness in oceanic lithosphere. A sample file, based on Mueller et al. [1997; J. Geophys. Res., 102(B2)] is available in age_1p5.grd.
(2) Consider running my newer program OrbData5, which makes use of seismic constraints on crust and lithosphere thickness from two additional .grd files. See further comments in OrbData_READ_ME.txt.