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 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)