Why the implementation of the solid-fluid coupling in the specfem2d and specfem3d_globe is different?

[lyuchao]

In the specfem2D/compute_coupling_viscoelastic_ac.f90, we can find that the definition of pressure is minus potential_dot_dot_acoustic. This can be understood due to the definition relation between pressure (P) and the potential of displacement(phi): P = - /partial_{tt} potential. We have benchmarked the specfem2d with our current DFDM 2D global, the waveforms are the same.

! compute pressure on the fluid/solid edge
  pressure = - potential_dot_dot_acoustic(iglob)

However, in the specfem3D/compute_coupling.f90, the pressure = - RHO_TOP_OC * accel_outer_core.
My understanding is that: accel_outer_core = /partial_{tt} phi (second derivatives with respect to time). This implementation is different from the specfem2d case. And it is also different from the previous velocity potential formulation (-rho* partial_t velocity potential).

pressure = - RHO_TOP_OC * accel_outer_core(iglob_oc)

So I feel confused.

After checking the source codes, it seems that RHO_TOP_OC is a ratio between RHO_TOP_OC and average density. That makes sense.

[danielpeter]

right, these different potential formulations were confusing me as well :)

there are many ways how to define the potentials, the most obvious choices are to either have a displacement potential as in

$${\bf u} = \nabla \chi \tag{1}$$

with pressure

$$p = - \kappa \nabla \cdot {\bf u} = - \rho \, \partial_{tt} \chi$$

or

$${\bf u} = \frac{1}{\rho} \nabla \chi \tag{2}$$

with pressure

$$p = - \kappa \nabla \cdot {\bf u} = - \partial_{tt} \chi$$

here, ${\bf u}$ is displacement.

so, for the globe version SPECFEM3D_GLOBE, the formulation follows this statement in the code file specfem3D.f90:

! Formulation in the fluid (acoustic) outer core:
! -----------------------------------------------
!
! In case of an acoustic medium, a displacement potential Chi is used
! as in Chaljub and Valette, Geophysical Journal International, vol. 158,
! p. 131-141 (2004) and *NOT* a velocity potential as in Komatitsch and Tromp,
! Geophysical Journal International, vol. 150, p. 303-318 (2002).
! This permits acoustic-elastic coupling based on a non-iterative time scheme.
! Displacement if we ignore gravity is then: u = grad(Chi)
! (In the context of the Cowling approximation displacement is
! u = grad(rho * Chi) / rho, *not* u = grad(Chi).)
! Velocity is then: v = grad(Chi_dot)       (Chi_dot being the time derivative of Chi)
! and pressure is: p = - rho * Chi_dot_dot  (Chi_dot_dot being the time second derivative of Chi).
! The source in an acoustic element is a pressure source.
! The potential in the outer core is called displ_outer_core for simplicity.
! Its first time derivative is called veloc_outer_core.
! Its second time derivative is called accel_outer_core.

that is the globe version follows formulation (1) as this turns out to be easier when gravity is included in the equations of motions. and as you mentioned, RHO_TOP_OC is density at the top of the outer core, and normalized (by an average density RHOAV) since in globe, all units and dimensions get normalized.

for the 2D version, the definition follows this statement in specfem2D.F90:

! In case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
!
! This permits acoustic-elastic coupling based on a non-iterative time scheme.
! Displacement is then:
!     u = grad(Chi) / rho
! Velocity is then:
!     v = grad(Chi_dot) / rho
! (Chi_dot being the time derivative of Chi)
! and pressure is:
!     p = - Chi_dot_dot
! (Chi_dot_dot being the time second derivative of Chi).
!
! The source in an acoustic element is a pressure source.
!
! First-order acoustic-acoustic discontinuities are also handled automatically
! because pressure is continuous at such an interface, therefore Chi_dot_dot
! is continuous, therefore Chi is also continuous, which is consistent with
! the spectral-element basis functions and with the assembling process.
! This is the reason why a simple displacement potential u = grad(Chi) would
! not work because it would be discontinuous at such an interface and would
! therefore not be consistent with the basis functions.

that is the 2D version (similar as the 3D_Cartesian version) uses formulation (2) for the potential definition, mostly for the reasons of pressure being continuous across discontinuities and therefore the potential is continuous, which lends itself well to the smooth GLL basis functions used in the spectral-element method together with the definition of the wavefields being on global nodes (instead of local nodes as in discontinuous Galerkin methods).

one last comment: the 2D version has a similar idea in mind like often done in exploration seismics, where the subsurface is approximated by an acoustic media. In such cases, we assume that the acoustic medium is allowed to have first-order discontinuities, i.e., velocity jumps from one subsurface layer to another one. using formulation (2) can easily account for this and still have pressure and potential being continuous. on the other hand, the acoustic domain in the globe version is implemented for the outer core. there we assume that the outer core is fluid with no internal discontinuities. in that case, formulation (1) is still fine to have a continuous potential as well, and thus the spectral-element discetization remains happy.

hope that helps out a bit to get this confusing thing sorted out... :)

[lyuchao]

Dear Daniel,
Thanks so much for your constructive summaries of the different displacement potentials used in specfem2d and specfem3d_globe.

I have checked that for models with a smooth density like the outer core of the Earth, both these two displacement potentials $\mathbf{u}=\frac{1}{\rho}\nabla(\varphi)$ and $\mathbf{u}=\frac{1}{\rho}\nabla(\rho\phi)$ ($\varphi$ and $\phi$ are two different displacement potentials) have exactly the same output waveforms, which are also the same as the reference waveforms calculated by the handy Distributional Finite Different Method.

Kind Regards
Chao Lyu