Gar6more2D
 All Classes Files Functions Variables Pages
Gar6more2D Documentation
Author
Julien Diaz and Abdelaaziz Ezziani

Description

This code compute the quasi-analytical solution of several wave equation in two layered media, using the Cagniard de Hoop method[1] ,[2] ,[9] ,[8] ,[7] ,[5] ,[4] ,[3] . It produces seismograms at given points.

The equations can be written in the general form.

\begin{eqnarray}\newcommand{\setR}{\hbox{I}\!\hbox{R}}\newcommand{\derp}[2]{\frac{\partial #1}{\partial #2}}\newcommand{\dsp}{\displaystyle} A(y)\derp {^2 U}{t^2}-B(y) U=\delta({\bf x}-{\bf x_s})\,f(t), \quad x\in\setR,\,y\in\setR \end{eqnarray}

where \(A\) and \(B\) are operators satisfying

\(A(y)=A^+, \, B(y)=B^+, \quad y>0\),
\(A(y)=A^-, \, B(y)=B^-, \quad y<0.\)

The code analytically compute the Green function \(u\) of the problem

\begin{eqnarray} A(y)\derp {^2 u}{t^2}-B(y) u=\delta({\mathbf x}-\mathbf{x_s})\,\delta(t), \quad x\in\setR, y\in\setR \end{eqnarray}

and convolves it with the source function \(f\). You can modify this function in the subroutine lib/libgeneral/source.F90. The convolution is done by a numerical integration, that is why the solution is only ``quasi-analytical''. You can improve the accuracy of the solution by increasing the number of intervals used for the integration (the variable Nint in the data file Gar6more2D.dat)

Acoustic

The code computes a seismogram at points \((x_i,y)_{i=1,Nx}\) of the pressure solution of the equations (in the following \(\mathbf{x_s}=(0,h)\))

  • Infinite Medium

    \begin{eqnarray*} \derp {^2 P^+}{t^2}-{c^+}^2 \Delta P^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad (x,y)\in\setR^2. \end{eqnarray*}

  • Free Boundary Condition

    \begin{eqnarray*} \derp {^2 P^+}{t^2}-{c^+}^2 \Delta P^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad x\in\setR, y>0, \end{eqnarray*}

    with the boundary condition

    \begin{eqnarray*} P^+=0 \quad x\in\setR, y=0. \end{eqnarray*}

  • Wall Boundary Condition

    \begin{eqnarray*} \derp {^2 P^+}{t^2}-{c^+}^2 \Delta P^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad x\in\setR, y>0, \end{eqnarray*}

    with the boundary condition

    \begin{eqnarray*} \frac{\partial P^+}{\partial y}=0 \quad x\in\setR, y=0. \end{eqnarray*}

  • Bilayered Medium

    \begin{eqnarray*} \derp {^2 P^+}{t^2}-{c^+}^2 \Delta P^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad x\in\setR, y>0,\\[10pt] \derp {^2 P^-}{t^2}-{c^-}^2 \Delta P^-=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad x\in\setR, y<0, \end{eqnarray*}

    with the transmission conditions
    \(P^+=P^-\),
    \(\dsp \rho^+ \frac{\partial P^-}{\partial y}=\rho^- \frac{\partial P^+}{\partial y},\quad x\in\setR, y=0.\)

The code also computes the velocity given by the relation~:

\[\derp {\mathbf V^\pm} t =-\frac{1}{\rho^\pm}{\mathbf \nabla} P^\pm.\]

If you want to compute the displacement \(U\), it can be easily computed by replacing \(f(t)\) by the primitive of the source function you are using. For instance, if you are using a Rickert, you'll have to consider a first derivative of a Gaussian for \(f\).

Acoustic/elastodynamic (isotropic)

The code computes a seismogram at point \((x_i,y)_{i=1,Nx}\) of the pressure (in the fluid) and the velocity (in the solid) solution of the equations

\begin{eqnarray} \derp {^2 P^+}{t^2}-{c^+}^2 \Delta P^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\quad x\in\setR, y>0,\\[10pt] \derp {^2 \mathbf{V}^-}{t^2}-(\lambda^-+2\mu^-) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}})+\mu^-\nabla\times(\nabla \times {\mathbf{V}^-})=0,\quad x\in\setR, y<0, \end{eqnarray}

with \(\mathbf{x_s}=(0,h)\) and the transmission conditions

\(\dsp \derp {V^-_y}t=-\frac{1}{\rho^+}\derp {P^+}{y}\)\(\quad y=0\),
\(\dsp (\lambda^-+2\mu^-) \derp {V_y^-}{y}+\lambda^- \derp{V_x^-}{x}=\derp{P^+}{t}\)
\(\dsp \derp{V_x^-}y+\derp {V_y^-} x=0.\)

on the interface \(y=0\). The code also computes the velocity in the fluid by using the relation

\[\derp {V^+} t =-\frac{1}{\rho^+}\nabla P^+.\]

Once again, if you want to compute the displacement \(U\), it can be easily computed by replacing \(f(t)\) by the primitive of the source function you are using.

Acoustic/poroelastic.

(see[4]) The code computes a seismogram at point \((x_i,y)_{i=1,Nx}\) of the pressure \(P^+\) and the displacement \(U^+\) (in the fluid) and the solid displacement \(U_s^-\) (in the poroelastic medium) solution of the equations

\[ \begin{array}{lr} \left| \begin{array}{l} \dsp \derp {^2 P^+}{t^2}-{c^+}^2 \Delta \chi^+=\delta({\mathbf x}-\mathbf{x_s})f(t),\\[10pt] \dsp \derp {^2U^+} {t} =-\frac{1}{\rho^+}\nabla P^+. \end{array} \right. \end{array}\]

for \(x\in\setR, y>0\) and

\[\begin{array}{lr} \left| \begin{array}{l} \dsp (1-\phi^-)\rho_s^- \derp {^2 \mathbf U_s^-}{t^2}+\phi\rho_f^-\derp {^2 \mathbf U_f^- }{t^2}-(\lambda^-+2\mu^-) \nabla (\nabla \cdot \mathbf U_s^-)+\mu^-\nabla\times(\nabla \times \mathbf U_s^-)+\beta\,\nabla P^-=0 ,\\[10pt] \dsp (1-a^-)\rho_f^- \derp {^2 \mathbf U_s^-}{t^2}+a^-\rho_f^-\derp {^2 \mathbf U_f^-}{t^2}+\nabla P^-=0\\[10pt] \dsp \frac 1 {m^-} P^- +(\beta^--\phi^-) \nabla \cdot \mathbf U_s^-+ \phi^- \nabla \cdot \mathbf U_f^-=0 \end{array} \right. \end{array}\]

for \(x\in\setR, y<0\), either with the open pore transmission conditions (if parameter open is set to 1 in Gar6more2D.dat))

\(\dsp \phi^-(U_{fy}^--U_{sy}^-)=U^+_y-U_{sy}^-,\)
\(\dsp P^-=P^+,\)
\(\dsp \big(\lambda^-+m^-\beta^-(\beta^--\phi^-)\big)\,\nabla\cdot \mathbf U_s^-+2\mu^-\derp{U_{sy}^-}{y}+m^-\beta^-\phi^-\nabla\cdot \mathbf U_f^- =-P^+,\)
\(\dsp\derp{U_{sx}^-}{y}+\derp{U_{sy}^-}{x}=0\),

or with the sealed pore transmission conditions (if parameter open is set to 0 in Gar6more2D.dat))

\(\dsp \phi^-(U_{fy}^--U_{sy}^-)=U^+_y-U_{sy}^-,\)
\(U_{fy}^-=U_{sy}^-,\)
\(\dsp \big(\lambda^-+m^-\beta^-(\beta^--\phi^-)\big)\,\nabla\cdot \mathbf U_s^-+2\mu^-\derp{U_{sy}^-}{y}+m^-\beta^-\phi^-\nabla\cdot \mathbf U_f^- =-P^+,\)
\(\dsp\derp{U_{sx}^-}{y}+\derp{U_{sy}^-}{x}=0\),

on the interface \(y=0\). The code does not compute the fluid displacement and the pressure in the poroelastic medium, but there is no particular difficulty to do that.

Elastodynamic

The code computes a seismogram at points \((x_i,y)_{i=1,Nx}\) of the velocity solution of the equations (in the following \(\mathbf{x_s}=(0,h)\))

  • (Infinite Medium)

    \begin{eqnarray*} \derp {^2 \mathbf{V}^+}{t^2}-(\lambda^++2\mu^+) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}^+})+\mu^+\nabla\times(\nabla \times {\mathbf{V}^+})=F(\mathbf x,t),\quad \mathbf x \in\setR^2. \end{eqnarray*}

    with

    \[F(\mathbf x,t)=\nabla \delta({\mathbf x}-\mathbf{x_s})f_P(t)+\nabla\times \delta({\mathbf x}-\mathbf{x_s})f_S(t)\]

    The first term of \(F\) represents a \(P-\)source while the second one represents a \(S-\)source.
  • (Free Boundary Condition)

    \begin{eqnarray*} \derp {^2 \mathbf{V}^+}{t^2}-(\lambda^++2\mu^+) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}^+})+\mu^+\nabla\times(\nabla \times {\mathbf{V}^+})=F(\mathbf x,t),\quad x\in\setR, y>0. \end{eqnarray*}

    \ with the boundary conditions

    \begin{eqnarray*} \dsp\derp{V_{y}^+}{y}=0 \hbox{ and } \dsp \derp{ V_{x}^+}{y} + \derp {V_{y}^+}{x} =0 \quad x\in\setR, y=0. \end{eqnarray*}

  • (Wall Boundary Condition)

    \begin{eqnarray*} \derp {^2 \mathbf{V}^+}{t^2}-(\lambda^++2\mu^+) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}^+})+\mu^+\nabla\times(\nabla \times {\mathbf{V}^+})=F(\mathbf x,t),\quad x\in\setR, y>0. \end{eqnarray*}

    with the boundary condition

    \begin{eqnarray*} {\mathbf{V}^+}=0\quad x\in\setR, y=0. \end{eqnarray*}

  • (Bilayered Medium)

    \begin{eqnarray*} \derp {^2 \mathbf{V}^+}{t^2}-(\lambda^++2\mu^+) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}^+})+\mu^+\nabla\times(\nabla \times {\mathbf{V}^+})=F(\mathbf x,t),\quad x\in\setR, y>0,\\[10pt] \derp {^2 \mathbf{V}^-}{t^2}-(\lambda^-+2\mu^-) {\mathbf \nabla} (\nabla \cdot {\mathbf{V}^-})+\mu^-\nabla\times(\nabla \times {\mathbf{V}^-})=0,\quad x\in\setR, y<0, \end{eqnarray*}

    with the transmission conditions\[10pt]
    \(\dsp V^+=V^-\),
    \(\dsp(\lambda^++2\mu^+)\derp{V_{y}^+}{y}=(\lambda^-+2\mu^-)\derp{V_{y}^-}{y}\),
    \(\dsp \mu^+\left(\derp{ V_{x}^+}{y} + \derp {V_{y}^+}{x}\right) =\mu^-\left(\derp{ V_{x}^-}{y} + \derp {V_{y}^-}{x}\right) ,\quad x\in\setR, y=0.\)
    Once again, if you want to compute the displacement \(U\), it can be easily computed by replacing \(f(t)\) by the primitive of the source function you are using.

Poroelastic.

(see[3]) The code computes a seismogram at point \((x_i,y)_{i=1,Nx}\) of the solid displacement \(U_s\) solution of the equations (in the following \(\mathbf{x_s}=(0,h)\))

  • (Infinite Medium)

    \[\hspace{-1cm} \begin{array}{lr} \!\!\!\!\left| \begin{array}{l} \dsp (1-\phi^+)\rho_s^+ \derp {^2 \mathbf U_s^+}{t^2}+\phi^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}-\!(\lambda^++2\mu^+) \nabla (\nabla \cdot \mathbf U_s^+)\!+\!\mu^+\nabla\times(\nabla \times \mathbf U_s^+)+\beta^+\nabla P^+\!=\!\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t) ,\\[10pt] \dsp (1-a^+)\rho_f^+ \derp {^2 \mathbf U_s^+}{t^2}+a^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}+\nabla P^+=\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t)\\[10pt] \dsp \frac 1 {m^+} P^+ +(\beta^+-\phi^+) \nabla \cdot \mathbf U_s^++ \phi^+ \nabla \cdot \mathbf U_f^+=\delta({\mathbf x}-\mathbf{x_s}) F_p(t) \end{array} \right. \end{array}\]

    for \((x,y)\in\setR^2.\) Actually the code computes the solution for each source \(F_s\) and \(F_p\) separately. If you want a bulk source ( \(F_s\)), set the parameter \(\it type\_source\) to 1, if you want a pressure source ( \(F_p\)), set the parameter \(\it type\_source\) to 2.
  • (Wall Boundary Condition)

    \[\hspace{-1cm} \begin{array}{lr} \!\!\!\!\left| \begin{array}{l} \dsp (1-\phi^+)\rho_s^+ \derp {^2 \mathbf U_s^+}{t^2}+\phi^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}-\!(\lambda^++2\mu^+) \nabla (\nabla \cdot \mathbf U_s^+)\!+\!\mu^+\nabla\times(\nabla \times \mathbf U_s^+)+\beta^+\nabla P^+\!=\!\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t) ,\\[10pt] \dsp (1-a^+)\rho_f^+ \derp {^2 \mathbf U_s^+}{t^2}+a^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}+\nabla P^+=\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t)\\[10pt] \dsp \frac 1 {m^+} P^+ +(\beta^+-\phi^+) \nabla \cdot \mathbf U_s^++ \phi^+ \nabla \cdot \mathbf U_f^+=\delta({\mathbf x}-\mathbf{x_s}) F_p(t) \end{array} \right. \end{array}\]

    for \(x\in\setR, y>0\), with the boundary conditions

    \[U_{fy}^+-U_{sy}^+=0 \hbox{ and } U_{s}^+=0, \quad x\in\setR, y=0.\]

  • (Free Boundary Condition)

    \[\hspace{-1cm} \begin{array}{lr} \!\!\!\!\left| \begin{array}{l} \dsp (1-\phi^+)\rho_s^+ \derp {^2 \mathbf U_s^+}{t^2}+\phi^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}-\!(\lambda^++2\mu^+) \nabla (\nabla \cdot \mathbf U_s^+)\!+\!\mu^+\nabla\times(\nabla \times \mathbf U_s^+)+\beta^+\nabla P^+\!=\!\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t) ,\\[10pt] \dsp (1-a^+)\rho_f^+ \derp {^2 \mathbf U_s^+}{t^2}+a^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}+\nabla P^+=\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t)\\[10pt] \dsp \frac 1 {m^+} P^+ +(\beta^+-\phi^+) \nabla \cdot \mathbf U_s^++ \phi^+ \nabla \cdot \mathbf U_f^+=\delta({\mathbf x}-\mathbf{x_s}) F_p(t) \end{array} \right. \end{array}\]

    for \(x\in\setR, y>0\), with the boundary conditions \[10pt]
    \(\dsp \alpha^+\nabla\cdot \mathbf U_s^+ +2\mu^+\derp{U_{sy}^+}{y} + m^+\beta^+\phi^+\,\nabla\cdot \mathbf U_f^+=0,\)
    \(\dsp \derp{ U_{sx}^+}{y} + \derp {U_{sy}^+}{x}=0,\) \( P^+=0,\)
    for \(x\in\setR, y=0\), with

    \[ \alpha^\pm=\lambda^\pm+m^\pm\beta^\pm(\beta^\pm-\phi^\pm). \]

  • (Bilayered Medium)

    \[\hspace{-1cm} \begin{array}{lr} \!\!\!\!\left| \begin{array}{l} \dsp (1-\phi^+)\rho_s^+ \derp {^2 \mathbf U_s^+}{t^2}+\phi^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}-\!(\lambda^++2\mu^+) \nabla (\nabla \cdot \mathbf U_s^+)\!+\!\mu^+\nabla\times(\nabla \times \mathbf U_s^+)+\beta^+\nabla P^+\!=\!\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t) ,\\[10pt] \dsp (1-a^+)\rho_f^+ \derp {^2 \mathbf U_s^+}{t^2}+a^+\rho_f^+\derp {^2 \mathbf U_f^+}{t^2}+\nabla P^+=\nabla\delta({\mathbf x}-\mathbf{x_s}) F_s(t)\\[10pt] \dsp \frac 1 {m^+} P^+ +(\beta^+-\phi^+) \nabla \cdot \mathbf U_s^++ \phi^+ \nabla \cdot \mathbf U_f^+=\delta({\mathbf x}-\mathbf{x_s}) F_p(t) \end{array} \right. \end{array}\]

    for \(x\in\setR, y>0\) and

    \[\hspace{-1cm} \begin{array}{lr} \!\!\!\!\left| \begin{array}{l} \dsp (1-\phi^-)\rho_s^- \derp {^2 \mathbf U_s^-}{t^2}+\phi^-\rho_f^-\derp {^2 \mathbf U_f^-}{t^2}-(\lambda^-+2\mu^-) \nabla (\nabla \cdot \mathbf U_s^-)+\mu^-\nabla\times(\nabla \times \mathbf U_s^-)+\beta^-\nabla P^-=0 ,\\[10pt] \dsp (1-a^-)\rho_f^- \derp {^2 \mathbf U_s^-}{t^2}+a^-\rho_f^-\derp {^2 \mathbf U_f^-}{t^2}+\nabla P^-=0\\[10pt] \dsp \frac 1 {m^-} P^- +(\beta^--\phi^-) \nabla \cdot \mathbf U_s^-+ \phi^- \nabla \cdot \mathbf U_f^-=0 \end{array} \right. \end{array}\]

    for \(x\in\setR, y<0\), with the transmission conditions on the interface \(y=0\)
    \(\phi^-\,(U_{fy}^--U_{sy}^-)=\phi^+\,(U_{fy}^+-U_{sy}^+),\)
    \(\dsp \alpha^-\nabla\cdot \mathbf U_s^- +2\mu^-\derp{ U_{sy}^-}{y}+m^-\beta^-\phi^-\,\nabla\cdot \mathbf U_f^-= \alpha^+\nabla\cdot \mathbf U_s^+ +2\mu^+\derp{U_{sy}^+}{y} + m^+\beta^+\phi^+\,\nabla\cdot \mathbf U_f^+,\)
    \(\dsp \mu^-(\derp {U_{sx}^-}{y}+\derp{U_{sy}}{x}^-)= \mu^+(\derp{ U_{sx}^+}{y} + \derp {U_{sy}^+}{x}),\)
    \(U_{sx}^-=U_{sx}^+,\quad U_{sy}^-=U_{sy}^+,\quad P^-=P^+,\)
    with

    \[ \alpha^\pm=\lambda^\pm+m^\pm\beta^\pm(\beta^\pm-\phi^\pm). \]

Remark

The code does not really compute the displacement, but its derivative (for some reasons related to the Cagniard-de Hoop method, see~[6][3]). Therefore, you have to replace \(f(t)\) by the primitive of the source function you are using to compute the displacement.

Installation and Usage

1) Please read the makefile and modify it according to your environment. In particular adapt the path to the lapack libraries. 2) Type make all 3) Modify the file Gar6more2D.dat (you may find examples in the directory EXAMPLES) 4) Run Gar6more2D.out 5) The results are stored in the files Ux.dat, Uy.dat and P.dat.

Description

I) The main program is lib/bin/Gar6more2D.F90.

II) The directory mod/ contains all files for the declaration of the common variables.

III) The directory lib/libgeneral contains functions used by most of the subroutines : the source function source.F90, the computation of the path in the complex plane path.F90, the computation of the derivate dp/dt derivee.F90.

IV) The directory includefic contains the interface files.

V) The directory lib/acousacous contains the subroutines for the computation of the incident wave (sub_incid_aa.F90), the reflected wave (sub_reflex_*_aa.F90) and the transmitted wave (sub_transmit_*_aa.F90) in the acoustic/acoustic case.

VI) The directory lib/acouselasto contains the subroutines for the computation of the incident wave (sub_incid_ae.F90), the reflected wave (sub_reflex_ae.F90), the transmitted P wave (sub_transmitp_ae.F90) and the transmitted S wave (sub_transmits_ae.F90) in the acoustic/elastodynamic case.

VII) The directory lib/acousporo contains the subroutines for the computation of the incident wave (sub_incid_ap.F90), the reflected wave (sub_reflex_ap.F90), the transmitted fast wave (sub_transmitf_ap.F90), the transmitted psi wave (sub_transmitpsi_ap.F90) and the transmitted slow wave (sub_transmits_ap.F90) in the acoustic/poroelastic case. The computation of the reflection and transmission coefficients is performed by calccoeff_ap.F90

VIII) The directory lib/elastoelasto contains the subroutines for the computation of the incident waves (sub_incid*_ee.F90), the reflected waves (sub_reflex*_*_ee.F90) and the transmitted waves (sub_transmit*_ee.F90) in the elastodynamic/elastodynamic case. The computation of the reflection and transmission coefficients (for P and S waves) is performed by calccoeff*_ee.F90.

IX) The directory lib/poroporo contains the subroutines for the computation of the incident waves (sub_incid_pp_*.F90), the reflected waves (sub_reflex*_*_pp.F90) and the transmitted waves (sub_transmit*_*_pp.F90) in the poroelastic/poroelastic case.