Static displacements off the charts?

[martijnende]

Hi Andrew,

I just played around with pyprop (v1.1.2) a bit to get some static displacements from a large earthquake (here taking Tohoku as an example), but the values I'm getting are many orders of magnitude higher than I'd expect. A quick example (mostly based on the documentation):

import numpy as np
import pyprop8 as pp
from pyprop8.utils import make_moment_tensor, rtf2xyz

layers = [(3.0,    1.8,   0, 1.02),
          (5.0,    4.5, 2.4, 2.5),
          (np.inf, 8.0, 4.5, 3.0)]

model = pp.LayeredStructureModel(layers)

stations = np.array([[ 5.3,  6.2],
                     [-2.1,  0.3],
                     [ 1.5, -3.1]]) # 3 pairs of (x,y) coordinates
depth = 3 # Seafloor instruments
receivers = pp.ListOfReceivers(stations[:,0], stations[:,1], depth)

event_x, event_y, event_dep = (0.1, 0.1, 30.) # Spatial location

# Tohoku EQ: https://rctwg.humboldt.edu/sites/default/files/duputel_et_al_2011_w-phase.pdf
M0 = 4.26e22
M = make_moment_tensor(strike=196, dip=12., rake=85., M0=M0)
M = rtf2xyz(M)

F = np.array([[ 0],
              [ 0],
              [ 0]]) # Force vector, expressed in Cartesian system

event_time = 0
source = pp.PointSource(event_x, event_y, event_dep, M, F, event_time)

static = pp.compute_static(model, source, receivers)
print(static)

Output:

[[ 1.86455305e+16  4.31039789e+16  1.24503706e+17]
 [-2.18813217e+16  6.16597481e+15  7.05473827e+16]
 [-4.38857015e+15 -1.76272979e+16  1.47671364e+17]]

These static offsets would be enough to get halfway to Alpha Centauri, so I suspect something is wrong (either with my input script or with the pyprop source).

[valentineap]

Hi Martijn,

Apologies, I think this is something that isn’t very well explained in the documentation/the examples are misleading.

There’s a pdf note explaining how the units of the output is related to the units of the various inputs here.

I think at present your output is in units of 10^-15 m. I believe you have:

• Length unit: km (i.e. station locations in km, velocities in km/s; alpha = 10^3 in the notation of that pdf)
• Densities: g/cc (beta = 10^3)
• Moment: Nm (epsilon = 1) Therefore output unit is given by epsilon/(alpha^4 beta) = 10^-15.

If you give all inputs in SI units, the output will be SI. In the examples, the seismic moment is scaled to be in GNm so that the output is in mm.

Hope that resolves the issue – let me know if not, or if anything is unclear. At some point I will find time to try and make things clearer in the docs…

Hope you're well!

[martijnende]

Thanks, that makes sense. But when I put all numbers in SI (essentially multiplying layers and all distances by 1000), I get:

[[ 2.56018345e-14  2.31090371e-14 -1.59758939e-14]
 [-1.84187340e-14 -1.52195682e-14  1.19634581e-15]
 [-6.89331328e-15 -1.03660986e-14  2.95758955e-15]]

Also, I receive the following warning:

pyprop8/_core.py:382: RuntimeWarning: Source-receiver distances exceed 200 km. Flat-earth approximation may not be appropriate. 

Exact code for reference:

import numpy as np
import pyprop8 as pp
from pyprop8.utils import make_moment_tensor, rtf2xyz

layers = np.array([(3.0,    1.8,   0, 1.02),
          (5.0,    4.5, 2.4, 2.5),
          (np.inf, 8.0, 4.5, 3.0)]) * 1e3

model = pp.LayeredStructureModel(layers)

stations = np.array([[ 5.3,  6.2],
                     [-2.1,  0.3],
                     [ 1.5, -3.1]]) * 1e3 # 3 pairs of (x,y) coordinates
depth = 3e3 # Seafloor instruments
receivers = pp.ListOfReceivers(stations[:,0], stations[:,1], depth)

event_x, event_y, event_dep = np.array((10, 10, 30.)) * 1e3 # Spatial location

# Tohoku EQ: https://rctwg.humboldt.edu/sites/default/files/duputel_et_al_2011_w-phase.pdf
M0 = 4.26e22
M = make_moment_tensor(strike=196, dip=12., rake=85., M0=M0)
M = rtf2xyz(M)

F = np.array([[ 0],
              [ 0],
              [ 0]]) # Force vector, expressed in Cartesian system

event_time = 0
source = pp.PointSource(event_x, event_y, event_dep, M, F, event_time)

static = pp.compute_static(model, source, receivers)
print(static)

Things are well over here, I hope you are too :)

[valentineap]

Ah. It looks like the issue is coming from pp.ListOfReceivers(), which does a cartesian-to-polar transformation and implicitly assumes the inputs are in km (because it uses PLANETARY_RADIUS=6371).

Clearly that will need to be fixed, probably by adding a 'units' option to `pp.ListOfReceivers'. However for the time being the easiest route is to work in km and rescale the seismic moment based on the note.

Choosing:

• Distance: km / Velocity: km/s (alpha=10^3)
• Density: g/cc (beta=10^3)

and seeking output in m (zeta = 1) implies that we need to specify moments in PN m, (epsilon = 10^15) i.e. we need to multiply your M0 value by 1e-15.

This seems to work sensibly:

import numpy as np
import pyprop8 as pp
from pyprop8.utils import make_moment_tensor, rtf2xyz

layers = np.array([(3.0,    1.8,   0, 1.02),
          (5.0,    4.5, 2.4, 2.5),
          (np.inf, 8.0, 4.5, 3.0)]) 

model = pp.LayeredStructureModel(layers)

stations = np.array([[ 5.3,  6.2],
                     [-2.1,  0.3],
                     [ 1.5, -3.1]]) # 3 pairs of (x,y) coordinates
depth = 3 # Seafloor instruments
receivers = pp.ListOfReceivers(stations[:,0], stations[:,1], depth)

event_x, event_y, event_dep = np.array((10, 10, 30.)) # Spatial location

# Tohoku EQ: https://rctwg.humboldt.edu/sites/default/files/duputel_et_al_2011_w-phase.pdf
M0 = 4.26e22 * 1e-15 # Rescale moment to give output in metres
M = make_moment_tensor(strike=196, dip=12., rake=85., M0=M0)
M = rtf2xyz(M)

F = np.array([[ 0],
              [ 0],
              [ 0]]) # Force vector, expressed in Cartesian system

event_time = 0
source = pp.PointSource(event_x, event_y, event_dep, M, F, event_time)

static = pp.compute_static(model, source, receivers)
print(static)

gives

array([[-25.19207871,  -4.78483244,  48.41510146],
       [ -1.23205709,   9.67687778, -17.56813853],
       [-21.09470066, -14.06958254,  23.69563566]])

which is the right order of magnitude (10s of metres) of slip for this event. Vertical motions seem a bit large compared to observations, but we are using a point-source representation in the very near field...

[martijnende]

Yes, that works for me. I suppose that the spatial derivatives are then expressed in units of m / km (scaling factor 1000)?

[valentineap]

Yes, I think that's correct for the spatial derivatives. Source parameter derivatives would be in m / PN, and the time derivative in m / s.

[martijnende]

Got it. Thanks!