OPEN RESEARCH
Files
for Research Paper:
Peter Bird, Jon Bryan May, and Michele M. C. Carafa [2026?]
“Fault Friction, Plate Rheology, and Mantle Torques from a Global Dynamic Model
of Neotectonics,” in
Journal of Geophysical Research/Solid Earth
These files, together with others previously published on this website and elsewhere, should allow researchers to reproduce all model-building, experimentation, and map-production in the above paper. Please also see USAGE NOTES at the end of this page for helpful advice.
I. Files used by OrbData2023 to compute nodal data, as described in Supporting Information: Updated Lithospheric Structure
ETOPO2v2g.grd
This file is based on a global DEM (NOAA National Geophysical Data Center,
2006) with grid-spacing of 2’ in both latitude and longitude. The filename
extension “.grd” indicates that the information has been transformed into Peter
Bird’s .grd format, which is described here:
http://peterbird.name/guide/grd_format.htm
The Fortran 90 conversion program was ETOPO2v2g_2_GRD.f90.
This file is quite large (300 MB), so for download purposes it has been
transformed into ETOPO2v2g.grd.zip (84
MB). An alternative would be to use the older ETOPO5.grd, which has 5’ spatial
resolution; this file is only 15 MB in .zip format: ETOPO5.zip. Either is
adequate for this project.
HFgrid14.grd
The global heat-flow prediction model of Lucazeau (2019), specifically version
HFgrid14 which incorporates other datasets (topography, seismic velocity,
tectonic setting, gravity) as well as traditional borehole and probe
measurements. Upper and lower bounds of 300 and 30 mW/m2 were
applied to avoid computational problems in ShellSet. Again, Lucazeau’s data
has been reformatted in Peter Bird’s .grd format (see above), and the Fortran
90 conversion program was HFgrid14_2_GRD.f90.
CRUST1_CrThick.grd
The CRUST1.0 crustal-thickness model of Laske et al. (2013), which has 1° resolution.
LAB_LithoRef2018.grd
The LithoRef2018.xyz model of the depth to the lithosphere/asthenosphere
boundary (LAB) by Afonso et al. (2019), which merged 6 tomographic models to
determine upper-mantle velocities, with gravity and geoid constraints to
estimate compositional density variations.
age_1p5.grd
Ages of seafloor lithosphere, interpolated from marine magnetic anomaly bands.
Based on the dataset of Mueller et al. (1997). Special integer codes are used
to indicate “continental” points and also points of “unknown age.”
OrbData2023.f90 and
OrbData2023-Win64seq.exe
These are the latest versions of utility program OrbData, which uses the global
data files above to compute nodal data for each node of a pre-existing F-E grid
(.feg) file, and then writes a modified version of that .feg file including the
nodal data. Version .f90 is the Fortran 90 source code; version -Win64seq.exe
is a 64-bit sequential executable for Windows. In this project, OrbData2023
was applied to existing .feg file Earth5R.feg
(from Paper I) to produce file Earth5-2023.feg.
II. Input Files Required to Recreate (or Plot) Preferred Model Earth5-271 “Asia”:
PB2002_plates.dig
Outlines of each of the 52 largest plates, according to the PB2002 model of
Bird (2003). This information will be used to assign a plate-affinity to each
node of the F-E grid.
PB2002_boundaries.dig
All plate-boundary segments in the PB2002 model of Bird (2003). This will be
used to decide the polarity of subduction zone megathrusts, and the two plates
involved in each.
Earth5-2023d.feg
For documentation of the structure of any .feg file, please see OrbWin_manual.pdf.
The node locations and element topology of this F-E grid file are inherited
from file Earth5R.feg of Paper I (Bird et al., 2008, JGR). The nodal data have
all been replaced, as explained in section I of this page, above. The suffix
“d” in the filename is a reminder of 2 special features specific to the “Asia”
experiment:
(1) The dips of all fault elements representing the East Anatolia fault zone
have been corrected from 45° to 90°; and
(2) Lithospheric Rheology Indices (LRIs) have been set to 0 for all elements inside
the Persia-Tibet-Burma orogen, and set to 1 for all elements outside. This
operation was performed with utility program Assign_LRs; its Fortran 90 source
code is Assign_LRs.f90. Its compiled
executable (for Windows) is Assign_LRs-Win64seq.exe.
The outline of the P-T-B orogen is in 1_selected_PB2002_orogen.dig.
This is an excerpt from the larger file PB2002_orogens.dig
published by Bird (2003) as part of the PB2002 plate-boundary model.
iEarth5-271.in
Input parameter file. The first line is a title. In each following line, the
parameter name (as it appears in Shells_v5.f90
source code) appears to the right of the parameter value(s). All values are in
the SI metric system: m, kg, s, Pa, J, W… All temperatures are on the absolute
Kelvin (K) scale. Specifically, integrated megathrust resistance (tauMax) is
in N m-1, all densities (rho___) are in kg m-3, thermal
conductivity (conduc) is in W m-1 K-1; volumetric
radioactive heat production (radio) is in W m-3. For definitions
and units of dislocation-creep flow-law parameters aCreep, bCreep, cCreep,
dCreep, & eCreep, see equation (1) in the paper.
LR_set_for_Earth5-2023d.dat
A table of Lithospheric Rheologies with LRI > 0. In this case, the only
entry is for LRI = 1, which indicates (per Earth5-2023d.feg) all areas outside
the Persia-Tibet-Burma orogen, and these are assigned the “Diorite2” rheology.
Earth5R-type4AplusA.bcs
One of 2 boundary-conditions (.bcs) files required for any model with parameter
iConve = 6. Type-4 velocity boundary conditions mean that Shells will
automatically calculate the correct nodal velocity from the plate-affinity of
the node and the PB2002 model of Bird (2003). In this file, velocity boundary
conditions apply to the leading edges of subducting plates (just below the
megathrust), and also to 2 or 3 relatively strong nodes in the interior of each
slabless plate. This file must be used in all iterations of a Shells_v5
solution EXCEPT the last one.
Earth5R-type4A.bcs
The second of 2 boundary-conditions (.bcs) files required for any model with
parameter iConve = 6. Type-4 velocity boundary conditions mean that Shells
will automatically calculate the correct nodal velocity from the plate-affinity
of the node and the PB2002 model of Bird (2003). In this file, velocity
boundary conditions apply ONLY to the leading edges of subducting plates (just
below the megathrust). This file must be used in the FINAL iteration of a
Shells_v5 solution.
III. Output Files Produced by Preferred Model Earth5-271 “Asia”:
tEarth5-271.txt
Human-friendly text file recording the progress of a Shells_v5 solution (or multiple
solution iterations, if iConve = 6). There will be at least one table (when
running Shells_v5) that shows how well the inner iteration on effective
viscosity of each integration point of each element has converged. When running
ShellSet, the output includes multiple tables of this kind, one for each
iteration. At the end of a t____.txt file produced by Shells_v5, there are
also long tables showing stress, strain-rate, and brittle/ductile transition
depth in each element. However, these per-element tables are omitted in the
ShellSet version of the file.
fEarth5-271.out
External horizontal forces exerted on each node of the grid, in Newtons. For
each node, the South component is listed first, followed by the East component.
The fact that 5 values (i.e., 2.5 nodes) appear in each line makes it
hard to find any particular node manually; instead, this information is best
visualized in a FiniteMap plot. Large values are forces needed to impose
velocity boundary conditions; small values are just numerical noise reflecting
slight imperfections in the solution of the linear system(s) that were supposed
to enforce vertically-integrated stress-equilibrium on a sphere. (In fact,
these small vectors are typically SO small that they disappear from FiniteMap
plots.)
qEarth5-271.out
Torque report. For each of the 52 plates of the PB2002 model (Bird, 2003) the
balance of 3 torques (lithostatic-pressure, side-strength, and basal-strength)
is reported, each in N m. Also, each torque is converted to a
statically-equivalent point force in Newtons applied at the planet’s surface.
If you are interested in only 1 or 2 plates, this format may be convenient.
However, to get a global overview, it is essential to plot this file with
FiniteMap (as in Figure 5 of the paper).
vEarth5-271.out
Horizontal long-term-average velocity at each node of the grid, in m s-1.
For each node, the South component is listed first, followed by the East
component. The fact that 4 values (i.e., 2 nodes) appear in each line
makes it hard to find any particular node manually; instead, this information
is best visualized in a FiniteMap plot.
IV. Scoring Datasets used by OrbScore3 or ShellSet, as described in Supporting Information: Updated Scoring Datasets:
GPS2014_selected_subset.gps
Geodetic Velocity (“GV”) scoring dataset. Obtained from Global Strain Rate Map
version 2.1 of Kreemer et al. (2014), which provides interseismic GPS
velocities at 18,356 locations.
robust_interpolated_stress_2023.dat
Stress Direction (“SD”) scoring dataset. We interpolate stress directions from
the World Stress Map Project (Reinecker et al., 2004; Heidbach et al., 2016;
2018) to a uniform global grid of points by the clustered-data algorithm of
Bird and Li (1996). The result is a set of interpolated stress directions at
764 grid points (Figure SI-3) for this project.
aggregated_offset_rates.dig
Fault Slip Rate (“FSR”) scoring dataset. We use selected fault traces (Figure
SI-2) and slip rates from the Global Active Fault Database (GAFD) of Styron and
Pagani (2020). Specifically, we only use rates which are reported with upper
and lower limits, because we note that most rates without uncertainties are
from rigid plate models (e.g., Bird, 2003), not from geologic data. We
also require the fault dip to be specified, for use in converting heave-rates
to slip-rates. After these two restrictions, the number of faults potentially
useful for scoring is reduced from 13,696 to 2,487.
magnetic_PB2002.dat
Seafloor Spreading Rate (“SSR”) scoring dataset. The same dataset as in Paper
I: 277 rates from the NUVEL-1 collection in Table 3 of DeMets et al. (1990),
which includes citations. We multiplied each of these rates by 0.9562 to
correct for the adjustment in the paleomagnetic reversal timescale that DeMets
et al. (1994) made in their revised NUVEL-1A model. As in Paper I, we added 35
rates which were published more recently; most are from back-arc spreading
systems not considered in the NUVEL models, but added in the PB2002 model of
Bird (2003).
Becker_SKS_2-degree_20220807-selected.dat
Seismic Anisotropy (“SA”) scoring dataset. We use 1,528 selected 2°-trapezoid averages from file
Becker_SKS_2-degree_20220807.dat, obtained from Thorsten Becker’s archive at
http://www-udc.ig.utexas.edu/external/becker/sksdata.html (last accessed
2022.10.24). The selection rules were exactly the same as in Paper I.
Note: OrbScore3 and ShellSet both allow the option to use a catalog of shallow earthquakes (in Bird’s .eqc format) as an additional scoring dataset, but that option was not used in this project.
V. USAGE NOTES
1. Downloading files
To download any file, right-click on its link, and choose “Save As …”, then
select the target directory on your own computer. As always (on this website),
ASCII data flat-files (such as .grd files) have an extra filename extension of
“.txt” added to make downloading go smoothly. Similarly, any Windows
executable file (.exe) will have an extra filename extension of “.bin” added to
make downloading go smoothly. After you have obtained each file, please delete
this extra filename extension by renaming the file.
2. Manual editing of F-E grid files
In order to change the topology of any F-E grid (.feg) file (i.e.,
adding or deleting elements; moving nodes; etc.), you need to use
interactive utility program OrbWin-Win64seq.exe
on a Windows computer equipped with a mouse or trackball. Please read the
electronic manual OrbWin_manual.pdf
first! A few very-useful “basemap” sources are WorldMap.dig
(coastlines) and/or PB2002_plates.dig
(outlines of 52 plates). Finally, be aware that after any editing, you will
have to perform 3 more necessary steps:
(a) Renumber the nodes with utility program OrbNumber.
(b) Manually modify the 2 boundary-conditions files (.bcs files) to reflect the
new node numbers. {Frankly, this will be tedious!}
(c) Recompute nodal data using OrbData2023, as explained in section I above.
3. Computing a single F-E model with Shells_v5
If you only want to reproduce the preferred model {Earth5-271 “Asia”}, then it
will be much easier to do this manually with Shells_v5 than it would be to
install and compile ShellSet (within a Unix partition, equipped with the Ubuntu
GUI, using specific free libraries etc.). Before starting, please skim
Paper I (Bird et al., 2008, JGR) to understand how a set of N runs of
Shells_v5 is needed to converge on a solution, where N @ 5 to 10. Each run (except the first) will
use the torque-report file produced by the previous run, and output an improved
estimate of the torque-report file. Also, all iterations except the last (1,
2, 3, …, N-1) will use boundary-conditions file Earth5R-type4AplusA.bcs,
while iteration N will use boundary-conditions file Earth5R-type4A.bcs.
4. Scoring a single F-E model with OrbScore3
If you perform step #3 above with Shells_v5, then it will be appropriate to use
the scoring files (presented above in section IV) to compute errors in model
predictions, using OrbScore3.
5. Running parameter-optimization experiments with ShellSet:
All files and information about ShellSet (source codes, installation and
compilation instructions, sample input files, etc.) was published online
by:
May, J. B., P. Bird, and M. M. C. Carafa [2024] ShellSet v1.1.0 parallel
dynamic neotectonic modelling: A case study using Earth4-049, Geosci.
Model Dev., 17(5), 6153-6171, https://doi.org/10.5194/gmd-17-6153-2024.
6. Creating maps to represent (almost) all input and output
datasets:
The common mapping utility for many of these files (those in .grd, .feg, .gps, .dig,
and .out formats) is FiniteMap.
This program runs interactively using a text-mode interface, and produces one
map per run. The user first selects the size of the virtual paper, the sizes
of the margins, the geographic center point of the map, the map scale, and the
map projection. Then, they can choose to plot one opaque “mosaic” layer (e.g.,
any .grd file, or v___.out velocity magnitudes, with automatic interpolation
between nodes), overlain by any number of semi-transparent “overlay” layers.
These different layers can come from different datasets, so it is possible to
visually compare model predictions to scoring datasets. Depending on what you choose
to plot, FiniteMap will require that you input some mixture of the files in
sections II, III, and/or IV above. Please use only these files, and do
not introduce any new versions, which would misrepresent the preferred model.
Some very-useful overlay sources are WorldMap.dig
(coastlines) and/or PB2002_plates.dig
(outlines of 52 plates). The output will be a single-page document in Adobe
Illustrator (.ai) file format. Thus, you will need a copy of Adobe Illustrator
(not necessarily the latest one!) to open the map, modify it, and then either
print it or save it in another electronic file format (such as object-oriented
.EPS, or bitmap in .JPG or .GIF).
Peter Bird
Professor Emeritus
Department of Earth, Planetary, and Space Sciences
UCLA
pbird@epss.ucla.edu
August 2025