<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">%\documentclass[paper]{geophysics}
\documentclass[manuscript,noblind]{geophysics}  %,revised]
\usepackage{subeqnarray,amsmath,amsfonts,latexsym,epsfig,color}
\usepackage{graphicx}
\usepackage{lineno}
\usepackage{url}
\usepackage{amsthm}
\newtheorem{theorem}{Theorem}

% An example of defining macros
\newcommand{\rs}[1]{\mathstrut\mbox{\scriptsize\rm #1}}
\newcommand{\rr}[1]{\mbox{\rm #1}}


%\usepackage{subeqnarray,amsmath,amsfonts,graphicx}
%\usepackage{oldgerm}
%\usepackage{rotating}
\bibliographystyle{seg} 
\oddsidemargin -.5cm
\newcommand{\BM}[1]{\mbox{\boldmath $ #1 $}}
\newcommand{\myfrac}[2]{\displaystyle \frac{#1}{#2}}
\newcommand{\ii}{{\rm i}}
\newcommand{\fo}{\rm}
\def\bv{{\mathbf v}}
\def\bu{{\mathbf u}}
\def\bw{{\mathbf w}}
\newcommand{\bsig}{{\boldsymbol \sigma}}
\newcommand{\btau}{{\boldsymbol \tau}}
\newcommand{\bvarep}{{\boldsymbol \varepsilon}}
\newcommand{\bnu}{{\boldsymbol \nu}}
\renewcommand\tilde{\widetilde}
\renewcommand\a{\alpha}
\def\b{\beta}
\newcommand\bA{\mathbf{A}}
\def\bmu{\boldsymbol\mu}
\renewcommand\d{\delta}
\newcommand\D{\Delta}
\newcommand{\fu}{\rm}
\def\dfrac{\displaystyle\frac}
\def\Sum{\displaystyle\sum}
\def\ai{\left\langle\langle}
\def\ad{\right\rangle\rangle}
\def\b{\beta}
\def\k{\kappa}
\def\opdiv{\operatorname{div}}
\def\l{\ell}
\def\D{\Delta}
\def\d{\delta}
\def\e{\varepsilon}
\def\dint{\displaystyle\int}
\def\g{\gamma}
\def\G{\Gamma}
\def\h{{\mathbf h}}
\def\hf{ f}
\def\hR{\widehat R}
\def\hu{ u}
\def\hz{\widehat z}
\def\hZ{\widehat Z}
\def\Im{\operatorname{Im}}
\def\v{{\mathbf v}}
\def\H{{\mathbf H}}
\def\n{{\mathbf n}}
\def\R{{\mathbb R}}
\def\b{\beta}
\def\lla{\langle\langle}
\def\Lam{\Lambda}
\def\lam{\lambda}
\def\some{{x,y,z}}
\def\Lip{\operatorname{Lip}}
\def\tr{\operatorname{tr}}
\def\curl{{\nabla\times}}
\def\div{{\nabla\cdot}}
\def\opcurl{\operatorname{curl}}
\def\opdiv{\operatorname{div}}
\def\nab{{\nabla}}
\def\grad{{\nabla}}
\def\g{\gamma}
\def\q{\quad}
\def\2q{\quad\quad}
\def\3q{\quad\quad\quad}
\def\s{\sigma}
\def\sG{\mathcal{G}}
\def\sL{\mathcal{L}}
\def\sB{\mathcal{B}}
\def\sA{\mathcal{A}}
\def\cV{\mathcal{V}}
\def\sm{\setminus}
\def\NC{\mathcal{NC}}
\def\sG{\mathcal{G}}
\def\sm{\setminus}
\def\t{\theta}
\def\T{\Theta}
\def\Tau{\mathcal T}
\def\talpha{\widetilde\alpha}
\def\tA{\tilde A}
\def\tB{\tilde B}
\def\tC{\tilde C}
\def\V{\Vert}
\def\z{\zeta}
\def\xo{(x,\omega)}
\def\&lt;{\left\langle}
\def\&gt;{\right\rangle}
\def\div{\nabla\cdot}
\def\grad{\nabla}
\def\lma{\left\langle}
\def\la{\left\langle}
\def\lla{\langle\langle}
\def\rra{\rangle\rangle} 
\def\rma{\right\rangle} 
\def\ra{\right\rangle}
\def\qq{\quad\quad}
\def\dt{{\varDelta t}}
\def\dx{{\varDelta x}}
\def\dy{{\varDelta y}}
\def\dz{{\varDelta z}}
\def\BbbR{\mathbb R}
\def\O{\Omega}
\def\o{\omega}
\def\Oh{{\mathcal O}}
\def\mD{{\mathcal D}}
\def\mS{{\mathcal S}}
\def\mT{{\mathcal T}}
\def\mC{{\mathcal C}}
\def\mN{{\mathcal N}}
\def\p{\partial}
\def\for{\,\forall}
\def\a{\alpha}
\def\nb{\nonumber}
\def\var{\varphi}
\def\w{\omega}
\def\G{\Gamma}
\def\R{{\mathbb R}}
\def\:{{\,:\,}}
\def\Set#1{{\left\{#1\right\}}}
\def\norm#1{{\left\|#1\right\|}}
\def\qnorm#1{{\left|\left|\left|#1\right|\right|\right|}}
\def\abs#1{{\left|#1\right|}}
\def\th#1{\theta_{jk}^{#1}}
\def\fis{S_1}
\def\fii{S_3}
\def\tot1{\sigma_{ij}^{(1,T)}}
\def\tot2{\sigma_{ij}^{2,T}}
\def\tot{\sigma_{ij}}
\renewcommand{\hat}{\widehat}
\newcommand{\pp}[2]{\frac{\partial {#1}}{\partial {\nu_{#2}}}} 
\newcommand{\pu}[2]{P_{\tau}u_{#1}^{#2}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\bchi}{{\boldsymbol \chi}}
\newcommand{\bpsi}{{\boldsymbol \psi}}
\newcommand{\bvarphi}{{\boldsymbol \varphi}}
\newcommand{\bgamma}{{\boldsymbol \gamma}}
\newcommand{\bTheta}{{\boldsymbol \Theta}}
\def\sA{\mathcal{A}}
\def\mV{\mathcal{V}}
\def\cV{\mathcal{V}}
\def\bU{{\mathbf U}}
\def\bE{{\mathbf E}} 
\def\bu{{\mathbf u}}
\def\be{{\mathbf e}}
\def\bv{{\mathbf v}}
\def\bx{{\mathbf x}}

\begin{document}

\title{Physics-Based Simulation of Poroelastic Anisotropy in Organic-Rich Mudrocks from the Vaca Muerta Formation}

\renewcommand{\thefootnote}{\fnsymbol{footnote}} 

\ms{} % manuscript number

\address{
\footnotemark[1]School of Earth Sciences and Engineering, \\
Hohai University, \\
Nanjing, 211100, China, jingba@188.com \\
\footnotemark[2]Universidad de Buenos Aires, 
Facultad de Ingenier\'\i  a, \\
Instituto del Gas y del Petr\'oleo, \\
Av. Las Heras 2214 Piso 3 \\
C1127AAR Buenos Aires, Argentina, gsavioli@fi.uba.ar \\
\footnotemark[3]Department of Mathematics,  Purdue University, \\
150 N. University  Street, West Lafayette, \\ 
Indiana, 47907-2067, USA, jesantos48@gmail.com \\
\footnotemark[4] Facultad de Ciencias Astron\'omicas y Geof\'\i sicas, Universidad Nacional de La Plata \\
, pgauzellino@gmail.com  \\
\footnotemark[5]Solaer Ingenier\'\i  a  , arcamus@gmail.com  \\}
\author{Juan Enrique Santos\footnotemark[2]\footnotemark[1]\footnotemark[3], Ariel S\'anchez Camus\footnotemark[4]\footnotemark[5], \\ Gabriela Beatriz Savioli\footnotemark[2], \\ Patricia Mercedes Gauzellino\footnotemark[4]and 
Jing Ba\footnotemark[1]}

\footer{Example}
\lefthead{Santos et al.}
\righthead{Anisotropy in Organic-Rich Mudrocks}

\maketitle

\begin{abstract}
This work presents a methodology for the advanced characterization of ultra-low permeability,
organic-rich mudrock reservoirs. The quiescent depositional environments of these formations
promote the development of a stratified microstructure, characterized by the preferential alignment
of clay platelets and organic matter parallel or sub-parallel to bedding. This finely layered,
multimineral architecture results in pronounced vertical transverse isotropy (VTI), reflecting a
strong mechanical anisotropy inherent to the rock fabric.
Using a theoretical rock physics model that integrates mechanical, mineralogical, and petrophysical
data, we estimate the stiffness tensor of a core sample extracted from the lower section of the Vaca
Muerta Formation, located within the oil generation window at a depth of 3100 meters in the
NeuquÃ©n Basin, Argentina. Frequency-dependent stiffness coefficients of the equivalent VTI
medium are derived from numerical harmonic compressibility and shear experiments conducted on
a representative 2D synthetic sample, modeled as a periodic sequence of two materials, each with
6 percent  porosity. Material 1 is a complex multimineral assemblage consisting of seven solid phasesâ€”
including 23 percent  kerogen treated as a mineral constituentâ€”whereas Material 2 consists of pure
kerogen.
The dry-rock stiffness tensor is initially computed and validated against laboratory phase velocity
measurements at 1 MHz, showing discrepancies below 10 percent. The synthetic sample is subsequently
saturated with hydrocarbon fluids to evaluate attenuation and dispersion effects associated with
wave-induced fluid flow (WIFF) and the resulting mode conversions. Fluid-saturated simulations
reveal that phase velocities decrease markedly with increasing gas content,  with
frequency variations depending on the fluid mixture, particularly above 100 Hz. Patchy gas-oil saturation
yields higher velocities than uniformly mixed fluids, while exhibiting significantly greater
attenuation due to enhanced WIFF. The frequency-dependent peak attenuation shifts toward lower
frequencies in patchy cases, in agreement with mesoscopic pressure diffusion theory. These results
emphasize the critical role of fluid composition and spatial distribution in seismic velocity and
attenuation, supporting the need for explicit modeling of patch-scale heterogeneity in
unconventional reservoir characterization.
\end{abstract}

\section{Introduction}

Organic-rich mudrock reservoirs are anisotropic formations typically characterized by elevated pore pressures, primarily due to the thermal cracking of organic matter into hydrocarbons
\citep{pinna2011,alduha2014} 
%(Pinna et al., 2011; Al Duhailan, 2014). 
These formations exhibit vertical transverse isotropy (VTI), defined by a vertical axis of symmetry and described by five independent stiffness coefficients 
â€”$(p_{11},
p_{13}, p_{33}, p_{55}$ and  $p_{66})$â€”
which govern the stressâ€“strain relationships of the equivalent VTI medium (Vernik and Liu, 1997). This equivalent medium captures the effective elastic behavior of the reservoir at seismic wavelengths significantly larger than the average lamination thickness \citep{vernik92}.


The petrophysical complexity of these rocksâ€”characterized by  micro- to nano-scale pore systems, ultra-low permeability, and pronounced permeability anisotropy 
%(Nelson, 2009; 
\citep{nelson2009,loucks2012,mokhtari2015,zoback2019}
%\cite{loucks2012}
%\cite{mokhtari2015}
%\cite{zoback2019}
%Loucks et al., 2012; 
%Mokhtari and Tutuncu, 2015; 
%Zoback and Kohli, 2019)
â€”is further amplified by the effects of thermal maturity. As organic matter evolves across the hydrocarbon generation window, significant transformations occur in pore structure, fluid composition and rheology, and in the mechanical properties of the solid rock frame \citep{vernik2011,yang2019,allan2016}. 
%( Vernik and Milovac, 2011;  
%Yang, 2019; 
%Allan et al., 2016). 
These coupled processes collectively exert a strong influence on the mechanical and seismic response of organic-rich mudrocks.


The objective of this study is to apply the two-dimensional numerical model developed in \cite{santos2015} 
â€”a physically grounded, non-phenomenological formulation based on Biotâ€™s theory of poroelasticity \cite{biot62} to investigate the anisotropic poroelastic response of the Vaca Muerta Formation (VMF), an organic-rich marlstone situated in the Neuqu\'en Basin, Argentina, under variable fluid saturation conditions and seismic range of  frequencies.

Rather than aiming for an exact match with laboratory-measured velocities, the primary goal is to evaluate whether the computed stiffness matrix provides a reliable approximation of the formationâ€™s effective mechanical behavior. Discrepancies between modeled and experimental results are  expected, largely due to the idealization of the medium as perfectly laminatedâ€”a valid assumption  for clay-rich shales, but less representative of marly facies, where compositional heterogeneities  such as dispersed carbonate grains and silt-sized particles reduce the intrinsic anisotropy. Additionally, the two-dimensional nature of the numerical model imposes further simplifications by not fully capturing the formationâ€™s three-dimensional microstructure.

As a basis for the numerical modeling, we employ the rock physics theory developed by \cite{carcione2005}, which integrates mineralogical and petrophysical data from a core sample extracted  from the VMF. This approach enables the estimation of the dry bulk modulus K$_m$ and shear modulus $\mu_m$ of the multi-mineral porous matrix. Then, the VMF is assumed to consist of a  sequence of very thin horizontal layers composed of kerogen (as the primary organic constituent)  and mineral aggregates, within which Biotâ€™s theory in the diffusive frequency range is applicable.

To determine the p${ij}$ stiffness coefficients of the equivalent VTI medium, a set of five boundary  value problems (BVPs) based on Biotâ€™s diffusive equations are formulated in the spaceâ€“frequency  domain. Each BVP defines a compressibility or shear test, which is solved numerically using the  Finite Element (FE) method across a range of relevant frequencies (see \cite{santos2015}
%(Santos and Carcione, 2015) 
for details). These tests are first applied to a dry synthetic sample to compute the stiffness tensor and  derive the corresponding dry-rock phase velocities. The simulated phase velocities exhibit good  agreement with laboratory measurements, with discrepancies remaining below 10 percent. These  results indicate that the stiffness matrix obtained through the numerical model provides a representative characterization of the elastic behavior of the lower section of the VMF formation at  the studied well location.


Following the dry-rock analysis, the synthetic sample is saturated with gas, oil, and water to  
investigate attenuation and dispersion effects associated with
mode conversions caused by wave-induced fluid flow (WIFF) \citep{white75,KM_11}.
Additionally, the case  of patchy gasâ€“oil saturation is examined and compared against scenarios 
of uniform saturation, in  order to evaluate the effects of fluid distribution on seismic velocity and attenuation.

\section*{The Biot Model for a fluid saturated  poroelastic solid}

Let us denote by $\bu(\bx) =  (\bu_s(\bx), \bu_f(\bx))$ the time Fourier transforms of the solid 
and relative fluid displacements. 
With $\omega = 2 \pi f$ denoting the angular frequency, it is assumed that on each thin horizontal layer the diffusive  Biot's equations hold:

\begin{eqnarray}\label{def_Lu}
(\nabla\cdot \bsig (\bu), {\rm i} \omega {\frac{\eta}{\kappa} \bu_f} -\nabla p_f(\bu)) = {\mathcal C}, \qquad 
%\end{eqnarray}
{\rm with} \quad 
%\begin{eqnarray}
{\mathcal C}=
\begin{pmatrix}\label{def_C}
0 I_2&amp; 0 I_2 \\
0 I_2 &amp;\dfrac{\eta}{\kappa} I_2
\end{pmatrix}, 
\end{eqnarray}
where $I_2$ is the $2\times 2$ identity matrix, $\eta$ is the fluid viscosity 
and  $\kappa$ is the frame permeability.
In \eqref{def_Lu} $\bsig$ and $p_f$   are the bulk stress and fluid pressure, respectively,  satisfying the  constitutive relations
\begin{eqnarray}
&amp;\sigma_{kl}(\bu) =  2 \mu\, e_{kl}(\bu_s)+\delta_{kl}
\left( \lambda_{u} \,\nabla\cdot \bu_s  + \alpha {\rm M}  \nabla\cdot \bu_f \right),\label{st_6}\\
&amp;p_f(\bu) =  - \alpha {\rm M} \nabla\cdot \bu_s -  {\rm M}  \nabla\cdot \bu_f.\label{st_7}
\end{eqnarray} 
The coefficient $\mu$ is the shear modulus of the dry frame. Let $\phi$ denote the porosity and  K$_s$, K$_m$, 
K$_f$, K$_u$ denote   the bulk modulus of the solid grain, the dry matrix, the  (single-phase) fluid,  and K$_u$ the undrained bulk modulus, respectively. Then    the other coefficients in \eqref{st_6}-\eqref{st_7} 
are given by the relations

\begin{eqnarray}\label{constitutive}
 {\rm K}_u = {\rm K}_m + \alpha^2 {\rm M}, \quad  \alpha = 1 - \dfrac{{\rm K}_m}{{\rm K}_s},\quad, \lambda_u= {\rm K}_u - \frac{2}{3} \mu, \quad 
  {\rm M} = \left(\dfrac{\alpha- \phi}{{\rm K}_s} 
+ \dfrac{\phi}{{\rm K}_f}\right)^{-1}. 
\end{eqnarray}
%\begin{eqnarray}\label{constitutive}
%{\rm K}{_u} = {\rm K}_m} + \alpha^2 {\rm M}, \quad  \alpha = 1 - \dfrac{{\rm K}_m}{{\rm K}_s},
% \quad \lambda_u= {\rm K}_u - \frac{2}{3} \mu,  \quad 
%{\rm M} = \left(\dfrac{\alpha- \phi}{{\rm K}_s} 
%+ \dfrac{\phi}{K_f}\right)^{-1}. 
%\end{eqnarray}



\section{Determination of the dry bulk and shear moduli for multimineral solids}
A core sample retrieved from a well at a depth of 3094.68 m within the VMF was dried under
laboratory conditions, and ultrasonic phase velocity measurements were conducted at a frequency
of 1 MHz. The results exhibit a clear VTI behavior in the formation. Figure 1 illustrates the
mineralogical composition along the well-bore, as well as  with the specific depth interval
from which the core was extracted.


To characterize the frequency-dependent elastic response of the equivalent VTI medium, numerical
harmonic simulationsâ€”both compressional and shearâ€”were performed on a representative
synthetic rock model. This model consists of a periodic alternation of two porous materials, each
with a porosity of 6 percent. Material 1 is a complex, multimineral assemblage composed of seven solid
phases, including 23 percent of  kerogen, which is treated as a solid mineral constituent. Material 2, in
contrast, is composed entirely of pure kerogen.

The fluid properties used in the simulations are listed in Table 1. Table 2
  provides the physical properties of the minerals at the effective pressure where the laboratory
data were obtained, together with the mineralogical composition and corresponding volume
fractions for the numerical model.
Based on these data, porosity  is 6 percent, and the average density of Material 1
is calculated to be 2486.69 kg/m$^3$.

The procedure described in \cite{carcione2005} is employed to estimate the dry bulk modulus
K$_m$ and shear modulus $\mu_m$ of the multimineral  solid phases  constituting Material 1 within
the VMF. 
For this purpose,  generalized version of Kriefâ€™s model \citep{krief_90,mavko2020} 
is applied to
estimate the dry bulk and shear moduli  of the seven dry mineral constituents. These individual moduli are then combined, accounting for their volumetric fractions, to
obtain the effective dry bulk and shear moduli of Material 1:
\begin{eqnarray}
 {\rm K}_m = 29.195 \,  {\rm GPa}, \,\,\mu_m = 17.782 \, {\rm GPa}
\end{eqnarray}


It is worth noting that the computed values of K$_m$ and $\mu_m$  lie within the Hashinâ€“Shtrikman (HS)
bounds  \citep{HS63,mavko2020}, i.e., between the theoretical upper and
lower limits,${\rm K}^{+}$  ${\rm K}^{-}$, $\mu^{+}$  $\mu{^-}$: 

\begin{eqnarray}\label{HSbounds}
{\rm K}^+ = 33.367 \, {\rm GPa}, \quad {\rm K}^- =  19.52 \,{\rm GPa},\quad \quad \mu^+ = 19.53 \, {\rm GPa},\quad   \mu^- = 7.88 \, {\rm GPa}.
\end{eqnarray}


Further analysis of the  core sample, allowed us to  estimate that the value of the 
Biot's  coefficient in \eqref{st_6} is  $\alpha = 0.48$  
which in turn  yields  the value K$_s$ = 56.14 GPa. for the solid grain  bulk modulus of 
the composite solid.

Material 2, composed entirely of kerogen, is characterized by a porosity  of  6  percent, 
a grain modulus of  7 GPa, a density of  1400 kg/m$^3$,  and 
dry matrix  bulk modulus  and shear moduli of  
K$_m$ = 1.29 GPa and 0.36 GPa, respectively.

For both materials  in the periodic sequence  the permeability, denoted  $\kappa$ 
takes the value 2.75 10$^{-18}$ m$^{2}$.

%{\bf quizas haya que ponerlo en nanodarcy (nd ?}


After having determined the dry bulk and shear moduli, porosity, and permeability of both constituent
materials, we applied time-harmonic compressibility and shear numerical harmonic experiments \cite{santos2015} 
to a representative two-
dimensional model of the core. These numerical experiments were designed to assess the
anisotropic (VTI) elastic behavior of the VMF.


\section{Computed VTI phase  velocities  using the dry-core data.  }

The harmonic simulations were performed on a dry, square sample with a side length of 2 mm,
consisting of four repetitions of a periodic sequence composed of 49 layers of Material 1 and one
layer of Material 2 per period. Each layer has a uniform thickness of  1. 10$^{-5}$ m. The computational domain was discretized using a 200$\times$200 mesh.

Table 3 summarizes the results of the VTI analysis, showing a good agreement between the
simulated and laboratory-measured phase velocities, with relative errors remaining below 10 percent.
These discrepancies are primarily attributed to the assumption of a perfectly laminated, finely
stratified medium, whichâ€”while appropriate for clay-rich shalesâ€”becomes less accurate in the
presence of more compositionally heterogeneous mudrocks such as marly facies.
As noted by \cite{sone13}, 
%Sone and Zoback (2013),
marls such as VMF often contain a significant proportion of
silt-sized grains and carbonate minerals dispersed in a less aligned matrix, leading to a reduction in
intrinsic anisotropy relative to well-laminated shales. Consequently, part of the mismatch between
measured and modeled velocities can be ascribed to the microstructural variability inherent in such
mixed-facies systems. Table 4 presents the corresponding numerically derived stiffness coefficients
for the equivalent VTI medium.








Figures 2 and 3 present polar plots of the phase and energy velocities, respectively, for the qP, qSV,
and SH waves in the dry sample, computed using the FE method at a frequency of 1 MHz. A clear
and strong VTI behavior is observed, as evidenced by the pronounced directional dependence of
wave velocities. In Figure 2, the qP wave reaches significantly higher velocities along the horizontal
(x) direction compared to the vertical (z) axis, which is consistent with the expected behavior of
VTI media. The SH wave pattern is governed by the stiffness coefficient p$_{66}$, capturing the
mediumâ€™s horizontal shear response. Figure 3 reveals a notable triplication in the qSV-wave energy
velocity pattern, which is a hallmark of strong elastic anisotropy. This phenomenon highlights the
anisotropic coupling between the qP and qSV modes, further evidencing the complexity of wave
propagation in the medium.





\section*{Analysis of fluid saturated VMF}

A fluid-saturated porous rock exhibits poroelastic behavior, wherein the effective stiffness of the rock
depends strongly on the frequency range. At low frequenciesâ€”typical of surface
seismic surveysâ€”pore fluids have sufficient time to redistribute in response to wave-induced
deformation. Under these conditions, the rock behaves in a relaxed regime, where pore pressure
equilibrates as seismic waves propagate through the medium and contributes minimally to the
effective stiffness.

Conversely, at higher frequenciesâ€”such as those encountered in sonic well logs or ultrasonic
laboratory measurementsâ€”the deformation occurs too rapidly for pore fluid pressure to equilibrate.
This corresponds to an unrelaxed regime, in which the fluid is effectively trapped, resulting in
localized pore pressure build-up that supports part of the applied stress. Consequently, the rock
exhibits higher stiffness and elevated elastic moduli.

In organic-rich mudrock reservoirsâ€”characterized by extremely low permeability and complex
pore-fluid distributionsâ€”the transition between relaxed and unrelaxed regimes becomes both
frequency-dependent and anisotropic. The corresponding transition frequency is governed by the
interplay between rock permeability and fluid viscosity \citep{santos2019}.

In this section, we investigate the effects of attenuation and dispersion induced by the presence of
gas, oil, and water in the pore space, including the case of patchy gasâ€“oil saturation. The effective
pressure employed in the numerical simulationsâ€”17.23 MPaâ€”is close to the in-situ value reported
by \cite{camus2024} for the well from which the core sample analyzed in this study was
extracted.

Figures 4 and 5 present the phase velocity and attenuation behavior of qP wave waves traveling 
normal to the bedding, denoted as '33' waves' ( p$_{33}$ mode)
under four fluid saturation scenarios: 100 percent  oil, 100 percent  gas, 10 percent gas plus 90 percent  oil, and a ternary fluid
mixture consisting of 40 percent gas, 52 percent oil, and 8 percent  water, obtained from the well's fluid saturation log
at the depth corresponding to the core sample.


Figure 4 shows that for '33' waves, the fully oil-saturated medium exhibits the highest phase
velocities, while the fully gas-saturated case yields the lowest. Partially saturated scenarios result in
intermediate velocities, with values decreasing progressively as gas content increases. Notably, even
small amounts of gas (e.g., 10 percent) significantly reduce wave velocity, underscoring the strong
sensitivity of the effective stiffness to the compressibility of the saturating fluid.
A key observation from the plot is that, although the velocities of the partially saturated casesâ€”10 percent 
gas with 90 percent  oil and the ternary mixture 
%(51 percent gas, 39 percent oil, 0.08  percent water)
 appear similar at lower
frequencies, their trajectories begin to diverge markedly above approximately 100 Hz. This
frequency-dependent separation suggests that the dynamic response of the medium becomes
increasingly sensitive to fluid composition at higher frequencies, revealing differences in fluid
mobility and interaction with the solid matrix under wave-induced pore pressure gradients.

Except for the fully oil-saturated case, the results indicate that within the seismic frequency range
(10â€“100 Hz), there are no significant differences in P-wave phase velocity associated with fluid
saturation variations across the volatile oil to dry gas window. However, in the sonic log frequency
range (approximately 1â€“20 kHz), velocity differences begin to emergeâ€”particularly between
partially saturated scenariosâ€”suggesting that the dynamic response becomes progressively
sensitive to fluid composition as frequency increases. It is worth mentioning that the P-wave
velocity observed for the ternary fluid mixture at well-log frequencies is in close agreement (about 13 percent discrepancy)  
with the
velocity measured in the wellâ€™s sonic log, further supporting the validity of the numerical model.
%{\bf  La nota la pongo para recordar que hay que cambiar la figura que estamos haciendo a esas saturaciones
% y el texto correspondiente
 
%Nota: Las saturaciones reportadas en el registro de pozo son Sw=8.5879999999999956E-002,
%Sg=0.51700999999999997, So=0.39711000000000002 y la Vp sonica (p33) es de 3450 m/s.
% FALTA DECIR A QUE FREQUENCIA ES LA VELOCIDAD SONICA}

Figure 5 presents the corresponding attenuation factor  curves (1000/Q) for â€™33â€™ waves as a function of
frequency, where Q denotes the quality factor. The fully oil-saturated case yields the highest
attenuation, with a prominent peak around 1 kHz, attributed to enhanced WIFF resulting from
higher viscosity and lower compressibility contrast. The 10 percent gas plus  90 percent  oil scenario also exhibits
strong attenuation, exceeding that of the ternary fluid mixture. In contrast, the 100 percent gas case shows
the lowest attenuation overall, with a peak shifted toward higher frequencies (100â€“200 kHz),
consistent with minimal fluidâ€“solid interaction and reduced energy dissipation.


Figure 6 examines the phase velocities of qP waves travelling parallel to the bedding,
denoted as '11'  waves (p$_{11}$ mode). As in Figure 4,
the fully oil-saturated case exhibits the highest velocities, while the ternary fluid mixture yields the
lowest across all frequencies. The 10 percent gas + 90 percent oil case follows a similar trend to the oil-
saturated scenario but displays consistently lower velocities throughout the entire frequency range.

The fully gas-saturated scenario remains nearly constant, with only a slight increase observed
between 10 kHz and 100 kHz. This limited sensitivity of the '11' phase velocity to full gas saturation
can be attributed to the direction of wave propagation along the bedding plane, where the effect of
WIFF is not significant. Consequently, the effective stiffness in this direction is practically
independent of frequency.

Figure 7 shows the attenuation curves for '11' waves across the same four saturation cases. The
overall behavior resembles that observed for '33' waves in Figure 5, with oil-saturated and 10 percent gas plus  90 percent  oil scenarios exhibiting the highest attenuation, peaking in the low- to mid-frequency range
(approximately 1â€“2 kHz). In contrast, the 100 percent  gas case displays the lowest attenuation overall,
with a broad, subdued peak shifted toward higher frequencies (~100â€“200 kHz). These results
further confirm that WIFF-driven dissipation is enhanced by fluid viscosity and heterogeneity in
compressibility, and is most pronounced when the fluid phases interact dynamically under pressure
gradients induced by seismic waves.


Comparing Figures 5 and 7 highlights a clear anisotropic attenuation behavior. For '33' waves
(vertical propagation), attenuation magnitudes are substantially higher, and the effects of fluid
saturation are more pronouncedâ€”particularly in the presence of gas. This behavior reflects the
strong coupling between the solid frame and pore fluids in the direction normal to the bedding,
where WIFF mechanisms are activated by pressure gradients across layers. These gradients arise
when waves propagate through alternating materials or saturation zones, leading to local fluid
motion, viscous losses, and strong frequency-dependent attenuation.

In contrast, '11' waves (horizontal propagation) exhibit significantly lower attenuation, and the
influence of fluid composition is much less evident. This is consistent with a scenario where the
wave travels parallel to the layers, encountering minimal heterogeneity and inducing negligible
pressure gradients between pores. As a result, WIFF is greatly reduced, and the rock behaves more
elastically, with attenuation controlled mainly by intrinsic frame properties rather than fluid
interaction.



% Figures 4 y 5  se generaron en 

%~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/23_abril_runs_2mm_sample/freq_loop_p33_4periods_dsize_2mm_100oil_100gas10gas_90oil_40_gas_53_oil_7water$ xmgrace p33_oil_nodo33/inv_qp33_fe.freq p33_gas_nodo44/inv_qp33_fe.freq p33_90_oil_wetting_10_gas_non_wetting_ex_biot/inv_qp33_fe.freq p33_water_oil_gas_nodo22/inv_qp33_fe.freq 


 
\section{Patchy saturation}
To consider fractal variations of gas and oil  saturations  we  use the von Karman self-similar
correlation function  for which the spectral density is \citep{frankel86},  

\begin{equation}\label{vonKarman}
S_d(r_x,r_z)= N_0 (1 + R^2 a ^2)^{-(H + E/2)} ,
\end{equation}
where $R=\sqrt{r_x^2+r_z^2}$ is the radial wavenumber, $a$ the correlation length, 
$H$ is a self-similarity Hurst  coefficient ($0 &lt; H &lt; 1$), $N_0$ is a normalization constant and
 $E$ is the Euclidean dimension.
 The von Karman correlation  \eqref{vonKarman}
 describes a self-affine, fractal processes of fractal dimension
$D=E+1-H$. In the experiments we chose E=2, D= 2.4 and $a$ to be much 
smaller than  the domain size (2 mm).

Figure 8 illustrates the binary patchy saturation map used to model patchy fluid distributions, where white
and black regions correspond to fully gas- and oil-saturated zones, respectively. This spatial
heterogeneity introduces local contrasts in fluid properties, particularly in compressibility and
viscosity, which significantly affect wave propagation through mesoscopic-scale WIFF
mechanisms.

Figure 9 displays the phase velocities of qP waves (p$_{33}$ mode) under two fluid distribution scenarios
â€”uniform and patchy saturationâ€”with identical bulk fluid composition (10 percent gas, 90 percent oil). The
uniform case assumes perfect mixing at the pore scale, producing a single-phase effective fluid with
increased compressibility, resulting in lower velocities. In contrast, the patchy case consists of
spatially segregated gas- and oil-filled regions. In this configuration, wave propagation is
preferentially governed by the stiffer oil-saturated zones, leading to systematically higher velocities
across the frequency spectrum. The divergence between the two models becomes particularly
evident above 100 Hz, consistent with mesoscopic WIFF activation.

Figure 10 displays the attenuation factor for '33' waves under both uniform and patchy saturation
conditions. Patchy saturation results in markedly higher attenuation due to enhanced WIFF at gasâ€“
oil interfaces. These effects promote localized pressure gradients and viscous fluid movement,
leading to increased energy dissipation. Additionally, the attenuation peak is shifted to lower
frequencies under patchy saturation, reflecting the longer diffusion times associated with larger-
scale heterogeneities and limited pressure communication between fluid patches.

Figure 11 compares the induced fluid pressure fields for the uniform (Figure 11 (a)) and patchy 
(Figure 11 (b))
saturation scenarios. In the uniform case, the pressure distribution is symmetric and layered, with
periodic horizontal bands corresponding to the propagation of compressional waves through a
homogeneous fluid mixture. The absence of spatial heterogeneity results in smooth pressure
gradients and relatively low pressure peak values (up to ~0.25 Pa), indicating limited WIFF.

In contrast, the patchy saturation case exhibits a markedly irregular pressure field, with high-
pressure concentrations developing at the interfaces between gas- and oil-filled regions. The
complex geometry of the fluid domains leads to heterogeneous pressure diffusion, higher amplitude
variations (up to ~0.45 Pa), and stronger WIFF effects. These results underscore the essential role of
fluid distribution geometry in modulating pore pressure responses and energy dissipation
mechanisms during wave propagation.




%Using a threshold values $S_o^*$,for each 
%cell with oil  saturation below or  above $S_b^*$ 
%we assign to that cell  either   gas or  oil saturation.  In this fashion we obtained 
%patchy gas-oil spatial distribution for 10 \% gas saturation displayed in  %Figure 8. 

\vskip20cm
% Estas figuras patchy se hicieron con  2 periods en
%~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/23_abril_runs_2mm_sample/%patchy_saturation/100x100_fractal_runs_2_periods/freq_loop_p33_patchy_10_gas_version1


\section*{Conclusions}

This study presented a physically grounded, non-phenomenological methodology for estimating the
effective anisotropic behavior of organic-rich mudrock reservoirs. The simulations were performed
on a two-dimensional, finely layered synthetic medium representative of a core sample from the
Vaca Muerta Formation. The procedure  yielded frequency-dependent stiffness tensors that showed good
agreement with laboratory-measured phase velocities, capturing the strong VTI anisotropy of the
formation.

Importantly, the model does not rely on empirical calibration but is derived from first principles,
integrating mineralogical, petrophysical, and mechanical properties in a physically consistent
manner. The inclusion of both patchy and uniformly mixed fluid saturation scenarios further
emphasized the critical influence of mesoscale fluid heterogeneities on seismic wave dispersion and
attenuation. The results obtained show excellent agreement with well-log data, further validating the
model's applicability under in-situ reservoir conditions.

By explicitly integrating the detailed multimineral composition of the solid matrix with realistic
fluid distributions, this work highlights the value of high-fidelity, physics-based numerical models
for investigating and understanding the poroelastic seismic behavior of unconventional reservoirsâ€”
particularly in the Vaca Muerta Formation, where lithological heterogeneity, spatial variations in the
stress field, and pore pressure fluctuations play a key role in the rockâ€™s mechanical response. 

This study reinforces the importance of explicitly accounting for anisotropic, frequency-dependent
poroelastic effects in organic-rich formations, demonstrating that rock physics models grounded in
first-principles can achieve remarkable agreement with field and laboratory dataâ€”without relying
on empirical corrections.

Future research will extend this framework to three-dimensional geometries and incorporate viscoelastic
and thermo-poroelastic couplings, which are critical for modeling reservoir evolution during
production.

 
\begin{thebibliography}{}
\itemsep0pt

\bibitem[\protect\citename{Al Duhailan, } 2014]{alduha2014} 
Al Duhailan, M. Petroleum-Expulsion Fracturing in Organic-Rich Shales: Genesis and 
Impact on Unconventional Pervasive Petroleum Systems. Doctor of Philosophyâ€™s thesis,
Colorado School of Mines, Golden, Colorado, 227 p. (2014). https://hdl.handle.net/
11124/17003.


\bibitem[\protect\citename{Allan et al., } 2016]{allan2016} 
Allan, A. M. Clark, A. C. Vanorio, T. Kanitpanyacharoen, W. and Wenk. H. R. On the
evolution of the elastic properties of organic-rich shale upon pyrolysis-induced
thermal maturation. Geophysics 2016; 81 (3): D263â€“D281. (2016). doi: https://doi.org
10.1190 /geo2015-0514.1

\bibitem[\protect\citename{Biot, }1962]{biot62} 
Biot, M. A.
Mechanics of deformation and acoustic propagation in porous media,
J. Appl. Phys. 33, 1482 -- 1498 (1962).


\bibitem[\protect\citename{Carcione et al, }2005] {carcione2005} 
Carcione, J.  M.,  Helle, H. B., Santos, J. E. and Ravazzoli, C. L. A constitutive equation and generalized Gassmann modulus for multimineral porous media. Geophysics, 70 (2), N17 â€“ N26 (2005). doi10.1190/1.1897035


\bibitem[\protect\citename{Frankel and Clayton, }1986]{frankel86}
Frankel, A. and  Clayton, R.~W. Finite difference simulation of seismic wave
  scattering: implications for the propagation of short period seismic waves in
  the crust and models of crustal heterogeneity, J. Geophys. Res. 91 
  6465 -- 6489 (1986).
  
 \bibitem[\protect\citename{Hashin and Shtrikman, }1963]{HS63}
 Hashin, Z. and Shtrikman, S. A variational approach to the theory of the elastic behavior 
 of multiphase materials. J. Mech. Phys. Solids, 11, 127-140 (1963).


\bibitem[\protect\citename{Krief et al., }1990]{krief_90}
Krief, M.,  Garat, J.,   Stellingwerff, J.  and Ventre, J.  
A petrophysical interpretation using the velocities 
of P and S waves (full waveform sonic), The Log Analyst 31, 355-369 (1990).


\bibitem[\protect\citename{Krzikalla  and M\"uller, }2011]{KM_11} 
Krzikalla, F.  and M\"uller, T. M.  
Anisotropic P-SV-wave dispersion and attenuation due to interlayer 
flow in thinly layered porous rocks, Geophysics 76 W 135 (2011); doi:10.1190/1.3555077.


\bibitem[\protect\citename{Loucks et al., } 2012]{loucks2012}
Loucks, R. G. Reed, R. M. Ruppel, S. C. and Hammes, U. Spectrum of pore types and networks in mudrocks and a descriptive classification for matrix-related mudrock pores: American Association of Petroleum Geologists Bulletin, v. 96, p. 1071â€“1098 (2012). doi: 10.1306/08171111061.

\bibitem[\protect\citename{Mavko et al, }2020]{mavko2020}
Mavko,G., Mujerki, T. and Dvorkin, H. The  Rock physics Handbook,  Cambridge University Press, Third Edition (2020).

\bibitem[\protect\citename{Mokhtari and Tutuncu, } 2015]{mokhtari2015}
Mokhtari, M. and Tutuncu A. N.  Characterization of anisotropy in the permeability of organic-rich shales. Journal of Petroleum Science and Engineering, Volume 133, 2015, Pages 496-506, ISSN 0920-4105 (2015). https://doi.org/10.1016/j.petrol.2015.05.024.

\bibitem[\protect\citename{Nelson, } 2009]{nelson2009} 
Nelson, P. H.  Pore-throat sizes in sandstones, tight sandstones, and shales. AAPG Bulletin 2009; 93 (3): 329â€“340  (2009). doi:https://doi.org/10.1306/10240808059.

\bibitem[\protect\citename{Pinna et al., }2011]{pinna2011} 
Pinna, G. Carcione, J. M. and Poletto, F. Kerogen to oil conversion in source rocks: Pore-pressure build-up and effects on seismic velocities. Journal of Applied Geophysics, 74, 229â€“235 (2011). doi: 10.1016/j.jappgeo.2011.05.006.

\bibitem[\protect\citename{Sanchez Camus et al., }2024]{camus2024}
S\'anchez Camus, A., Gauzellino M.P., and Ramos R.  Application of the Vernik method for
pore-pressure estimation in the Vaca Muerta Formation: Insights for geomechanical analysis and
quantitative seismic interpretation in organic-rich mudrocks. The Leading Edge 43: 774â€“781.
https://doi.org/10.1190/tle43110774.1.


\bibitem[\protect\citename{Santos and Carcione, }2015] {santos2015} 
Santos, J. E. and  Carcione, J. M. Finite-element harmonic experiments to model fractured induced
anisotropy in poroelastic media, Comp. Methods Appl. Mech. Engr. 283 (2015). 


\bibitem[\protect\citename{Santos et al., }2019]{santos2019} 
Santos, J. E., Savioli, G.B., Carcione, J. M.  and Ba, J. 
Effect of capillarity and relative permeability on Q anisotropy of
hydrocarbon source rocks, Geophys. J. Int.  218   1199-1209 (2019). 


\bibitem[\protect\citename{Sone and  Zoback, }2013]{sone13}
Sone, H  and Mark D. Zoback, M. D. Mechanical properties of shale-gas reservoir rocks â€” 
Part 1: Static and dynamic elastic properties and anisotropy, Geophysics  78 (5) 
https://doi.org/10.1190/geo2013-0050.1


\bibitem[\protect\citename{Vernik and Liu, } 1997]{vernik1997}
Vernik, L. and Liu, X.  Velocity anisotropy in shales: A petrophysical study.
Geophysics, 62, no. 2, 521â€“532 (1997). doi: 10.1190/1.1444162.
  
\bibitem[\protect\citename{Vernik  and Milovac, }2011]{vernik2011}
Vernik, L.  and Milovac, J. Rock physics of organic shales, 30 (3) 
318-323 (2011).

\bibitem[\protect\citename{Vernik  and Nur, }1992] {vernik92}
Vernik, L. and Nur, A. 1992, Ultrasonic velocity and anisotropy
of hydrocarbon source rocks, Geophysics  57 (5) 727â€“735 (1992),
doi:10.1190/1.1443286.


\bibitem[\protect\citename{White et al., }1975]{white75} 
White, J. E.   Mikhaylova, N. G.  and Lyakhovitskiy, F. M.   
Low-frequency seismic waves in fluid saturated layered rocks, 
Physics of the Solid Earth  11  654-659 (1975).

\bibitem[\protect\citename{Yang, } 2019]{yang2019} 
Yang Y. Rocks micro-structural and elastic changes during thermal maturation. Doctor of
Philosophy thesis, Department of Geophysics, Stanford University (2019). 

\bibitem[\protect\citename{Zoback and Kohli, } 2019]{zoback2019}  
Zoback, M., and Kohli, A. 2019. Unconventional Reservoir Geomechanics: Shale Gas, Tight Oil, and Induced Seismicity (2019). doi: 10.1017/9781316091869.


\end{thebibliography}  

\newpage

%\begin{center}
%{\bf Table 1. Properties of the sandstone}    
%\[  
%\begin{tabular}{cc}
%\hline
%\hline
%Grain  bulk modulus, $K_{s}$   &amp;    33.4 GPa \\
%\  density, \ \  $\rho_{s}$   &amp;   2650 \ kg/m$^3$ \\   
%\hline
%Dry-matrix bulk modulus, $K_{m}$ \ \   &amp; 1.3 GPa \\
% \ \ \ \ \ \ \ \ \ \ \ \ \ \         shear modulus, $\mu$  \ \  &amp; 1.4 \ GPa \\
%\ \ \ \ \ \  porosity, $\phi$   \ \ &amp;       0.3  \\  
%\ \ \ \ \ \ \ \ \ \ \ permeability,  $\kappa$  \ \ &amp;   10$^{-12}$ \ m$^2$  \\
%\hline
%\hline
%\end{tabular}
%\]
%\end{center}


\section*{TABLES}

%\begin{center}

\begin{center}
{\bf Table 1. Fluid properties.}
%\end{center}
\[
\begin{tabular}{lccccc}
\hline
\hline
 Fluid &amp;  Bulk modulus (Pa) &amp; Density (kg/m$^3$) &amp; Viscosity (Pa . s)  \\
\hline 
Air  &amp;  1.01325 10$^5$  &amp; 1.225            &amp;  1.805 10$^-5$  \\
%
Water  &amp;   2.25 10$^{9}$   &amp; 1000   &amp;  0.001  \\
Oil  &amp; 0.57 10$^9$    &amp; 700   &amp;  0.01   \\
\end{tabular}
\]
\end{center}


\vspace{0.4cm}

\begin{center}
{\bf Table 2. Mineralogical properties obtained from X-ray diffraction (XRD) analysis of core sample,
used for Material 1 at an effective pressure of 17.23 MPa }
%\end{center}
\[
%\begin{tabular}{cccccc}
\begin{tabular}{lccccc}
\hline
\hline
Mineral  &amp; Bulk modulus (GPa) &amp; Shear modulus (GPa) &amp; Density (kg/m$^3$) &amp; Proportions (\%) \\
\hline 
Kerogen  &amp;     7.0   &amp; 2.         &amp;  1400  &amp; 23  \\
Clay     &amp;     25    &amp; 9          &amp;  2700  &amp; 0.3727 \\
Quartz   &amp;     45    &amp; 55         &amp;  2700  &amp; 0.1461\\
Calcite  &amp;     80    &amp; 40         &amp;  2800  &amp; 0.1068\\
Plagioclase&amp;   80    &amp; 40         &amp; 2800   &amp; 0.0257 \\
Dolomite  &amp;   100    &amp; 50         &amp; 2900   &amp; 0.0237 \\
Pyrite   &amp;    170    &amp; 110        &amp; 5000   &amp;0.035 \\
\end{tabular}
\]
\end{center}

\vspace{0.4cm}

\begin{center}
{\bf Table 3. Phase velocities computed, measured and error  percentage.}
\[
\begin{tabular}{lccccc}
\hline
\hline
Phase velocity vp (m/s)&amp; Computed  &amp; Measured&amp; Percentage error \\
\hline 
v11  &amp; 4644.378    &amp; 4331       &amp;  7.2  \% \\
v33  &amp; 3804.510    &amp; 4217.47    &amp;  9.8 \% \\
v55  &amp; 1974.760    &amp; 2193.61    &amp;  9.9   \% \\
v66  &amp; 2742.641    &amp; 2581 &amp;     6.2 \%   \\
\end{tabular}
\]
\end{center}


\vspace{0.4cm}

\begin{center}
{\bf Table 4.  Numerical pij values (Pa) from the VTI experiments. Frequency is 1 MHz, 
%{\bf valores definitivos  para 4 periodos}
}
%\end{center}
%\label{table1}       
%---------
\[
\begin{tabular}{lccccc}
\hline
\hline
p11 = 
%\hline 
(49981147009.91, 1315.04)\\
\hline 
 p33 = 
(33538888987.67, 11716.97) \\
\hline
 p55 = 
(9036079273.1829395, 0.)\\
\hline
 p66 = 
(17429665052.36, 504048.90) \\
 \hline
p13 =  (11164071155.78, 3909.98)\\
\hline
\end{tabular}
\]

\end{center}

\newpage
\section*{LIST OF FIGURES}


\noindent
Figure 1. Minerals and its proportions in the well from where the core sample was obtained.

\noindent
Figure 2. Polar representation  of  phase  velocities of qP,
qSV  and SH waves  computed using the FE method. Frequency is 1 MHz.

\noindent
Figure 3. Polar representation  of energy velocities of qP,  qSV and SH 
waves  computed using the FE method. Frequency is 1 MHz.

\noindent
Figure 4. '33'  waves phase velocity as function of frequency.

{\fo
\noindent
Figure 5. '33' waves  attenuation factor as function of frequency.

\noindent
Figure 6. '11'  waves phase velocity as function of frequency.

\noindent
Figure 7. '11' waves  attenuation factor  as function of frequency. 

 \noindent
Figure 8. Fractal patchy gas-oil spatial distribution. White zones correspond to full gas saturation and black regions to full oil saturation, with overall gas saturation 10 percent. Fractal dimension is 2.4, correlation length is 0.09  mm. 

\noindent
Figure 9. Phase velocity of '33' waves as function of frequency for  10 \% gas, 90 \% oil and uniform and patchy saturation.

\noindent
Figure 10. Attenuation factor of '33' waves  as function of frequency for  10 \% gas, 90 \% oil and uniform and patchy saturation.

\noindent
Figure 11. Fluid pressure  for the case of 10 \% gas, 90 \% oil  and uniform (a) and patchy (b)  saturation.

}


 

\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=16cm,height=14cm,angle=0]{figure1.eps}
\caption {Minerals and its proportions in the well from where the core sample was obtained.
}
% \vspace{0.3cm}

\label{Fig1} 
\end{figure}


\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=-90]{figure2.eps}
\caption {Polar representation  of  phase  velocities of qP,
qSV  and SH waves  computed using the FE method. Frequency is 1 MHz.
}
% \vspace{0.3cm}
\label{Fig2} 
\end{figure}
\vskip20cm

\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=-90]{figure3.eps}
\caption {Polar representation  of energy velocities of qP,  qSV and SH 
waves  computed using the FE method. Frequency is 1 MHz.  
}
% \vspace{0.3cm}
\label{Fig3} 
\end{figure}
\vskip20cm

 \begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure4.eps}
\caption {'33'  waves phase velocity as function of frequency.} 
% \vspace{0.3cm}
\label{Fig4} 
\end{figure}
\vskip20cm



\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure5.eps}
\caption {'33' waves  attenuation factor as function of frequency. 
}
% \vspace{0.3cm}
\label{Fig5} 
\end{figure}
\vskip20cm

% Figuras 6 y 7 se hicieron en 
%~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/23_abril_runs_2mm_sample/%freq_loop_p11_dsize_2mm_100oil_100gas10gas_90oil_40_gas_53_oil_7water

\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure6.eps}
\caption {'11'  waves phase velocity as function of frequency. 
} 

% \vspace{0.3cm}
\label{Fig6} 
\end{figure}
\vskip20cm



\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure7.eps}
\caption {'11' waves  attenuation factor  as function of frequency. 
}
% \vspace{0.3cm}
\label{Fig7} 
\end{figure}
\vskip20cm



\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure8.eps}
\caption {Fractal patchy gas-oil spatial distribution. White zones correspond to full gas saturation and black regions to full oil saturation, with overall gas saturation 10 percent. Fractal dimension is 2.4, 
correlation length is 0.09  mm. 
}
% \vspace{0.3cm}
\label{Fig8} 
\end{figure}
\vskip20cm 

%Las figuras patchy se hicieron con este path   

%~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/23_abril_runs_2mm_sample/freq_loop_p33_4periods_dsize_2mm_100oil_100gas10gas_90oil_40_gas_53_oil_7water/p33_90_oil_wetting_10_gas_non_wetting_ex_biot$ xmgrace inv_qp33_fe.freq ~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/patchy_saturation/2025_patchy_saturation_runs/200x200_fractal_runs_4_periods/200x200_fractal_runs_4_periods/freq_loop_p33_patchy_10_gas_version1/inv_qp33_fe.freq



\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure9.eps}
\caption {Phase velocity of '33' waves as function of frequency for  10 \% gas, 90 \% oil 
and uniform and patchy saturation.
}
% \vspace{0.3cm}
\label{Fig9} 
\end{figure}
\vskip20cm


%Las figuras patchy se hicieron con este path
%~/santos2024/papers/paper_con_camus/first_draft_paper_april_9/23_abril_runs_2mm_sample/freq_loop_p33_dsize_2mm_100oil_100gas10gas_90oil_40_gas_53_oil_7water/90_oil_wetting_10_gas_non_wetting xmgrace vp33_fe.freq  ../../patchy_saturation/100x100_fractal_runs_2_periods/freq_loop_p33_patchy_10_gas_version1/vp33_fe.freq

\begin{figure} 
\vspace{0.5cm}
\includegraphics[width=12cm,angle=0]{figure10.eps}
\caption {Attenuation factor of '33' waves  as function of frequency for  10 \% gas, 90 \% oil 
and uniform and patchy saturation.
} 
% \vspace{0.3cm}
\label{Fig10} 
\end{figure}
\vskip20cm


\begin{figure} 
%\vspace{0.5cm}
(a) \includegraphics[width=8cm,angle=0]{figure11a.eps}
(b) \includegraphics[width=8cm,angle=0]{figure11b.eps}
\caption {Fluid pressure  for the case of 10 \% gas, 90 \% oil  and uniform (a) and patchy (b)  saturation.} 
% \vspace{0.3cm}
\label{Fig11} 
\end{figure}
\vskip20cm




\end{document}

\section{Acknowledgements}
The authors wish to thank the Centro de SimulaciÃ³n Computacional para Aplicaciones TecnolÃ³gicas (CSC) CONICET, Argentina and the Research Center for Advanced Computation (RCAC) Purdue University for useful support and collaboration. This research was partially supported by Universidad de Buenos Aires (UBACyT 20020190100236BA).


EJEMPLOS DE TABLES en un paper nuestro anterior

\begin{center}
{\bf Table 1. Properties of the sandstone}    
\[  
\begin{tabular}{cc}
\hline
\hline
Grain  bulk modulus, $K_{s}$   &amp;    33.4 GPa \\
\  density, \ \  $\rho_{s}$   &amp;   2650 \ kg/m$^3$ \\   
\hline
Dry-matrix bulk modulus, $K_{m}$ \ \   &amp; 1.3 GPa \\
 \ \ \ \ \ \ \ \ \ \ \ \ \ \         shear modulus, $\mu$  \ \  &amp; 1.4 \ GPa \\
\ \ \ \ \ \  porosity, $\phi$   \ \ &amp;       0.3  \\  
\ \ \ \ \ \ \ \ \ \ \ permeability,  $\kappa$  \ \ &amp;   10$^{-12}$ \ m$^2$  \\
\hline
\hline
\end{tabular}
\]
\end{center}

\vspace{0.4cm}

\begin{center}
{\bf Table 2. Properties of the saturant fluids}
\[
\begin{tabular}{cc}
\hline
\hline
Brine  bulk modulus, $K_w$   &amp;  2.2  GPa  \\ 
  density, $\rho_w$   &amp;  975  kg/cm$^3$ \\
\   viscosity,  $\eta_w$ &amp;   0.001 Pa $\cdot$ s  \\
\hline
Oil bulk modulus, $K_o$   &amp;  2  GPa  \\ 
 density, $\rho_o$   &amp;      870  kg/cm$^3$ \\
	\ \   viscosity, $\eta_o$ &amp;         0.3 Pa $\cdot$ s  \\		
\hline
Gas  bulk modulus, $K_g$    &amp;   0.0044515 GPa  \\ 
density,  $\rho_g$   &amp;   42.316  kg/m$^3$ \\
\ \ viscosity, \ $\eta_g$ &amp;     1.1186 $\times 10^{-5}$ Pa $\cdot$ s  \\		
\hline
\hline
\end{tabular}
\]
\end{center}





\bibitem[Adams and Fournier, 2003]{Adams2003}
 Adams, R. A. and Fournier, J. F., 2003, 
 Sobolev Spaces, 2nd ed., Elsevier, UK. 
  

  \bibitem[Armstrong, 1984]{Armstrong1984}
 Armstrong, B. H., 1984, Models for thermoelastic attenuation of waves in
  heterogeneous solids: Geophysics, {\bf 49}, 849--1121,
  \url{https://doi.org/10.1190/1.1441718}.

  
\bibitem[Biot, 1956]{Biot1956}
 Biot, M. A., 1956, Thermoelasticity and irreversible thermodynamics:
  Journal of Applied Physics, {\bf 27}, 240--253,
  \url{https://doi.org/10.1063/1.1722351}.
 
 
\bibitem[Carcione, 2022]{Carcione2022}
 Carcione, J. M., 2022,  Wave Fields in Real Media. Theory and numerical
  simulation of wave propagation in anisotropic, anelastic, porous and
  electromagnetic media,  4th ed., revised and extended: Elsevier, Amsterdam. 
  
 
   
\bibitem[Carcione et al., 2019a]{Carcione2019a}
Carcione, J.M., F. Cavallini, E. Wang, J. Ba, and L. Y. Fu, 2019a, Physics
  and simulation of wave propagation in linear thermoporoelastic media:
  Journal of Geophysical Research: Solid Earth, {\bf 124}, 8147--8166,
  \url{https://doi.org/10.1029/2019JB017851}.
  
   
\bibitem[Carcione et al., 2019b]{Carcione2019b}
Carcione, J. M.,  Z.W. Wang, W. Ling, E. Salusti, J. Ba, and L.Y. Fu, 2019b, 
   Simulation of wave propagation in linear thermoelastic media:
  Geophysics, {\bf 84}, T1--T11,
  \url{https://doi.org/10.1190/geo2018-0448.1}.
  
  
\bibitem[Carcione et al., 2020]{Carcione2020}
Carcione, J. M., D. Gei, J. E. Santos, L. Y. Fu, and J. Ba, 2020, Canonical
  analytical solutions of wave-induced thermoelastic attenuation: Geophysical
  Journal International, {\bf 221}, 835--842,
  \url{https://doi.org/10.1093/gji/ggaa033}.
  
  
\bibitem[Duvaut and Lions, 1976]{Duvaut1976}
Duvaut, G. and J. L. Lions, 1976, Inequalities in mechanics and physics,
  Springer-Verlag, Berlin, 1976. 

 
\bibitem[Girault and Raviart, 1981]{Girault1981}
Girault, V. and P. A. Raviart, 1981, Finite element approximation of the
  Navier Stokes equations, Springer-Verlag, Berlin.

\bibitem[Lifshitz et al., 2000]{Lifshitz2000}
 Lifshitz, R., and M. L. Roukes, 2000, Thermoelastic damping in micro- and
  nanomechanical systems: Physical Review B, {\bf 61},
  \url{https://doi.org/10.1103/PhysRevB.61.5600}.

  
\bibitem[Lord and Shulman, 1967]{Lord1967}
Lord, H., and Y. S. Shulman, 1967, A generalized dynamical theory of
  thermoelasticity: Journal of the Mechanics and Physics of Solids, {\bf 15},
  299--309, \url{https://doi.org/10.1016/0022-5096(67)90024-5}.

 
\bibitem[Nedelec, 1980]{Nedelec1980}
Nedelec, J. C., 1980, Mixed finite elements in ${R}^3$: Numerische
  Mathematik, {\bf 35}, 315--341,
  \url{https://doi.org/10.1007/BF01396415}.
    
\bibitem[Raviart and Thomas, 1975]{Raviart1975}
 Raviart, P. A. and J. M. Thomas, 1975, Mixed finite element method for
  2$^{nd}$ order elliptic problems: Mathematical Aspects of the Finite Element
  Methods, Lecture Notes of Mathematics, Springer, {\bf 606}, 292--315,
  \url{https://doi.org/10.1007/BFb0064470}.
  
  
\bibitem[Rudgers, 1990]{Rudgers1990}
Rudgers,  A. J., 1990, Analysis of thermoacoustic wave propagation in
  elastic media: The Journal of the Acoustical Society of America, {\bf 88}, 
  \url{https://doi.org/10.1121/1.399856}.


\bibitem[Santos et al., 2021]{Santos2021}
 Santos, J. E., J. M. Carcione, and J. Ba, 2021, Existence and uniqueness of
  solutions of thermo-poroelasticity: Journal of Mathematical Analysis and
  Applications, {\bf 499}, \url{https://doi.org/10.1016/j.jmaa.2020.124907}.
  
\bibitem[Santos et al., 1988]{santos88}
 Santos, J. E., J. Douglas, Jr., M. E. Morley, and O. M. Lovera, 1988,  
 Finite element methods for a model for full waveform acoustic logging; IMA Journal of Numerical Analysis,
 {\bf 8}, 415--433.
 
 
\bibitem[Savage, 1966]{Savage1966}
 Savage, J. C., 1966, Thermoelastic attenuation of elastic waves by cracks:
  Journal of Geophysical Research (1896-1977), Correction: 1967, JGR, 72,
  6387.,  {\bf 71}, 3929--3938,
  \url{https://doi.org/10.1029/JZ071i016p03929}.
  

\bibitem[Sharma, 2008]{Sharma2008}
Sharma, M. D., 2008, Wave propagation in thermoelastic saturated porous
  medium: Journal of Earth System Science, {\bf 117}, 951--958,
  \url{https://doi.org/10.1007/s12040-008-0080-4}.


\bibitem[Treitel, 1959]{Treitel1959}
Treitel, S, 1959, On the attenuation of small-amplitude plane stress waves
  in a thermoelastic solid: Journal of Geophysical Research (1896-1977), {\bf 64}, 
  661--665, \url{https://doi.org/10.1029/JZ064i006p00661}.  
  
\bibitem[Wang et al., 2020]{Wang2020}
Wang, Z.W.,  L. Y. Fu, J. Wei, W. Hou, J. Ba, and J. M. Carcione, 2020, On
  the green function of the Lordâ€“Shulman thermoelasticity equations;
  Geophysical Journal International, {\bf 220}, 393--403,
  \url{https://doi.org/10.1093/gji/ggz453}.

\bibitem[Wang et al., 2021]{Wang2021}
Wang, E., J. M. Carcione, F. Cavallini, M. Botelho, and J. Ba, 2021, 
  Generalized thermo-poroelasticity equations and wave simulation; Surveys in
  Geophysics, {\bf 42}, 133--157,
  \url{https://doi.org/10.1007/s10712-020-09619-z}.
  
  
\bibitem[Wei et al., 2020]{Wei2020}
Wei, J., L. Y. Fu, Z. W. Wang, J. Ba, and J. M. Carcione, 2020, Green
  function of the Lord-Shulman thermo-poroelasticity theory: Geophysical
  Journal International, {\bf 221}, 1765--1776,
  \url{https://doi.org/10.1093/gji/ggaa100}.

\bibitem[White et al., 1975]{White1975}
White, J.E., N. Mikhaylova, and F. M. Lyakhovitskiy, 1975, Low-frequency
  seismic waves in fluid saturated layered rocks: Izvestija Acad. Sci. USSR,
  Phys. Solid Earth, {\bf 11}, 654--659.

  
\bibitem[Zener, 1938]{Zener1938}
Zener, C., 1938, Internal friction in solids ii. general theory of
  thermoelastic internal friction: Physical Review, {\bf 53}, 90--99.
 
\bibitem[Zener, 1946]{Zener46}
Zener, C., 1946, Anelasticity of Metals, American Institute of Mining and
Metallurgical Engineers, Pub. no. 1992.
  
</pre></body></html>