[CIG-SHORT] long-term simulation of subduction seismic cycles with rate and state friction law

[Brad Aagaard]

On 9/9/13 9:56 AM, Shaoyang Li wrote:
> Dear Brad,
>
> I am trying to apply RS friction law on subduction fault based on the 2d
> subduction mesh in Pylith/examples. As a test, I apply the
> velocity-weakening properties on seimogenic part of slab top interface
> and rest of the interface a creeping rate. I have problems:
>
> 1, the two faults on subduction interface have no common node. If the
> two faults beside each other, pylith can not read them and creat
> cohesive nodes. If I leave one node between the two faults, pylith can
> read and calculate. Why this can happen?

There are two issues here. (1) Creating faults that "meet" to form a 
single fault is not currently feasible in PyLith. In rare circumstances 
you can get lucky and the cohesive cells would be created in just the 
right way that it would work, but this is not guaranteed. See the 
diagram of how we form cohesive in our JGR paper 
(http://onlinelibrary.wiley.com/doi/10.1002/jgrb.50217/abstract). We 
continue to think about better ways to form cohesive cells to make this 
robust and eliminate nonuniqueness. (2) This means it is not feasible to 
drive a shallow fault with friction by deep constant prescribed creep. 
The solution is to use a single fault with rate-state friction and 
stable and unstable regions driven by boundary deformation. This should 
be possible but is more difficult to setup.

> 2, Going on with one node left between two faults, I apply creeping
> velocity on deeper fault (the rest nodes of coseimic nodes on slab top
> interface) and slab bottom interface, and Columb friction on coseimic
> fault on slab top. It works. When I change the friction law into RS
> friction. Pylith calculates a while and gives errors:
>
> WARNING! Fault opening with nonzero traction., v_fault: 38, opening:
> 1.11594e-11, normal traction: -1.00349e-06
> mpinemesis: faults/FaultCohesiveDyn.cc:375: virtual void
> pylith::faults::FaultCohesiveDyn::integrateResidual(const
> pylith::topology::Field<pylith::topology::Mesh>&, PylithScalar,
> pylith::topology::SolutionFields*): Assertion `fabs(tractionNormal) <
> _zeroTolerance' failed.
> [0]0:Return code = 0, signaled with Aborted
>
> I use the mesh from the given example, it should have good quality. So I
> am very confused why this errors. I try to use gravity to constrain it
> and no help.

There are two things to notice here. (1) The normal traction is 
compressive but very small. This can be alleviated by using 
gravitational body forces, but you may also want to add a little 
additional compressive traction right at the surface if you want it to 
remain closed. Alternatively, you can allow opening at the surface. (2) 
The fault opening is very small. You need to increase the 
friction.zero_tolerance value so that such small values force the fault 
to remain closed. This is due to the fact that quasi-static solutions 
use iterative solvers. I recommend reviewing Session V from the June 
workshop 
(http://www.geodynamics.org/cig/community/workinggroups/short/workshops/cdm2013/agenda), 
noting the discussion of tolerances.

Regards,
Brad