Very long period wave simulation

[lucela]

Hi,

I am currently using SPECFEM Globe to produce low frequency synthetics. It seems like the resolution period cannot go above 68s. Is it possible to simulate very long period seismograms, for example 100s or 200s period ?

Thanks a lot !

[lsawade]

Hi @lucela,

I regularly compute synthetics in the range of 20, 400 seconds and have yet to hit the limit. specfem3D_globe should be able to handle seismograms to 400-500s IIRC (I forget which specific specfem paper mentions where exactly the cowling approximation breaks down).

Could you elaborate on your setup? That is, how have you set these settings?

# forward or adjoint simulation
SIMULATION_TYPE                 = ???
NOISE_TOMOGRAPHY                = ??? 
SAVE_FORWARD                    = .false.

# number of chunks (1,2,3 or 6)
NCHUNKS                         = ???

# angular width of the first chunk (not used if full sphere with six chunks)
ANGULAR_WIDTH_XI_IN_DEGREES     = ??.d0  
ANGULAR_WIDTH_ETA_IN_DEGREES    = ??.d0
CENTER_LATITUDE_IN_DEGREES      = ??.d0
CENTER_LONGITUDE_IN_DEGREES     = ??.d0
GAMMA_ROTATION_AZIMUTH          = ??.d0

# number of elements at the surface along the two sides of the first chunk
# (must be multiple of 16 and 8 * multiple of NPROC below)
NEX_XI                          = ???
NEX_ETA                         = ???

# number of MPI processors along the two sides of the first chunk
NPROC_XI                        = ???
NPROC_ETA                       = ???


model                           = ???

[lucela]

Hi @lsawade ,

Thanks for your answer, that's good news !
Yes I couldn't find what was the maximum period of my simulations, I just found that the minimum period resolved cannot be bigger than 68s for my settings.

`# forward or adjoint simulation
SIMULATION_TYPE = 1 # set to 1 for forward simulations, 2 for adjoint simulations for sources, and 3 for kernel simulations
NOISE_TOMOGRAPHY = 0 # flag of noise tomography, three steps (1,2,3). If earthquake simulation, set it to 0.
SAVE_FORWARD = .false. # save last frame of forward simulation or not

`# number of chunks (1,2,3 or 6)
NCHUNKS = 6

`# angular width of the first chunk (not used if full sphere with six chunks)
ANGULAR_WIDTH_XI_IN_DEGREES = 90.d0 # angular size of a chunk
ANGULAR_WIDTH_ETA_IN_DEGREES = 90.d0
CENTER_LATITUDE_IN_DEGREES = -20.2d0
CENTER_LONGITUDE_IN_DEGREES = 179.4d0
GAMMA_ROTATION_AZIMUTH = 0.d0

# number of elements at the surface along the two sides of the first chunk # (must be multiple of 16 and 8 * multiple of NPROC below)
NEX_XI = 64
NEX_ETA = 64

`# number of MPI processors along the two sides of the first chunk
NPROC_XI = 2
NPROC_ETA = 2

MODEL =transversely_isotropic_prem_plus_3D_crust_2.0

[lucela]

I also did it with one chunk, and with NPROC_XI = 8 and NPROC_ETA = 8

[lsawade]

Ok, your settings don't look crazy!

First, a couple of comments

• If you do a 1 chunk simulation you will get spurious arrivals from the boundaries even if you use stacey boundary conditions. Stacey isn't perfect so be aware. Those data are probably not great to use for periods this long since you would want to compare them to long-period, observed data, which are free from any boundary effects.
• Using this type of model transversely_isotropic_prem_plus_3D_crust_2.0, then neither is perfect. The process here stretches and deforms the crust, and wherever crust2.0's moho is shallower than PREM's, PREM's mantle values are somewhat "extended". So again beware here. Note that this is probably not a problem for you.

Finally, a couple of questions:

• when you say "resolution period cannot go above 68s" what do you mean by that? How does this manifest? Could it be a problem with the way you process the seismograms?
• Are you setting PARTIAL_PHYS_DISPERSION_ONLY = .false.?
• Could you send a plot of the raw seismogram you are producing?

[lucela]

Thanks for all of your comments, it's very useful.
My main goal is to use specfem synthetics in a W phase inversion algorithm to see if the W phase is in the synthetics. The W phase is very long period and low amplitude, which arrives between the P and S first arrival.

• Earth model : I am wondering which velocity model I should use and if for such long periods, there will be a noticeable difference, on the waveforms and mainly for the W phase. I am thinking of trying several of the others global models available for specfem globe.

• "resolution period cannot go above 68s", I mean for one chunk, the minimum number of processors per side of chunk is 8, so for a chunk of 90degrees wide (the maximum size it can be), the minimum period resolved is 68s, according to the equation given in the manual.

• I have PARTIAL_PHYS_DISPERSION_ONLY = .true. Should I set it to false ?

• Here are some seismograms for the simulation I sent you the inputs. WEL is 2000km away from the source, RAO 1300km, MSVF 900km. The event is a magnitude 6.9.

8bd1f4a75.png)

Thanks a lot

[lucela]

Sorry, for the second bullet point I meant number of elements by side of chunk, not processors.

[lsawade]

Thanks for all of your comments, it's very useful. My main goal is to use specfem synthetics in a W phase inversion algorithm to see if the W phase is in the synthetics. The W phase is very long period and low amplitude, which arrives between the P and S first arrival.

I don't know whether anyone is looking into this, could be very cool! I forgot about this paper: Morales-Yáñez et al. 2020

Earth model: [...]

A good starting point, even with a 3D solver is PREM. Because you have simple reference points. Especially with the literature on the W-phase. Once, that is figured out you can move to 3D models.

"resolution period cannot go above 68s", [minimum period resolved is 68s]

Indeed the minimum period, i.e., max frequency is 1/68 Hz. The equation in the manual is an estimate, so take it with a grain of salt. The true resolution may lie a little over or a little under the estimated value. To be on the safe side I tend to filter a little more conservatively. Say, keep periods between [70s,300s]. Now, there is a problem here. If you want to see 300 s periods in your spectrogram clearly you need a long seismogram (or you pad it with zeros). The seismograms that you are showing are ~30min, into which 300s fits a few times, but a few times only. So, keep that in mind. For seismograms with long periods like this we usually do global simulations.

I have PARTIAL_PHYS_DISPERSION_ONLY = .true. Should I set it to false?

Yes, right now you are not using "full physical dispersion". It's computationally or storage-wise slightly more expensive (not sure which one), but the seismograms are more accurate.

Here are some seismograms for the simulation I sent you the inputs. WEL is 2000km away from the source, RAO 1300km, MSVF 900km. The event is a magnitude 6.9.

These two seismograms seem to be unfiltered. 68s warrants only one wavelength per 68 seconds. There may be some numerical noise related to the length of the half duration. Try to filter your streams between 68s and 250s.

stream.filter('bandpass', freqmin=1.0/250.0, freqmax=1.0/68.0)