search

ydluo / qdyn

A Quasi-DYNamic earthquake simulator
56 stars 28 forks source link

Reproducing the results from Kaneko and Ampuero (GRL, 2011) #1

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago 
Hello,

I am trying to reproduce the simulations documented in the paper by Kaneko and 
Ampuero (GRL, 2011) titled: "A mechanism for preseismic steady rupture fronts 
observed in laboratory experiments".  

I believe that I have become relatively acquainted with the newest version of 
QDYN.  I have setup the problem using the parameters provided in the paper as 
follows:

p = qdyn('set');

%fromulation
p.MESHDIM=2;%dimesion (3D)
p.THETA_LAW=2;%slip law = 2,aging law = 1, aging (no heal) = 0
p.RNS_LAW=1;%1=rate-and-state with cutoff
p.V1=100; %cutoff velocities
p.V2=1e-7;%cutoff velocities

%geometry/nodes
p.L=150e-2;
p.W=50e-2;
p.NX=512*4;
p.NW=1;
p.Z_CORNER=-150e-2;
p.N=p.NX*p.NW;
p.DW(1:p.NW)=p.W/p.NW;
p.DIP_W(1:p.NW)=0;

%Initial condition
p.SIGMA=Normal_Stress; %input via function in Matlab
p.V_SS=10e-6;
p.TH_SS=p.DC/p.V_SS;
p.IC=ceil(p.N/2);
p.V_0 =1.01*p.V_SS ;

%BC's
P.FINITE=1;%0=periodic, 1=surrounded by steady slip

%Constitutive parameters
    %I have just discretized the a and b constitutive parameters here so
    %that the left and right hand side of the fault is VS (37.5 cm, 
    %a-b=0.01) and the middle section is VW (75 cm, a-b=-0.0008)

edge_effect=512; 
b_1=0.00*ones(1,edge_effect);
b_2=0.0108.*ones(1,p.N-2*edge_effect);
p.B=[b_1,b_2,b_1];
p.A=ones(1,p.N).*0.01;
p.DC=0.1e-6;
p.MU=32e9;

%time stepping
twm=1.5*0.297E-09;
year = 3600*24*365;
p.TMAX=twm*year;
p.ACC = 1e-14;

%output
p.OX_SEQ=0; %type of fortran file output (0 = only fort.19)
p.NTOUT=1;%temporal interval for snapshot outputs 
p.NXOUT=1;%spatial interval
p.NSTOP=1; %0=stop @TMAX, 1=end of localization , 2=first slip rate peak

My question is:

1. In the paper the boundary conditions are specified as follows:
 "slip rate delta_dot_Load = 5 cm/year is prescribed on the two 75 cm long outer portions of the fault".  How is this specified in the code?  Is this periodic boundary conditions?

I am performing similar laboratory experiments to Nieslen et al. (2010) and 
observe slow propagating fronts prior to nucleation.  I am working with 
Professor Glaser at the University of California in Berkeley. Please feel free 
to email me at pa.selvadurai@berkeley.edu also.

Regards,

Paul Selvadurai

Original issue reported on code.google.com by pa.selva...@gmail.com on 9 Apr 2013 at 6:06

GoogleCodeExporter commented 9 years ago 
Hi Paul,

Kaneko and Ampuero (2011) used Nadia Lapusta's code, which assumes periodic 
boundary conditions and loading by steady slip on a fault section. QDYN also 
has an option for periodic boundaries but it is instead loaded by remote 
displacements on a boundary parallel to the fault, at a distance W. The 
boundary conditions of Kaneko and Ampuero (2011) cannot be exactly reproduced 
with QDYN. However, by setting MESHDIM=1 and FINITE=1 you can apply a steady 
slip velocity over the semi-infinite region of the fault outside the 
rate-and-state region. That might be sufficiently similar.

Also note that Lapusta's code includes complete elastodynamics, but QDYN is 
only quasi-dynamic. You should expect some differences in the results due to 
that too.

I also found the following errors in your input parameters:  
+ the simulations are on a 1D fault embedded in a 2D medium, hence you must set 
MESHDIM=1.
+ velocity cut-off was not used, you must set RNS_LAW=0 (sorry this parameter 
is not documented in the manual).

Let me know if that helps.

Pablo Ampuero

Original comment by ampu...@caltech.edu on 11 Apr 2013 at 5:24

GoogleCodeExporter commented 9 years ago 
Just too clarify I have attached a figure with my understanding of the code 
(QDYN_example_2D.jpg).  If so, is it possible to formulate the problem as 
depicted in the separately attached figure (Possible_formulation.jpg).

Original comment by pa.selva...@gmail.com on 11 Apr 2013 at 10:56

Attachments:

GoogleCodeExporter commented 9 years ago 
Sorry I just noticed that this figure was hard to read due to missing dashed 
lines.

Original comment by pa.selva...@gmail.com on 12 Apr 2013 at 1:15

Attachments:

GoogleCodeExporter commented 9 years ago 
I am afraid there is some confusion there. The boundary condition FINITE=1 is 
only implemented for 2D media, but your left figures are 3D. MESHDIM=1 means 
that the fault is a one-dimensional line contained in a two-dimensional elastic 
medium.
With MESHDIM=1, the parameter NW is ignored, only NX is considered. Slip 
parallel to X (mode II) or normal to X (mode III) follow similar equations, if 
you set MU=(shear modulus)/(1 - Poisson's ratio) in mode II.

You can think of MESHDIM=1 as the special case of a 3D problem where slip is 
invariant along one of the axis, as if fault cells were infinitely elongated in 
one direction.

Original comment by ampu...@caltech.edu on 12 Apr 2013 at 6:00

GoogleCodeExporter commented 9 years ago 
Sorry for the confusion.  It is completely clear now.  Thank you kindly.  I 
removed the figures to avoid further confusion.

Original comment by pa.selva...@gmail.com on 12 Apr 2013 at 6:27

GoogleCodeExporter commented 9 years ago

Original comment by ampu...@caltech.edu on 2 Nov 2014 at 11:00