GEOS icon indicating copy to clipboard operation
GEOS copied to clipboard

Viscoelastic wave propagator

Open sframba opened this issue 10 months ago • 1 comments

This PR introduces attenuation for the isotropic Elastic Wave Propagator (SEM, 2nd order formulation). Specifically, this PR follows the approach of [1], which is based on the satndard Fichtner approach [2] for frequency-indpendent quality factors.

This formulation is based on the standard-linear-solid (SLS) description. Compared to the article, the method has been adapted to the SEM formulation with a Leap-Frog 2nd order timestepping scheme. Since the stress field $\sigma$ is not explicitly stored, this allows reducing the memory footprint of the auxiliary variables from $6 L$ to $3L + 3$, where $L$ si the number of auxiliary SLSs.

With the present formulation, the standard elastic wave equation is augmented as follows:

$\partial_t^2 u - \nabla\cdot C\varepsilon(u) + \sum\limits_{\nu=1}^L\psi^a_\nu = f$, $\partial_t\psi^a_\nu + \omega_\nu\psi^a_\nu = \omega_\nu y_\nu\nabla\cdot C^a\varepsilon(u)$,

where $\varepsilon$ is the strain tensor, $C^a$ is the isotropic attenuative elasticity tensor (see [1]), and $(\omega_\nu, y_\nu)_{\nu=1}^L$ are a set of $L$ reference angular frequencies and adimensional anelasticity coefficients, respectively.

The attenuative matrix $C^a$ is computed from the attenuative Lamé coefficients $\lambda_a, \mu_a$, obtained from the original $\lambda, \mu$ Lamé fields using $\lambda_a+2\mu_a=(\lambda+2\mu)/Q_p$, $\mu_a=\mu/Q_s$ (see [1]). $Q_p$ and $Q_s$ are the pressure and shear elastic quality factor coefficients, which are cell-based fields required as an input and can be either defined on the mesh or given as constant fields via a FieldSpecification. A higher value of $Q$ indicates less attenuation.

Note that if $\sum\limits_{\nu=1}^L y_\nu > \mathrm{\min}(Q)$, some artifacts may appear in the solution, usually in the form of a long-term stable zero-velocity wave focused on the source point(s). If this happens, a warning is given by GEOS.

These coefficients $\omega_\nu$ and $y\nu$ are needed to describe the behavior of $Q$ as a function of frequency, and can be computed as described in [2] given a minimal Q requirement as well as a frequency range. For example, for a frequency range $[9-110Hz]$ and $Q_{\mathrm{min}}=30$, a relatively flat $Q$ (independent of frequency) can be obtained using $\omega_1\approx 69 s^{-1}$, $\omega_2\approx 592 s^{-1}$, $y_1\approx 1.64$ and $y_2\approx 1.75$.

The SLS attenuation mechanism is activated by adding the option attenuationType="sls" to the ElasticWaveSEM solver. SLS parameters can be entered via the slsReferenceAngularFrequencies and slsAnelasticityCoefficients fields. If no SLS parameters are specified, the following default values are used: $L=1$ $\omega_1=2\pi f$, $y_1=\frac{2 Q_0}{Q_0-1}$,

where $f$ is the peak frequency of the source and $Q_0$ is the minimal value of $Q_p$ and $Q_s$ on the entire domain.

note This mechanism can be easily added to the upcoming VTI elastic wave propagator. The kernel introduced here is modular and only the non-attenuative stiffness vector needs to be changed (by calling the appropriate underlying kernel) in order to apply attenuation to the VTI case. The isotropic attenuation mechanism does work even with an anisotropic elasticity kernel, and in fact has some advantages as discussed in [1].

[1] P.-T. Trinh et al, "Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible Cartesian-based mesh", Geophysics, Vol. 84, No. 1.

[2] A. Fichtner, M. van Driel, "Models and Fréchet kernels for frequency-(in)dependent Q", Geophysical Journal International, Volume 198, Issue 3, September, 2014, Pages 1878–1889.

sframba avatar Apr 12 '24 17:04 sframba

Codecov Report

Attention: Patch coverage is 97.04142% with 5 lines in your changes missing coverage. Please review.

Project coverage is 53.77%. Comparing base (7a18ecf) to head (82001dc). Report is 102 commits behind head on develop.

Files with missing lines Patch % Lines
...sSolvers/wavePropagation/shared/WaveSolverBase.cpp 81.81% 4 Missing :warning:
...econdOrderEqn/isotropic/ElasticWaveEquationSEM.cpp 98.27% 1 Missing :warning:
Additional details and impacted files
@@             Coverage Diff             @@
##           develop    #3080      +/-   ##
===========================================
+ Coverage    53.69%   53.77%   +0.07%     
===========================================
  Files         1010     1011       +1     
  Lines        85663    85816     +153     
===========================================
+ Hits         45998    46146     +148     
- Misses       39665    39670       +5     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

codecov[bot] avatar Apr 26 '24 14:04 codecov[bot]