significant difference between prem_a_5s and prem_a_2s

[Thomas-Ulrich]

Hi,

I've been using Instaseis to generate teleseismic from a dynamic rupture model and compare them with observations. I recently realized that using a database including anisotropy significantly improved the fit to the surface waves. For that I used prem_a_2s on http://ds.iris.edu/ds/products/syngine/.
Then, I've downloaded the prem_a_5s on a local server for quicker retrieval of teleseismic data. But great was my surprise when I compare the synthetics obtained with 2s:

teleseismic_comparison_Sulawesi_whole_signala2s.pdf

With those obtained with the 5s database:

teleseismic_comparison_Sulawesi_whole_signala5s.pdf

Note that I use a 66-450 s band-pass filter for comparing the data, so I would expect similar results.
[myst.filter('bandpass', freqmin=0.00222222, freqmax=0.015, corners=4, zerophase=True)]

Any idea about what could explain the differences? Which database should I trust (I guess the 2s)? I could share my scripts for plotting these files if needed. Note that I m using a finite source (40 points sources).

Thanks in advance,

Thomas Ulrich.

Just for reference, here are the database that I m using:

ReciprocalMergedInstaseisDB reciprocal Green's function Database (v10) generated with these parameters:
    components           : vertical and horizontal
    velocity model       : prem_ani
    attenuation          : True
    dominant period      : 4.950 s
    dump type            : displ_only
    excitation type      : dipole
    time step            : 1.214 s
    sampling rate        : 0.824 Hz
    number of samples    : 3214
    seismogram length    : 3900.0 s
    source time function : gauss_0
    source shift         : 8.497 s
    spatial order        : 4
    min/max radius       : 5620.0 - 6371.0 km
    Planet radius        : 6371.0 km
    min/max distance     : 0.0 - 180.0 deg
    time stepping scheme : symplec4
    compiler/user        : ifort         1600 by ri32xag5 on login22
    directory/url        : ../../../../../import/freenas-m-05-seissol/ulrich/prem_a_5s_merge_compress2
    size of netCDF files : 123.1 GB
    generated by AxiSEM version GIT_VERSION at 2017-03-27T19:03:08.000000Z

Syngine service version:  1.0.4
    components           : vertical and horizontal
    velocity model       : prem_ani
    attenuation          : True
    dominant period      : 2.100 s
    dump type            : displ_only
    excitation type      : dipole
    time step            : 0.512 s
    sampling rate        : 1.952 Hz
    number of samples    : 7048
    seismogram length    : 3609.9 s
    source time function : gauss_0
    source shift         : 3.586 s
    spatial order        : 4
    min/max radius       : 5621.0 - 6371.0 km
    Planet radius        : 6371.0 km
    min/max distance     : 0.0 - 180.0 deg
    time stepping scheme : symplec4
    compiler/user        : ifort         1500 by ri32xag4 on login21
    directory/url        : http://service.iris.edu/irisws/syngine/1
    size of netCDF files : 1.4 TB
    generated by AxiSEM version SVN_VERSION at 2015-12-10T18:42:23.000000Z

[rmodrak]

Hi Thomas,

I am not an AxiSEM developer, and I'm not sure how relevant the following is.

I tried to reproduce your comparison, except rather than downloading synthetics directly from syngine, I downloaded Green's functions from syngine and used these to generate the synthetics myself. This is slightly different from your procedure probably, but it is the only way I am readily set up to do it.

For me, everything looks fine. I am not seeing the same discrepancy. I am attaching my script and output figure for reference.

-Ryan

EDIT: clarified comments in code
compare_syngine_models.txt
compare_syngine_models.pdf

[rmodrak]

Some context--I was coincidentally doing an extremely similar experiment, so I went ahead and posted above. At best, I hope it might narrow the scope of debugging. At worst, it might not be relevant.

[Thomas-Ulrich]

Hi,
Thx for your answer and your script (as it depends on mtuq, which itself depends on python 2.7, and as my install is based on python 3.6 I cannot run the script straightforwardly, but in any case the plot look as I would have expected).

I will try to simplify my script and share it, so that any problem could be more easily spotted.

Thomas.

[Thomas-Ulrich]

Hi,
Here is a script to reproduce the discrepancies I m experiencing between instaseis database:
generateTeleseismicAndCompareSulawesi_simple.txt
and here is the ouput:
comparison_prem_ani_obs.pdf
Most probable is that I m doing something wrong, but I cannot spot the problem.

Thomas.

[martinvandriel]

Hi Thomas,
the stf used in the computation of the two databases is different, so you need to make sure to use the same stf for comparison. From your script I understand you try to do that, but I guess you may need to set reconvolve_stf=True in the calls to get_seismogram() for the stf to be actually used.

Cheers,
Martin

[Thomas-Ulrich]

Hi,
you are right, I was wrongly using get_seismogram().
But that does unfortunately not solve the problem as I m using a finite source in my more complex setup.
Here is the corrected script featuring the same discrepancies with a finite source (of 1 point source):
generateTeleseismicAndCompareSulawesi_simple_FiniteSource.txt
And the output:
comparison_prem_ani_obs.pdf
Thomas.

[Thomas-Ulrich]

Hi, Just following up on this problem.
I did some comparison with prem_a_2s, prem_a_5s and prem_a_10s for a different earthquake model (Sumatra, 2004), and got matching signals at Periods 100-450 for prem_a_2s and prem_a_10s but different signals for prem_a_5s.
So I think prem_a_2s have to be trusted over prem_a_5s.
(or probably there are different ingredients in prem_a_5s compared to the other 2 data bases).

[martinvandriel]

Hi,
I think you are right, there is some problem with the 5s model. Here is what I got.

import matplotlib.pyplot as plt
import numpy as np


import instaseis
dbs = []

dbs.append(instaseis.open_db("syngine://prem_a_10s"))
dbs.append(instaseis.open_db("syngine://prem_a_5s"))
dbs.append(instaseis.open_db("syngine://prem_a_2s"))



receiver = instaseis.Receiver(latitude=42.6390, longitude=74.4940)
source = instaseis.Source(
latitude=89.91, longitude=0.0, depth_in_m=12000,
    m_rr = 4.710000e+24 / 1E7,
    m_tt = 3.810000e+22 / 1E7,
    m_pp =-4.740000e+24 / 1E7,
    m_rt = 3.990000e+23 / 1E7,
    m_rp =-8.050000e+23 / 1E7,
    m_tp =-1.230000e+24 / 1E7)


for db in dbs:
    print(db.info.period)
    source.set_sliprate_lp(db.info.dt, db.info.npts, 0.01)
    tr = db.get_seismograms(source=source, receiver=receiver, components='Z',
                            dt=0.1, reconvolve_stf=True,
                            remove_source_shift=False)[0]
    print(tr)
    tr.filter('lowpass', freq=0.01, corners=8)
    tr.filter('highpass', freq=0.002, corners=8)
    plt.plot(tr.times(), tr.data, label=f'{db.info.period}s')

plt.legend()

plt.show()

@sstaehler I assume you computed these, do you have any idea what could be the reason? Was there any change in the end with the Q models between the 2 runs?

[sstaehler]

Long time ago...
(and I do not have access to the machine anymore)
The only commit with changes in background_models.f90 that is close to this in time is
geodynamics/axisem@474758a, which added a Q=80 layer in the lithosphere.
Is it possible that one of the models was calculated using the an older version? I'm traveling right now, but it should be possible to do a quick 10s run with the latest version and compare.

[krischer]

@sstaehler Did you ever to that run?