Heterogeneity of rock's fabric can induce heterogeneous distribution of immiscible fluids in natural reservoirs, since the lithological variations (mainly permeability) may affect fluid migration in geological time scales, resulting in patchy saturation of fluids. Therefore, fabric and saturation inhomogeneities both affect wave propagation. To model the wave effects (attenuation and velocity dispersion), we introduce a double double‐porosity model, where pores saturated with two different fluids overlap with pores having dissimilar compressibilities. The governing equations are derived by using Hamilton's principle based on the potential energy, kinetic energy, and dissipation functions, and the stiffness coefficients are determined by gedanken experiments, yielding one fast P wave and four slow Biot waves. Three examples are given, namely, muddy siltstones, clean dolomites, and tight sandstones, where fabric heterogeneities at three different spatial scales are analyzed in comparison with experimental data. In muddy siltstones, where intrapore clay and intergranular pores constitute a submicroscopic double‐porosity structure, wave anelasticity mainly occurs in the frequency range (10⁴–10⁷ Hz), while in pure dolomites with microscopic heterogeneity of grain contacts and tight sandstones with mesoscopic heterogeneity of less consolidated sands, it occurs at 10³–10⁷ Hz and 10¹–10³ Hz (seismic band), respectively. The predicted maximum quality factor of the fast compressional wave for the sandstone is the lowest (approximately 8), and that of the dolomite is the highest. The results of the diffusive slow waves are affected by the strong friction effects between solids and fluids. The model describes wave propagation in patchy‐saturated rocks with fabric heterogeneity at different scales, and the relevant theoretical predictions agree well with the experimental data in fully and partially saturated rocks.

The standard Biot‐Gassmann theory of poroelasticity fails to explain strong compressional wave velocity dispersion experimentally observed in 12 tight siltstone with clay‐filled pores. In order to analyze and understand the results, we developed a new double‐porosity model of clay squirt flow where wave‐induced local fluid flow occurs between the micropores in clay aggregates and intergranular macropores. The model is validated based on the combined study of ultrasonic experiments on specimens at different saturation conditions and theoretical predictions. The presence of a sub‐pore‐scale structure of clay micropores contained in intergranular macropores, where the fluid does not have enough time to achieve mechanical equilibrium at ultrasonic frequencies and thus stiffens the rock matrix, provides a suitable explanation of the experimental data. Moreover, the model provides a new bound for estimating the compressional wave velocity of tight rocks saturated with two immiscible liquids. The theoretical predictions indicate that the velocity variation between gas‐ and liquid‐saturated specimens is predominantly induced by the clay squirt stiffening effect on the rock matrix and not by fluid substitution. The effect contributes more than 90% to the variation in the porosity range of 0–5%. Thus, clay squirt flow dominates the relationships between compressional wave velocity and pore fluid in tight rocks.

We develop a numerical algorithm for simulation of wave propagation in linear nonisothermal poroelastic media, based on Biot theory and a generalized Fourier law of heat transport in analogy with Maxwell model of viscoelasticity. A plane wave analysis indicates the presence of the classical P and S waves and two slow waves, namely, the Biot and the thermal slow modes of propagation, which present diffusive behavior under certain conditions, depending on viscosity, frequency, and the thermoelastic constants. The wavefield is computed with a direct meshing method using the Fourier differential operator to calculate the spatial derivatives. We propose two alternative time‐stepping algorithms, namely, a first‐order explicit Crank‐Nicolson method and a second‐order splitting method. The Fourier differential operator provides spectral accuracy in the calculation of the spatial derivatives. Modeling the thermal diffusive mode is relevant for high‐temperature high‐pressure fields and since it leads to mesoscopic attenuation by mode conversion of the fast waves to the thermal waves.

The fractal texture (or fabric) of porous media, which supports fluid flow at different scales, is the main cause of wave anelasticity (dispersion and attenuation) on a wide range of frequencies. To model this phenomenon, we develop a theory of wave propagation in a fluid saturated infinituple‐porosity media containing inclusions at multiple scales, based on the differential effective medium (DEM) theory of solid composites and Biot‐Rayleigh theory for double‐porosity media. The dynamical equations are derived from first principles, that is, based on the strain (potential), kinetic, and dissipation energies, leading to generalized stiffness and density coefficients. The scale of the inclusions can be characterized by different distributions. The theory shows that the anelasticity depends on the size (radius) of the inclusions, parameter θ (exponential distribution), mean radius r₀ and variance σ_r² (Gaussian distribution) and the fractal dimension D_f (self‐similar distribution). When D_f = 2, θ = 1 and σ_r² = 4, the three distributions give the same P‐wave velocities and attenuation, since each added inclusion phase has nearly the same volume fraction. For the modeling results, the range of anelasticity of D_f = 2.99/θ = 1/σ_r² = 4 is broader than that of D_f < 2.99/θ < 1/σ_r² < 4. To confirm the validity of the model, we compare the results with laboratory measurements on tight sandstone and carbonate samples in the range 1 Hz–1 kHz, Fox Hill sandstone (5 Hz–800 kHz) and field measurements of marine sediments (50 Hz–400 kHz). This comparison shows that the model successfully describes the observed anelasticity.

We have developed a numerical algorithm for simulation of wave propagation in linear thermoelastic media, based on a generalized Fourier law of heat transport in analogy with a Maxwell model of viscoelasticity. The wavefield is computed by using a grid method based on the Fourier differential operator and two time-integration algorithms to cross-check solutions. Because the presence of a slow quasistatic mode (the thermal mode) makes the differential equations stiff and unstable for explicit time-stepping methods, first, a second-order time-splitting algorithm solves the unstable part analytically and a Runge-Kutta method the regular equations. Alternatively, a first-order explicit Crank-Nicolson algorithm yields more stable solutions for low values of the thermal conductivity. These time-stepping methods are second- and first-order accurate, respectively. The Fourier differential provides spectral accuracy in the calculation of the spatial derivatives. The model predicts three propagation modes, namely, a fast compressional or (elastic) P-wave, a slow thermal P diffusion/wave (the T-wave), having similar characteristics to the fast and slow P-waves of poroelasticity, respectively, and an S-wave. The thermal mode is diffusive for low values of the thermal conductivity and wave-like for high values of this property. Three velocities define the wavefront of the fast P-wave, i.e., the isothermal velocity in the uncoupled case, the adiabatic velocity at low frequencies, and a higher velocity at high frequencies.