Attenuation with and without read_velocities_at_f0 flag

[kenenig]

Description

I wanted to run identical simulations showing the difference in seismograms for both elastic and anelastic cases. I chose to use the read_velocities_at_f0 flag since the velocity profile I wanted to input had much slower relaxed velocities than specfem was calculating when I input infinite frequency velocities. The reference frequency flag ended up creating a more realistic amplitude reduction but the anelastic waveform arrived earlier than the identical elastic run. The only thing changing between the elastic and anelastic run is setting viscoelastic attenuation to be false or true. Does the flag revert to using non-causal attenuation anywhere in the workflow?

For context, I am inputting velocities from a mid-ocean ridge (MOR) geodynamic model. The relaxation spectrum used to calculate the relaxed velocities are from Yamauchi and Takei 2016. In specfem, I am using a line of point sources to simulate a plane wave from a teleseism. There are 101 receivers across the 300 km in the center above the MOR model. There is another 200 km on either side to avoid the absorbing boundary response in the receivers.

I attached a zip file with all my inputs as well as figures showing the normalized waveforms from select stations.

Specfem2D_inputs.zip

Affected SPECFEM2D version

SPECFEM2D v8.1.0

Your software and hardware environment

Macbook Pro (Monterey, Version 12.5.1), M1 Processor

Reproduction steps

1. Download zip file. 
2. Open MOR1_ref_f0_off or MOR1_ref_f0_on (depending on if you want to see results with or without reference frequency flag)
2. Run MOR_driving_script.sh

Screenshots

No response

Logs

No response

OS

Mac

[danielpeter]

looking at your zip-file, both Par_files in your example specify:

ATTENUATION_f0_REFERENCE        = 1  
READ_VELOCITIES_AT_f0           = .false.

thus, the READ_VELOCITIES_AT_f0 and ATTENUATION_f0_REFERENCE won't be considered and don't affect your simulation. the main difference in your Par_files are the Qkappa and Qmu values of the provided model materials.

in your case, the MOR1_ref_f0_off/ example specifies Qkappa == Qmu == 9999.d0. this is basically no attenuation effect, even if you turn on ATTENUATION_VISCOELASTIC = .true.. I'm not sure though why in your figures you see an amplitude change as these Q values shouldn't produce noticeable attenuation.

as explained in the documentation, turning on attenuation will assume that the velocities are given at unrelaxed state, unless the READ_VELOCITIES_AT_f0 is on: https://github.com/SPECFEM/specfem2d/wiki/03_mesh_generation

the code however also has a slightly different way then how it shifts velocities to the frequency range of the simulation. using READ_VELOCITIES_AT_f0 should specify an f0 frequency that is close the simulation range, otherwise the computation of attenuation with standard linear solids won't be a flat frequency spectrum anymore (this depends also on the number of SLS specified in the Par_file).

[kenenig]

Sorry I should have specified. In the shell script MOR_driving_code.sh, it run two simulations: one where attenuation is turned off, then another where attenuation is turned on and there is a sempar command that changes ATTENUATION_VISCOELASTIC accordingly. For the shell script in the MOR_ref_f0_on folder, there is an additional sempar command that turns on the read_velocities_at_f0 command. I also made sure the frequency I calculated the velocities at was the same as the source frequency (1Hz).

Also, I used a lot of materials to define the geodynamic model. The materials and regions toward the base of my input model has much higher quality factor values than the rest of the model. I attached one of the output snapshots from when I ran MOR_ref_f0_on with attenuation turned on and off.

Having the frequency spectrum to be not flat would be ideal actually. My intention was to see the difference in waveforms using different anelastic scaling relationships and those frequency spectrums tend to have a slope. Is there a way to do that?

Anelastic: Elastic:

[danielpeter]

regarding non-flat attenuation spectrums, that's not implemented - and I would suggest against shifting and using fewer SLS to create a drop-off of the spectrum as that will lead to rather unpredictable spectrum behaviors.

what is the console output for your runs above?

[kenenig]

Okay sounds good. Here are the output logs for the elastic and anelastic case with the reference frequency flag on.

output_elastic.log.zip output_anelastic.log.zip

[danielpeter]

thanks for the output files - that's quite a bit of output for that many material definitions, sources and time step outputs :)

first, I notice that you write out seismograms at every time step. that will slow down the simulation quite a bit. consider outputting the seismograms only every 1000 (or higher) steps if needed - or as usually done only at the very end of the simulation, by setting in your Par_file, say

..
# interval in time steps for writing of seismograms
# every how many time steps we save the seismograms
# (costly, do not use a very small value; if you use a very large value that is larger than the total number
#  of time steps of the run, the seismograms will automatically be saved once at the end of the run anyway)
NTSTEP_BETWEEN_OUTPUT_SEISMOS   = 8000
..

the other work-around would be to specify a tomography file and use the tomo model instead of defining all these materials per elements, but that requires a bit of preparation and is more for having a convenient input and won't affect the simulation duration much.

regarding the attenuation part, you do get the correct output for your settings - to explain:

• all sources use a Ricker wavelet with a 1 Hz dominant frequency. Thus, the spectrum of each source is dominated by a 1 Hz signal.
• using the Ricker wavelet and the READ_VELOCITIES_AT_f0 flag, the code assumes that your material definitions are all set for your source reference frequency. (note that the ATTENUATION_f0_REFERENCE parameter in the Par_file has only an effect for Dirac or Heaviside sources)

attenuation has two main effects, one being the shift of the apparent signal velocity, i.e., a phase shift, and the other the decrease of the signal energy. the shift of signal velocity depends on its frequency difference to the reference frequency f0.

when using the READ_VELOCITIES_AT_f0 flag, the model velocities are defined at 1 Hz, i.e., there is no phase shift on the dominant signal due to attenuation - no matter how strong it is. for the source frequencies away from the 1 Hz, there would be a shift, but given your STF is dominated by the 1 Hz signal, you hardly see a "pulse broadening". the main effect you get and see is the amplitude decrease. if your source would excite a broader range of dominant frequencies, then the frequencies further away from 1 Hz would experience a stronger phase shift, i.e., you see the "pulse broadening" and wave group separation.

and in case the READ_VELOCITIES_AT_f0 flag is .false., the model velocities are assumed to be given in the unrelaxed state, i.e., at infinite frequency. the attenuation phase shift in this case has the effect of slowing down your signals.

[kenenig]

Yeah, I was having issues with the external tomo file at first, so I ended up using the the default meshing internal mesher.

All this sounds good, I'll try calculating the velocities at a different frequency than the source frequency next. Thank you for your help !

[kenenig]

Unfortunately, I tried using velocities calculated at 0.083 Hz instead while keeping the source frequency at 1 Hz but there was still no phase delay. It looks identical to when the frequencies were the same. Do you know why it could be doing this?

[danielpeter]

updating the velocity model will just speed up or slow down the whole signal according to your changes. the READ_VELOCITIES_AT_f0 will always assume your velocity model is given at the source frequency in your setting (using a Ricker wavelet). it is thus unaware if your velocity model is supposed to be at 1 Hz or at another frequency, it just takes the source frequency as reference frequency.

to highlight the differences between elastic and viscoelastic seismograms, you would have two options: A. you shift your velocities to the source frequency by hand and run the simulation with the READ_VELOCITIES_AT_f0 flag turned on and the compare this result against the elastic simulation with the velocity model at infinite frequency.

B. you use your velocity model for the unrelaxed state (infinite frequency) and don't use the READ_VELOCITIES_AT_f0 flag to let attenuation figure out the phase shift required for your source signals, and compare the results with the ATTENUATION_VISCOELASTIC flag turned on and off.

in principle, the viscoelastic seismograms for A. and B. should be identical, unless the attenuation model isn't flat and covering the whole simulated frequency range. in that case, you might need to increase the number of standard linear solids (N_SLS parameter).