[CIG-SHORT] significant discrepancy between coseismic displacement field of homogeneous FEM model and Okada analytical solution

[Brad Aagaard]

Yu Zhou,

The main issue in comparing PyLith prescribed slip results with Okada is 
that Okada uses uniform slip patches while PyLith uses a continuous slip 
field that is C0 continuous. At the edges of a large uniform slip patch, 
there is a stress concentration, which would require very fine 
discretization with PyLith. If you want to compare PyLith against an 
Okada solution, the preferred approach is to construct a slip 
distribution that tapers linearly to zero at the buried edges of the 
fault. This is more difficult to represent in Okada. We have used many 
overlapping uniform slip patches of varying sizes stacked together to 
form a slip distribution with very fine stair stepping to approximate 
the linear variation represented exactly in PyLith. Our results show 
nice convergence using the error metric Matt mentioned.

Regards,
Brad


On 07/18/2013 03:52 AM, ZHOU Yu wrote:
> Dear Pylith developers:
>
> I build up a homogeneous FEM model to compare coseismic ground displacement field output by Pylith with Okada analytical solution,
> but I find unneglectable discrepancy between the two.  I've examined possible reasons and tried hard, but I cannot eliminate it through all my efforts.
> I would be very grateful if you can help me with this discouraging problem.
>
> My FEM model settings are as the following:
> Model dimension: 1800kmx2000kmx600km
> Fault geometry: strike=0, dip=90, 40km long, 12km wide
> Tetrahedron element size: 1km on or near the fault, gradually larger to 100km on the boundaries.
> Coseismic fault slip: 10m strike slip, homogeneous on the fault
> Material type: homogeneous Poisson solid
> Boundary conditions: fixed along the normal direction for 4 lateral surfaces and bottom surface, free on ground surface
>
>
> In the attached tarball please find mesh_shortfaullt.jou, which is the journal file for mesh generating in CUBIT. Its output data file mesh_shortfault.exo is not included due
> to its big size(>50Mbyte). All other files that Pylith needs can also be found in the tarball, please just run:
>   pylith eqsim.cfg
>
> Output hdf5 files will be stored in ./output. Then Please run the matlab script homoslip_diff.m.
> This Matlab script reads displacement at every node on the ground surface in FEM mesh, then interpolates disp field at a group of gridded points.
> Theoretical disp field is also calculated with Okada solution at those gridded points. Differences between Okada disp and FEM disp are shown as a vector map
>   and 3 grid images. Discrepancy for disp-Y component, which is along strike direction, is larger than 6% in most area, increasing with distance from the fault.
>
> At the beginning, the interval size was set to 2km on fault interface, so I thought perhaps the mesh was not dense enough near the fault. After I changed it to 1km,
> the result didn't improve as much.
>
> Then I thought about interpolation method. Since variables in a tetrahedron are linear, it should be right to employ bilinear interpolation.
>
> Later, I thought that strain component exy is very large due to pure strike slip in the neighbourhood of the fault, I set interval size to 0.5km along x direction
>   near the fault(see line 41 in mesh_shortfaullt.jou). The number of nodes jumped to bigger than 500,000, but such an excessively dense mesh didn't help either.
>
> I've also considered that slip at a single node in FEM means linear slip distribution tapering into zero at the neighbouring nodes. However, if the mesh is sufficiently
> dense, and slip distribution is uniform on the entire fault surface, this effect will be neglectable.
>
> After all these computation and consideration, now I wonder whether this discrepancy is kind of intrinsic error of Pylith software and cannot be eliminated by users.
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>