I3ELVIS User Guide

User Manual:

Open the PDF directly: View PDF .
Page Count: 29

Download
Open PDF In Browser	View PDF

A
User's Guide
to
I3ELVIS
in Subduction and Collision setup

by Ria Fischer
Geophysical Fluid Dynamics
Department of Earth Sciences
ETH Zurich

Zurich
October 2012

Contents
1

Physics

1.1

1.2
1.3
2

. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
reactions
. . . . . .
. . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

t3c files

2.1
2.2

2.3

3

2

Basic physical principals . . . . . . . . . . . . .
1.1.1 The continuity equation . . . . . . . . .
1.1.2 The Navier-Stokes Euqation . . . . . . .
1.1.3 Heat conservation equation . . . . . . .
1.1.4 Rheology . . . . . . . . . . . . . . . . .
1.1.5 Impact treatment . . . . . . . . . . . . .
1.1.6 Computation of crust . . . . . . . . . .
1.1.7 Phase transitions, melting and hydration
I2ELVIS . . . . . . . . . . . . . . . . . . . . . .
I3ELVIS . . . . . . . . . . . . . . . . . . . . . .

9

le.t3c . . . . . . . . . . . . . . . . . . . . . . . . .
init.t3c . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Grid parameter description . . . . . . . . .
2.2.2 Rock type description . . . . . . . . . . . .
2.2.3 Boundary conditions . . . . . . . . . . . . .
2.2.4 Box description . . . . . . . . . . . . . . . .
2.2.5 Temperature box description . . . . . . . .
mode.t3c . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Timestepping description . . . . . . . . . .
2.3.2 General parameters . . . . . . . . . . . . . .
2.3.3 Erosion and Sedimentation parameters . . .
2.3.4 Velocity- and Pressure-iterations parameters
2.3.5 Temperature-iterations parameters . . . . .
2.3.6 Hydration and melting parameters . . . . .
2.3.7 Collision velocity parameters . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

Raw output files: .prn

3.1
3.2
3.3
3.4
3.5
3.6

Part
Part
Part
Part
Part
Part

2
2
2
3
4
5
6
6
7
8

I: General Information . . . . . . .
II: Rock type information . . . . .
III: Nodes information . . . . . . .
IV: Gridline positions . . . . . . . .
V: Boundary Condition Equations
VI: Markers . . . . . . . . . . . . .

9
9
9
10
13
14
15
17
17
18
18
19
20
21
21
22

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

22
23
23
24
24
25

4

Numerical problems and solutions

26

5

Usage of paraview

27

5.1
5.2

27
27
27

Conversion from raw output to Paraview les . . . . . . . . . . . . . . . . . .
Visualisation with Paraview . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 Visualize composition . . . . . . . . . . . . . . . . . . . . . . . . . . .

References

28

1

1 Physics

1 Physics
1.1 Basic physical principals
1.1.1 The continuity equation

The continuity equation describes the conservation of mass, while it is displaced in a continuous medium. In its Lagrangian form it reads the following,

Dρ
+ ρ∇ · ~v = 0,
Dt
where ρ denotes material density, ~v denotes displacement velocity and
grangian time derivative.

(1.1)

D
denotes the LaDt

For many geological media like the crust or the mantle, where temperature and pressure are
not too large and no phase changes occur, which would lead to larger volume changes, one
can assume the following incompressibility condition,

Dρ
∂ρ
=
+ ~v ∇ρ = 0,
Dt
∂t

(1.2)

which means that the density of material points does not change in time.
This leads us to the following incompressible continuity equation, which is the same in its
Eulerian and Lagrangian form,

∇ · ~v = 0,

(1.3)

The incompressible continuity equation very widely used in numerical geodynamic modelling,
although it is often a rather strong simplication.

1.1.2 The Navier-Stokes Euqation

The Navier-Stokes equation of motion in its full form reads the following,
0
∂σij
∂P
Dvi
−
+ ρgi = ρ
,
∂xj
∂xi
Dt

where σij is the strain-rate and ~g = (gx , gy , gz ) is the gravity vector.

2

(1.4)

1 Physics
Dvi
, is much smaller
Dt
compared to the gravitational force and therefore, be neglected. This leads to the Stokes
In highly viscous ows the right-hand side of (1.4), the inertial forces ρ

equation for creeping ow,

0
∂σij
∂P
−
+ ρgi = 0.
∂xj
∂xi

(1.5)

Under Boussinesq approximation the density is assumed to be constant, except in the buoyancy force term, where temperature and volatile content play an important role (?). Taking
into account the Boussinesq approximation, density ρ(T, P, c) in the buoyancy term ρgi may
vary locally as a function of temperature T , pressure P and composition c,
0
∂σij
∂P
−
= −ρ(T, P, c)gi .
∂xj
∂xi

(1.6)

1.1.3 Heat conservation equation

The heat conservation equation, also called temperature equation, describes the heat balance
in a convective medium, taking into account changes due to internal heat generation, advection
and conduction. The Lagrangian heat conservation equation reads as follows,


ρCp

DT
Dt


= −∇ · ~q + Hr + Ha + Hs + HL ,

(1.7)

with ~
q = −k(T, p, c)∇T , where thermal conductivity k(T, P, c) depends on temperature, pressure and rock composition c. Hr , Ha , Hs , HL denote radioactive, adiabatic, shear and latent
heating.
Adiabatic and shear heating have shown to be important in many tectonic situations, which
is why they are not taken as constant.

Hr = const.

(1.8a)

Ha = T α~v ∇P

(1.8b)

0
Hs = σij
˙0ij − ˙ij(elastic)

HL = const.




(1.8c)
(1.8d)

The resulting set of the above equations, together with equations (1.6) and (1.8), is called the
extended Boussinesq approximations.

3

1 Physics
1.1.4 Rheology
Equation of State

ρ = ρr [1 + β (P − Pr )] × [1 − α (T − Tr )]

(1.9)

where ρr is the reference density given separately for each rock type, Pr = 1.0 bar is the
reference pressure, Tr = 298.15 K is the reference temperature, α is the thermal expansion
and β is the compressibility.

Viscosity

Plastic yield strength

σyield = C + sin(φdry )(1 − λ)P

(1.10)

where σyield denotes the shear stress limit after which plastic yielding occurs, C is the cohesion,
Pf luid
φdry is the eective internal friction angle in dry rock, λ = 1 −
is the pore uid pressure
Psolid
factor and P = Psolid is the mean stress of the solid.

Peirl’s creep (Katayamo & Karato, 2008)

"

k #q )
E
+
P
V
σ
a
a
II
2
= AP eirl σII
exp −
1−
RT
σP eirl
(

˙II

(TODO: No!

Visco-plastic!)

deviatoric strain-rate

˙0ij

(1.11)

A visco-elasto-plastic rheology is employed, with the

being composed of the following components,

˙0ij = ˙0ij(viscous) + ˙0ij(elastic) + ˙0ij(plastic ,

(1.12)

where

1 0
σ ,
2η ij

(1.13a)

0
1 Dσij
,
2µ Dt

(1.13b)

˙0ij(viscous) =
˙0ij(elastic) =
˙0ij(plastic) = χ

4

0
σij
∂G
=
χ
0
∂σij
2σII

for G = σII = σyield .

(1.13c)

1 Physics
0 denotes the deviatoric stress tensor, µ denotes the shear moduwhere η denotes viscosity, σij
q
0
lus, G is the plastic potential, σyield is yield strength, σII = 12 σij2 is second deviatoric stress

invariant and χ is plastic potential.
Generally the strain tensor ij can be dened as a function of displacement ~u = (ux , uy , uz ),


∂uj
∂ui
+
.
∂xj
xi
q
0
and the second strain rate invariant is given by ˙II = 12 ˙ij2
ij =

1
2



(1.14)

The viscosity η is dened as follows,


η=

2
σII

(n−1)



Fn
E + PV
exp
,
AD
RT

(1.15)

where AD , E, V and n are experimentally dened ow parameters, R is the gas constant and

F is a dimensionless factor depending on the type of experiment (triaxial compression, simple
shear).
The viscous constitutive relationship relates stress σij with strain ij ,
0
σij
= 2η ˙0ij + δij ηbulk ˙kk ,

(1.16)

0 is the deviatoric stress, ˙0 is the deviatoric strain-rate, ˙
where σij
kk is the bulk strain-rate,
ij

and η and ηbulk are shear and bulk viscosity.

1.1.5 Impact treatment

The actual impact is not part of the model. The model only starts after the intrusion of the
impactor into the parent body. Processes like crater excavation, redistribution of impactor
and parent body material around the planet or decompression melting are not considered.
A simplied model takes into account the thermal anomaly created by the impactor. A
region called the isobaric core, of uniform temperature increase and shock pressure around
the impactor can be found (Senshu et al., 2002).

1
Ric = 3 3 rimp

(1.17)

where Ric is the radius of the isobaric core and rimp is the radius of the impactor.

5

1 Physics
Thermal anomaly in the isobaric core has been approximated by

∆T =

?

in the following way,

4π ψ ρP GRP2
9 F cP

(1.18)

where ψ is the eciency of conversion of kinetic energy to thermal energy and in this thesis
assumed to be 0.3.
Outside the isobaric core, for r > Ric , the thermal anomaly ∆T is decaying exponentially,
according to the following rule (Senshu et al., 2002; ?),


T (r) = ∆T

Ric
r

4.4

(1.19)

1.1.6 Computation of crust

This algorithm is only implemented in the 3D code.
Silicate melt within a certain depth is positively buoyant (ddepthmelt = 2e5) and rises up to the
surface (Golabek et al., 2011). Only markers with a melt fraction between 1% and 20% are
considered, as this corresponds roughly to the pyroxene fraction in a fertile mantle (Golabek
et al., 2011). Silicate melt on markers fullling these criteria now percolates upwards through
the mantle and at the surface crust is formed by freezing of the silicate melt.

1.1.7 Phase transitions, melting and hydration reactions
Melting

+20

Melting occurs a soon as the melt fraction is larger than 0. +20 is added to the type rock type
(not for 11/13/31/33 and 19/20/21). A rock type >20 signies that the material is partially
molten.

Freezing

-20

If the melt fraction becomes 0, -20 is subtracted from the type (not for 11/13/31/33 and
19/20/21). A rock type <20 signies that the material is fully solid.

Hydrated Peridotite

11 -> 34 (-> 14) (and 31 -> 14)

Hydrated wet mantle peridotite (type 11) can be molten (type 34) and resolidied again to
quenched dry mantle peridotite (type 14).

6

1 Physics
Hydrated, wet mantle peridotite (11) -> Molten Peridotite (34) -> Resolidied dry quenched
mantle peridotite (14)
(31 -> 14)

Mantle hydration

[9,10,12,14] -> 11

If water is present, lithospheric (type 9) and asthenospheric (type 10) peridotite as well as dry
peridotite from the shear zone (type 12) and the resolidied quenched (type 14) peridotite
are hydrated (type 11).

Crust hydration

5/6 -> 17/18

If water is present, continental upper (type 5) and lower (type 6) crust are hydrated (type
17+18).

Layered Sedimentation

Sequence of 3,4,3,4,3,4,....

Formation of new crust

-> 16

Antigorite weakening

13 -> 11 (11 -> 13)

Serpentinization of hydrated peridotite depending on Antigorite pressure and temperature
eld following (??)

Eclogitization

No type change.

Happens for type 7/8 (upper/lower oc. crust) and 27/28 (molten upper/lower oc. crust).

1.2 I2ELVIS

To model two-dimensional creeping ow under extended Boussinesq approximation, with both
thermal and chemical buoyancy, the conservative nite-dierence code I2ELVIS (?Gerya and
Yuen, 2007) is used, which operates on a staggered grid and uses the moving marker technique.
See also part ??, chapter ?? for a more detailed discussion. Silicate material is assumed to have
temperature-, pressure-, strain-rate and melt fraction-dependant visco-elasto-plastic rheology.
Furthermore, impact heat, batch melting of silicates and phase changes have all been taken
into account.

7

1 Physics
1.3 I3ELVIS

The 3D models have been carried out by the 3D numerical I3ELVIS (Gerya and Yuen, 2007)
code which is also based on a conservative nite dierence method with a marker-in-cell
technique and multigrid solver (?Gerya and Yuen, 2007). See also part

??,

chapter

??

for a

more detailed discussion. Additionally, the 3D code also features impact heat, batch melting
of silicates and phase changes as discussed in Golabek et al. (2011) for the 2D case. The
initial thermal-chemical model setup (including initial conditions, boundary condition and
uid/melt transport mechanism) and numerical approach are kept as similar as possible to
the 2D models. Furthermore, the 3D code also features computation of the primordial crust
from silicate melt.

Parameter
Radius of planetary body
Radius of impactor core
Radius of nal core rel. to Rplanet
Temperature of impactor core
Temperature of protocore
Temperature of diapirs
Mean temperature of nal core
Mean temperature of silicate mantle
Mean temperature of planetary body
Mean density of nal core
Mean density of silicate mantle
Mean density of planetary body
Volume fraction of iron (3D)
Mass fraction of iron (3D)
Gravitational acceleration
Table 1:

8

Symbol

Value

Unit

RM ars
ric
rcore
Tic
Tp
Td
T̄c
T̄m
T̄tot
ρ̄c
ρ̄m
ρ̄tot
fF e,vol
fF e,mass
g

3389
232 − 500
0.5
1300 − 2300
1300 − 2500
1300 − 2300
−
−
−
−
−
−
0.1
0.2
3.73

km
km
%
K
K
K
K
K
K
kg m−3
kg m−3
kg m−3
%
%
m s−2

TODO: complete table

List of parameters

2 t3c les

2 t3c files

2.1 file.t3c

The le 'le.t3c' contains only one number, which gives the number of the NEXT le to write.
The number given here is the actual number of the le, if all les would be numbered from the
beginning (including the initial le), starting with 0. The le name though can be dierent
(and is given in the 'mode.t3c'-le) and can contain any dierent number.

2.2 init.t3c

Null-point is in the frontal upper left corner of the grid.

2.2.1 Grid parameter description

The rst part of 'init.t3c' describes some basic grid parameters:

9

2 t3c les
Parameter
xnumx

Description
Number of nodes in x

ynumy

Number of nodes in y

znumz

Number of nodes in z

mnumx
mnumy
mnumz
xsize
ysize
zsize
pxinit
pyinit
pzinit
pinit
GXKOEFF
GYKOEFF
GZKOEFF
timesum
nonstab
xnonstab

Number markers per cell in x
Number markers per cell in y
Number markers per cell in z
Dimension of model in x
Dimension of model in y
Dimension of model in z
x-coordinate of initial pressure cell
y-coordinate of initial pressure cell
z-coordinate of initial pressure cell
Initial pressure in initial pressure cell
Gravitational acceleration in x
Gravitational acceleration in y
Gravitational acceleration in z
Starting time
number of random number generated
maximum random displacement of
markers in x direction
maximum random displacement of
markers in y direction
maximum random displacement of
markers in z direction
Number of initial output le

ynonstab
znonstab
markers types le
name Y(Name) N(0)
data output le name

TYPE

Name of initial output le
type of data output

Initial value
16N + 5 (4 multigrid levels)
16N + 5 (4 multigrid levels)
16N + 5 (4 multigrid levels)
(xnumx − 1)/2
(ynumy − 1)/2

-

Unit
[m]
[m]
[m]

-

[P a]
[m/s2 ]
[m/s2 ]
[m/s2 ]
[years]

-

-

name_0.prn
b (= binary), no
other supported

-

Table 2: Grid parameter description

2.2.2 Rock type description

In the next part a table lists all possible types of material compositions together with there
rheological properties.
(TODO: Add ductile rheology equation.
Chp.1) (TODO: Add density equation)

10

But maybe not here but in

2 t3c les
Parameter
rocknum
markn0
markn1
marks0
marks1
Nu or marknu
DE or markdh
DV or markdv
SS or markss
MM or markmm
LL or markll
a0 or marka0
a1 or marka1
b0 or markb0
b1 or markb1
e0 or marke0
e1 or marke1
RO or markro
bb or bRo or markbb
aa or aRo or markaa
CP or markcp
Kt or markkt
Kf or markkf
Kp or markkp
Ht

Description
Rock number
Individual lower viscosity limit
Individual upper viscosity limit
Individual lower stress limit (TODO: ???)
Individual upper stress limit
Newtonian viscosity
Activation energy Ea
Activation volume Va
dislocation-diusion transition stress σcrit
Stress exponent for creep law
Pore uid pressure factor 1 − λ (see eq. 1.10)
Cohesion (see eq. 1.10)
unused
Sinus of eective internal friction angle in dry rock
(see eq. 1.10)
unused
Koef in ductile rheology weakening (TODO: ???)
Koef in ductile rheology weakening (TODO: ???)
Density
Thermal expansion α (see eq.1.9)
Compressibility β (see eq.1.9)
Heat capacity
Thermal conductivity
Temperature dependency coecient of conductivity
Pressure dependency coecient of conductivity
heat generation

Unit
-

[P aM M ∗ s]
[J]
[J/bar]
[P a]
(P ower)
(koef )
[P a]

[deg°]/[rad]
-

(TODO: ???)
(TODO: ???)
[kg/M 3 ]
[1/K]
[1/kbar]
[J/kg]
[W t/(m ∗ K)]
[W t/m]
[1/bar]
[W t/kg]

Table 3: Rock type parameters

11

2 t3c les
Nr.
00
01
02
03
04
05
06
07
08
09
10
11
12
13
14

Description
Air
Water

15
16
17

(TODO: ???)

(TODO: Sediments 1)

Sediments 2
Sediments 3
Upper Continental Crust
Lower Continental Crust
Upper Oceanic Crust
Lower Oceanic Crust
Lithospheric Mantle
Asthenospheric Mantle
Hydrated Mantle
Shear Zone
Serpentinized Mantle
Resolidied peridotite/ quenched mantle
Newly formed crust
Hydrated Upper Crust

18

Hydrated Lower Crust

19
20
21
22
23
24
25
26
27
28
29
30
31

-

32
33

Molten Shear Zone

34
35
36
37
38
+50
+100
NAN

(TODO: Molten Sediments 1)

Molten Sediments 2
Molten Sediments 3
Molten Upper Continental Crust
Molten Lower Continental Crust
Molten Upper Oceanic Crust
Molten Lower Oceanic Crust
Molten Lithospheric Mantle
Molten Asthenospheric Mantle
(TODO: Molten
Mantle 1)

(TODO:
Molten
tinized Mantle)

Hydrated
Serpen-

Rheology
WET QUARTZITE RANALLI 1995
WET QUARTZITE RANALLI 1995
WET QUARTZITE RANALLI 1995
dry,felsic - WET QUARTZITE RANALLI 1995
dry,felsic - WET QUARTZITE RANALLI 1995
Basalts - WET QUARTZITE RANALLI 1995
Gabbro - An75 Ranalli1995
dry peridotite - DRY OL Ranalli1995
dry peridotite - DRY OL Ranalli1995
wet peridotite - WET OL Ranalli1995
dry mantle peridotite - WET OL Ranalli1995
wet peridotite
dry peridotite
Basalt
felsic - Hydrated WET QUARTZITE RANALLI
1995
felsic - Hydrated WET QUARTZITE RANALLI
1995

felsic
molten gabbro
molten basalt
molten gabbro
dry peridotite
dry peridotite
dry peridotite

Molten Peridotite

(TODO: dry or wet?)

Molten Newly Formed Crust
Molten Hydrated Upper Crust
Molten Hydrated Lower Crust
Water markers
External Composition
Undened composition

black

(TODO: ???)

light yellow
molten basalt
molten gabbro
molten gabbro

Table 4: Rock types

12

2 t3c les
Compositions 0-19 are all solids (except air and water). 20-39 are the equivalent melts. 50-89
are uid markers. 100-139 are external compositions of the equivalent type, waiting outside
the boundary to come into the model. They have the same properties as their internal
equivalents.

2.2.3 Boundary conditions

Boundary conditions are prescribed separately for each of the six faces of the cube. Boundary
conditions on each face don't have to be homogeneous. Special 'open boundary conditions'
can be prescribed, where e.g. the last marker gets 99% of the velocity of the second to last
and so on.

Parameter
Val
m10
m11
m20
m21
m30
m31
Const
Koef
nshiftx
nshifty
nshiftz
Koef1
nshiftx1
nshifty1
nshiftz1
Koef2
nshiftx2
nshifty2
nshiftz2

Description
BC type
starting node position in x
ending node position in x
starting node position in y
ending node position in y
starting node position in z.
ending node position in z

Range of values
P, Vx, Vy, Vz, T, X, Y, Z, M

Ignored if Koef1 is 0
Ignored if Koef1 is 0
Ignored if Koef1 is 0
Ignored if Koef2 is 0
Ignored if Koef2 is 0
Ignored if Koef2 is 0
Table 5: Boundary condition parameters

P, Vx, Vy, Vz, T boundary conditions

For all nodes in the given range [m10, m11] ×

[m20, m21] × [m30, m31] the following boundary condition is applied:
A(x, y, z) = Const + Koef · A(x + nshif tx, y + nshif ty, z + nshif tz)

(2.1)

E.g. the following part in the init.t3c le:
/ Left_Bondary

13

2 t3c les
/Val___m10_m11___m20_m21___m30_m31___Const___Koef_dm1_dm2_dm3
Vx____0___0_____0___y−1___0___z−1___3e−10___0____0___0___0

translates into the following boundary condition:

Vx (x, y, z) = 3 × 10−10

X, Y, Z coordinates definition

| x = 0, y ∈ [0, y − 1], z ∈ [0, z − 1]

(2.2)

For all nodes in the given range [m10, m11] × [m20, m21] ×

[m30, m31] the following equation is applied, shown at the example of the x-grid:

X(x, y, z) = X(x − 1, y, z) + Const + (Koef − Const) ·

(x − m10)
(m11 − m10)

| Koef 1 = 0 (2.3a)




Koef 1
(x − m10)
X(x, y, z) = X(x − 1, y, z) + exp log(Const) + log
·
Const
(m11 − m10)

M marker grid set to cell

(2.3b)

For all nodes in the given range [m10, m11] × [m20, m21] ×

[m30, m31] additional markers are added. All parameters (Koef , Koef 1, Koef 2, nshif tx,
nshif tx1,..,nshif tz1,nshif tz2) need to be set. nshif tx gives the shift in x, nshif ty1 gives
the shift in y , nshif tz2 gives the shift in z and all other shift parameters are ignored. If Koef
is set, random nonstability is set on the new marker eld in X in the following way:

markx = x+

rand()% (bConstc ∗ 2 + 1) − bConstc X(x + 1, y, z) − X(x, y, z)
·
·Koef (2.4)
Const
nshif tx

Random nonstability is set in Y for Koef 1 > 0 and in Z for Koef 2 > 0 in the same way.

2.2.4 Box description

Arbitrary shapes of cubes are prescribed by setting the (x,y,z) coordinate of each of the 8
corners.

14

2 t3c les
Parameter

Description

Range of value

Type
[x0 , y0 , z0 ]
[x1 , y1 , z1 ]
[x2 , y2 , z2 ]
[x3 , y3 , z3 ]
[x4 , y4 , z4 ]
[x5 , y5 , z5 ]
[x6 , y6 , z6 ]
[x7 , y7 , z7 ]

Rock type
front upper left
front lower left
front upper right
front lower right
back upper left
back lower left
back upper right
back lower right

0 − 140

Table 6: Box parameters

Coordinates can be given either relative [0,1], or absolute [m0, mMAXSIZE], where MAXSIZE
cannot be larger than the maximum domain size in this direction. With the leading 'm', the
value is interpreted as an absolute value in meters.
e.g.
/Type=RockType
/Type__x0__y0__z0___x1__y1__z1_____x2__y2__z2_____x3__y3__z3
/_____x4__y4__z4___x5__y5__z5_____x6__y6__z6______x7__y7__z7
/ Asthenosphere
10____0_m15000_0___0__1. 1__0_____1. 1_m15000_0_____1. 1 _1 . 1__0
______0_m15000_1___0__1. 1__1_____1. 1_m15000_1_____1. 1 _1 . 1__1

2.2.5 Temperature box description

Arbitrary shapes of cubes are prescribed by setting the (x,y,z) coordinate of each of the 8
corners.

15

2 t3c les
Parameter

Description

Type
[x0 , y0 , z0 ]
[x1 , y1 , z1 ]
[x2 , y2 , z2 ]
[x3 , y3 , z3 ]
[x4 , y4 , z4 ]
[x5 , y5 , z5 ]
[x6 , y6 , z6 ]
[x7 , y7 , z7 ]
t0
t1
t2
t3
t4
t5
t6
t7

Box type: 0 - simple box; 1&4 - age box; 5&6 - transitional
front upper left
front lower left with x1 = x0 , z1 = z0
front upper right with z2 = z0
front lower right with x3 = x2 , z3 = z0
back upper left with x4 = x5 , z4 = z7
back lower left with z4 = z7
back upper right with x6 = x7 , z6 = z7
back lower right
simple box: temperature P0 ; age box: surface temperature
simple box: temperature P1 ; age box: initial temperature
simple box: temperature P2 ; age box: thermal diusivity κ in P0
simple box: temperature P3 ; age box: thermal diusivity κ in P2
simple box: temperature P4 ; age box: thermal diusivity κ in P4
simple box: temperature P5 ; age box: thermal diusivity κ in P6
simple box: temperature P6 ; age box: characteristic diusion time
simple box: temperature P7 ; age box: overprinted linear geotherm [/m]
Table 7: Temperature box parameters

Coordinates can be given either relative [0,1], or absolute [m0, mMAXSIZE], where MAXSIZE
cannot be larger than the maximum domain size in this direction. With the leading 'm', the
value is interpreted as an absolute value in meters.

Box type 0: simple box

The given temperatures are set to the given coordinates and tem-

perature is interpolated linearly within the box.

Box type 1&4: age box

Within the given coordinates the temperature is given by conductive

cooling in the following way:

y
T (y, t) = T1 − (1 − erf ( √ ))(T1 − T0 ),
2 κτ

(2.5)

k
l2
is the thermal diusivity interpolated over the corner of the box and τ =
ρcp
κ
is the characteristic diusion time.
where κ =

Box type 5&6: transitional box

No new temperatures are set but existing temperatures at

the coordinates of the box corners are read. Those are then interpolated linearly through out
the box.

16

2 t3c les
E.g. the continental crust and lithosphere is modelled in a simple box by a geothermal gradient
of ∼ 17 K/km spanning the whole width of the model and starting at 12 km depth with a
surface temperature of 273 K , going down to a depth of 90 km:

/T_BOXES_DESCRIPTION
/Typ__x0________x2_______y0________y1_______y2_______y3_______z0_______
_x5_______x7_______y4_______y5_______y6_______y7_______z7_______
t0___t1___t2___t3___t4___t5___t6___t7
/ Oceanic_and_continental_geotherms
0
− 0.000001 m350001 m12000
m90000
m12000
m90000 − 0.000001
− 0.000001 m350001 m12000
m90000 m12000
m90000 1.000001
273 1595 273 1595 273 1595 273 1595

2.3 mode.t3c

2.3.1 Timestepping description

The rst/second line gives the lename and type of the initial le.
In the next very large block one timestep corresponds with one line and the data is only valid
for one timestep. Seperate timestep parameters are given in table 8.

Parameter
SAVEFILE
TYPE
cyc0max
maxxystep
maxtkstep
maxtmstep
nubeg
nuend
p0koef
p1koef
p2koef
unused
multinum
stp1
multinum2

Description
name of the le (if this name is non-unique, the le will be overwritten without
warning)
le type: b = binary (no other supported)
number of iterations
maximum step size in dimension (TODO: ???)
maximum step size in temperature (TODO: ???)
maximum step size in time (TODO: ???)
global lower viscosity cut-o
global upper viscosity cut-o
Pressure penalty factor 1, Relaxation parameter in Continuity equation
Pressure update interpolation from coarser multigrid levels
Pressure penalty factor 2 (TODO: ???)
number of multigrid levels
number of iterations/ repeats for whole timestep
number of multigrid levels 2 (TODO: ???)
Table 8: Timestep parameters

17

2 t3c les
2.3.2 General parameters

The following parameters are mainly used to control the general behaviour of the program.

Parameter
loadmod
printmod
crustmod
dynamod
0num
movemod
tempmod
markmod
gridmod
outgrid
densimod
stp100
CTreset
smeltextract
sthdatabase
p2vmod
lestop

Description
load from data le (1) or set initial conditions (0)
print information on the monitor, Yes (1)/ No (2)
print information on crustthickness in les for each nth timestep; 0 =
disable
do dynamo calculations for each nth timestep; 0 = disabled
number of otput le Names
do velocity-pressure iterations, solve continuity and Stokes equation
do temperature iterations, solve heat transfer equation
move markers Y(1-simple,2-Runge-Kutta4)/N(0)
recalculate density and viscosity
marker move out of grid Y(0)/N(1) Orthogonal Only (2)
mode of density calculation: 0-constant, 1-PT-dependent, 2-TDbase 3PT-dependent+WaterTDbase
(TODO: ???)

composition/temperature reset for water/air at 100 km above surface
Y(1)/N(0)
extract melt when moving markers Y(1)/N(0)
Use of Mars thermodynamic database Y(1) or standard database N(0)
convert each nth prn to vtr /N(0)
number of timesteps to execute before exiting

Default
1
1
1
1
2
1
2
3
9000
1
1
1
10
50

Table 9: General parameters

2.3.3 Erosion and Sedimentation parameters

Parameter
eroslev
sedilev
waterlev

Description
Erosion level: markers above this depth which are neither sticky
air nor water get converted to sticky air
Sedimentation level: sticky air or water below this depth gets
converted to sediments
Water level: sticky air markers below water level are converted to
water and water markers above water level are converted to sticky
air
Table 10: Erosion and Sedimentation parameters

18

Initial value
8000
20000
12000

2 t3c les
2.3.4 Velocity- and Pressure-iterations parameters

The following parameters are mainly needed to solve the Continuity and Stokes equation with
multigrid. Some these parameters already appeared in table 8 and their global value given
here will be overwritten with the new value for each timestep.

Parameter
cyc1max
DIVVMIN
STOKSMIN
DIVVMAX
STOKSMAX
multinum
multicyc
multinnn
p0koef
p1koef
p2koef
v0koef
v1koef
nubeg
nuend
nukoef
viscmod
viscoutermod
spheryn

Description
unused
Continuity equation lower error bound
Stokes equation lower error bound
Continuity equation upper error bound
Stokes equation upper error bound
number of multigrid levels; overwritten by separate
timestep value
Number of whole multigrid V-cycle iterations
V-cycle structure: number of GS-iterations for each multigrid level in upcycle and downcycle.
Global pressure penalty factor 1, Relaxation parameter
in Continuity equation; overwritten by separate timestep
value
Global pressure update interpolation from coarser multigrid
level; overwritten by separate timestep value
Global pressure penalty factor 2 (TODO: ???); overwritten by separate timestep value
Velocity penalty facto. Relaxation parameter in vx-,vy-,vzStokes equation
(TODO: ???)

global lower viscosity cut-o; overwritten by separate
timestep value
global upper viscosity cut-o; overwritten by separate
timestep value
Average ν for pressure optimisation
Eective viscosity mode. 0 = lin. interp; 1 = exp. interp;
2 = inverse interp;
viscosity in space/air/water; 1-gradual increase in space,
2-gradual increase in water/air
Spherical gravity. 0 = o; 1 = on;

Initial value
3000
3e−03
5e+01
0e−03
3e−03
4
1
4 16 16 32 0;
4 16 16 32 64
3.0e−01

1.0e−00
0.0e−00
1.0e−00
1.0e−00
1e+18
1e+25
0.0
0
2
0

Table 11: V and P iterations parameters

viscmod: Effective viscosity interpolation

Type of interpolation done to obtain eective vis-

cosity. viscmod = 0 uses linear interpolation (eq. 2.6a), viscmod = 1 uses exponential interpolation (eq. 2.6b) and viscmod = 2 uses inverse interpolation (eq. 2.6c).

ηef f =

1X
ηi
8

(2.6a)

i

19

2 t3c les

ηef f

!
1X
= exp
(log (ηi ))
8

(2.6b)

i

ηef f =

1
8

1
P
i

(2.6c)

1
ηi

2.3.5 Temperature-iterations parameters

The following parameters are mainly needed to solve the temperature equation with multigrid.

Parameter
cyc2max
HEATMIN
multinumt
multicyct
multittt
t0koef
t1koef
heatdif
frictyn
adiabyn

Description
unused
Temperature equation lower error bound
Number of multigrid levels
Number of whole multigrid V-cycle iterations
V-cycle structure: number of GS-iterations for each multigrid level in upcycle and downcycle
(TODO:
(TODO:
(TODO:
(TODO:
(TODO:

???)
???)
???)
???)
???)

Table 12: T iterations parameters

20

Initial value
2500
1e−4
0
1
1; 0
1.0e−00
1.0e−00
1.0
1
1

2 t3c les
2.3.6 Hydration and melting parameters

Parameter
tkpor
zmpor
vyuid
vymelt
dmwamin
tdeep
dtdeep
drdeep
zdeep
vdeep
nudeep
dxwater
dywater
dzwater
maxwater
minmelt
maxmelt

Description

Initial value

(TODO: ???)
(TODO: ???)

97300000.0
75000
−3e−09
−3e−09
1e−1
1880.0
100.0
000.0
660000.0
670000.0
1e+21
2e+3
2e+3
2e+3
5e−1
1e−2
1e−2

Initial uid velocity
Initial melt velocity
Minimum water release dierence
(TODO:
(TODO:
(TODO:
(TODO:
(TODO:
(TODO:

???)
???)
???)
???)
???)
???)

Fluid extension in x
Fluid extension in y
Fluid extension in z

(TODO: ???)
(TODO: ???)
(TODO: ???)

Table 13: Hydration and melting parameters

2.3.7 Collision velocity parameters

Linearly change collision velocity. Use initial constant velocity set from boundary conditions. Between timebeg and timeend linearly change the collision velocity to the nal velocity

velocitykf.

Parameter
timebeg
timeend
velocitykf

Description
Begin velocity change
End velocity change
Final collision velocity

Initial value
20e+6
25e+6
0

Table 14: Collision velocity parameters

21

3 Raw output les: .prn

3 Raw output files: .prn

Valid for subduction/collision setup

3.1 Part I: General Information

Parameter
xnumx

Description
Number of nodes in x

ynumy

Number of nodes in y

znumz

Number of nodes in z

mnumx
mnumy
mnumz
xsize
ysize
zsize
pxinit
pyinit
pzinit
pinit
GXKOEFF
GYKOEFF
GZKOEFF
rocknum
bondnum
marknum
n1
timesum
ival1
gridcur
gridtot

Number markers per cell in x
Number markers per cell in y
Number markers per cell in z
Dimension of model in x
Dimension of model in y
Dimension of model in z
x-coordinate of initial pressure cell
y-coordinate of initial pressure cell
z-coordinate of initial pressure cell
Initial pressure in initial pressure cell
Gravitational acceleration in x
Gravitational acceleration in y
Gravitational acceleration in z

Starting time

Initial value
16N + 5 (4 multigrid levels)
16N + 5 (4 multigrid levels)
16N + 5 (4 multigrid levels)
(xnumx − 1)/2
(ynumy − 1)/2

[m]
[m]
[m]

-

-

[P a]
[m/s2 ]
[m/s2 ]
[m/s2 ]

-

[years]

Table 15: prn-le general information block

22

Unit
-

3 Raw output les: .prn
3.2 Part II: Rock type information

Parameter
markn0
markn1
marks0
marks1
marknu
markdh
markdv
markss
markmm
markll
marka0
marka1
markb0
markb1
marke0
marke1
markro
markbb
markaa
markcp
markkt
markkf
markkp
markht

Description
Individual viscosity limit (TODO: ???)
Individual viscosity limit (TODO: ???)
Individual stress limit (TODO: ???)
Individual stress limit (TODO: ???)
Coecients in ductile rheology (TODO: ???)
Powerlaw stress
Koef in ductile rheology (TODO: ???) Activation
energy
Koef in ductile rheology (TODO: ???) Activation
volume
Koef in ductile rheology (TODO: ???) Max. Stess
up to dislocation creep/diusion creep
Koef in ductile rheology (TODO: ???)
Koef in ductile rheology (TODO: ???)
Koef in ductile rheology Coh (TODO: ???)
Koef in ductile rheology Coh (TODO: ???)
Koef in ductile rheology sinus/ friction angle

Unit
-

Koef in ductile rheology sinus/ friction angle

[deg°]/[rad] (TODO: ???)

Koef in ductile rheology weakening (TODO: ???)
Koef in ductile rheology weakening (TODO: ???)
Density
Density koef b
Density koef a
Heat capacity
Thermal conductivity
Temperature dependeny koef in conductivity

(TODO: ???)
(TODO: ???)
[kg/M 3 ]
[1/K]
[1/kbar]
[J/kg]
[W t/(m ∗ K)]
(TODO: ???)

Pressure dependeny koef in conductivity (TODO:

(TODO: ???)

heat generation

(W t/kg)

(TODO: ???)
(TODO: ???)

(TODO: ???)
???)

[P aM M ∗ s]
[J]
[J/bar]
[P a]

(P ower)
(koef )
[P a] (TODO: ???)
[P a] (TODO: ???)
[deg°]/[rad] (TODO: ???)

Table 16: prn-le Rock type information block, see also tbl 3.

3.3 Part III: Nodes information

For each node the following parameters are given. The order is the following, where n is the
number of nodes nodenum:

pr0 , vx0 , ...ht0 , pr1 , .., ht1 , .., prn , .., htn

23

3 Raw output les: .prn
Parameter
pr
vx
vy
vz
ro
nu
tk
cp
et
kt
ht

Description
Pressure
Velocity in x
Velocity in y
Velocity in z
Density
Viscosity
Temperature
Heat capacity
Thermal expansivity
Thermal conductivity
Heat sources

Unit
[P a]
[m/s]
[m/s]
[m/s]
[kg/m3 ]
[P a ∗ s]
[K]
[J/kg]
[1/K]
[W t/m/K]
[W t/kg]

Table 17: prn-le Marker information block

3.4 Part IV: Gridline positions
Position of gridlines in x direction

xnumx numbers are given

Position of gridlines in y direction

ynumy numbers are given

Position of gridlines in z direction

znumz numbers are given

3.5 Part V: Boundary Condition Equations

For each boundary condition there 5 values neede. The total number of boundary conditions
is bondnum. The boundary conditions are all of the following general form:

CU RP AR = CON ST + KOEF 1 ∗ P AR1

Parameter
m2
m3
bondv1[m3][0]
bondv1[m3][1]
bondn1[m3]

Description
Index in Indexmatrix bondm
Index in BC-matrices, saved in bondm
CON ST value in BC equ.
KOEF 1 value in BC equ.
P AR1+1 value in BC equ., P AR1 = 0 means no boundary

Table 18: prn-le Marker information block

24

3 Raw output les: .prn
3.6 Part VI: Markers

The total number of markers is marknum. For each marker the following parameters is given:

Parameter
markx
marky
markz
markk
markw
markd
markex
marktm
markc1
markc2
markt

Description
Momentary marker position in x
Momentary marker position in y
Momentary marker position in z
Temperature in K
Water percentage
Density
Meltfraction?
creation time
Melt composition Granite part ??
Melt composition Dacite part ??
Rock type of marker

Table 19: prn-le Marker information block

25

4 Numerical problems and solutions

4 Numerical problems and solutions
- Program stops because pressure gets too high.
Solution: This might indicate an awkward or inappropriate model setup. Try to lower the
two pressure penalty factors from table 8 to [0.5, 0.3, 0.1, 0.001]. If that doesn't help, try
another model setup.
- Timestep doesn't converge within the time-limit of the Brutus queue.
Solution: Use a longer queue, e.g. 24h-queue instead of the 8h-queue. If the problem still
persists try one of the following:
a) Increase continuity equation error limit DIVMIN in tbl. 11 to 3e−3
b) Decrease multinum1 and multinum2 in tbl. 8 to 3.
c) Increase both pressure penalty factors in tbl. 8.

Choose something along

[0.05, 0.11, 0.31, 0.41].
d) Decrease continuity equation error limit DIVMIN in tbl. 11 to 1e−3.
e) Try

the

following

congurations

for

the

pressure

penalty

[p0koef,p1koef,p2koef,p3koef] = [0.31, 1.00, 0, 0] , [0.41, 1.00, 0, 0], [0.11, 0.50, 0, 0]
f) Make 4 v-cycles series. In tbl. 11 set multicyc = 4.
- Unstable model, model goes into very fast convection
Solution: Try one of the following:
a) Decrease maximum timestep maxtkstep from 500 to 100.
b) Decrease maximum timestep maxtmstep from 5e3 to 1e3.
c) Increase lower viscosity limit from 1018 to 1019 .

26

factors

5 Usage of paraview

5 Usage of paraview
5.1 Conversion from raw output to Paraview files

prn2vtr

5.2 Visualisation with Paraview
5.2.1 Visualize composition

Use the following pipeline:
- Threshold lter:
set upper and lower threshold to the same composition number. Don't set the 'All Scalars'
option.
- Contour lter:
Contour by composition. Choose x.001 as value for the isosurface.
- GenerateSurfaceNormals lter:
This generates a polygonal mesh, which has a smoothing eect.

27

References

References
Gerya, T. V. and Yuen, D. A. (2007). Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems. Physics of the Earth and Planetary

Interiors, 163(1-4):83105.

Golabek, G. J., Keller, T., Gerya, T. V., Zhu, G., Tackley, P. J., and Connolly, J. A. D.
(2011). Origin of the martian dichotomy and tharsis from a giant impact causing massive
magmatism. Icarus, 215(1):346357.
Senshu, H., Kuramoto, K., and Matsui, T. (2002). Thermal evolution of a growing mars. J.

Geophys. Res., 107(E12):5118.

28



---

Source Exif Data:

File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 29
Producer                        : MiKTeX pdfTeX-1.40.12
Creator                         : TeX
Create Date                     : 2014:06:05 09:27:41+02:00
Modify Date                     : 2014:06:05 09:27:41+02:00
Trapped                         : False
PTEX Fullbanner                 : This is MiKTeX-pdfTeX 2.9.4487 (1.40.12)

EXIF Metadata provided by EXIF.tools