%\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\<{\left\langle} \def\>{\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& 0 I_2 \\ 0 I_2 &\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} &\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}\\ &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 < H < 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}$ & 33.4 GPa \\ %\ density, \ \ $\rho_{s}$ & 2650 \ kg/m$^3$ \\ %\hline %Dry-matrix bulk modulus, $K_{m}$ \ \ & 1.3 GPa \\ % \ \ \ \ \ \ \ \ \ \ \ \ \ \ shear modulus, $\mu$ \ \ & 1.4 \ GPa \\ %\ \ \ \ \ \ porosity, $\phi$ \ \ & 0.3 \\ %\ \ \ \ \ \ \ \ \ \ \ permeability, $\kappa$ \ \ & 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 & Bulk modulus (Pa) & Density (kg/m$^3$) & Viscosity (Pa . s) \\ \hline Air & 1.01325 10$^5$ & 1.225 & 1.805 10$^-5$ \\ % Water & 2.25 10$^{9}$ & 1000 & 0.001 \\ Oil & 0.57 10$^9$ & 700 & 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 & Bulk modulus (GPa) & Shear modulus (GPa) & Density (kg/m$^3$) & Proportions (\%) \\ \hline Kerogen & 7.0 & 2. & 1400 & 23 \\ Clay & 25 & 9 & 2700 & 0.3727 \\ Quartz & 45 & 55 & 2700 & 0.1461\\ Calcite & 80 & 40 & 2800 & 0.1068\\ Plagioclase& 80 & 40 & 2800 & 0.0257 \\ Dolomite & 100 & 50 & 2900 & 0.0237 \\ Pyrite & 170 & 110 & 5000 &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)& Computed & Measured& Percentage error \\ \hline v11 & 4644.378 & 4331 & 7.2 \% \\ v33 & 3804.510 & 4217.47 & 9.8 \% \\ v55 & 1974.760 & 2193.61 & 9.9 \% \\ v66 & 2742.641 & 2581 & 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}$ & 33.4 GPa \\ \ density, \ \ $\rho_{s}$ & 2650 \ kg/m$^3$ \\ \hline Dry-matrix bulk modulus, $K_{m}$ \ \ & 1.3 GPa \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ shear modulus, $\mu$ \ \ & 1.4 \ GPa \\ \ \ \ \ \ \ porosity, $\phi$ \ \ & 0.3 \\ \ \ \ \ \ \ \ \ \ \ \ permeability, $\kappa$ \ \ & 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$ & 2.2 GPa \\ density, $\rho_w$ & 975 kg/cm$^3$ \\ \ viscosity, $\eta_w$ & 0.001 Pa $\cdot$ s \\ \hline Oil bulk modulus, $K_o$ & 2 GPa \\ density, $\rho_o$ & 870 kg/cm$^3$ \\ \ \ viscosity, $\eta_o$ & 0.3 Pa $\cdot$ s \\ \hline Gas bulk modulus, $K_g$ & 0.0044515 GPa \\ density, $\rho_g$ & 42.316 kg/m$^3$ \\ \ \ viscosity, \ $\eta_g$ & 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.