man(1) relax man page man(1)
NAME
relax - Evaluates the deformation due to fault slip, surface loading or
inflation and the time series of postseismic relaxation that follows
due to fault creep or viscoelastic flow.
SYNOPSIS
relax [-h] [--dry-run] [--help] [--no-grd-output] [--no-proj-output]
[--no-relax-output] [--no-stress-output] [--no-txt-output] [--no-vtk-
output] [--no-xyz-output]
DESCRIPTION
relax computes nonlinear time-dependent viscoelastic deformation with
powerlaw rheology and rate-strengthening friction in a half space due
to coseismic stress changes, initial stress, surface loads, and/or mov-
ing faults.
OPTIONS
-h print a short message and abort calculation
--dry-run
write lightweigh geometry files and abort calculation
--help print a short message and abort calculation
--no-grd-output
cancel output in GMT grd binary format
--no-proj-output
cancel output in geographic projection
--no-relax-output
cancel output of the postseismic contribution
--no-stress-output
cancel output of stress tensor in any format
--no-txt-output
cancel output in text format
--no-vtk-output
cancel output in Paraview VTK format
--no-xyz-output
cancel output in GMT xyz format
--with-stress-output
export stress tensor
--with-vtk-output
export output in Paraview VTK format
--with-vtk-relax-output
export relaxation to VTK format
ENVIRONMENT
The environment variable OMP_NUM_THREADS controls the number of threads
used by OpenMP at execution. For example, the calls
OMP_NUM_THREADS=4 ./relax
and
export OMP_NUM_THREADS=4
relax
produce the same output.
INPUT PARAMETERS
grid size (sx1,sx2,sx3)
Defines the number of samples in the three directions. sx1, sx2
and sx3 must be powers of two for optimal efficiency of the Fast
Fourier Transform. sx3 must be even for internal memory manage-
ment considerations.
sx2
+------------+
1 / /|
x / / | x1 (north)
s / / | /
/ / | /
+------------+ + +-------> x2 (east)
s | | / |
x | | / |
3 | | / x3 (depth)
| |/
+------------+
sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
dx1,dx2,dx3 control the sampling size of the model and the size
of the computational domain.
beta is a roll-off parameter (0 <= beta <= 0.5) that controls
the tapering of the slip distribution on faults. beta=0.2-0.3
are recommended values.
nyquist controls how the slip distribution is evaluated: is a
fault patch dimension is smaller than nyquist*dx, then it is
evaluated using the Okada's equations; otherwise, fault slip is
represented by equivalent body forces, which leads to a much
faster calculation.
The extent of the computational grid is exported to
wdir/cgrid.vtp for visualization in Paraview.
origin position & rotation (xo,yo,rot)
Indicates the coordinates (xo,yo) and orientation of the compu-
tational coordinate system. Use 0 0 0.
observation depths for displacements and stress
Indicates depths at which map views of displacement and stress
are exported in GMT.
output directory (wdir)
All output files are written to the specified directory.
elastic parameters and the buoyancy wavelength (lambda, mu, gamma)
The uniform Lame parameter (lambda), shear modulus (mu) and the
buoyancy wavelength (gamma=(1-nu)*rho*g/mu), where nu is Pois-
son's ration, rho is the density of the crust, g is the gravity
acceleration. For the Earth, typical values are lambda=mu=30
GPa, and gamma = 8.33e-7 /m = 8.33e-4 /km
integration time and output time parameters (T,dt,a)
Integration time (T) refers to the duration of the calculation
in physical units. dt is the time step of output files (dt<T).
Negative integer values for dt indicates output tied to inter-
nally optimized time steps: -1 corresponds to output every time
step, -2 to output every other time steps. Scaling parameter a
modifies the internally evaluated time step. 0 < a < 1 improves
the accuracy of the time evolution. a > 1 reduces the accuracy
of the explicit time integration but speeds up the calculation.
number of observation planes (nop)
Observation planes are planar surface of arbitrary orientation
where displacement and stress are exported in ASCII and GMT .grd
format for visualization. Integer nop indicates the number of
such planes. If nop > 0 then relax asks for the geometry of each
planes, with one line per plane, as follows:
# nb x1 x2 x3 length width strike dip
where nb is an index running from 1 to nop, x1, x2 and x3 are
the reference coordinates, length and width, and strike and dip
are the dimension and the orientation of the observation plane.
These parameters are defined in Section FAULT GEOMETRY.
number of observation points (np)
Observation points are locations where displacement and stress
are exported as time series in ASCII. Integer no indicates the
number of such points. If np > 0 then relax asks for the name
and location of each point, with one line per point, as follows:
# nb NAME x1 x2 x3
where nb is an index running from 1 to np, NAME is a four-char-
acter name used to identify the output file, x1, x2 and x3 are
the point coordinates. Time series of displacement and stress at
these points are written to file NAME.txt, where NAME is the
user-provided name.
number of stress observation segments (nsp)
Stress observation segments are fault patches where stress
(shear, normal, dip-shear, strike-shear, Coulomb stress) evalu-
ated and exported in GMT and VTK formats. This is how Coulomb
and other time-dependent stress calculations are carried out in
relax. Integer nsp indicates the number of such patches. If nsp
> 0 then relax asks for the definition of each fault patch, with
one line per patch, as follows:
# nb x1 x2 x3 length width strike dip friction
where nb is an index running from 1 to nsp, x1, x2, x3, length,
width, strike and dip are the position, dimension and orienta-
tion of the fault patches and friction is the friction coeffi-
cient (usually chosen at 0.6) used to compute Coulomb stress.
The geometry parameters are defined in section FAULT GEOMETRY.
All receiver faults for Coulomb stress calculations are exported
in wdir/rfaults-dsigma-0000.vtp for visualization in Paraview.
number of pre-stress interface (npsi)
Pre-stress interfaces specify at what depth and how pre stress
changes. If npsi > 0, then relax requires the depths and stress
values at each interface, one line per interface, as follows:
# nb depth sigma11 sigma12 sigma13 sigma22 sigma23 sigma33
where nb is an index running from 1 to npsi, depth is the depth
where pre-stress changes, and sigma11, 12, 13, 22, 23, and 33
and the components of the symmatric stress tensor.
number of linear viscous interfaces (nlvi)
Viscous interfaces specify at what depth and how the viscosity
changes in the Earth, and define the background 1-D viscosity
model that can be subsequently modified using ductile zones. If
nlvi > 0, then relax requires the depths and viscosity and cohe-
sion values at each interface, one line per interface, as fol-
lows:
# nb depth gammadot0 cohesion
where nb is an index running from 1 to nlvi, depth is the depth
where cohesion and gammadot0 change, gammadot0 is the fluidity
(defined as gammadot0 = mu / eta, where eta is the viscosity),
the reciprocal of the Maxwell relaxation time, and cohesion is
the minimum value of stress to drive viscoelastic flow. The def-
inition of the 1-D model is explained in Section DEPTH-DEPENDENT
STRUCTURE.
All viscous interface are exported to wdir/linearlayer-nb.vtp ,
where nb is the interface index, for visualization in Paraview.
The definition of the 1-D depth-dependent model is followed by:
number of linear ductile zones (nldz)
Ductile zones are volumes where the background viscosity is
ammended. If nldz > 0, then relax requires the list of ductile
zones, defined as
# nb dgammadot0 x1 x2 x3 length width thickness strike dip
where nb is an index running from 1 to nldz, dgammadot0 is the
modifier to the background fluidity, x1, x2, x3, length, width,
thickness, strike and dip are the position, dimension and orien-
tation of the rectangular volume. The fluidity used to drive
viscoelastic flow is gammadot0+dgammadot0. If gammadot0+dgam-
madot0<=0, no flow occurs. Therefore, setting large negative
values of dgammadot0 makes the region elastic. The geometric
parameters are defined in Section LATERAL VARIATIONS OF VISCOUS
PROPERTIES.
All ductile zones are exported to wdir/weakzones-linear.vtp for
visualization in Paraview, including when computation is aborted
with the --dry-run option.
number of nonlinear viscous interfaces (nnlvi)
Nonlinear viscous interfaces specify at what depth and how the
power-law rheology parameters change in the Earth, and define
the background 1-D viscosity model that can be subsequently mod-
ified using ductile zones. Viscoelastic relaxation in relax can
have ontributions from both linear and nonlinear rheologies. If
nnlvi > 0, then relax requires the depths, viscosity, power and
cohesion at each interface, one line per interface, as follows:
# nb depth gammadot0 power cohesion
where nb is an index running from 1 to nnlvi, depth is the depth
where cohesion and gammadot0 change, gammadot0 is the reference
fluidity, power is the power-law rheology power exponent (strain
rate = gammadot0 ( tau / mu ) ^ power, where tau is the coseis-
mic stress change plus the prestress), and cohesion is the mini-
mum value of stress to drive viscoelastic flow.
The definition of the 1-D depth-dependent power-law model is
followed by:
number of nonlinear ductile zones (nnldz)
Nonlinear ductile zones are volumes where the background nonlin-
ear viscosity is ammended. If nnldz > 0, then relax requires the
list of nonlinear ductile zones, defined as
# nb dgammadot0 x1 x2 x3 length width thickness strike dip
where nb is an index running from 1 to nnldz, dgammadot0 is the
modifier to the background fluidity, x1, x2, x3, length, width,
thickness, strike and dip are the position, dimension and orien-
tation of the rectangular volume. The power exponent of the duc-
tile zone is the same as in the background model.
All ductile zones are exported to wdir/weakzones-nonlinear.vtp
for visualization in Paraview, including when computation is
aborted with the --dry-run option.
number of friction interfaces (nfi)
Friction interfaces define the variations of fault friction
properties with depth, using the framework of rate-strengthening
friction. If nfi < 0, relax requires the depth, reference veloc-
ity, strengthening parameter and cohesion at each depth, one
line per interface, as follows:
# nb depth gamma0 (a-b)sigma friction cohesion
where nb is an index running from 1 to nfi, depth is the depth
where friction properties change, (a-b)sigma is the reference
stress (typically of the order of 1 MPa), friction is the fric-
tion coefficient (usually 0.6) and cohesion is the stress
enveloppe. If nfi > 0 the list of interface is followed by a
definition of faults where stress-driven slip occurs:
number of afterslip planes (nap)
Afterslip planes are rectangular surfaces where stress-driven
slip occurs. If nap > 0, relax requires the list of afterslip
planes, as follows:
# nb x1 x2 x3 length width strike dip rake
where nb is a index running from 1 to nap, x1, x2, x3, length,
width, strike and dip are the position, dimension and orienta-
tion of the fault plane and rake is a +-90 constrain on the rake
of afterslip. If |rake| > 360, the constraint is ignored. Some
of these parameters are defined in Section FAULT GEOMETRY.
All afterslip planes are exported in wdir/aplane-nb.vtp , where
nb in the patch index, for visualization in Paraview.
number of interseismic loading shear faults (nisf)
Interseismic shear faults are faults that move at a user-defined
constant rate. If nisf > 0, relax requires the list of faults.
number of interseismic loading opening dykes (niod)
Interseismic opening dykes are intrusions that open at a user-
defined constant rate. If niod > 0, relax requires the list of
dykes.
number of events (ne)
Events are moments in time when new internal or external forces
act of the system (ne >= 1). If ne = 1, then a list of shear
faults, opening dyke and surface tractions are required and the
change occurs at t = 0. If ne > 1, then a list of shear faults,
opening dyke and surface tractions are required for each event.
The first event occurs at time 0 and each new event is pre-
scribed a time of occurrence. Having multiple events allows the
user to model the effect of a sequence of earthquakes, or to
prescribe time-dependent loads.
number of shear dislocations (strike-slip and dip-slip faults) (nsd)
Shear dislocations are rectangular slip patches. If nsd > 0,
relax expects a list of such slip patches, as follows
# nb slip x1 x2 x3 length width strike dip rake
where nb is an index running from 1 to nsd, x1, x2, x3, length,
width, strike dip are the position, dimension and orientation of
the slip patch; slip and rake are the slip amplitude and rake.
For positive slip, rake = 0 indicates left-lateral slip, and for
positive slip and shallow dip (dip <= 90), rake = 90 indicate
thrust motion. These parameters are defined in Section FAULT
GEOMETRY.
All faults are exported to wdir/rfaults-e.vtp , where e is the
event number, for visualization in Paraview. Export to
wdir/rfaults-e.xy allows visualization with GMT.
number of tensile cracks (nts)
Tensile cracks are dykes with opening or closure of the elastic
walls. If nts > 0, relax expects a list of cracks:
# nb opening x1 x2 x3 length width strike dip
where nb is an index running from 1 to nts, opening is the nor-
mal motion of the walls, and the other parameters define the
position, orientation and dimension of the cracks.
number of surface loads (nsl)
Surface loads are surface tractions in the vertical direction
coming from the loading and unloading of lakes, dams or the
freezing or melting of ice. If nsl > 0, relax expects a list of
surface loads, defined with their geometry and weight, as fol-
lows:
# nb x1 x2 length width t3 T phi
where nb is an index running from 1 to nsl, x1, x2, length and
width define the position and dimension of the load, t3 is in
units of stress (force/surface), positive down, and T can be a
period (T > 0 implies stress=t3*sin(2 pi/T + phi) or not (T <= 0
implies stress = t3 H(t), with H(t) the Heaviside function).
time of next event (te)
If the computation includes several events (ne > 0), the second
and subsequent events are preceded by their time of occurrence.
EXAMPLE INPUTS
The line starting with the '#' symbol are comments.
CALLING SEQUENCE
relax < input.dat
or
relax <<EOF
# this line is a comment
`cat input.dat`
EOF
COSEISMIC DISPLACEMENT
Computes coseismic displacements due to uniform fault slip:
relax --no-proj-output <<EOF
# grid size (sx1,sx2,sx3)
256 256 256
# sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
0.05 0.05 0.05 0.2 0
# origin position & rotation
0 0 0
# observation depths for displacements and stress
0 0.5
# output directory
output_dir
# elastic parameters and gamma = (1-nu) rho g / mu = 8.33e-7 /m = 8.33e-4 /km
30 30 8.33e-4
# integration time (t1)
0 -1 1
# number of observation planes
0
# number of observation points
0
# number of stress observation segments
0
# number of prestress interfaces
0
# number of linear viscous interfaces
0
# number of powerlaw viscous interfaces
0
# number of friction interfaces
0
# number of interseismic loading strike-slip and opening
0
0
# number of coseismic events
1
# number of shear dislocations (strike-slip and dip-slip faults)
1
# index slip x1 x2 x3 length width strike dip rake
1 1 -1 0 0 2 1 0 90 0
# number of tensile cracks
0
# number of dilatation sources (Mogi source)
0
# number of surface loads
0
EOF
POSTSEISMIC VISCOELASTIC DEFORMATION
Computes time-dependent postseismic viscoelastic deformation
driven by stress induced by fault slip:
relax --no-proj-output <<EOF
# grid size (sx1,sx2,sx3)
512 512 512
# sampling size, smoothing & nyquist (dx1,dx2,dx3,beta,nq)
0.5 0.5 0.5 0.2 0
# origin position & rotation
0 0 0
# observation depths for displacements and stress
0 10
# output directory
viscoelastic
# elastic parameters and gamma = (1-nu) rho g / mu = 8.33e-7 /m = 8.33e-4 /km
30 30 8.33e-4
# integration time (t1)
10 -1 0.5
# number of observation planes
0
# number of observation points
0
# number of stress observation segments
0
# number of prestress interfaces
0
# number of linear viscous interfaces
1
# nb depth gammadot0 cohesion
1 20 1 0
# number of linear ductile zones
0
# number of powerlaw viscous interfaces
0
# number of friction interfaces
0
# number of interseismic loading strike-slip and opening
0
0
# number of coseismic events
1
# number of shear dislocations
1
# index slip x1 x2 x3 length width strike dip rake
1 1 -10 0 0 20 10 0 90 0
# number of tensile cracks
0
# number of dilatation sources
0
# number of surface loads
0
EOF
FAULT GEOMETRY
Static dislocation sources are discretized into a series of planar seg-
ments. Slip patches are defined in terms of position, orientation, and
slip, as illustrated in the following figure. For positive slip, a zero
rake corresponds to left-lateral strike-slip motion and a 90 degree
rake corresponds to a thrust motion (when dip is smaller than 90
degrees).
N (x1)
/
/| strike
x1,x2,x3 ->@-------------------------- E (x2)
|\ p . \ w
:-\ i . \ i
| \ l . \ d
:90 \ s . \ t
|-dip\ . \ h
: \. | Rake \
| --------------------------
: l e n g t h
Z (x3)
Slip distributions are defined as a list of slip on individual patches,
for example:
# number of shear dislocations
4
# nb slip x1 x2 x3 length width strike dip rake
1 0.4 0 0 0 1.3 2.3 18 57 0
2 1.1 0 1 0 1.3 2.3 18 57 0
3 2.7 0 0 2 1.3 2.3 18 57 0
4 0.2 0 1 2 1.3 2.3 18 57 0
DEPTH-DEPENDENT STRUCTURE
Depth-dependent variations of properties is obtained from the interpo-
lation of a series of tie points, following the method employed in the
PREM model. For example, the 1-D model below
@------------------------> (modulus)
|.
| .
| .
z1 | + v1
| .
| v3 .
z2,z3 | + - - + v2
| |
| |
| | v4
z4,z5 | + - - - - - + v5
| |
| :
| |
| :
|
Z (x3)
is specified as follows:
# number of interfaces
6
# nb depth value
1 0 0
2 z1 v1
3 z2 v2
4 z3 v3
5 z4 v4
6 z5 v5
and the last value v5 is continued down to the bottom extension of the
computational grid.
LATERAL VARIATIONS OF VISCOUS PROPERTIES
Lateral variations of viscous properties can occur in rectangular vol-
umes of arbitrary orientation and dimension. The geometry of the anoma-
lous ductile zones is defined with the reference position (x1,x2,x3),
length, width, thickness, strike and dip, as illustrated below. The
final value of the fluidity that controls viscoelastic flow is the sum
of the background value defined in the depth-dependent model and the
value in the ductile zones.
N (x1)
/
/| strike
x1,x2,x3 ->@-------------------------- E (x2)
|\ \ w +
:-\ \ i /
| \ \ d / s
:90 \ \ t / s
|-dip\ \ h / e
: \ \ / n
| -------------------------- k
: l e n g t h / c
| / i
: / h
| / t
: /
| +
Z (x3)
The input is defined as follows:
# number of ductile zones
1
# nb dgammadot0 x1 x2 x3 length width thickness strike dip
1 -1 0 0 0 1 1 1 0 90
SEE ALSO
Rousset B., S. Barbot, J.-P. Avouac and Y.-J. Hsu, "Postseismic Defor-
mation Following the 1999 Chi-Chi Earthquake, Taiwan: Implication for
Lower-Crust Rheology", J. Geophys. Res., 2012
Bruhat L., S. Barbot and J.-P. Avouac, "Contributions of Afterslip and
Viscoelastic Flow Following the 2004 Parkfield Earthquake", J. Geophys.
Res., v. 116, B08401, 11 PP., 2011, doi:10.1029/2010JB008073
Barbot S. and Y. Fialko, "A Unified Continuum Representation of Post-
seismic Relaxation Mechanisms: Semi-Analytic Models of Afterslip, Poro-
elastic Rebound and Viscoelastic Flow", Geophys. J. Int., v. 182, 3, p.
1124-1140, 2010, doi:10.1111/j.1365-246X.2010.04678.x
Barbot S. and Y. Fialko, "Fourier-Domain Green Function for an Elastic
Semi-Infinite Solid under Gravity, with Applications to Earthquake and
Volcano Deformation", Geophys. J. Int., v. 182, no. 2, pp. 568-582,
2010, doi:10.1111/j.1365-246X.2010.04655.x
Barbot S., Y. Fialko, Y. Bock, "Postseismic Deformation due to the
Mw6.0 2004 Parkfield Earthquake: Stress-Driven Creep on a Fault with
Spatially Variable Rate-and-State Friction Parameters", J. Geophys.
Res., vol. 114, B07405, 2009, doi:10.1029/2008JB005748
BUGS
No known bugs.
AUTHOR
Sylvain Barbot (sbarbot@ntu.edu.sg)
COPYRIGHT
RELAX is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your
option) any later version.
RELAX is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License along
with RELAX. If not, see <http://www.gnu.org/licenses/>.
1.0.3 02 Nov 2012 man(1)