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

Originally these parameters were all surface-related. New sections have been added for processes which aren't surface processes.

The parameters for each process are defined within a "section". The start of a section is marked with a keyword. Within a section, input parameters must be in a particular order, but sections may be in any order within the file.

Section keywords:


ADIABATIC section

This section defines parameters for adiabatic heating.

ADIABATIC

adiabatic_heat

adiabatic_heat
= 1 for adiabatic heating
= 0 for no adiabatic heating

ARC section

The ARC module simulates a volcanic arc using multiple, stacked, arc regions.

ARC

Indicates begin of ARC section

An arc region of the model is defined by lagrangian particles whose mechanical property number is ‘color1’. This is mat1, the original material in the arc region. More than one arc region may be defined, each with its own properties.

Arc regions are specified by material and thermal ‘boxes’ defined in the normal manner in the section of input where these boxes are defined. These boxes will be considered normal Sopale material and thermal boxes unless they contain arc lagrangian particles with mechanical material property ‘color1’ (mat 1). Within the boxes only eulerian elements that contain ‘color1’ lagrangian particles are treated as arc eulerian elements. In cases where eulerian elements contain a mixture of arc (color1) and non-arc lagrangian particles a fraction of the element is calculated based on the fraction of arc lagrangian particles, and arc processes are applied to this fraction of the eulerian element and the arc lagrangian particles. Other materials are unaffected by arc processes

Arc regions evolve by the progressive addition of material whose mechanical property number is ‘color2’. This is mat2 and usually takes the form of an infiltrating melt (e.g. basaltic melt) to the arc region or a metsomatic fluid to other regions.

mat2 is added using a fractional rate of change function (time vs dvdt). That is, dvdt is the fractional volume of mat2 added (delta_volume per unit volume of the arc region/per unit time) (units 1/s). The rate of mass additions is dvdt x mat2_density. The dvdt(t) time function for each arc region is defined in the input.

As mat2 is added the density and volume of mat1 evolve. mat1 therefore becomes a ‘composite material’ reflecting the addition of mat2 to mat1

Conservation of heat

mat2 is added at a high temperature and is assumed to cool to the ambient temperature of the arc region during the timestep. During cooling it releases heat, both specific and latent heat (if it solidifies). The heat released is used to heat the arc region comprising the evolving product of the addition of mat2 to mat1. This consumes specific heat and latent heat if mat1 melts.

The temperature of the 4 corner nodes of each eulerian element are increased by ¼ of the elemental temperature increase calculated by the heat balance.

The properties of mat2 and mat1(including the way it evolves as the composite) are taken to be constant. The user must select a consistent set of properties that represent the average temperature dependence of mat2 and a similar set of properties for the evolving mat1 that represent both its average temperature dependence and the average way its properties change as mat2 is added.

Mixing and Reaction Models

The evolution of mat1 (the composite) as mat2 is added can be either that of a Mixing model or a Reaction model.

Mixing model: volume of mat1 increases by the volume of mat2 that is added and the chemical density evolves as the average density of mat1 (the composite) as mat2 is added. Mass is conserved.

Reaction model: mat2 is assumed to react with mat1 to produce the evolving composite. The evolving chemical density of the reaction product (composite) is specified as input as a function of the fractional amount of mat2 that has been added. Mass is conserved. Therefore the volume of the composite evolves in a manner that conserves mass but it is not the sum of the mat1 and mat2 constituent volumes.

The changing chemical density of the lagrangian particles of mat1 (the composite) is stored as the updated chemical density of the material set that contains color1 (mat1). Owing to the changing density, material sets containing arc color1 (mat1) mechanical materials should be defined separately from other material sets and have only one member (color1 number) in the set.

In both models the volume of each eulerian element (cell) is modified each timestep to reflect the evolving volume of the composite. The volume change is enforced using nodal forces to expand or contract the element.

Restarts:

  1. All material densities are reset using the material densities in the input file.
  2. The density of any material containing the color1 of an arc region is reset from the value in the arc restart file.

The rate of change of volume is specified by the time and dvdt values. The function is a step function; it keeps its previous value until the current time exceeds one of the t(i) values, and then the value changes to v(i).

There should be dvdt value for time = 0.0

Initial values of dvdt can be calculated:

  • choose a percentage volume change, e.g. volchge = 1% = 1e-2
  • dvdt = volchge / dt

Each arc region is terminated by [end]. The ARC section also needs a terminating [end]

color1 color2 dens2 temp2

color1
color1 - the mechanical material set color number of the arc material (initially mat1, but evolving as mat3)
color2
color2 - the mechanical color number of the added-injected material (mat 2). If non-positive, no injection will take place. Non-positive values allowed since early April 2015.
dens2
dens2 - density of mat2, the added material. Used in the thermal calculations and in calculating the mixing density of the composite (mat3). Note mat2 is commonly added as a melt and therefore dens2 would be the melt density. However, to simplify the calculation, mat2 is assumed to be added at its solid density (but at its melt temperature). This is justified because the injected melt cools and solidifies rapidly and, therefore, the mixing or reaction between mat2 and mat3 is considered to take place in the solid state.
temp2
temp2 - temperature at which mat2 is added (Kelvin). Used in the thermal calculations.

Properties of the two materials. Only used in the thermal calculations. Thermal properties of mat1 (the composite) used in the finite element calculation are specified in the normal input

Cp1 L1 T1l T1i T1s X1h X1l

properties of mat1 (the composite) used to calculate heat absorbed by the composite and the associated increase in temperature.

Cp1
specific heat. Use an average that represents both the solid and melt, if necessary.
L1
latent heat. Use an average that represents both the solid and the melt, if necessary.
T1l
liquidus temperature
T1i
intermediate transition temperature
T1s
solidus temperature
X1h
fraction of melt that solidifies during cooling from Tl to Ti
X1l
fraction of melt that solidifies during cooling from Ti to Ts

Cp2 L2 T2l T2i T2s X2h X2l

properties of mat2 (the added material) used to calculate heat added to the system

Cp2
specific heat. Use an average that represents both the solid and melt, if necessary.
L2
latent heat. Use an average that represents both the solid and the melt, if necessary.
T2l
liquidus temperature
T2i
intermediate transition temperature
T2s
solidus temperature
X2h
fraction of melt that solidifies during cooling from Tl to Ti
X2l
fraction of melt that solidifies during cooling from Ti to Ts

Specification of addition of mat2 The rate of addition of mat2 is specified by the time and dvdt values. The function is a step function; it keeps its previous value until the current time exceeds one of the t(i) values, and then the value changes to v(i).

There should be dvdt value for time = 0.0

Initial values of dvdt can be calculated:

  • choose a percentage volume change, e.g. volchge = 1% = 1.00d-2
  • dvdt = volchge / dt

time(1) time(2) time(3) ...

Times (in seconds)

dvdt(1) dvdt(2) dvdt(3) ...

rate of volume change per second per unit volume for cells of color1. units are m3 per second per m3. These values correspond to those at time(1), time(2), time(3) above.
Volumes are calculated by area (m2) x unit length of strike direction (m). The rate of addition is scaled by the initial volume of the arc region, the volume containing color1 lagrangian particles. This means that the total volume rate of addition of mat2 will remain constant while dvdt is constant. Note, however, the arc region volume may expand owing to the addition of mat2, which means that each eulerian arc element receives progressively less mat2 as the volume of arc region increases. This approach is consistent with a known rate of addition of mat2 which is distributed uniformly through the evolving arc region. Units of dvdt are volume units (m3) per second per m3.

[reaction_model]

Indicates the begin of reaction_model inputs. Optional. If not present, the mixing model will be used.

this is the reaction model In the reaction case the added mat2 reacts with mat1 to produce a composite (mat3) with a specified progressively changing density. The associated change in the volume of the composite is not known (and is not the volume of mat1 plus mat2), Instead, it must be calculated from the mass balance before and after the reaction: Mass balance equation is:

rho3_final x vol3_initial ( 1 + delta_vol3) = rho3_initial x vol3_initial + rho2 x vol2 (eq. 1)

initial and final are the initial and final values at the beginning and end of the timestep rho3_final is the final density of mat3 (the composite) at the end of the timestep from equation above

vol3_initial is the intial volume of mat3 (the composite) at the beginning of the timestep

delta_vol3 is the fractional change in volume of mat3 (the composite) during the timestep, that is the actual volume change divided by vol3_initial

delta_rho3 is the change in the density of the composite from the beginning to end of the timestep and is specified in terms of the volume of mat2 that is added, delta_vol2 (see eq.4 below)

rho2 is the density of mat2

vol2 is the absolute volume of mat 2 added this timestep

delta_rho3 = rho3_final - rho3_initial

Solving for delta_vol3

delta_vol3 = (delta_rho3/ rho3_final) + (rho2/ rho3_final) x (vol2 / vol3_initial) (eq. 2)

where delta_vol3 is the fractional volume change of the composite. This can also be written as

delta_vol3 = (delta_rho3/ rho3_final) + (rho2/ rho3_final) x delta_vol2 (eq. 3)

where delta_vol2 = (vol2 / vol3_initial) the fractional volume of mat 2 added this timestep with respect to vol3_initial. Note eq.3 takes into account that the effect of adding vol2 of mat2 in causing a change delta_vol3 in the reactant (composite) progressively declines as vol3 increases. This is because vol2 becomes a smaller proportion of the composite.

delta_rho3 must be calculated from the input which specifies the total change in rho3 (rho3_increase) that results from adding a total fractional volume of mat2 (f_tot).

Specify how the density of mat1 (the composite) evolves as mat2 is added

f_tot rho3_increase

f_tot
This is the total fractional volume of mat2 that must be added so that the final density of the composite changes from dens1 by an amount rho3_increase. (For example 0.3d0 means that 30% of mat2 by volume must be added to mat 1 in order that the composite (mat3) changes density by rho3_increase.)
rho3_increase
rho3_increase This is the total change in density of the composite (mat3) (SI units) as a result of reacting with f_tot of added mat2. Increasing density is positive.

It is assumed that the density of the composite changes linearly as mat3 is added and reacts

For each timestep the change in composite density (delta_rho3) is given by

delta_rho3 = rho3_increase x f / f_tot; where f = vol2 / vol3_initial (eq. 4)

That is adding vol2 of mat2 initially changes the density in proportion to the initial volume of mat1 (vol1) but this volume evolves as vol3_initial at the beginning of each timestep. Therefore adding vol2 of mat2 each timestep has progressively less effect on the density change during the reaction because this volume of mat2 represents a progressively smaller fraction of vol3, the volume of mat3.

The process of adding mat2 to the arc region is the same for both the Mixing and Reaction Models. That is, the volume of mat2 that is added per timestep is specified. This volume is distributed uniformly within the arc region initially occupied by mat1 and later occupied by mat3 (the composite). Only the consequences of Mixing vs Reaction of mat2 with mat3 are different.

[end]

Indicates end of reaction_model.

[end]

Indicates end of an arc region.

[end]

Indicates end of the ARC section.

Since: sopale_nested version 2013-11-28_1.0
Changed: sopale_nested version 2013-11-28_2.0
Changed: sopale_nested version 2014-01-26_1.0 added thermal properties
Changed: May 2015. Added reaction model.


ARC_DELV section

The cell volume changes that result from density changes in an arc can be applied either as forces when solving for velocities in the FE mechanical calculation, or added as vertical velocities after the FE mechanical calculation.

ARC_DELV

Indicates start of the ARC_DELV section.

value

Indicates where the volume changes should be applied.
force : apply the cell volume change as nodal forces acting in the FE mechanical calculation. This is previous behaviour.
velocity : apply the cell volume change as a vertical velocity after the FE mechanical calculation. Cell volume changes accumulate from the bottom upward.

Default: force

Since: June 2015


COMPACTION section

Since sopale_nested version 2009-07-26_14

COMPACTION

Indicates start of the COMPACTION section.

n_compact_mat

number of materials to compact

compaction_method

0 => no compaction
2 => compaction by using volume change to modify the y components of the velocity solution.

compaction_rho_fluid

density of fluid to use in compaction computations

For each material to be compacted.....

compact_matset_num(i)

material number

compact_or_contract(i)

0 => COMPACT
1 => CONTRACT

rho_grain(i) surface_porosity(i) compaction_coeff(i)


EROSION section

This section starts with the word EROSION (on a line by itself).

EROSION

erosion_model

erosion_model
type of erosion model you want to use
= 0 no erosion. No other inputs in the section will be read.
= 1 combinations of total erosion and total deposition depending on nerode_type
= 2 slope-dependent erosion/deposition
= 3 diffusive erosion (added by Susanne Buiter)
= 4 combines climate-modulated slope-dependent erosion above reference height (refh) and partial/total filling of sediment accommodation space below the reference height. Added in code-may27-04-ero4.

erosion_rate OR erode_p OR dummy

erosion_rate (If erosion_model=2)
erosion rate per second when the surface slope is 1

erode_p (If erosion_model=3)
diffusion coefficient k (in m2/s)

dummy (If erosion_model=1,4)
not used

nerode_type   nerode_ends

Each nodal position, absolute height h = y(i) along the top surface of the model is potentially subject to erosion (removal of material) or deposition (adding of material). Erosion occurs when y(i) .gt. refh. Deposition occurs when y(i). lt .refh (refh is defined in the file sopale1_plus_i). nerode_type specifies when erosion and deposition will occur.

sopale uses the value of refh to determine whether erosion or deposition is done. If the surface is above refh and erosion is turned on, then the surface will be eroded. Similarly, if the surface is below refh and deposition is turned on, then material will be deposited.

nerode_type
= 1 erosion above refh, no erosion or deposition below refh.
= 2 erosion and deposition: both erosion above refh and deposition below refh.
= 3 deposition below refh and no erosion above refh.
nerode_ends
= 0: No erosion of surface endpoints.
  • This is the recommended value when deformation in the model will not propagate to the boundaries.
  • removed version 2013-06-28_1.0. Use '3' or '5'.
= 1: Erodes points 1 and nx by same amount as points 2 and nx-1, respectively.
  • This gives no-slope boundaries at the edges for surface processes.
  • removed version 2013-06-28_1.0. Use '4' or '6'.
= 3:
  • In LS, no erosion of surface endpoints.
  • In any nest, the endpoints are set to values (same position) from LS.
= 4:
  • In LS, erode points 1 and nx by same amount as points 2 and nx1-1, respectively.
  • In any nest, the endpoints are set to values (same position) from LS.
= 5:
  • In LS, no erosion of surface endpoints.
  • In any nest, an extra value (interpolated from the LS surface) on each side is used, and the concept of an 'endpoint' is gone.
  • Since version 2013-09-17_1.0
= 6:
  • In LS, erode points 1 and nx by same amount as points 2 and nx1-1, respectively.
  • In any nest, an extra value (interpolated from the LS surface) on each side is used, and the concept of an 'endpoint' is gone.
  • Since version 2013-09-17_1.0

er_stablax   method_slope_avg

only for erosion_model = 2, 4

er_stablax
= Lax stabilization factor
The larger this value is, the more diffusive stabilization is applied.
method_slope_avg
= 1 calculate average slope by taking the absolute values of the slopes and then averaging
= 2 calculate average slope by taking the slope between points 1 and 3.
Do not use method_slope_average = 2 because it underestimates local slopes and consequently the amount of erosion.

In erosion models 2 and 4, slope dependent erosion may lead to node by node (n.b.n) oscillations in the surface (development of a saw tooth surface) which may further lead to 'failure and slumping' of the positive saw teeth.

To reduce the tendency for saw teeth to grow (ie an instability to develop) a small component of numerical diffusion is added to the surface processes in erosion models 2 and 4. The theory behind this Lax stabilization is in W. H. Press et al., Numerical Recipes, p.623-629 (Edition 1) and is implemented in the erosion models using equation 17.1.18. The stabilization allows the erosion to satisfy the Courant condition for stability of the surface.

er_stablax is the diffusion coefficient which determines the amount of diffusion and therefore the numerical smoothing of the surface.

The diffusion of the surface used to stabilize slope-dependent erosion is calculated by the explicit forward differencing finite difference method.

y(i,n+1) - y(i,n) = (dt/dx2) * y(i+1,n) - 2y(i,n) + y(i-1,n)

Where y(i) are the vertical positions of the surface nodes...either the nodes of the FE grid (traditional method) or the surface string (string method) and n is the timestep number.

The left side is the amount of diffusive erosion at node y(i) between timestep n and n+1.

dt is the timestep length (secs)

dx is the horizontal spacing between nodes (metres) (FE grid for traditional and particle spacing on surface string for string method)

The above equation is only convergent/stable for dt/ dx2 < ½

This means that for most choices of dt and dx the above formula has to be multiplied by a factor, er_stablax to ensure stability. This is Lax stabilization. That is

dt/dx2 * er_stablax < ½

or

er_stablax > (dx2) / 2 dt

Example:

1000 yr timestep => dt = 3.15 x 1010 sec
dx = 1km = 106m
gives er_stablax ~ 10-5

Note this is the maximum value for stability based on nominal values of dt and dx which the user chooses.

However, in cases where dt varies to satisfy the CFL criterion the maximum value of dt should be used in calculating er_stablax. Smaller values will be stable.

Similarly if dx varies the minimum value of dx should be used to calculate er_stablax.

It is important to remember that numerical diffusive stabilization is only a secondary surface process. The physical erosion is the slope-dependent part. The best choice of er_stablax is that it satisfies the above criterion (the maximum value possible for stability) but that the user value should be as small as possible (factor of 10-2 to 10-3 less if possible) so that diffusion is minimized but the slope-dependent erosion is stabilized and no ``saw teeth`` develop. Essentially, the numerical diffusion is designed to keep the growth of the saw teeth under control

Values if er_stablax that are too large may lead to excessive surface smoothing.

It is IMPORTANT to NOTE that Lax stabilization is only turned on where erosion and deposition are turned on. This means that when using nerode_type 1 or 3 the parts of the surface respectively below refh and above refh will not be subject to Lax stabilization and may develop n.b.n oscillations.

It is therefore recommended to use nerode_type =2, unless users have a special need otherwise.

When users only want surface diffusion and no slope-dependent erosion in erosion models 2 and 4, set erosion_rate = 0 and choose and then TEST a suitably small value of er_stablax. Use nerode_type =2 if diffusion is required everywhere.

In erosion model 4 the climate modulation does not affect the numerical diffusion of the surface. It only modulates the slope-dependent erosion and the sediment deposition.

Erosion Model 4

These inputs are for erosion model 4.

nslope

nslope
= number of points in the definition of the piecewise linear time-varying function for slope-dependent erosion rate

time_slope_dep_rate(i)   slope_dep_rate(i)

time_slope_dep_rate(i)
= time in seconds
slope_dep_rate(i)
= erosion rate per second when the surface slope is 1 at the specified time

nsedfil

nsedfil
= number of points in the definition of the piecewise linear time-varying function for sediment fill factor.

time_sedfil_rate(i)   sedfil_rate(i)

time_sedfil_rate(i)
= time in seconds
sedfil_rate(i)
= the rate of sedimentation, in m/s, at time time_sedfil_rate(i), in the accommodation space available for sedimentation below refh. The actual sedimentation rate will be scaled by the accomodation space available.
In the code:
      get_erosion = scalefactor * curr_sedfil_rate * dt_e

where

      scalefactor = (dabs(h_e - refh)/1.0d0)

(the 1.0d0 is in units of metres, to make the scale factor dimensionless)

This means that sedimentation rates will be greater in basins where the accommodation space is large and less where there is little accommodation. The sedimentation rate, sedfil_rate, varies according to the time function time_sedfil_rate.

NOTE: 1 cm/yr = 3.17097d-10 m/s

nclimate climate_flag

climate_flag < 10: climate function type 1

In climate function type 1 the spatial variation of the climate function is entirely specified by the user using the inputs listed below.

nclimate
= number of points in the definition of the piecewise linear spatial function for climate modulation of erosion rate.
climate_flag
= 0 no climate modulation
= 1 climate modulation is only applied to erosion rate
= 2 climate modulation is only applied to sedimentation rate (sedfil_rate)
= 3 climate modulation is applied to both erosion rate and sedimentation rate

climate_position(i)   climate_value(i)

climate_position(i)
= x position
climate_value(i)
= climate modulation factor at x=climate_position(i)
Depending on the value of climate_flag (see above), the erosion rate, or the sedimentation rate, or both, is/are multiplied by climate_value.

climate_flag ≥ 10: climate function type 2

Since: April, 2016

Climate function type 2 is an orographic climate function designed for orogen-type models where there is one high elevation orogen. The climate function decreases across the left facing orogenic front to simulate the effect of decreasing orographic precipitation between the orogenic foreland and the elevated platueau which is considered to be in a rain shadow, as is the rest of the model to the right of this location.

The climate function is tuned such that starting at left Eulerian grid boundary the climate_value(1) is 0. It then increases linearly to climate_value(2) (recommended 1.0d0) at a distance delta_climate_position1 to the left of the orogenic front. It then decreases linearly across the orogenic front (region of high elevation) to a position within the orogen at distance delta_climate_position2 to the right of the orogenic front. At this position the climate function has value climate_value(3) (recommended 0.0d0- 0.05d0) and remains at this value to climate-position(4), the right boundary of the Eulerian grid.

In Climate function type 2 Sopale recalculates some input values of the climate function to take into account the evolving model surface configuration, in particular the advection of the position of the orogenic front (x_position_half_surface_max)

Best to use Lagrangian string definition of surface to ensure selection of a single point for x_position_half_surface_max for both LS and SS regions.

nclimate
= ignored
climate_flag
= 10 no climate modulation
= 11 climate modulation is only applied to erosion rate
= 12 climate modulation is only applied to sedimentation rate (sedfil_rate)
= 13 climate modulation is applied to both erosion rate and sedimentation rate

climate_position(1)   climate_value(1)

climate_position(1)
= 0.0d3 (or x-position of left boundary of model Eulerian grid)
climate_value(1)
= 0.0d0

climate_position(2)   climate_value(2)

climate_position(2)
= user choice (Recommend 200.0d3, recalculated by Sopale, see below)
climate_value(2)
= user choice (Recommend 1.0d0)

climate_position(3)   climate_value(3)

climate_position(3)
= user choice (Recommend position in first high elevation area, recalculated by Sopale, see below)
climate_value(3)
= user choice (Recommend 0.0d0- 0.05d0)

climate_position(4)   climate_value(4)

climate_position(4)
= position of right boundary of model Eulerian grid
climate_value(4)
= user choice (Recommend 0.0d0)

delta_climate_position1 delta_climate_position2 ntime_step_climate2

delta_climate_position1
= distance left of the orogenic front for location of climate_position(2) (positive number)
delta_climate_position2
= distance right of the orogenic front for the location of climate_position(3)
ntime_step_climate2
= number of timesteps between calls to subroutine_climate_function _type_2 to recalculate values of climate_position(2) and climate_position(3)

er_max_n

Since: September 6, 2014

The inputs here apply a limit to the calculated erosion value at each point on the eulerian grid. The limits are applied only when the top of the eulerian grid is above refh. The calculated erosion value is usually negative (unless overcome by the er_stablax term (stab * er_stablax * dt), so the er_max_v(i) values should be negative in order to be effective. The limiting erosion rate er_max_v(i) will be multiplied by dt and then compared to the calculated erosion value.

The dt value used will the value in effect for the current ALE cycle.

The points below define a time varying maximum erosion rate. The maximum erosion rate used at any time will be interpolated from these values.

er_max_n
number of points in the time varying limit to effective erosion rate

er_max_t(i) er_max_v(i)

er_max_t(i)
time of this point
er_max_v(i)
maximum effective erosion rate

EROSION_OUTPUT section

EROSION_OUTPUT

Indicates start of section.

X, the surface, and the average erosion will be written to a file at time steps which are multiples of ts_incr. Averaging of erosion first starts at time step ts_start; subsequently it starts after file output. Both X and the surface values are from the time step at which the file is written.

The 'dat' output files can be be plotted using gnuplot; the first plot command below will generate a plot of the surface. The second will generate a plot of the erosion. gnuplot can be daunting, inscrutable, immensely capable, and rewarding. For details, see
/usr/share/doc/gnuplot-doc/html/index.html

dguptill@scram:-sopale_nested$ gnuplot
gnuplot> plot './A19_10DEPCML_erosion_p00_ts000010.dat' using 1:2 with linespoints
gnuplot> plot './A19_10DEPCML_erosion_p00_ts000010.dat' using 1:3 with linespoints

Since sopale_nested version 2011-12-01_7.0

ts_start ts_incr ts_end file_type

The output file name is <nameout1>erosion_p<nn>_ts<tttttt>.<file_type>, where:

nn is the process number (00 for LS, 01 for the first nest, etc)
tttttt is the time step

ts_start
Erosion averaging will be started at this time step.
ts_incr
Files will be written at multiples of this time step.
ts_end
Erosion averaging will terminate at this time step.
file_type
'dat' (without the quotes)

ISOSTASY section

This section defines parameters for isostasy. The implementation of flexural isostasy in SOPALE is based on the "beam code" in microfem (an older model code). Explanation of beam, from Juliet Hamilton's notes on microfem.

isost, n_beam_elem, load_extended_beam

  isost 
   = 1  use Airy (local) isostasy
   = 11 use Airy (local) isostasy with tangential boundary conditions
   = 2  use flexural isostasy
   = 21 use flexural isostasy with tangential boundary conditions

  n_beam_elem  (only needed for versions code-apr11-05-extbeam and later)
     = the number of elements in the beam (for flexural isostasy)
       This allows the user to make the beam extend beyond the Eulerian
       grid on the right-hand side.  For example, if your Eulerian grid
       is 100 elements wide, then setting n_beam_elem to 200 makes the
       beam twice as wide as your Eulerian grid.  The extra beam length
       extends to the right of your Eulerian grid.

       NOTE: if you are using an extended beam and erosion, it is 
       recommended that you set nerode_ends to 0. (nerode_ends is found 
       in the EROSION section of surfaceproc_i.)
       You should also check to make sure that your definition of 
       the flexural rigidity (see below) covers the entire length of the
       extended beam.

  load_extended_beam (only needed for versions code-apr13-05-extbeamtest and later)
     = 1  Load the extended part of the beam with the same loads as on
          beam element nx-1 (i.e. the rightmost element in the Eulerian grid).
     = 0  Set the load on the extended part of the beam to 0 (zero).

     NOTE: The extended beam is always loaded for the initial isostatic 
     adjustment.  For any isostatic calculations during the model run, the
     extended beam is loaded according to the value of load_extended_beam.


rho_below
   Density of the material below the model (i.e. density of material in
   which the model "floats")

If isost=2 or isost=21 (indicating flexural isostasy) then you need the
following additional lines:

nbreak
   Column number for node which defines a break in the beam.
   Beam 1 is the right-hand beam extending from columns nbreak+1 to nx1.
   Beam 2 is the left-hand beam extending from columns 1 to nbreak.
   If you just want one beam, you can use nbreak = 0; this will make
   the entire beam Beam 1.

ibbc(1,k), k=1,4
   Type of boundary conditions for Beam 1 (to the right of nbreak, 
   excluding nbreak).
      ibbc(i,1) = deflection or force for left end of beam i
                  0 = imposed deflection  (Dirichlet)
                  1 = imposed force       (Neumann)
      ibbc(i,2) = rotation or moment/torque, left end of beam i
                  0 = imposed rotation       (Dirichlet)
                  1 = imposed moment/torque  (Neumann)
      ibbc(i,3) = deflection or force for right end of beam i
                  0 = imposed deflection   (Dirichlet)
                  1 = imposed force        (Neumann)
      ibbc(i,4) = rotation or moment/torque right end of beam i
                  0 = imposed rotation       (Dirichlet)
                  1 = imposed moment/torque  (Neumann)

bbc(1,k), k=1,4
   Boundary conditions corresponding to the BC flags in the line above,
   e.g. if ibbc(1,1)=1, corresponding to an imposed deflection, then
   bbc(1,1) would be the value of the imposed deflection.

ibbc(2,k), k=1,4
   Boundary conditions for Beam 2 (to the left of nbreak, including nbreak).
   These are defined the same way as for Beam 1 (see above).
   Beam 2 has an additional possible value for ibbc(2,3) corresponding
   to the right-hand end of Beam 2.
   
   ibbc(2,3) = -1 means that Beam 2 is slaved to Beam 1; the deflection
                  for this end of Beam 2 must have the same deflection
                  as the left-hand end of Beam 1.
                  Note that this imposed value will result in an incorrect
                  overall isostatic balance of loads on Beam 2 in the 
                  region of nbreak because the slaved deflection
                  overrides the isostatic balance of Beam 2.
   
bbc(2,k),k=1,4  
   Boundary conditions corresponding to the BC flags in the line above.
   These are defined in the same way as for Beam 1.

nflex
    Number of intervals in the piecewise linear definition of
    the lateral variation of the flexural rigidity of the beams.
    These intervals are defined along the beam, left to right.

For k= 1, nflex, we have:

xflxl(k),xflxr(k),tmpflx(k)
    xflxl(k)  = x coordinate of left end of interval
    xflxr(k) =  x coordinate of right end of interval
    tmpflx(k) = value of flexural rigidity within the interval
    
    Sopale will use linear interpolation to fill in the flexural
    rigidity values for any sections of the beam which are not covered
    by the intervals defined above.
    Note that the flexural rigidity is constant within any of the
    above defined intervals.  You will need nflex .ge. 2 if you want 
    linear interpolation of the flexural rigidity, i.e. not a constant
    value of the whole interval.

ncfb
    Number of lines defining additional loads
    These loads are in addition to the weight of the model, which is 
    calculated and applied as equivalent nodal loads on the finite 
    element beam. The additional nodal loads represent additional forces 
    other than isostatic (buoyancy) forces.

For i=1,ncfb, we have:

node,cfb(2*node-1),cfb(2*node)
     node = column number where the additional beam load is placed
     cfb(2*node-1)  = vertical (z) component of additional load at this node
     cfb(2*node)    = additional torque load at this node



Example:

ISOSTASY
21      ! isost (1=Airy,11=Airy with tangential bcs, 2=flexural, 21=flexural with tangential bcs)
3.3d3    ! rho_below  
0        ! nbreak
1 1 1 1         ! (ibbc(1,k),k=1,4)
0.d0 0.d0 0.d0 0.d0    ! (bbc(1,k),k=1,4)
0 0 1 1     ! (ibbc(2,k),k=1,4)
0.d0 0.d0 0.d0 0.d0    ! (bbc(2,k),k=1,4)
1   ! nflex
0.d0  500.d3  1.d22   ! xflxl(k),xflxr(k),tmpflx(k)
1    ! ncfb
1 0.d0 0.d0   ! node,cfb(2*node-1),cfb(2*node)


MATCHANGE section

This section defines parameters for material phase changes.

The traditional (prior to version 2009-07-26_28.0) matchange code uses cell (Eulerian elemental pressures and temperatures interpolated from the nodal values) pressure and temperature when applying matchange rules to a Lagrangian particle.

When sopale_nested is running in nested mode, cell sizes are usually quite different (LS vs SS); hence the cell pressures and temperatures are different. Hence the result of any matchange rule can be different when cell pressures and temperature are used.

Aside from the undesirable model results, this breaks the restart code. For details, see the (closed) ticket describing different temperatures on a restart.

In version 2009-07-26_28.0 (and later), when applying a matchange rule to a Lagrangian particle, there is the option of using a pressure and temperature that are calculated for that particle. When this option is used, the result of the rule should be the same, no matter whether the rule is applied in LS or SS.

Each Lagrangian particle already has a temperature (calculated by interpolating Eulerian nodal values - t1), so we can use that.

Lagrangian Particle Pressures associated with MATCHANGE

It only remains to select, calculate, and save a pressure for each Lagrangian particle. This is a new array.

When comparing with nature, materials at depth experience the dynamical pressure (mean stress). Therefore the preferred pressure for matchanges that are functions of P,T conditions is mean stress (nodpress contains the Eulerian nodal values)

However, Sopale pressures are sometimes subject to a checkerboard mode which may make nodpress values unreliable and subject to checkerboard errors and flickering between timesteps. Users should use nodpress for the Lagrangian particle pressure unless they experience the checkerboard mode errors. In the latter case it is acceptable to use plitho which is the lithostatic pressures.

MATCHANGE

Indicates start of section.

num_rules_in_file [ PT ]

num_rules_in_file
The number of phase-change rules which the user wants to define. A phase-change rule must specify two material numbers and the conditions under which one material changes to the other. (The materials must have their properties defined in SOPALE1_i.) The conditions for phase change are defined with a polygon in Pressure vs Temperature.
PT
Determines which pressure and temperature are used when applying matchange rules to a Lagrangian particle.

'dynamical'
Uses dynamical pressure, interpolated to each Lagrangian particle. Uses Lagrangian particle temperature. The dynamical pressures can have significant errors, for example from 'checkerboard mode'. Therefore dynamical should only be used with CAUTION and the user needs to check the results carefully and plot dynamical pressures.
'plitho'
use lithostatic pressure, interpolated to each Lagrangian particle. Use Lagrangian particle temperature.
'traditional'
use the cell (Eulerian elemental) pressure and temperatures. This option is the code prior to version 2009-06-26_28.0, and causes errors. It is included for those who want to re-verify previous results...USE ONLY IF YOU MUST, AND EXPECT ERRORS!

Results vary depending on which is used. When using 'dynamical', pressure differences between LS and SS large enough to have significant effects on matchanges have been observed. Consider these results from model CollExTM_H1_p1, time step 2, lag number 509905:
dynamical
LS: 0.1204496E+10
SS: 0.1191629E+10
plitho
LS: 0.1022932E+10
SS: 0.1022928E+10
traditional
LS: 0.9760302E+09
SS: 0.9760303E+09

To conclude, based on this one example,
  1. dynamical pressures agree to 2 significant digits
  2. plitho pressures agree to 5 significant digits
  3. traditional pressures agree to 6 significant digits

Since: version 2009-07-26_28.0
Default, version 2009-07-26_28.0: 'dynamical'
Default, version 2011-12-01_5.0 and later: 'plitho'

For each phase-change rule i, we have the following lines:

color1(i) color2(i) n_vertices(i) set_strain_to_zero(i) [ mech_flag(i) ther_mat1(i) ther_mat2(i) ther_flag(i) ]

If the rule is triggered, any lagrangian particles whose material is mat1(i) will be transformed to material mat2(i). The polygon describes the P-T conditions for triggering the rule.

color1(i)
current mechanical color number
color2(i)
new mechanical color number
n_vertices(i)
number of points in phase-change polygon
set_strain_to_zero(i)
0 => no change to accumulated strain during the matchange
1 => strain reinitialized to zero during the matchange
set_strain_to_zero is used to reinitialize the second invariant of the total strain for a Lagrangian particle that undergoes a particular matchange that transforms color1(i) to color2(i).
NOTE: set_strain_to_zero must be used with care. Each matchange operation works independently of the others. This means, for example that if a particular color2 must always have zero strain accumulated strain, ALL matchanges that create this color must have have set_strain_to_zero =1. Failure to do this will produce results that will set the strain to zero for some matchanges and not for others. This will mean that this color will have a mixture of Lagrangian particles some with strain reinitialized and some not.
Since version 2009-07-26_25.0

Since version 2011-12-01_9.3 and 2012-12-21_1.0
The switch matchange_code_version affects thermal matchanges.

mech_flag(i)
= 'yes' enables this rule for mechanical colors. This is the default.
= 'no' disables this rule for mechanical colors.
ther_mat1(i)
current thermal material number
ther_mat2(i)
new thermal material number
ther_flag(i)
= 'yes' enables this rule for thermal materials.
= 'no' disables this rule for thermal materials.

For each vertex j, we have:

phase_change_poly_x(i,j) phase_change_poly_y(i,j)

phase_change_poly_x(i,j)
temperature in degrees K
phase_change_poly_y(i,j)
pressure in Pa

color_inside_polygon reversible [ ther_mat_inside_polygon ]

The final line of input for each matchange rule

color_inside_polygon
the mechanical color number corresponding to the inside of the phase-change polygon
reversible
0 => phase change is reversible (i.e. color2 can change back to color1)
Reversability disallowed Since version 2008-10-08
1 => phase change is not reversible (i.e. once color1 has changed into color2, it can't change back)

Since version 2012-12-21_1.0

ther_mat_inside_polygon
the thermal material number corresponding to the inside of the phase-change polygon

NOTE: The polygon vertices must be specified in counter-clockwise order, and the P-T axes must have the origin at the lower left corner:

Example:

MATCHANGE
1    ! number of phase-change rules
! 1st mech color, 2nd mech color, # points in polygon, reset_strain
           5               7                5              0 
  823.  1.2d9    ! temperature (degrees K), pressure (in Pa)
 1473.  2.5d9    ! temperature (degrees K), pressure (in Pa)
 1473   3.3d10   ! temperature (degrees K), pressure (in Pa)
  623.  3.3d10   ! temperature (degrees K), pressure (in Pa)
  623   3.3d9    ! temperature (degrees K), pressure (in Pa)
7  1    ! mechanical color inside polygon, 1=reversible,0=non-reversible

In this case, the mech color inside the phase-change polygon is 7.


MATCHANGE_DELV section

The cell volume changes that result from color (material) changes can be applied either as forces when solving for velocities in the mechanical FE calculation, or added as vertical velocities after the mechanical FE calculation.

MATCHANGE_DELV

Indicates start of the MATCHANGE_DELV section.

value

Indicates where the volume changes should be applied.
force : apply the cell volume change as a force acting in the FE mechanical calculation. This is previous behaviour.
velocity : apply the cell volume change as vertical velocity after the FE mechanical calculation. ell volume changes accumulate from the bottom upward.

Default: force

Since: June 2015


PROGRADE section

This section starts with the word PROGRADE (on a line by itself).

PROGRADE

progradation_model

progradation_model
= 0 no progradation
= 3 Lykke Gemmer's fill-and-spill
= 4 half-Gaussian progradation + uniform aggradation
= 5 hybrid progradation, fill-and-spill sedimentation
= 7 half-Gaussian progradation + local aggradation
= 9 pure aggradation with user-defined geometry
= 10 notes (pdf)
= 11 notes (pdf)
= 12 notes (pdf)
= 13 notes (pdf)
= 14 level-slope...slope...slope-level function

progradation model 3

no documentation available

progradation models 4 and 7

numprograde

numprograde
= number of times at which the progradation and aggradation parameters are specified.

For each time, we need:

tprograde(i), h0(i), hinf(i), x0(i), prograde_length(i), vprograde(i), haggrade(i), vaggrade(i)

tprograde(i)
= model time (in seconds) at which the profile values are set to these values
h0(i)
= surface height at x.le.x0(i) (metres)
hinf(i)
= limit of surface height at x.eq.infinity (metres)
x0(i)
= initial x position where the half-Gaussian part of the profile starts (metres)
prograde_length(i)
= Gaussian width of the profile as defined above (metres)
vprograde(i)
= progradation velocity (m/s) velocity at which progradation function moves in positive x direction (0.317097d-9 = 0.01m/year). Negative values give retrogradation of the profile but do not result in erosion of sediments above the profile.
haggrade(i)
= base level for aggradation (metres)
vaggrade(i)
= aggradation velocity (m/s) which is uniform for all x velocities at which sediments aggrade (0.317097d-9 = 0.01m/year). Negative values give degradation of the profile but do not result in erosion of sediments above the profile.

NOTE:

  • The meaning of haggrade and vaggrade depends on the progradation model. Models 4 and 7 use these number differently. See below.
  • All heights and positions are with respect to the origin (x=0, y=0)
  • All parameter values are the values for the time interval i and are reset at the beginning of time interval i+1.
  • The form of the profile function is quite flexible but it is necessary to the user to calculate the initial values for each time interval carefully if steps/jumps in the positions of the profile are to be avoided. For example, x0(i+1) = x0(i) + vprograde * (t(i+1)-t(i)) correctly positions the profile at the beginning of interval(i+1) to the position it had at the end of interval(i). Similarly, setting haggrade(i+1) to haggradet at the end of the i'th tme interval ensures continuity in the position of haggradet.
  • Parameter values specified for tprograde(numprograde) apply for all times larger than this value.

Notes On Progradation Model 4

This is a half-Gaussian progradation, plus uniform aggradation, that vary with time according to a time- stepping function. In progradation model 4 haggrade and haggradet are increments in metres that are added/subtracted from the half gaussian profile at all positions. The general form of the sediment profile function is

  h(x) = (h0(i)-hinf(i)) * exp(-((x-x0(i))**2)/prograde_length(i)**2))
          + hinf(i) + haggrade(i) 

The values of h0(i), hinf(i), x0(i), prograde_length(i), and haggrade(i) are reset at the beginning of each time interval(i) and do not vary during the time interval.

During the time interval x0(t) and haggrade(t) vary with time since the beginning of the time interval according to

  x0(i, deltat) = x0(t=t(i)) + vprograde(i) * deltat
  haggrade(i,deltat) = haggrade(t=t(i)) + vaggrade(i)*deltat

  where 
    deltat is the time elapsed since the beginning of this time interval
    t(i)   is the absolute model time at the beginning of the time interval

Notes On Progradation Model 7

This is a half-Gaussian progradation, plus local aggradation, that vary with time according to a time- stepping function. In progradation Model 7 haggrade and haggradet are absolute values in metres of the progradation/aggradation profile such that the profile is set to haggradet whereever it is less than haggradet.

haggrade is the initial height of haggradet for each time interval and haggradet is raised or lowered at vaggrade with time. This means the effect of aggradation only occurs where haggradet is greater than the progradation profile.

For example, the seaward end of the progradation profile can be independently determined by aggradation if hinf is set to a value that is less than haggrade. Variations in haggradet will then specify the profile at the seaward end and not hinf provided haggradet remains greater than hinf.

The general form of the sediment profile function is

   h(x) = (h0(i)-hinf(i)) * exp(-((x-x0(i))**2)/prograde_length(i)**2))
          + hinf(i) 

   if h(x) is greater than haggradet. 
      h(x) is set to haggradet where h(x) is less than haggradet.

The values of h0(i), hinf(i), x0(i), prograde_length(i), and haggrade(i) are reset at the beginning of each time interval(i) and do not vary during the time interval. During the time interval x0(t) and haggrade(t) vary with time since the beginning of the time interval according to

  x0(i, deltat) = x0(t=t(i)) + vprograde(i) * deltat
  haggrade(i,deltat) = haggrade(t=t(i)) + vaggrade(i)*deltat 

  where
    deltat is the time elapsed since the beginning of this time interval
    t(i) is the absolute model time at the beginning of the time interval

progradation model 5

no documentation available

progradation model 9

no documentation available

progradation model 10, 11, 12, 13

notes (pdf)

progradation model 14

numprograde

numprograde
= number of intervals at which the progradation and aggradation parameters are specified.

For each prograde interval:

tprograde(i), vprograde(i), haggrade(i), vaggrade(i)

tprograde(i)
= model time (in seconds) at which the profile values are set to these values
vprograde(i)
= progradation velocity (m/s) velocity at which progradation function moves in positive x direction.
haggrade(i)
= base level for aggradation (metres)
vaggrade(i)
= aggradation velocity (m/s) vertical velocity of haggrade

x(1) x(2) ...

x(1) x(2) ...
points on the x axis. maximum 20 points.

level(1) level(2) ...

level(1) level(2) ...
level at the 'x' points.

PTT_TRACK section

This section defines the Lagrangian particles to be tracked during the model run. When a Lagrangian particle is "tracked", its pressure and temperature are periodically saved to a file, for use in PTt (Pressure-Temperature-time) plots.

PTT_TRACK

Indicates start of section.

ptt_track_flag

ptt_track_flag
0 => no particle tracking
1 => track Lagrangian particles

num_tracked_ptcls save_ptt_interval dens_ptt

num_tracked_ptcls
the number of Lagrangian particles you want to track.
save_ptt_interval
the timestep interval to use for saving the PT data.
(e.g. save_ptt_interval=100 means save the tracked data every 100 timesteps)
dens_ptt
density value used to calculate static pressure.
Used for consistency with Microfem output.

For i = 1 to num_tracked_ptcls, there should be a line specifying each Lagrangian particle to be tracked:

ptcl_col(i) ptcl_row(i)

ptcl_col(i)
the column number (in the Lagrangian grid) of the Lagrangian particle to be tracked
ptcl_row(i)
the row number (in the Lagrangian grid) of the Lagrangian particle to be tracked

When tracking is on, Sopale will create files called ptt_p_c#_r# (where c#_r# is the row number, column number), one file for each tracked particle. Within the ptt_p_c#_r# file, each line has the following values:

Column Value
1timestep
2time in millions of years
3temperature in Celsius
4 dynamic pressure (epress*) divided by 10**8, => it is in kbars
5x position of particle (in km)
6y position of particles (in km)
7depth of particle (in km)
8"static pressure" (in kbar), calculated from dens_ptt, assuming the overburden is uniformly of this density.
(this is a very approximate pressure estimate).
9dynamic pressure (epress*) at the centre of the Eulerian element where this particle is found (in Pa)
10lithostatic pressure (plitho*) at the centre of the Eulerian element where this particle is found (in Pa)
11thermal material
12density of the thermal material
13**maximum temperature (Celcius)
14**maximum dynamic pressure
15**maximum static pressure

* epress and plitho are the names of the variables within Sopale; epress corresponds to dynamic pressure; plitho corresponds to lithostatic pressure. They will be good estimates of the dynamic pressure and weight of the overburden (including topography) experienced by the Lagrangian particle if the vertical spacing of the Eulerian grid is small.

** Since version 2009-07-26_22.7


MAXPTOUTPUT section

MAXPTOUTPUT

Records maximum pressure and temperature events for lagrangian particles.

ts_start ts_mod ts_end

The internal data will be updated at time steps which are multiples of ts_mod. The updating will start at ts_start and terminate at ts_end.

(output interval specifications)

There are two ways of specifying when the output files will be written

ts1 ts2 ts3 ...

list of time steps at which the output file will be written. File names will be:
<nameout1>maxptoutput_ts<timestep>.dat
where timestep is 6 digits.

SameAsNsave2

output files will be written at the same time as type 2 files. File names will be:
<nameout1>maxptoutput_f<nn>_o
where nn is the same as the frame number in type 2 files.
This is recommended for files that will, or may, be used in a restart.

xul yul xbl ybl

xul yul
the x,y coordinates of the upper left corner (meters)
xbl ybl
the x,y coordinates of the lower left corner (meters)

xbr ybr xur yur

xbr ybr
the x,y coordinates of the lower right corner (meters)
xur yur
the x,y coordinates of the upper right corner (meters)

screen_factor

The dynamical pressure is known to have spikes that coincide with material changes. These spikes contaminate the record of maximum dynamical pressure. screen_factor can be used to avoid that contamination.
If (dynamical_pressure > screen_factor x lithostatic_pressure) then
the value of dynamical_pressure will be ignored.

All lagrangian particles inside the box defined by the 4 points will be tracked. As particles enter the box, they will be added to the tracked list. As particles leave the box, they will be removed from the tracked list.

For each lagrangian particle within the box defined by the user, these data items will be recorded:

Data Item Accompanying Data
number of the lagrangian particle
model time of output file
x, y coordinates at 'model time'
maximum lithostatic pressure time, temperature
maximum dynamical pressure time, temperature
maximum temperature time, lithostatic pressure, dynamical pressure

For example:

#  lag number   model time    lag x pos    lag y pos   plitho_max         time  temperature   pdynam_max
       843590   0.1250E-01   0.9000E+03   0.6204E+03   0.1095E+01   0.1250E-01   0.5442E+03   0.0000E+00
       843591   0.1250E-01   0.9000E+03   0.6198E+03   0.1117E+01   0.1250E-01   0.5511E+03   0.0000E+00

         time  temperature   temper_max         time       plitho       pdynam
   0.0000E+00   0.0000E+00   0.5442E+03   0.1250E-01   0.1095E+01   0.1435E+01
   0.0000E+00   0.0000E+00   0.5511E+03   0.1250E-01   0.1117E+01   0.1453E+01

Units
pressure GPa (sopale_nested uses Pa)
temperature C (sopale_nested uses K)
time Myr (sopale_nested uses s)
x km (sopale_nested uses m)
y km (sopale_nested uses m)

Since version 2011-12-01_8.0


SEDIMENT section

This section starts with the word SEDIMENT (on a line by itself).

SEDIMENT

nsedcolor

nsedcolor
= number of time/sediment material data lines

For each time/sediment material data line:

tchron_strat(i) sedcolor(i) [ sedcolort(i) ]

tchron_strat(i)
= time (in seconds)
sedcolor(i)
= the mechanical material number for the sediment that will be deposited when the model time is greater than or equal to tchron_strat(i)

Since sopale_nested version 2009-07-26_18,
sedcolort(i) is an optional parameter.

sedcolort(i)
= the thermal material number for the sediment that will be deposited when the model time is greater than or equal to tchron_strat(i)
Unless added_sed_colort is negative, it overrides this.

WATERLOAD section

This section starts with the word WATERLOAD.


waterloading  = 1 use nodal weights for water loading
              = 0  don't use water loading

refh_water, rho_water, stab_waterweight
    refh_water  = sea level (for code-sep17-03, assumed to be parallel 
		  with the x axis)
    rho_water   = density of water
    stab_waterweight = stabilization factor for the water weight
                       (versions code-dec04-03 and later)

NOTE: code-sep17-03 version has modified water loading to act as
      a pressure normal to the (land) surface.  This only works
      if the model is constructed in the x-z coordinate system
      such that the surface of the water (refh_water) is horizontal
      (i.e. parallel to the x coordinate axis) and that gravity acts
      in the minus z direction. 
      [ Will only work if ang_gravity is -90.]


OUTPUTS section

Since sopale_nested version 2011-09-06_1.0

OUTPUTS

x_begin x_end flag frequency

Defines a range in x within which lagrangian densities are printed to a file. The filename will be

<nameout1><process_number>.lag_density

where 'process_number' is 0 for LS, and a positive integer for any nest.

Each line of the file, for each time step, for each lagrangian particle, will contain:

  1. lag number
  2. x coordinate
  3. y coordinate
  4. depth (below surface)
  5. mechanical color
  6. density
x_begin
begin of range
x_end
end of range
flag
= 'all': the density of all lags will be printed.
= anything-else: the density of lags whose material is compactible will be printed.
frequency
n => output files will be written at time steps which are multiples of this number. This input flag was not present in version 2011-09-06_1.0

[surface_strings]

Since sopale_nested version 2014-04-24_1.0

Discontinued sopale_nested version 2016-01-28

Surface process calculations - erosion, sedimentation, progradation, surface advection - can be performed on the top of the eulerian grid. This is the default.

If this section is present, those calculations will be performed on two strings of lagrangian-like particles.

Two strings are maintained; one for the surface, and one for the sediment base. The surface string is the surface of the eulerian grid. The sediment base string is the initial eulerian grid top, modified during time stepping by eulerian velocities - but not affected by surface processes.

These particle strings are initially formed from a copy of the top of the lagrangian grid. The initial spatial density is then increased by a user-specified factor. During time-stepping, the particles move with a velocity in both x and y interpolated from neighbouring points on the eulerian grid. A minimum spacing of particles is maintained by injecting new particles if required.

[surface_strings]

Indicates start of section.

sep_min sep_max

sep_min
a number; fraction of the cell width. If the distance between any adjacent particles in a surface string becomes less than this, one particle will be deleted.
sep_max
a number; fraction of the cell width. If the distance between any adjacent particles in a surface string becomes larger than this, one particle will be injected.
example:
0.1 0.3

[end]

Indicates end of section.

[surface_strings_v2]

Since sopale_nested version 2016-01-28

Surface process calculations - erosion, sedimentation, progradation, surface advection - can be performed on the top of the eulerian grid. This is the default.

If this section is present, those calculations will be performed on two strings of lagrangian-like particles.

Two strings are maintained; one for the surface, and one for the sediment base. The surface string is the surface of the eulerian grid. The sediment base string is the initial eulerian grid top, modified during time stepping by eulerian velocities - but not affected by surface processes.

These particle strings are initially formed from a copy of the top of the lagrangian grid. The initial spatial density is then increased by a user-specified factor. During time-stepping, the particles move with a velocity in both x and y interpolated from neighbouring points on the eulerian grid. A minimum spacing of particles is maintained by injecting new particles if required.

[surface_strings_v2]

Indicates start of section.

separation

spacing
spacing, in meters, of points on the surface string. If the distance between adjacent particles in a surface string becomes less than 0.5 of this number, particles will be deleted. If the distance between adjacent particles in a surface string becomes greater than 1.5 of this number, particles will be inserted.
example:
1000.0

[end]

Indicates end of section.


This page was last modified on Wednesday, 20-Apr-2016 11:02:42 ADT
Comments to geodynam at dal dot ca