Dalhousie Geodynamics Group


Department of Oceanography
Dalhousie University
1355 Oxford Street
PO Box 15000
Halifax, NS, B3H 4R2
Canada




SOPALE

Sopale Nested

Processing SOPALE Output

Printing

MOZART

TMM (Thermo-mechanical Model)

old SOPALE documentation

Geodynamics home page

sopale_nested input

Notation

In many of the lines below, there are items in square brackets []. The square brackets are meant to indicate that the item(s) inside are optional. This DOES NOT APPLY to the section header line [SOPALE1_i]. It must be entered as written here.

Lines whose first character is '!' are treated as comment lines; they are flushed from the input stream.

[SOPALE1_i]

Indicates start of section.

template_i

template_i
number of lines in the model description, which follows immediately.

comment line

comment line
There should be template_i of these.
sopale_nested will echo these lines to the standard output file. The modeller is encouraged to use these lines as a description of the model.

nmodels

nmodels
sopale_nested reads and ignores this.

modelsn(1) modelsn(2) ... modelsn(nmodels)

modelsn(i), i=1,nmodels
sopale_nested reads and ignores this.

modelkey

modelkey
sopale_nested reads this with an a300 format, echos it to the text output file, then forgets it.

amodel benchmark

amodel
sopale_nested reads and ignores amodel.
sopale_nested reads and ignores benchmark.

wres rres align_time [ update_viscosities [ cond_weight ]]

wres
= number of timesteps between writing of 2 restart files
rres
determines whether to read restart files.
= 0   means don't read restart files.
= 1   means read restart files.
align_time
(only has effect when rres = 1)
= 0   On restart, set the time to zero and timestep to 1. This is appropriate when you are using restart files which just establishes the state you want.
= 1   On restart, set the time and timestep to those found in the restart files. This is good for restarting an interrupted run.
Since version 2009-07-26_22.9:   When restarting with the pressure pump on (type_bc=300), sopale_nested will stop if align_time is not 1.
Since version 2009-07-26_26.0:   When restarting with the pressure pump on (type_bc=300), and pump_time greater than dt, sopale_nested will stop if align_time is not 1.
Said another way:
  1. When the pressure pump is on
  2. The requirement use aligntime relates to the value of pump-time. If pump_time is greater than the timestep length aligntime must be used....Rajesh-type models. If pump_time is less/equal to the timestep length or not defined the user has the option to use aligntime or not.
  3. If the pressure-pump is off the user has the option to use aligntime or not

Since version 2009-07-26_22.

update_viscosities
An optional parameter. Rationale
viscosities are normally updated on the second, and subsequent, iterations of every time step.
= 0 means do not update viscosities.
= 1 means update viscosities. This is the default.

Since version 2009-07-26_27.0

cond_weight
An optional parameter.
Expected value is between 0.0 and 1.0; default is 0.5
Previous versions of sopale_nested implement cond_weight=1.0

When a finite element calculation is done to compute the initial temperature field (type_bct0=3), and thermal conductivities vary with temperature, the calculation sometimes does not converge, and the error flips back and forth between two solutions, one positive, the other negative . In an attempt to prevent this "flipping", we use a weighted average of the thermal conductivities from the current and previous iterations. The formula is
cond_weight*current + (1-cond_weight)*previous

There are test results of a (previously) non-converging temperature field, where the convergence error alternated between 21.xxx and -21.yyy degrees.

While convergence for cond_weight=0.5 in some tests is very quick, there is a good possibility that this is chance, as opposed to good design.

Convergence or non-convergence of the finite element calculation will depend both on the initial temperature(s) and the boundary conditions. Your mileage will vary.

nsave1

nsave1
number of saves for Eulerian grid
(Maximum number of saves is 200)

isave_1(1) isave_1(2) ... isave_1(nsave1)

(isave_1(i), i = 1, nsave1)
[only needed if nsave1 > 0 ]
timesteps at which to save the Eulerian grid

nsave2

nsave2
number of saves for Lagrangian grid
(Maximum number of saves is 200)

isave_2(1) isave_2(2) ... isave_2(nsave2)

(isave_2(i), i =1, nsave2)
[only needed if nsave2 > 0 ]
timesteps at which to save the Lagrangian grid

nameout1 [ nameout1_r ]

nameout1
soaple_nested uses this as a run name. Output file names are prefixed with it. Copies of the input file(s) are made; the copies have this pre-pended to their original name. Restart file names may begin with it. The readability of file lists is improved if the last character is an underscore.

length limit is 200.
nameout1_r
An optional input. run name for the restart files. '_i' will be appended to this to create the input file names.
Default: nameout1
Since: implementation date not immediately at hand

time_nstep dt outunt [ dt_adjust dt_min_f dt_max_f dt_incr_f time_end [ dt_ref ]]

time_nstep
number of timesteps in this model run, unless dt_adjust is 'yes'

If dt_adjust is 'yes' the model timestep lengths will vary and the model will be run under cfl (Courant_Friedrichs_Lewy factor) control to limit the maximum advection per timestep.
dt
timestep length in seconds, or initial timestep length if dt_adjust is 'yes'
    Restarts with dt_adjust = no
  • dt (from this file) will be used.
    Restarts with dt_adjust = yes
  • If dt is found in the restart file, it will be used, provided it fits within the range defined by dt, dt_min_f, dt_max_f.
  • If dt is found in the restart file, and does not fit within the range defined by dt, dt_min_f, dt_max_f, then dt (from this file) will be used.
outunt
sopale_nested read and ignores this.
dt_adjust
yes => The time step length (dt) will be adjusted each time step
any other value => dt will not be adjusted.

First implementation version 2013-05-30_1.0 . frame numbers on restart may be inconsistent.
retrofitted to version 2013-03-21_1.1 . frame numbers on restart may be inconsistent.
Revised in version 2013-06-28_1.0 . frame numbers on restart may be inconsistent.
Revised in version 2013-09-17_1.0

In version 2013-09-17_1.0:
  • In each time step (after the thermal loop) the eulerian grid is scanned, horizontal and vertical directions separately, and the cells with the largest cfls identified. A dt is calculated for each direction. (See the ucompute input line below for details of the user supplied cfl.)

    dt_x = cfl * cell_width / x velocity_at_upper_left_corner
    dt_y = cfl * cell_height / y velocity_at_upper_left_corner


  • If both of the cfl values are less than the user supplied cfl, a new dt is calculated:

    dt_new = dt_current + dt_incr_f * (min(dt_x, dt_y) - dt_current)

  • If either of the cfl values are larger than the user supplied cfl, dt is set to the smallest of dt_x, dt_y. The rationale here is
    • when increasing dt, it is done gradually.
    • when decreasing dt, it is done all at once, since the measured cfl is larger than desired.
    dt_new = min (dt_x, dt_y)
  • Otherwise dt is not changed.
    dt_new = dt_current
  • The new dt is restrained to the user-defined limits. See dt_min_f and dt_max_f.
  • Once LS and any nests have calculated their new dt separately, the smallest is chosen for the next time step.
dt_min_f
dt will not be allowed to become smaller than dt * dt_min_f.
only used if dt_adjust is 'yes'.
Default: 0.1
Suggested: 0.1
dt_max_f
dt will not be allowed to become larger than dt * dt_max_f.
only used if dt_adjust is 'yes'.
Default: 10.0
Suggested: 10.0
dt_incr_f
dt_incr_f moderates the rate at which dt is increased.
only used if dt_adjust is 'yes'.
Default: 0.1
Suggested: 0.1
time_end
When dt is adjustable, using the number of time steps to determine the length of a model run is not useful. So this number (in million years) is used.
Default: no default
dt_ref
When dt_adjust='yes', dt_ref is used to map timestep numbers to time when writing output files. When doing restarts, setting this number to the dt used to calculate the time step numbers for each type of output file will ensure that frame numbers on the output files behave consistently.
Default: dt
added in Sept 2013

bigu bigt

bigu
multiplier for penalizing velocity boundary conditions
Recommended value: 1d7
bigt
multiplier for penalizing the temperature boundary conditions
Recommended value: 1d7

length height nx1 ny1 betaxg betayg length2 height2 nx2 ny2

The variables in this line are a partial definition of the eulerian and lagrangian grids.

length
length of Eulerian rectangle
height
height of Eulerian rectangle
nx1
number of columns of Eulerian grid
ny1
number of rows of Eulerian grid
betaxg
sopale_nested reads and ignores this.
betayg
sopale_nested reads and ignores this.
length2
length of Lagrangian rectangle
height2
height of Lagrangian rectangle
nx2
number of columns of Lagrangian grid
ny2
number of rows of Lagrangian grid

NOTE: surface is at "height", not 0
nodes are numbered starting at top surface
x=0 is at left edge

npert

npert
sopale_nested reads and ignores this.

nfpert

nfpert
sopale_nested reads and ignores this.

iformul

iformul
a parameter determining that plane strain calculation will used in Sopale
Required value = 2

hists(1) hists(2) hists(3) ... hists(11)

(hists(i),i=1,11)
sopale_nested reads and ignores this.

ang_gravity gravity pcohes visminmax(1) visminmax(2) psmooth p_bdry_layer splasmin vismin_bdry_layer edotmax

ang_gravity
angle of gravity in degrees
The convention is that ang_gravity is positive when measured anti-clockwise from the positive X direction. Gravity directed in positive X direction has ang_gravity=0.
gravity
value of gravity, in SI units (m/s2
assumed to act toward decreasing Z
pcohes
lower bound for pressure (pressures below this are set to pcohes)
visminmax(1)
minimum viscosity limit (any viscosities below this are set to visminmax(1))
visminmax(2)
maximum viscosity limit (any viscosities above this are set to visminmax(2))
psmooth
= 1 for pressure smoothing
= anything else, no pressure smoothing
NOTE: [Jun 12, 2006] In versions of the code other than the ones listed below, there are two bugs which prevent psmooth from having any effect (i.e. no pressure smoothing is done). The following versions of the code have both these bugs fixed:

  • code-feb07-06-adiabat, where the executable is dated later than Jun 8, 2006
  • code-mar01-06-pttpaths, where the executable is dated later than Jun 8, 2006

See the email here for information on these bugs if you need to fix it in your own code. The first email in the link deals with the first bug; the second email deals with the second bug, which was discovered and fixed later.

p_bdry_layer
minimum pressure for setting minimum splas and viscosity; i.e. if the pressure is less than p_bdry_layer, then the splas and viscosity are kept above the minimums specified in splasmin and vismin_bdry_layer.
splasmin
minimum splas limit for areas where epress(i) is less than p_bdry_layer
vismin_bdry_layer
minimum viscosity limit for areas where epress(i) is less than p_bdry_layer.
edotmax
maximum value of the second invariant of the strain rate

see NOTE on the use of edotmax

n_matsets

n_matsets
number of material sets (aka rock types) whose mechanical properties are to be defined

A material set is a collection of materials whose properties are all the same. This allows the user to define several materials which have the same mechanical properties but different material numbers. This can be useful in tracking particular areas during model evolution.

Said another way, a material set is a rock type with one set of mechanical properties, but with several mineralologies. The different mineralologies are represented by different "colors".
NOTE: FOR VERSIONS OF THE CODE BEFORE code-nov16-05-matsets, material sets were not used. Instead, the variable color_number was used to define the total number of materials.

mixhyp refsrate redphi nredphi plas_fluid_method

mixhyp
mixing model for viscous material properties (viscous rheologies)
= 2, meaning viscous rheologies act in parallel.
Each viscous rheology is assumed to have the same second invariant of strain rate (they act in parallel). Stresses are determined from their respective effective viscosities X strain rate. Total second invariant of stress is sum of individual stresses. The total viscous stress is compared with the plastic yield stress to determine whether flow is viscous or plastic.
Since version 2009-07-26_8:
= 1, meaning viscous rheologies act in series.
Each viscous rheology is assumed to have the same second invariant of stress (they act in series). The total second invariant of strain rate is the same as the strain rates of each viscous rheology acting alone. The total viscous stress is compared with the plastic yield stress to determine whether flow is viscous or plastic.
This variable used to be floating point. Values like '2.0d0' have been seen in older input files. The variable is now an integer, and values like "2.0d0" should be changed to "2".
refsrate
reference strain rate to start viscosity computation viscous creep flow models.
redphi
mechanical material number for all elements greater than nredphi
nredphi
element number. Eulerian elements greater that this number will use the properties of material redphi for mechanical calculations. Note: elements are numbered by row, so this can be used to set a basal boundary layer. The mechanical material properties overrides those specified through the normal input process (Quick and dirty method for setting all elements to material redphi - use redphi= 0)
NOTE: redphi and nredphi are part of a special method for setting mechanical material properties, used primarily to create a bottom boundary layer. This sets the material properties for the mechanical calculations only, it does NOT change the actual color (material) number of the Eulerian elements. Set nredphi > total number of Eulerian elements to suppress this method. See note on using redphi and nredphi.
plas_fluid_method
Method for calculating plastic yielding for a frictional plastic incompressible material using a Drucker-Prager yield, for plane-strain and including various options for the effects of pore fluid pressure (see notes below).
plas_fluid_method Description Comments
0 Defaults to old behaviour in which phi values (see below) are interpreted as effective phi's modified by pore fluid pressure. In subroutine viscosities1, yielding is given by
  splas = epressure*sin(phi) + coh
1
Note: coh is not multiplied by cos(phi)
Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress).
If lambda_pore_prs <= 10**-30, then
  splas = epressure*sin(phi) + coh
else
  splas = (epressure - pwater)*(1-lambda_pore_prs)*sin(phi)+ coh

NOTE: Methods 1 and 3 have errors in the cohesion term. These methods are available for comparison with models that used the old behaviour. Methods 2 and 4 have the correct cohesion term.

Models using plas_fluid_method = 3 and plas_fluid_method = 4 may have problems converging in the first timestep unless lambda_blend1 and lambda_blend2 are used. The problem may be more pronounced with values of lambda_pore_prs greater than 0.4.

2 Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress).
If lambda_pore_prs <= 10**-30, then
  splas = epressure*sin(phi) + coh*cos(phi)
else
  splas = (epressure - pwater)*(1-lambda_pore_prs)*sin(phi) + coh*cos(phi)
3
Note: coh is not multiplied by cos(phi)
Assume lambda_pore_prs = fluid pressure/overburden weight. In the code, overburden weight is called plitho.
If lambda_pore_prs <= 10**-30, then
  splas = epressure*sin(phi) + coh
else
  splas = (epressure - (pwater + lambda_pore_prs*plitho))*sin(phi) + coh
4 Assume lambda_pore_prs = fluid pressure/overburden weight. In the code, overburden weight is called plitho.
If lambda_pore_prs <= 10**-30, then
  splas = epressure*sin(phi) + coh*cos(phi)
else
  splas = (epressure - (pwater + lambda_pore_prs*plitho))*sin(phi) + coh*cos(phi)
5 The effect of hydrostatic fluid pressures on epressure is taken into account. Method assumes that the fluid pressure is equal to the hydrostatic fluid pressure everywhere. Because fluid pressure is hydrostatic, lambda_pore_prs defaults to 1. (in code, hydrostatic weight is called phydro)
  splas = (epressure - (pwater + lambda_pore_prs*phydro))*sin(phi) + coh*cos(phi)
When you use plas_fluid_method=5, lambda_pore_prs1 and lambda_pore_prs2 must be set to 1 in SOPALE1_i.

NOTE: The hydrostatic pressure (phydro) is applied to ALL materials in the model. User beware!

Definitions
splas
is the magnitude of the second invariant of the stress required for plastic yielding.
epressure
is the mean stress calculated by the finite element solution. It includes the effect of water loading because waterweight is added to the body force vector in subroutine "compute_body_forces".
pwater
is the weight of the water (rho_water*d*-gravity) applied as a vertical load in "compute_body_forces".
plitho
is the weight (lithostatic pressure) of the model computed at the centre of the Eulerian elements. It DOES NOT include the water load.
phydro
is the weight of the model assumed to be composed of water (rho=1000kg/m**3) computed at the centre of the Eulerian elements. It DOES NOT include the water load.

Material Sets

For each material set (i = 1, n_matsets) include these lines of input.

color_numbers_for_material_set

color_numbers_for_material_set
This is a list of material numbers, indicating the material numbers which belong to this material set. The items in the list are separated by commas. Each item is a material number or a range of material numbers. Example: 1,5,17-21 This indicates that materials 1, 5, 17,18,19,20,21 belong to this material set.

density [ maximum_yield_stress [ mixing_reaction_scaling ]]

density
density of materials in material set i
maximum_yield_stress
maximum yield stress
Default: 1.0d12
Since: June 2014
mixing_reaction_scaling
mixing_reaction_scaling =1 optional flag only required if this material set is subject to weakening that results from mixing or reaction (metasomatism) in an arc_region. In such cases the material density changes as a consequence of mixing or reaction (see ARC section of User Guide). In this case the variable density (above) changes as the model evolves and is interpreted as the current density.

If (mixing_reaction_scaling. eq. 1) the evolving density is used to vary the scaling factor nlcreep_sc in the nonlinear (power-law) creep flow law (see input below) in which case the scaling factor varies linearly from nlcreep_sc1 to nlcreep_sc2 over the range that the density changes. ( No scaling will occur if nlcreep_sc1 =nlcreep_sc2.). The density scaling replaces the strain-weakening scaling of nlcreep_sc.

The scaling for strain weakening is
nlcreep_sc = nlcreep_sc1 - ((straine – nlcreep_strain1) / (nlcreep-strain2 –nlcreep_strain1)) x (nlcreep_sc1 – nlcreep_sc2) Equation 1

In the case of density weakening ( mixing_reaction_scaling) the scaling becomes
f_current = f_initial - ((rho3_current – rho3_initial) / (rho3_final – rho3_initial)) x (f_initial – f_final), Equation 2
so that the input parameters listed below are reinterpreted. Note the current_strain (straine) is replaced by the current_denisity (density)

If the density of this material set varies as a consequence of mixing or reaction its density will be updated automatically by the ARC module (if this is a restart the density = current_density will be read from the restart file)

Default: 0
Since: May 2015

bulkv visc_on(1) visc_on(2) visc_on(3) visc_on(4) [ bulkv_factor ]

bulkv
"bulk viscosity"
This is not actually a realistic value of bulk viscosity, but rather a numerical parameter used to maintain incompressibility of the system during solution by the penalty method. This value is typically about 7 orders of magnitude larger than the lowest viscosity in the model.

"visc_on(i), i=1,4" determines which of the viscous flow laws are active during the calculation.

visc_on(1)
Flag for flow law number 1; grain size sensitive creep.
0 => do not use.
1 => use.

visc_on(2)
Flag for flow law number 2; nonlinear creep.
0 => do not use.
1 => use.

visc_on(3)
Flag for flow law number 3; Peierl's law creep.
0 => do not use.
1 => use.

visc_on(4)
Flag for flow law number 4; pressure solution creep.
0 => do not use.
1 => use.

bulkv_factor
bulkv_factor is a number less than or equal to 1.0d0 (the default) used to multiply and rescale the bulk viscosity (bulkv) in Eulerian elements with this mechanical material when there is an enforced volume change as part of the mechanical finite element solution (ie delv_vdt_phase and/or delv_vdt_arc are not equal to 0.0d0). The scaling scaling makes the element more compressible and also reduces the nodal forces required to enforce the volume change. It does not apply to delv_vdt_compact.

Default: 1.d0
Since: May 2015

phi1 phi2 phie1 phie2 lambda_pore_prs1 lambda_pore_prs2 lambda_pore_prse1 lambda_pore_prse2

phi1
angle of friction in degrees corresponding to strain level phie1
phi2
angle of friction in degrees corresponding to strain level phie2
phie1
strain level 1
phie2
strain level 2, which must be greater than phie1
	      
	      Diagram:
	      
	                              ________________ phi2
	                             /|
	                            / |
	                           /  |
	                          /   |
	      phi1 ______________/    phie2 (strain level 2)
	                         |
	                        phie1 (strain level 1)    
lambda_pore_prs1
pore pressure ratio corresponding to strain level lambda_pore_prse1
lambda_pore_prs2
pore pressure ratio corresponding to strain level lambda_pore_prse2
lambda_pore_prse1
first strain level
lambda_pore_prse2
second strain level, which must be greater than lambda_pore_prse1
      
	      
                                  ________________ lambda_pore_prs2
                                 /|
                                / |
                               /  |
                              /   |
  lambda_pore_prs1 __________/    lambda_pore_prse2
                             |
                            lambda_pore_prse1

More info: see note on use of phi and lambda_pore_prs

See note on plastic yield criterion and pore fluid pressure

coh1 coh2 cohe1 cohe2

coh1
cohesion corresponding to strain level cohe1
coh2
cohesion corresponding to strain level cohe2
cohe1
first strain level
cohe2
second strain level, which must be greater than cohe1
        
	   Diagram:
	      
	                              ________________ coh2
	                             /|
	                            / |
	                           /  |
	                          /   |
              coh1 ______________/    cohe2
	                         |
	                        cohe1

GS gscreep_mu gscreep_m gscreep_T gscreep_strain1 gscreep_strain2 gscreep_a0 gscreep_a1 gscreep_a2 gscreep_updown_T gscreep_VdivR

[GRAIN SIZE SENSITIVE CREEP]
gscreep_mu
viscosity at level of grainsize a0
To suppress gs_creep in parallel, take gs_creep_mu = 0
NOTE: See Philippe's note on viscosities
gscreep_m
[2] power of grainsize dependency [m* in formula below]
gscreep_T
[3] temperature at which mu is (1/e) x that viscosity at a0 and T [T* in formula below]
gscreep_strain1
[4] S1* (see fig below)
gscreep_strain2
[5] S2* (see fig below)
NOTE: Always use S1* < S2*
gscreep_a0
[6] parameter a0* [ a0* in formula below ]
gscreep_a1
[7] parameter a1* [see fig below ]
gscreep_a2
[8] parameter a2* [see fig below ]
gscreep_updown_T
For case exp (- Td /T*), gscreep_updown_T = +1
For case exp (- T*/Td ) gscreep_updown_T = -1
	                               [2]
	              d                m*         d
	          mu(T  ) = mu* (a/a* )   exp (- T  /T*)
	                            0                [3]  
	                           [6]                
                 d = dynamical



	    [8] a1* _ _ _ _ _ _ _ _ _ ______________
	                             /|
	                            / |
	                           /  |
	                          /   |
	    [7] a2* _____________/    |
	                         |    |
	                         |    |
	                        S1*   S2*
	                        [4]   [5]
gscreep_VdivR
Parameter for adding effect of activation volume V to the creep laws for computing viscosities.
VdivR = V*/R, where V* is the activation volume and R is the universal gas constant 8.3144 J/(K x mol).
Use 0.d0 if you don't want to have effect of activation volume.

PL nlcreep_eps nlcreep_sc1 nlcreep_sc2 nlcreep_n nlcreep_T nlcreep_strain1 nlcreep_strain2 nlcreep_updown_T nlcreep_VdivR [ peierls_threshold peierls_factor ]

[NONLINEAR, 'POWER-LAW', CREEP]

nlcreep_eps
[1] in the formula below = Aps = ½ x 3**((n+1)/2) x Auni
Where Aps (A invariant form) is the experimental value of A, Auni, converted to the plane strain or tensor invariant form, assuming that the experiment that determined A was uniaxial deformation. If experiment was not uniaxial the conversion formula must be recalculated

In the nonlinear (power-law) creep formulation used here the second invariant of the stress (sigma2) is scaled by nlcreep_sc. This scaling factor represents how the stress and viscosity are modified by comparison with the standard flow law (nlcreep_sc = constant) by the following factors

1) no scaling of the flow law nlcreep_sc = constant (nlcreep_sc1 = nlcreep_sc2 = constant);

2) strain weakening in which nlcreep_sc varies linearly from nlcreep_sc1 to nlcreep_sc2 over the range of second invariant of the strain nlcreep_strain1 to nlcreep_strain2;

1) or 2) are assumed unless the input has the optional input mixing_reaction_scaling = 1 (see density input)

3) weakening in which the above variables are interpreted in terms of the effects of mixing or reaction in an arc region controlled by mixing or reaction of this mechanical material (mat3) with an added (injected) second material (mat2) that progressively changes the density of this material as more of this material is added (see description of mixing and reaction in ARC section).

If (mixing_reaction_scaling. eq. 1) then this material set is weakened by the effect of mixing or reaction such that nlcreep_sc varies linearly from nlcreep_sc1 to nlcreep_sc2 over a range of density defined by the inputs nlcreep_strain1 = density (the density at which weakening begins) and nlcreep_strain2 is the density beyond which no further density scaling occurs.

Note the inputs for the mechanical material properties (listed here) may, or may not be consistent with the densities and value of rho3_increase in the reaction part of the ARC section. Consistency requires that nlcreep_strain2 = initial_density and nlcreep_strain1 = initial_density + rho3_increase. However, the user can define a range of density increase that is different from the range defined in the ARC section.
nlcreep_sc1
[2] in the formula below = initial strain weakening or if ( mixing_reaction_scaling . eq.1) initial scaling factor ( f_initial)
nlcreep_sc2
[3] in the formula below = final strain weakening or if ( mixing_reaction_scaling . eq.1) final scaling factor ( f_final)
(Note these can also be used to give strain, or mixing_reaction_scaling hardening)
nlcreep_n
[4] in the formula below = power-law exponent n (ie value of the power in the power law flow equation)
nlcreep_T
[5] in the formula below = T* = Q/R where Q is the activation energy (kJ/mol) and R the universal gas constant, 8.3144 J/(K x mol)
nlcreep_strain1
[6] in the formula below = value of strain for the onset of strain weakening, or if ( mixing_reaction_scaling . eq.1) value of density at which scaling starts (rho3_initial)
nlcreep_strain2
[7] in the formula below = value of strain at which strain weakening ends, or if ( mixing_reaction_scaling . eq.1) value of density at which scaling ends (rho3_final)
nlcreep_updown_T
[8] in the formula below = -1, 1, or 0 (see below) Use -1 for normal calculations
nlcreep_VdivR
[for versions later than code-jul26-04-testV]
Parameter for adding effect of activation volume V to the creep laws for computing viscosities.
[9] in equation above
VdivR = V*/R, where V* is the activation volume (SI units) and R is the universal gas constant 8.3144 J/(K x mol).
Use 0.d0 if you don't want to have the effect of activation volume.
peierls_threshold
If the calculated sigma2 is larger than this, it is multiplied by peierls_factor. After the multiplication, sigma2 is forced to be at least peierls_threshold.
Since: July 2014
Restrictions: must be positive.
Default: 1.d20
peierls_factor
Since: July 2014
Restrictions: must be positive, and less than or equal to 1.d0
Default: 1.d0

The formula below is the flow law written as strain rate (edot) calculated from stress (sigma) where both are the square root of the second invariants.


	               [1]             [4]                    [9]
	                |               n*                     |
	       edot =  Aps x   (sigma/sc)    exp –((T* + P x VdivR)/ T)
	                              |             |
                                   sc1=[2]         [5]    if updown= -1
                                   sc2=[3]


	       [2]  sc1* _____________ 
                        	      |\
                        	      | \
                        	      |  \
                        	      |   \
        	[3]  sc2* _ _ _ _ _ _ | _ _\__________________
                        	      |    |
                        	      |    |
                        	      S1*  S2*
 
                       	              [6]  [7]

P is the lithostatic pressure (Pa) and T is temperature (K) both calculated by Sopale.

The actual formulae used in Sopale are the inverse of the equation above, ie the stress Sigma2 is calculated from the strain rate (square root second invariant). See subroutine viscosities1_2


 if (nlcreep_updown_T.eq.-1) then
  sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n)) *
  dexp((nlcreep_T + plitho * nlcreep_vdivr) / (nlcreep_n * tlg(l)))
 endif
 if (nlcreep_updown_T.eq.1) then
  sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n)) *
  dexp((tlg(l))/(nlcreep_n * (nlcreep_T + plitho * nlcreep_vdivr))
 endif
 if (nlcreep_updown_T.eq.0) then
  sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n))
 endif

PL pei_eps pei_Tn pei_sc pei_T2 pei_strain1 pei_strain2

[PEIERL'S LAW CREEP]

pei_eps
[1] in diagram below
reference strain rate
pei_Tn
[2] in diagram below
[use 0 to suppress in parallel]
pei_sc
[3] in diagram below
pei_T2
[4] in diagram below
pei_strain1
not used
pei_strain2
not used
	 PEIERL'S LAW POWER-LAW BREAKDOWN 
	  NOTE: E = epsilon
	  
	                           [2]
	         .   .             Tn*/T          [4]
	         E = E* (sigma/sc*)       exp ( - T2* / T ) 
	             [1]       [3]

PS eps p d sc1 sc2 t updown_t strain1 strain2

PRESSURE SOLUTION CREEP

eps
is -
p
is -
d
is -
sc1
is -
sc2
is -
t
is -
updown_t
is -
strain1
is -
strain2
is -
         .   .  
       ( e / eps ) = (sigma/sc) * exp -(T*/ T  ) * (1/ (T . dp ) )

                  |
             sc1  |____________ 
                  |            |\
                  |            | \
                  |            |  \
                  |            |   \
             sc2  |_ _ _ _ _ _ | _ _\__________________
                  |            |    |
                  |            |    |
                  +------------+----+---------------------------
                           strain1  strain2

texpand tdensity melt_weaken temp_solidus temp_melt_weakened visc_melt_weakened [ visc_bingham ]

texpand
volume coefficient of thermal expansion.
See notes on Mantle Thermal Expansivity.
Since version 2013-02-08_1.0: If negative, the value will be ignored and 4 parameters for a variable coefficient of thermal expansion will be read from the next line.

temperature1 temperature2 alpha1 alpha2

temperature1
For temperatures <= temperature1, the coefficient of thermal expansion will be alpha1.
temperature2
For temperatures between temperature1 and temperature2, the coefficient of thermal expansion will be a linear interpolation between alpha1 and alpha2.
For temperatures => temperature2, the coefficient of thermal expansion will be alpha2.
tdensity
reference temperature for density (in absolute degrees)
density is
rho(T) = rho0 (1 - alpha(T - tdensity))
melt_weaken
'0' => no melt weakening
'1' => material melt weakens
'2' => material melt weakens with extra parameters
Since June 2014: if melt_weaken=2, the extra parameters will be read from the next line.
The pressure and temperature pairs (0, melt2_t0), (melt2_p1, melt2_t1), (melt2_p2, melt2_t2) define a line in P/T space. That line is the solidus of the material. The liquidus is solidus + melt2_t_int.

melt2_t0 melt2_t1 melt2_p1 melt2_t2 melt2_p2 melt2_t_int

temp_solidus
temperature (degrees Kelvin) at which melt weakening starts
temp_melt_weakened
temperature (degrees Kelvin) at which melt weakening is complete
visc_melt_weakened
linear viscosity of melt-weakened material
visc_bingham
bingham viscosity. Since version 2009-07-26_22.2. Optional. Default is 0.0


END OF INPUT FOR EACH MATERIAL'S MECHANICAL PROPERTIES

type_bc

type_bc
Defines the type of boundary condition to apply to the model.

The following table summarizes the type_bc values. For detailed descriptions of the type_bc input, please see the document on "SOPALE Boundary Condition Input"

type_bc
Input variables and summary description
300

vxleft_up vxleft_down yleft1 yleft2 rollleft vxright_up vxright_down yright1 yright2 rollright vytop0 vxbottom0 reference_pressure_column reference_material [ pump_time [ isos_time_scale [ delta_vel_factor [ delta_vel_toler [ start_time flag ]]]]]]

The recommended, updated version of type_bc = 3.

Defines horizontal velocities at side boundaries (can have inflow and outflow); can set top boundary to rigid top (no vertical movement); can set bottom boundary to no-slip (no horizontal movement)

Other features: side roller boundaries are more consistently implemented than in type_bc=3; "pressure pump" feature adjusts horizontal velocities at side boundaries to keep the pressure within the model constant. Details for type_bc = 300.

Since version 2009-07-26_21.2 an optional parameter: pump_time
Since version 2009-07-26_23.0 an optional parameter: isos_time_scale
Since version 2009-07-26_29.0 an optional parameter: delta_vel_factor
Since version 2012-06-08_2.0 an optional parameter: delta_vel_toler. See type_bc 301 for details.
Since July 2014; two optional parameters

  • start_time
  • flag
start_time
If present, sets the start time for the boundary condition definition on this line.
flag
'more' => indicates there will be another boundary condition. Spelling flexible provided the word is not 'last'.
'last' => indicates this is the last boundary condition.
301
(current)

(Since version 2013-01-29_1.0)

An extension of type 300.

In type 300, constant velocities are applied to the top and bottom of each side of the model. The velocities are choosen so that the mass (or volume, if density is assumed to be constant) in the model remains constant.

In type 301, the velocities applied to the top of each side of the model are adjusted every time step so that the forces applied approach a 'target_force'. The bottom velocities are then adjusted so that the mass (or volume, if density is assumed to be constant) in the model remains constant. Forces are measured from the top down to yleft1 on the left side, and from the top down to yright1 on the right side.

Details for type_bc=301

vxleft_up target_force_left yleft1 yleft2 rollleft

vxleft_up
Initial value of vxleft_up (like bc_type 300)
target_force_left
target force (from the top down to yleft1) for the left side.
'ft' in the calculation of vxleft_up
(others like bc_type 300)

vxright_up target_force_right yright1 yright2 rollright

vxright_up
Initial value of vxright_up (like bc_type 300)
target_force_right
target force (from the top down to yright1) for the right side.
'ft' in the calculation of vxright_up
(others like bc_type 300)

vytop0 vxbottom0

(like bc_type 300)

reference_pressure_column reference_material [ pump_time [ isos_time_scale [ delta_vel_factor [ delta_vel_toler ]]]]

The parameters in square brackets are optional after version 2013-01-29_1.0. In version 2013-01-29_1.0 they are mandatory.

reference_pressure_column
(like bc_type 300)
reference_material
(like bc_type 300)
pump_time
(like bc_type 300)
isos_time_scale
(like bc_type 300)
delta_vel_factor
(like bc_type 300)
delta_vel_toler
When the pressure pump calculates delta_vel, the calculation includes subtracting plithob in the reference column from the average plithob. If abs ((ref_plithob - ave_plithob) - 1.0) < delta_vel_toler, the calculation is skipped and delta_vel is set to zero.

mantle_wind_vel fraction_left

mantle_wind_vel
mantle_wind_vel (in the usual SI units). This velocity is added to the velocites in the lower region ....below the lithosphere, after all of the other velocities have been set using the rescaling velocity. Note that the addition of this velocity is with the same sign on both sides....in one side, out the other.

fraction_left
Determines how the bottom velocities are distributed between the left and right sides of the model.

'separate': Each side will be balanced separately. After vxleft_up is re-calculated, vxleft_down will be set so that volume_in = volume_out on the left side. The right side calculation will be similar.

<number>: The volume change due to the upper velocities will be distributed (in the lower region) between the left and right sides of the model.
  • fraction_left * volume_change will go into, or out of, the left side.
  • (1 - fraction_left) * volume_change will go into, or out of, the right side.
0.5 will give equal volumes on each side.

use_tectonic_forces force_column force_process [ filter_length filter_alpha ]

use_tectonic_forces
= 1: use tectonic forces.
force_tect = total horizontal driving force (force_dyna) - mean stress (force_epress).
See .out file for values of force_tect, force_dyna and force_epress
= 0: use total forces
force_column
column number (measured from both left and right side) to be used when calculating the forces. Column 1 may give very jittery results. Column 10 is suggested.
force_process
After the force type (use_tectonic_forces) and the column number (force_column) have been selected, there is a choice of how to process the result before using it in the velocity rescaling formula.
'raw': use the force calculated from the current time step.
'filtered': use the output of a low-pass time series filter, where the time series is the boundary force of the previous filter_length time steps. The filter parameter filter_alpha can be used to set the filter sensitivity.
filter_length
Default: 20
filter_alpha
Initial tests suggest that a value 0.2 gives some lag, and much smoothing, whereas 0.4 gives less lag and less smoothing.
Default: 0.4

vxleft_up_min vxleft_up_max vxleft_up_velocity_scale

vxleft_up_min, vxleft_up_max
In each timestep after vxleft_up is re-calculated, the result is compared to these minimum and maximum velocities and set to one of these values if 'out of range'.
vxleft_up_velocity_scale {in version 2013-01-29_1.0: vxleft_up_factor}
'Vd' in the calculation of vxleft_up.

vxright_up_min vxright_up_max vxright_up_velocity_scale

vxright_up_min, vxright_up_max
In each timestep after vxright_up is re-calculated, the result is are compared to these minimum and maximum velocities and set to one of these values if 'out of range'.
vxright_up_velocity_scale {in version 2013-01-29_1.0: vxright_up_factor}
'Vd' in the calculation of vxright_up.

ksplitx1 tsmooth ismooth leftsmooth rightsmooth

ksplitx1
number of subintervals per Eulerian delta x at which markers are embedded for top surface tracking (use 2, or < 4)
tsmooth
number of timesteps between two numerical smoothings of the top surface, i.e. smoothing is done every tsmooth timesteps. On a restart, if align_time=0, smoothing will be done if ismooth is > 0. NOTE: in the timestep when smoothing is done, smoothing is done ncyc2 times (see subroutine vylifts).
** CAUTION ***

Surface smoothing may not be consistent with all types of tectonic and surface processes. In earlier versions surface smoothing is applied to the vertically advected Eulerian tectonic surface and then surface and other processes are applied to the smoothed surface.

in the nov21-02 version of the code and later, surface smoothing is applied to the final surface after the effects of advection owing to the surface processes have been added.

Numerical surface smoothing is best avoided if possible. Instead, it is preferable to use the diffusion option in the surface processes - see subroutine erode. If smoothing is all that is required (not true surface processes) ensure that the material properties of any injected Lagrangian particles that correspond to sediment are the same as those of the material that comprises the surface that is being smoothed.

ismooth
number of smoothing passes per smoothing timestep
Recommended value = 1
leftsmooth
Eulerian column number defining the left edge of the portion of the surface to be smoothed.
The minimum value is 2.
rightsmooth
Eulerian column number defining the right edge of the portion of the surface to be smoothed.
The maximum value is nx - 1.

vscale verror verror1 miterv irescale kinitvis force_iter1 force_iter force_continue mech_trouble_iter

vscale
velocity scale for scaling velocity errors in Picard iterations.
Convergence is defined to be achieved when
maxdv/vscaleverror
where maxdv is the greatest change in velocity (in either x or y direction) since the previous iteration
verror
relative error for timesteps greater than 1 in Picard iteration
verror1
relative error for timestep 1 in Picard iterations
miterv
maximum number of iterations in Picard iterations. If convergence is not achieved by miterv iterations, the model will stop unless force_continue is set greater than 0.
Actually, the code always does one more iteration than miterv. This could be called a bug. However the behaviour (bug) has existed for so many years that it seems the lesser of two evils to leave the code as is.
irescale
scheme used for updating viscosities
0 = calculate viscosities according to material properties if stress above yield viscosities are rescaled
[viscosity rescaling created by strain rate (see earlier notes)]
kinitvis
maximum number of iterations per timestep before the viscosity is reinitialized and computation is started from scratch. Should always check by setting kinitvis = 1
The above is the original documentation. At present (2010-04-16) kinitvis is a one-shot parameter; viscosities will be initialized when the iteration number = kinitvis, in addition to when they would otherwise be intialized.
force_iter1
number of iterations to force for the first timestep
This number of iterations will always be performed even if convergence is achieved earlier. If convergence has not been reached after the forced iterations, the code will continue to iterate until convergence is reached or until the maximum iterations (miterv) is reached.
0 = don't force a minumum number of iterations (normal behaviour)
force_iter
number of iterations to force for timesteps beyond the first timestep in the Picard iterations.
This number of iterations will always be performed even if convergence is achieved earlier. If convergence has not been reached after the forced iterations, the code will continue to iterate until convergence is reached or until the maximum iterations (miterv) is reached.
0 = don't force a minumum number of iterations (normal behaviour)
force_continue
This is a flag to indicate whether to force continued timesteps when Picard divergence occurs.
0 = if Picard divergence occurs, stop the model (normal behaviour)
1 = if Picard divergence occurs, don't stop the model but continue to next timestep anyway. USE WITH CAUTION! ONLY TO BE USED IN SPECIAL CASES!
mech_trouble_iter
iteration number at which to write an output file. This output file is given frame number 99, which means that its name will be something like "g01_p00_f99_o" or "g02_p00_f99_o". This file is useful for debugging problem models where convergence is not reached on the last timestep and the restart file doesn't give any clues why.

ncyc2 color2op lagtest eulcolorrule eul_density_rule

ncyc2
number of cycles in Lag/Eul remeshing/tracking per step
Suggested value = 2
color2op
sopale_nested reads and ignores this.
lagtest
sopale_nested reads and ignores this.
eulcolorrule
rule to use for coloring Eulerian cells from Lagrangian particles
0 = use old method: color Eulerian cell the same as the last Lagrangian particle found in the cell. Note that this will bias the material colors in favour of Lagrangian particles with higher node number.
1 = use "majority rules": Eulerian cell is colored according to the majority (actually plurality) of Lagrangian particles; i.e. the color which is represented by the most number of Lagrangian particles in the Eulerian cell. If there is not a plurality, then the Eulerian cell's existing color is not changed.
eul_density_rule
1 = calculate cell mass using the densities of its Lagrangian particles.
Note that each Lagrangian particle has two densities; one used in the finite element calculations, the other calculated by the compaction module.
2 = calculate cell mass using eulerian cell material densities. Obsolete.
Must be 1 if there are any material phase change rules.
Since 2016-06-30: Must be 1.

ucompute tcompute ale_compute dt_compute tbottom ttop e_temperature cfl lambda_blend1 lambda_blend2

ucompute
sopale_nested stops if this is not 1.
tcompute
1 = compute a temperature field using an FE computation
0 = no FE computation
Use 0 if you are running a mechanical-only model (no thermal-mechanical coupling)
ale_compute
sopale_nested stops if this is not 1.
dtcompute
sopale_nested reads and ignores this.
tbottom
temperature at bottom boundary (in Kelvin)
ttop
temperature at surface (in Kelvin)
e_temperature
2014-02-03: disabled. sopale_nested reads and ignores this.
cfl
Courant–Friedrichs–Lewy factor; used to adjust timestep length (dt) when dt_adjust = 'yes'. suggest approx 0.2, no greater than 0.3
lambda_blend1
Iteration number during the first timestep where the calculation of plastic yielding starts to incorporate the value of lambda_pore_prs (fluid pressure ratio)
lambda_blend2
Iteration number during the first timestep where the calculation of plastic yielding has finished incorporating the full value of lambda_pore_prs (fluid pressure ratio).

lambda_blend1 and lambda_blend2 define an interval during the first timestep where you want to "blend in" the value of lambda_pore_prs (fluid pore pressure ratio) in the calculation of plastic yielding (splas).

In timestep 1, starting at iteration lambda_blend1, the value of lambda_pore_prs is gradually increased from zero until it reaches its full value at iteration number lambda_blend2.

This technique avoids problems with negative values of plastic yielding in the first timestep, when the model has not yet stabilized. Normally used in conjunction with plas_fluid_method (see above).

NOTE: You should make sure that force_iter1 is at least as large as lambda_blend2.

steady analy esteady

steady
sopale_nested reads and ignores this.
analy
sopale_nested reads and ignores this.
esteady
sopale_nested reads and ignores this.

Thermal section (only needed if tcompute = 1)

maxter dtempm ther_advection th_gamma_time th_cn tempsmooth

maxter
maximum error for thermal iterations
For steady=0, this is a relative error
For steady=2, this is an absolute error
dtempm
temperature error (units are temperature). If temperature error more than this, continue iterating
MAY NOT BE USED (value is set in code)
ther_advection
1 = there is advection of the temperature field
0 = there is no advection of the temperature field
th_gamma_time
numerical parameter for integration in time
Choose between .5 and 1, typically 1
Trapezoidal schemes for time integration
th_cn
Crank-Nicolson parameter
Required value = 1
tempsmooth
timestep interval between successive smoothings of temperature field

type_bct

type_bct
boundary condition type for thermal problem during the model run.
1 = convection - top and bottom.
2 = convection - left and right sides.
3 = for subduction and deep model thermo-mechanical (dmtm)

Only if type_bct is 3

bct3(1) bct3(2) bct3(3) bct3(4) bct3(5) bct3(6) bct3(7) bct3(8) bct3(9) bct3(10) bct3(11) [ start_time flag ]

bct3(1)
top temperature (degrees K) [see tshift in the code]
bct3(2)
0 => bottom is imposed temperature, value bct3(3)
1 => bottom is imposed flux, value bct3(3)
Since version 2009-07-26_10:
2 => bottom is imposed temperature that varies with space and time. Read extra input lines below.
3 => bottom is imposed flux that varies with space and time. Read extra input lines below.
bct3(3)
temperature or flux at base (depending on bct3(2)).
Since version 2009-07-26_10:
If bct3(2) = 2 or 3 the temperature or flux varies in space and time, and bct3(3) is ignored. Read input below.
bct3(4)
0 => On left boundary, impose initial geotherm above bct3(5) and bct3(7) for flux below.
1 => On left boundary, impose bct3(6) as temperature above bct3(5) and bct3(7) as flux below.
2 => On left boundary, impose bct3(6) as temperature above bct3(5) and initial geotherm below.
Since November 2014
3 => On left boundary, impose a depth-varying temperature from the surface down to the bottom of the range defined by the depth - temperature pairs. Below that, apply bct3(7) as flux. Read the depth - temperature pairs below.
bct3(5)
Position on left boundary about which the boundary condition changes (see bct3(4)) - distance in metres measured from origin upwards. However, when bct3(4) is 3, the point where the boundary condition changes is determined by the last of the depth - temperature pairs.
bct3(6)
Temperature on left boundary when bct3(4) is 1 or 2. Not used when bct3(4) is 0 or 3.
bct3(7)
Lateral flux imposed at left boundary below bct3(5) when bct3(4) is 0 or 1 or 3. Not used when bct3(4) is 2.
bct3(8)
0 => On right boundary, impose initial geotherm above bct3(9) and bct3(11) for flux below.
1 => On right boundary, impose bct3(10) as temperature above bct3(9) and bct3(11) as flux below.
2 => On right boundary, impose bct3(10) as temperature above bct3(9) and initial geotherm below.
Since November 2014
3 => On right boundary, impose a depth-varying temperature from the surface down to the bottom of the range defined by the depth - temperature pairs. Below that, apply bct3(11) as flux. Read the depth - temperature pairs below.
bct3(9)
Position on right boundary about which the boundary condition changes (see bct3(8)) - distance in metres measured from origin upwards. However, when bct3(8) is 3, the point where the boundary condition changes is determined by the last of the depth - temperature pairs.
bct3(10)
Temperature on right boundary when bct3(8) is 1 or 2. Not used when bct3(8) is 0 or 3.
bct3(11)
Lateral flux imposed at right boundary below bct3(9) when bct3(8) is 0 or 1 or 3. Not used when bct3(8) is 2.

Since July 2014; two optional parameters
  • start_time
  • flag
start_time
If present, sets the start time for the boundary condition definition on this line.
flag
'more' => indicates there will be another boundary condition. Spelling flexible provided the word is not 'last'.
'last' => indicates this is the last boundary condition.

Only if bct3(2) = 2 or 3
Since version 2009-07-26_10
Space-time variation of basal temperature (when bct3(2) = 2) or heat flux (when bct3(2) = 3)

num_space_points

= number of points in the piecewise function defining the space/time variation of the basal temperature or flux.

x_base_thermal_positions(i), i=1,num_space_points

= x positions of piecewise function defining space/time variation of the basal temperature or flux.

val_thermal_bc(i), i=1,num_space_points

= values of temperature or flux in the piecewise function.

num_time_points

= number of points in the piecewise function defining the time scaling variation.

time_values(i), i=1,num_time_points

= time values (in seconds) of the piecewise function defining the time scaling variation.

val_time_values(i), i=1,num_time_points

= values of the piecewise function defining the time scaling variation.

Only if bct3(4) = 3
Since November 2014

These are the upper left side thermal boundary conditions. A sequence of lines, each one containing a depth and temperature, terminationed by a line with [end].

depth(1) temperature(1)

depth(1)

The depths are measured from the current surface, (i.e start at zero) and increase downwards. Since the surface height will, in general, vary with time, the vertical distance over which the lower part of the boundary conditions - flux - is applied will vary.

temperature(1)

The temperatures are in degrees K.

depth(2) temperature(2)

...

[end]

Only if bct3(8) = 3
Since November 2014

These are the upper right side thermal boundary conditions. A sequence of lines, each one containing a depth and temperature, terminationed by a line with [end].

depth(1) temperature(1)

depth(1)

The depths are measured from the current surface, (i.e start at zero) and increase downwards. Since the surface height will, in general, vary with time, the vertical distance over which the lower part of the boundary conditions - flux - is applied will vary.

temperature(1)

The temperatures are in degrees K.

depth(2) temperature(2)

...

[end]


An example. All lines that begin with a '!' are treated as comments.

 273. 2 21.d-3 3 1210000. 0. 0. 3 1210000. 0. 0.   0.d0  last
!     ^
!     |  indicates the need for 
!  1) space-varying bottom temperature boundary conditions
!  2) time-varying factors for the bottom boundary conditions
!
!  1) the space varying bottom temperature boundary conditions
! --------------------------------------------------------------------
!  the number of points
5
! the x positions
0      100.d3   500.d3  1250.d3  1750.d3
! the temperatures at those positions
2000   2300      2500   2300      2000
!
!  2) the time-varying factors for those bottom temperatures
! --------------------------------------------------------------------
! the number of time points
3
! the time values, in seconds
0.0    3.1536d11    6.3072d11
! the factors
1.0     0.5         0.1
!
! left side depth - temperature profile
! --------------------------------------------------------------------
!  273. 2 21.d-3 3 1210000. 0. 0. 3 1210000. 0. 0.   3.1536d11  notlast
!                ^
!                |  indicates the need for left side depth - temperature pairs
! The left-side depth temperature pairs
 0.d3      273
 120.d3   1600
 250.d3   1700
 550.d3   1800
 800.d3   1900
 1050.d3  2000
 [end]
!
! right side depth - temperature profile
! ---------------------------------------------------------------------
!  273. 2 21.d-3 3 1210000. 0. 0. 3 1210000. 0. 0.   3.1536d11  notlast
!                                 ^
!                                 |  indicates the need for right side depth - temperature pairs
! The right-side depth temperature pairs
 0.d3      273
 120.d3   1600
 250.d3   1700
 550.d3   1800
 800.d3   1900
 1050.d3  2000
 [end]
! ----------------------------------------------------------------------

th_dt was_cfl vzgiven

th_dt
timestep for thermal computation
Is computed internally when steady=0 and dtcompute=1 (steady=0 means non steady state calculation)
NOTE: normally th_dt should be the same as dt. For code-feb07-06-adiabat and later, if th_dt is set to 0.d0 or a negative number, it is reset by the code to be the same as dt. If th_dt is a positive number, but different from dt, the code will print a warning message.
was_cfl
not used.
vzgiven
only used for benchmarks 20, 21, 22, 23
not used in deep thermo-mechanical models

ncolor2t added_sed_colort

ncolor2t
number of thermal materials
added_sed_colort
thermal material number to use for injected sediment particles
Since version 2009-07-26_18:
If present, and positive, this overrides sedcolort(i).

For each thermal material (i = 1, ncolor2t) we have the following 5 lines of input:

thermod1_conductivity(i) OR cond1(i) cond2(i) temp1(i) temp2(i)
thermod1_rau0(i)
thermod1_cp(i)
thermod1_fheating(i)
thermod1_power_permass(i)

Said another way...

Choose one column or the other
thermod1_conductivity(i) cond1(i) cond2(i) temp1(i) temp2(i)
thermod1_rau0(i) thermod1_rau0(i)
thermod1_cp(i) thermod1_cp(i)
thermod1_fheating(i) thermod1_fheating(i)
thermod1_power_permass(i) thermod1_power_permass(i)

thermod1_conductivity(i) OR cond1(i) cond2(i) temp1(i) temp2(i)

If there is only one number on this line, it will be read as:
thermod1_conductivity(i)
thermal conductivity (K) for thermal material i
units are W/(m K) or user units

Since version 2009-07-26_10:
If there is more than one number on this line, it will be read as:
cond1(i) cond2(i) temp1(i) temp2(i)
cond1(i) and cond2(i) are thermal conductivities (K). Units are W/(m K) or user units
temp1(i) and temp2(i) are temperatures. Units are degrees Kelvin.

These four numbers define a conductivity that varies with the element temperature. For temperatures less than temp1(i), the conductivity is cond1(i). For temperatures greater than temp2(i), the conductivity is cond2(i). For temperatures between temp1(i) and temp2(i), the conductivity is interpolated.
This triggers iteration of the initial steady state temperature calculation.

thermod1_rau0(i)

thermod1_rau0(i)
thermal density of thermal material i (density used by thermal code)
units kg/m3 or user units
POTENTIAL INCONSISTENCY: For accurate calculations, the thermal and mechanical material properties should be consistent. Make sure that the densities are the same.

thermod1_cp(i)

thermod1_cp(i)
Specific heat for thermal material i
units J/(kg K) or user units
POTENTIAL INCONSISTENCY: For accurate calculations, the thermal and mechanical material properties should be consistent. Choose the value of the specific heat Cp appropriately.

thermod1_fheating(i)

thermod1_fheating(i)
frictional (strain) heating coefficient for thermal material i
Fraction of the internal mechanical work that is converted to heat (can also be interpreted as the efficiency of conversion of mechanical energy to heat).
WARNING: Strain heating couples the mechanical and thermal parts of SOPALE. The internal mechanical rate of work is calculated per unit volume of the mechanical Eulerian elements and stored in array dissip(6,nel) where nel=Eulerian element number. This rate of internal dissipation (Watts/m**3) is added to the internal heating in the thermal calculation per unit volume of the thermal Eulerian elements. To be consistent, the mechanical and thermal Eulerian element meshes (and nodal grids) must be the same.
NOTE: Radioactive heating is specified as rate of energy generation/unit mass (Watts/kgm) in input and converted to rate of energy generated per unit volume (Watts/m3) by multiplying by the material thermal density (thermod1_rau0(i) for thermal material i). This conversion is not required for the internal dissipation because dissip(6,nel) is dissipation per unit volume. The array element dissip(7,nel) stores the dissipation in the mechanical Eulerian elements per unit mass, where the density used in the conversion is the mechanical material dnesity.

thermod1_power_permass(i)

thermod1_power_permass(i)
internal heat production (e.g. radiactive heat generation) for thermal material i
units watts/kg
NOTE: need to convert from A in units of watts/m3 to H in units of watts/kg. A/rau0 = H

OUTFIL

OUTFIL
sopale_nested reads and ignores this.

MATFIL

MATFIL
sopale_nested reads and ignores this.

ICASE

ICASE
sopale_nested reads and ignores this.

CACHSZ

CACHSZ
sopale_nested reads and ignores this.

LEVEL

LEVEL
sopale_nested reads and ignores this.

Thermal section (only needed if tcompute = 1)

MATFIL

MATFIL (for thermal systems)
sopale_nested reads and ignores this.

ICASE

ICASE (for thermal systems)
sopale_nested reads and ignores this.

CACHSZ

CACHSZ (for thermal systems)
sopale_nested reads and ignores this.

LEVEL

LEVEL (for thermal systems)
sopale_nested reads and ignores this.

ibinread

ibinread
sopale_nested reads and ignores this.

filename1

filename1
sopale_nested reads and ignores this.

filename2

filename2
sopale_nested reads and ignores this.

filenameg1

filenameg1
sopale_nested reads and ignores this.

filenameg3

filenameg3
sopale_nested reads and ignores this.

filenameg4

filenameg4
sopale_nested reads and ignores this.

mechanical color assignment

Here we assign a initial mechanical color (not material) to the lagrangian particles.

ncolor2 [ s1 s2 ... ]

ncolor2
number of Lagrangian mechanical color boxes defined on Lagrangian grid
s1 s2 ...
seed used to initialize the pseudo-random number generator. On the scrXm family and the stable, when using binaries compiled with Intel's ifort, this should be two integers, maximum 9 digits each. Other compilers may need a different number. gfortran's pseudo-random number uses a seed with eight integers. Eight is the maximum number allowed for at present.
Only used if strain2_type is 'random'.
Default: 123456789 362436069 521288629 316191069 987654321 458629013 582859209 438195021
Since version 2009-07-26_29.0

For each box (i = 1, ncolor2), there are 3 lines of input. Coordinates are in meters, time is in seconds.

x1 y1 x2 y2

x1 y1
coordinates of upper left corner
x2 y2
coordinates of lower left corner

x3 y3 x4 y4

x3 y3
coordinates of lower right corner
x4 y4
coordinates of upper right corner

newcol2 [ strain2_type strain2_magnitude [ runtime-flag time ]]

newcol2
the mechanical color number assigned to all Lagrangian particles in this box
strain2_type
'constant: The strain of the Lagrangian particles will be initialized strain2_magnitude.
'random': The strain of the Lagrangian particles will be initialized to a number from the pseudo-random number generator in the range -strain2_magnitude to strain2_magnitude.
Default: 'constant'
Since version 2009-07-26_29.0
strain2_magnitude
magnitude used when initializing the strain of Lagrangian particles.
Default: 0.d0
Since version 2009-07-26_29.0
runtime-flag
'F' => the box will be applied at model initialization.
'T' => the box will be applied at time='time'
Default: F
Since version 2014-03-??
time
Time (in seconds) at which the box will be applied. Only used if 'runtime-flag' is 'T'
Default: 0.0
Since version 2014-03-??

thermal material assignment

Here we assign a initial thermal material to the lagrangian particles.

You need to define at least one thermal material box, even if you're not doing thermal calculations. If you are doing a mechanical-only model, define a dummy thermal material box which covers your entire domain and consists of thermal material 1. (You don't have to define a thermal material).

ncolor2

ncolor2
number of thermal material boxes defined on Lagrangian grid

For each box (i = 1, ncolor2), there are 3 lines of input. Coordinates are in meters, time is in seconds.

x1 y1 x2 y2

x1 y1
coordinates of upper left corner
x2 y2
coordinates of lower left corner

x3 y3 x4 y4

x3 y3
coordinates of lower right corner
x4 y4
coordinates of upper right corner

newcol2 [ runtime-flag time ]

newcol2
the thermal material number assigned to all Lagrangian particles in this box
runtime-flag
'F' => the box will be applied at model initialization.
'T' => the box will be applied at time='time'
Default: F
Since version 2014-03-??
time
Time (in seconds) at which the box will be applied. Only used if 'runtime-flag' is 'T'
Default: 0.0
Since version 2014-03-??

dmtm n_skip_thermal

dmtm
should be 1
n_skip_thermal
number of timesteps where the thermal calculations are skipped (i.e. if n_skip_thermal = 50, then the code will not perform thermal calculations for the first 50 timesteps)

Only if dmtm = 1, then read the following:

type_bct0 [ temperature_initial maxter_i [ max_th_iters ] ]

type_bct0
Defines how the initial Eulerian temperature field is calculated.
0 => the initial temperature field is computed based on the inputs in the ETEMP section of sopale1_plus_i. This gives a 2D initial steady-state temperature variation with depth.
3 => use a finite element calculation to compute the initial steady-state temperature field.

Note that
  • the ETEMP values for temperature will be overwritten if type_bct0=3.
  • the type_bct0=3 temperatures will be overwritten if EPERTURB-T is present.
  • For more about temperature initialization, see here.

    Since version 2009-07-26_10

    temperature_initial
    initial temperature for the steady state thermal solution.
    Since version 2009-07-26_22: If this is negative, the ETEMP field values are used as the initial temperature for the steady state thermal calculation.
    The default is -1.0
    maxter_i
    The steady state thermal solution is declared converged when the maximum change in temperature is below this.
    Since version 2009-07-26_22:
    Default is 1.d0 .
    max_th_iters
    maximum number of thermal iterations. Used to limit both the initial steady state and time-stepping thermal solution. Default value is 100.

    Only If type_bct0=3, the next line defines parameters for the boundary conditions to use when solving for the initial steady-state temperature field.

    bct30(1) bct30(2) bct30(3) bct30(4) bct30(5) bct30(6) bct30(7) bct30(8) bct30(9) bct30(10) bct30(11)

    bct30(1)
    top temperature (degrees K) [see tshift for the code]
    bct30(2)
    = 0 bottom is imposed temperature, value bct30(3)
    = 1 bottom is imposed flux, value bct30(3)
    Since version 2009-07-26_10:
    = 2 bottom is imposed temperature that varies with space. Read extra input lines below.
    = 3 bottom is imposed flux that varies with space. Read extra input lines below.
    bct30(3)
    = temperature or flux at base (depending on value of bct30(2).
    Since version 2009-07-26_10:
    If bct30(2) = 2 or 3 the temperature or flux varies in space, and bct30(3) is ignored. Read input below.
    bct30(4)
    = 0 On left boundary, impose initial geotherm above bct30(5) and bct30(7) for flux below bct30(5)
    = 1 On left boundary, impose bct30(6) as temperature above bct30(5) and bct30(7) as flux below bct30(5)
    = 2 On left boundary, impose bct30(6) as temperature above bct30(5) and initial geotherm below bct30(5)
    bct30(5)
    = z Position on left boundary about which the boundary condition changes (see bct30(4)) - distance in metres or user units measured from origin upwards
    bct30(6)
    = T Temperature on left boundary when bct30(4) =1 or 2. Not used when bct30(4) = 0.
    bct30(7)
    = flux Lateral flux imposed at left boundary below bct30(5) when bct30(4) = 0 or 1. Not used when bct30(4) = 2.
    bct30(8)
    = 0 On right boundary, impose initial geotherm above bct30(9) and bct30(11) for flux below bct30(9)
    = 1 On right boundary, impose bct3(10) as temperature above bct30(9) and bct30(11) as flux below bct30(9)
    = 2 On right boundary, impose bct30(10) as temperature above bct30(9) and initial geotherm below bct30(9)
    bct30(9)
    = z Position on right boundary about which the boundary condition changes (see bct30(8)) - distance in metres or user units measured from origin upwards
    bct30(10)
    = T Temperature on right boundary when bct30(8) =1 or 2. Not used when bct30(8) = 0.
    bct30(11)
    = flux Lateral flux imposed at right boundary below bct30(9) when bct30(8) = 0 or 1. Not used when bct30(8) = 2.

    NOTE no geotherm
    Full calculation of diffusion/production thermal equilibrium
    Note need for consistent temperature input - use Kelvin

    Only if bct30(2) = 2 or 3
    Since version 2009-07-26_10
    Spacial variation of basal temperature (bct30(2) = 2) or heat flux (bct30(2) = 3)

    num_space_points

    = number of points in the piecewise function defining the spacial variation of the basal temperature or flux.

    x_base_thermal_positions(i), i=1,num_space_points

    = x positions of piecewise function defining spacial variation of the basal temperature or flux.

    val_thermal_bc(i), i=1,num_space_points

    = values of temperature or flux in the piecewise function.

    ntperts

    ntperts
    number of boxes for temperature perturbation
    if zero, no boxes (below) will be read.

    For each box, read in coordinates of corners and temperature perturbation value (tpert). The value of tpert is added to the initial temperature values. For more about temperature initialization, see here.

    xx1 xx2 xx3 xx4

    xx1
    x coordinate of top left corner of box (units: metres)
    xx2
    z coordinate of top left corner of box
    xx3
    x coordinate of bottom left corner
    xx4
    z coordinate of bottom left corner

    xx5 xx6 xx7 xx8

    xx5
    x coordinate of bottom right corner
    xx6
    z coordinate of bottom right corner
    xx7
    x coordinate of top right corner
    xx8
    z coordinate of top right corner

    tpert

    tpert
    Add this to the temperatures inside the box.


    This page was last modified on Friday, 01-Jul-2016 13:50:40 ADT
    Comments to geodynam at dal dot ca