# A Simulation Method for the Computation of the Effective P-Wave Velocity in Heterogeneous Rocks

Angel Javier A. J. Omella, Magdalena Grigoroscuta-Strugaru, Julen Alvarez-Aramberri, Vincent Darrigrand, David Pardo, Hector Gonzalez, Carlos Santos

► To cite this version:

Angel Javier A. J. Omella, Magdalena Grigoroscuta-Strugaru, Julen Alvarez-Aramberri, Vincent Darrigrand, David Pardo, et al.. A Simulation Method for the Computation of the Effective P-Wave Velocity in Heterogeneous Rocks. Computational Mechanics, 2021, 10.1007/s00466-020-01966-3 hal-03092551

HAL Id: hal-03092551

https://hal.science/hal-03092551

Submitted on 2 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Simulation Method for the Computation of the Effective P-Wave Velocity in Heterogeneous Rocks

Ángel Javier Omella · Julen Alvarez-Aramberri · Magdalena Strugaru · Vincent Darrigrand · David Pardo · Héctor González · Carlos Santos

Received: date / Accepted: date

**Abstract** We propose a set of numerical methods for the computation of the frequency-dependent effective primary wave (P-wave) velocity of heterogeneous rocks. We assume the rocks internal micro-structure is provided by given micro computed tomography (micro-CT) images. In the low/medium frequency regime, we propose to solve the acoustic equation in the frequency domain by a Finite Element Method (FEM). We employ a Perfectly Matched Layer (PML) to truncate the computational domain and we show the need to repeat the domain a sufficient number of times to obtain accurate results. To make this problem computationally tractable, we equip the FEM with non-fitting meshes and we precompute multiple blocks of the stiffness matrix. In the high-frequency range, we solve the eikonal equation with the Fast Marching Method (FMM). The results computed via this method match for vertical transversely isotropic media (VTI) with the effective compressional wave velocity calculated by the Wyllie-time-average formula. Extensive numerical

results confirm the validity of the proposed methods and illustrate the effect of density, porosity, and the size and distribution of the pores on the effective compressional wave velocity.

**Keywords** Effective P-wave Velocity · Finite Element Method · Fast Marching Method

**Mathematics Subject Classification (2010)** MSC 86-08 · MSC 65N30 · MSC 74Q05

## 1 Introduction

In geophysical applications, it is of paramount importance to characterize the effective velocity of an elastic wave traveling through the Earth's subsurface. This information enables to identify gas saturated rocks as well as the porosity and other key petrophysical properties of the subsurface, and it is critical for an initial assessment of a reservoir.

According to the frequency of operation, we distinguish three different measurement acquisition systems for the characterization of subsurface elastic properties: (a) seismic data (typically below 100 Hz), (b) logging sonic measurements (from 2 kHz up to 300 kHz), and (c) core samples (analyzed at 300 kHz-2 MHz). After proper interpretation of the results via advanced numerical methods (see [9, 39, 50]) we obtain an Earth map with its macro-scale velocities. Unfortunately, maps obtained with measurements acquired at dissimilar frequencies are essentially different since a given heterogeneous rock exhibits dispair effective velocities at different frequencies; see e.g., [11, 30, 32, 36, 37, 58, 59]. This undesired phenomena occurs due to the nature of wave propagation in heterogeneous rocks [42].

Biot's theory [8] describes the wave propagation phenomena over a porous rock. However, it ignores the

<footnote>Á. J. Omella (✉)
Universidad del País Vasco (UPV/EHU), Leioa, Spain.
E-mail: ajaome@gmail.com</footnote>

<footnote>M. Strugaru · J. Alvarez-Aramberri
Basque Center for Applied Mathematics (BCAM),
Bilbao, Spain.</footnote>

<footnote>V. Darrigrand
CERFACS. Toulouse, Midi-Pyrnes, France.</footnote>

<footnote>D. Pardo
Universidad del País Vasco (UPV/EHU), Leioa, Spain.
Basque Center for Applied Mathematics (BCAM),
Bilbao, Spain.
Ikerbasque, Bilbao, Spain.</footnote>

<footnote>H. González · C. Santos
Repsol Tech Lab, Madrid, Spain.</footnote>
microscopic level and assumes that all the minerals of the rock have the same bulk and shear moduli [40]. Moreover, the application of this theory using elastic or poro-elastic models is computationally expensive. In here, we focus on estimating the effective P-wave velocity using the acoustic equation and without incurring in a prohibitive computational cost. Hence, we are interested in the first arriving compressional (P) wave that propagates along heterogeneous rocks in the excitation direction. Specifically, we focus only in the first wavefront travel-time, neglecting errors in the amplitude values. Thus, although compressional-to-shear (P-to-S) and S-to-P energy conversions may occur, they may be ignored, and an acoustic approximation of elastic waves is valid in this case, overcoming the limitations exposed in [11].

We need to analyze the structure of a given porous rock and estimate its effective P-wave velocities at different frequencies. With X-ray micro-CT technology (see [19] for a review), it is possible to characterize the internal structure of rock samples at pore scale, identifying with precision the spatial distribution of its constituents. The current resolution of micro-CT available at laboratories of oil companies provides a set of around 1500 images of 1200 x 1200 pixels per core sample. This allows us to identify the material of each voxel in three dimensions (3D) with sub-pore resolution.

After producing a micro-CT scan, we need to analyze its effective velocities at different frequencies. The high-quality requirements for the cores destined to ultrasonic testing makes expensive the laboratory measurements and often can only be reliably estimated at high frequencies (above 500 kHz) to avoid undesired reflections. Some geometrical recommendations to properly analyze cores for ultrasonic testing are available in [2].

On the other hand, the use of relatively simple bounds and “averaging” formulas is extended in the industry for computing the stiffness tensor and the effective velocity from the single properties of each constituent. Some averaging formulas are the so-called Voigt [61] and Reus [49] bounds, Voigt-Reuss-Hill average [31], Hashin-Shtrikman bounds [28, 29], Kuster-Toksz formulation for low-porosity rocks [34], Backus-average [3, 6, 53], Gassmann’s Equation [7, 18, 25], Wyllie time average [20, 52, 62], Raymer’s [48] and Gardner’s [24] relations. These analytical and semi-analytical formulas omit the specific spatial distribution of the rock constituents. In a layered 1D media perpendicular to the direction of propagation, the Backus average is exact in the low-frequency limit, while the Wyllie time-average drives the high-frequency regime.

However, when it comes to realistic two-dimensional (2D) and three-dimensional (3D) heterogeneous rocks, all aforementioned theories fall short to provide accurate effective velocities at different frequencies.

There also exist multiple effective media theories, see [17, 23, 44]. For example, Differential Effective Medium (DEM) [43] considers different pore-type inclusions. Models for cracked medium are available in [13, 16, 22, 27, 33, 60]. There also exist a plethora of multi-scale methods for heterogeneous materials (see e.g., [10, 15, 21] or [26, 38] for a review). While some multi-scale methods are difficult to implement for industrial applications [1], others are unsuitable for high-frequency computations where it is crucial to determine the location of the micro-constituents in the rock at pore-scale. To the best of our knowledge, existing homogenization techniques often fail to solve the problem in the wide range of frequencies, as needed by our application.

Herin, we propose a set of simple-to-implement numerical methods that are capable of handling different frequencies and are suitable for the oil industry needs. These methods produce accurate approximations of the first P-wave arrival corresponding to heterogeneous rocks.

The numerical results presented throughout this work exhibit the need of repeating the rock several times in the low-frequency regime. This increases the size of the problem (information of thousands of millions of pixels in 3D problems with possibly different material properties) and the computational cost, becoming the problem computationally untractable via conventional numerical methods. Moreover, since this is a wave propagation problem over an open domain, truncating the computational domain becomes a great challenge that should be carefully treated to preserve the accuracy of the homogenized velocities. Finally, it is also challenging to estimate effective velocities from a given (simulated) wave propagation solution.

In the low-medium frequency range, we solve the frequency-dependent wave equation by a FEM. To avoid undesired reflections from the domain boundary, we implement a Perfectly Matched Layer (PML). At low frequencies, we observe that this PML adversely affects the results. Thus, we need to extend virtually the original domain by repeating it multiple times. This increases dramatically the problem size. Traditional FEM employ fitting meshes in which material properties are continuous within each element. The large number of pixels attached to the necessity of the rock repetition leads to a prohibitive computational cost when using this technique, especially in 2D and 3D. To overcome this challenge, we propose the use of non-fitting meshes
[14], which reduces dramatically the number of degrees of freedom when the domain is large. However, this technique needs to be carefully employed to preserve the accuracy of the results. The repetition of the original domain allows to precompute blocks in the stiffness matrix in order to reduce the integration and assembling time.

After simulating the complex-valued pressure field in the frequency domain, we apply two methods to estimate the effective compressional wave velocity: (a) Prony's method [46], which provides an accurate value of effective compressional wave velocity at low frequencies, but it diverges at large ones, and (b) counting the number of wavelengths, which provides valid (although sometimes exhibiting low accuracy) approximations of the effective velocities across the entire frequency spectrum.

In the high-frequency limit, we solve the eikonal equation to approximate the travel time of the wave by the FMM [55]. Results in 1D heterogeneous rocks match with those obtained by the application of the Wyllie Time-Average [62]. This confirms the accuracy of our proposed method. Thus, we apply it to estimate the high-frequency P-wave velocity limit in 2D and 3D problems, where analytical methods are unavailable for general porous rocks.

The remainder of the manuscript is organized as follows. Section 2 introduces the theoretical bounds in layered media for the P-wave velocity. Section 3 describes the acoustic formulation in time and frequency. Then, we explain the physical need to repeat the rock, and we introduce two techniques to accelerate the computations: non-fitting meshes and precomputed block matrix. After that, we describe two methods to obtain the effective compressional wave velocity. To conclude this section, we analyze the high-frequency range, in which we propose to solve the eikonal equation to estimate the effective velocity. Section 4 focuses on the numerical results. The first four considered experiments analyze numerically the techniques presented for low-medium frequencies; experiments 5-8 show the influence of the main physical quantities (density, porosity, the size of the pore, and the distribution in the sizes of the pores) in the effective velocities. Experiment 9 shows the effective velocity profile computed in non-periodic VTI rocks. Experiment 10 describes the necessity to extend the domain in the perpendicular dimension of the P-wave propagation direction. Experiment 11 illustrates the method scalability in higher dimensions. Experiment 12 shows the results for horizontal transversely isotropic (HTI) rocks. Experiment 13 checks the validity of the FMM in a synthetic rock. Finally, experiment 14 describes numerical results obtained from images of

real rocks acquired in the laboratory via micro-CT. The last section is devoted to conclusions.

## 2 Theoretical Bounds in Layered Media

We consider a heterogeneous rock composed of various materials. Backus [3] showed that a stratified medium composed of VTI layers behaves like a homogeneous VTI media in the long-wavelength limit. To determine it, we search for “a global” fourth-order stiffness tensor $C_{eff}$ that allows to treat the multi-material body as a homogenized effective one. The Backus homogenized $c_{ijkl}^*$ components of $C_{eff}$ are defined by relations among the components of $C$, and by using the volumetric weighted average defined as $\{\cdot\}$.

Attending to this homogenization, the vertically propagating P-wave velocity [35] is:

$$v_{P,v} = \sqrt{\frac{c_{3333}^*}{\rho^*}}, \qquad (1)$$

where $c_{3333}^* = \{c_{3333}^{-1}\}^{-1}$, and $\rho^* = \{\rho\}$ is the volumetric weighted average of the densities.

Given a VTI layered rock, we consider known some information about each i-constituent (or layer): the volume fraction ($\phi_i$), and some physical properties as density ($\rho_i$) and compressional wave velocity ($v_i$). Following (1), we use the Backus average to obtain the effective P-wave velocity $v_{eff}^B$:

$$v_{P,v} = v_{eff}^B = \sqrt{\frac{c_{3333}^*}{\sum_{i=1}^{n} \phi_i \rho_i}}, \quad \text{where} \quad \frac{1}{c_{3333}^*} = \sum_{i=1}^{n} \frac{\phi_i}{\rho_i v_i^2}. \quad (2)$$

Wyllie time-average [62] describes the effective velocity in high-frequency regimes $v_{eff}^{TA}$ as:

$$\frac{1}{v_{eff}^{TA}} = \sum_{i=1}^{n} \frac{\phi_i}{v_i}. \qquad (3)$$

Figure 1 shows a schematic representation of the effective velocity with respect to the frequency. Specifically, we represent (a) the Backus average zone, where the effective velocity tends to $v_{eff}^B$; (b) the time-average zone where the effective velocity approximates $v_{eff}^{TA}$; and (c), the transition zone that occurs between the two aforementioned zones. The frequency limits in this transition zone are unclear and an analytical expression for the effective velocity is unknown. Indeed, within this transition zone, the behavior of the effective velocity with respect to the frequency is often not even monotonic [59].
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Schematic representation of the three frequency regimes and the effective velocity profile with respect to frequency. $\lambda$ is the frequency, and $d$ represents the pore size average.</figcaption></figure>

## 3 Mathematical Model and Solution Method

In this section, we consider a rock with heterogeneous density and we first focus on the low/medium frequency range. We introduce the mathematical model, followed by the FEM adopted to solve it. Next, we discuss the need to repeat virtually the original domain several times until it becomes at least $n$-wavelengths long, where $n$ is typically between two and four. To make these problems computationally tractable, we consider a FEM with non-fitting meshes and precomputed stiffness matrices. Then, we use the solution obtained with FEM to evaluate the effective velocity via two different techniques, both based on the post-processing of the solution sampled at multiple points. To conclude this section, we focus on the high-frequency regime. Specifically, we introduce the eikonal equation and a Fast Marching Method (FMM) to solve it.

### 3.1 Acoustics

When the density varies with position [5], the strong formulation of the acoustic problem in a fluid domain $\Omega$ is governed by the time-dependent wave equation,

$$ \nabla \cdot \left( \frac{1}{\rho} \nabla \hat{u} \right) - \frac{1}{\rho c^2} \frac{\partial^2 \hat{u}}{\partial t^2} = \hat{f} \quad \text{in } \Omega, \quad (4) $$

where $\hat{u} \equiv \hat{u}(\mathbf{x}, t)$ is the sound pressure, $c \equiv c(\mathbf{x})$ is the propagation velocity of the sound wave in the fluid, $\rho \equiv \rho(\mathbf{x})$ is the material density, and $\hat{f} \equiv \hat{f}(\mathbf{x}, t)$ is the volumetric stationary source term. In order to have unicity in the solution, we consider a

Sommerfeld radiation condition [54, 57], which imposes proper decay conditions at infinity:

$$ \lim_{r \to \infty} \left( r \left( \frac{\partial \hat{u}}{\partial r} + \frac{1}{c} \frac{\partial t}{\partial \hat{u}} \right) \right) = 0, \quad (5) $$

where $r$ denotes the radial component in the spherical coordinate system.

### 3.2 Low-Medium Frequencies Analysis

#### 3.2.1 Mathematical Model

Applying a Fourier transform in time to equation (4), we obtain the following hyperbolic partial differential equation in the frequency domain:

$$ \nabla \cdot \left( \frac{1}{\rho(\mathbf{x})} \nabla u(\mathbf{x}, \omega) \right) + \frac{k^2(\mathbf{x})}{\rho(\mathbf{x})} u(\mathbf{x}, \omega) = f(\mathbf{x}, \omega) \quad \text{in } \Omega, \quad (6) $$

where $u$ and $f$ stand for the Fourier transforms in time of $\hat{u}$ and $\hat{f}$, respectively, and $k(\mathbf{x}) = \frac{\omega}{c(\mathbf{x})}$ is the medium wavenumber at angular frequency $\omega$. In addition, the Sommerfeld radiation condition in spherical coordinates becomes:

$$ \lim_{r \to \infty} \left( r \left( \frac{\partial u}{\partial r} - iku \right) \right) = 0. \quad (7) $$

#### 3.2.2 Finite Element Method (FEM)

To solve equation (6) using a FEM, we first introduce a polygonal domain $\Omega_h$, over which we generate the mesh $T_h = \{K\}$. Then, we define the space of finite elements $V_h := \{v_h \in H^1(\Omega_h) \mid v_h|_K \in P_p(K), \forall K \in T_h$, and $v_h|_{\partial\Omega_h} = 0\}$, where $P_p(K)$ is the space of polynomial functions of order $p$. We then rewrite (6) as a discrete variational problem:

$$ \left\{ \begin{array}{l} \text{Find } u_h \in V_h, \text{ such that:} \\ - \int_{\Omega_h} \frac{1}{\rho} \nabla u_h \cdot \nabla v_h + \int_{\Omega_h} \frac{k^2}{\rho} u_h v_h = \int_{\Omega_h} f v_h, \quad (8) \\ \text{for all } v_h \in V_h, \end{array} \right. $$

where $u_h$ is the approximated solution of $u(\mathbf{x}, \omega)$ that appears by solving a system of linear equations.

We employ a PML [4] to truncate the computational domain and impose homogeneous Dirichlet boundary conditions on its outer part. This ensures a proper Sommerfeld radiation condition and a lack of reflections from the boundary.
**Repetition Strategy.** We first execute simulations at low frequencies –inside the Backus regime– and we observe incorrect values of the effective velocity. This occurs because the PML avoids boundary reflections, while in this application we seek for the homogenized velocities when the rock is repeated multiple times in the subsurface. This problem also affects laboratory experiments and partially explains the difficulties encountered when analyzing rocks at low frequencies. To overcome it in numerical simulations, we repeat the rock domain sample several times until we obtain a computational domain in which the effective velocity is unaffected by a surrounding (homogeneous slower) media. Experiments 3, 4 and 10 in Section 4 illustrate this strategy.

**Non-fitting Meshes.** The standard FEM employs fitting meshes to assemble the bilinear form associated with the material. In particular, it assumes that density $\rho$ and wavenumber $k$ are sufficiently smooth within each cell $K$. However, this assumption makes the resulting simulator prohibitively expensive in our case since it would force us to consider thousands of millions of elements. For the sake of simplicity, we focus on the diffusive term of (8), and we write for a single element $S$:

$$ \int_S \frac{k^2}{\rho} \phi_j \phi_i = \text{area}(S) \frac{k_S^2}{\rho_S} \int_{\hat{S}} \hat{\phi}_j \hat{\phi}_i, \quad (9) $$

where $\phi_i$ is the $i$-th shape function associated with the element $S$, and $\hat{\phi}_i$ is the $i$-th shape function corresponding to the reference element $\hat{S}$.

The key ingredient in non-fitting meshes [14] is to group two or more cells, with possibly different material properties, into a single macro-element $K$ (see Figure 2). Let $s$ be the number of sub-cells (elements) grouped into a macro-element: $K = \bigcup_{l=1}^s S_l$, where $S_l$, $l = 1, ..., s$. Assuming that the material properties are constant inside each sub-cell $S_l$, we rewrite equation (9) for non-fitting meshes as:

$$ \int_K \frac{k^2}{\rho} \phi_j \phi_i = \sum_{l=1}^s \int_{S_l} \frac{k^2}{\rho} \phi_j \phi_i \\ = \sum_{l=1}^s \text{area}(S) \frac{k_{K,l}^2}{\rho_{K,l}} \int_{\hat{S}_l} \hat{\phi}_j \hat{\phi}_i, \quad (10) $$

where $\hat{S}_l$ is the sub-cell associated with the reference macro-element $\hat{K}$, and as before, $\phi_i$ is the $i$-th shape function associated with the reference macro-element $\hat{K}$.

By using non-fitting meshes, we select larger macro-elements and perform exact integration. This reduces

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Sketch of the shape functions (blue lines) and the numbering of elements over a formation with six pixels. Black and white pixels represent materials with a different density and velocity.</figcaption></figure>

the number of unknowns in the finite element model and the corresponding computational requirements (time and memory). For more details about this technique, we refer the reader to [14].

**Precomputed Matrix.** We consider an original master domain $\Omega_0$ discretized by the master mesh $\mathcal{T}_0 = \{K_1^0, K_2^0, ..., K_m^0\}$. We often repeat $\Omega_0$ in space $n$-times to generate the entire domain $\Omega = \bigcup_{i=1}^n \Omega_i$, along with the mesh $\mathcal{T}$ that discretizes the domain $\Omega$, $\mathcal{T} = \bigcup_{i=1}^n \mathcal{T}_i = \{K_1, K_2, ..., K_{m \cdot n}\}$, where $\mathcal{T}_i$ is the discretization of each $\Omega_i$.

We define the operator $\mathcal{F}: K_j \in \mathcal{T} \mapsto K_i^0 \in \mathcal{T}_0$ such that given an element $K_j \in \mathcal{T}$ it returns the corresponding element in the master mesh $\mathcal{T}_0$.

The matching of materials and geometry associated with multiple copies of the same original rock allows us to write for one element:

$$ \int_K \frac{k^2}{\rho} \phi_j \phi_i = \int_{\mathcal{F}(K)} \frac{k^2}{\rho} \phi_j \phi_i. \quad (11) $$

To extend this equality to all terms of the bilinear form, we compute the sub-matrices of each element of the master grid before assembling the global matrix. Then, in order to assemble the matrix associated with $\mathcal{T}$, we search and copy each element matrix into it.

In our application, we employ two sets of precomputed matrices: one for the original rock, that contains the transmitter and its surrounding refinements, and another one for the repeated rock without transmitter and refinements (see Figure 3).

### 3.2.3 From the FEM Solution to an Effective Velocity

At low frequencies, a wave traveling through a heterogeneous media behaves as a sum of plane waves. Figure 4 shows the real part of the pressure computed by Equation (8). In such scenarios, Prony's method [45, 46] accurately estimates the effective velocity. From a uniform sampling of a signal, Prony's method builds a linear combination of damped complex exponentials to
<figure><img src="image_4.png" /><figcaption>Fig. 3 2D model problem with the original rock in green repeated at one side in the propagation direction (N = 3) and at both sides in its perpendicular direction (M = 2). The transmitter, in red, covers the entire domain in the no-propagation direction. The red cross represents the center of the transmitter $T_x$. The blue and black dashed lines represent two lines of analysis. The blue crosses represent the first ($R_{x1}$) and last ($R_{xn}$) receiver of each analysis line. We precompute two different matrices: one for the blocks of the green striped area that contains the transmitter, and another one for the blocks in the brown striped area. The PML zone is filled with gray color.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5 Real part of the solution at 48.3 kHz in 1D rock with 125 mm as size of the pore. The white background color indicates a material with $v_{solid} = 4500$ m/s; the gray background represents a material with $v_{fluid} = 800$ m/s. The density of both materials is 1000 Kg/m³.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 Real part of the solution at 1.27 kHz in a 1D rock with a porous size equal to 125 mm. The white background color indicates a material with $v_{solid} = 4500$ m/s; the gray background represents a material with $v_{fluid} = 800$ m/s. The density of both materials is 1000 Kg/m³.</figcaption></figure>

reconstruct the signal, providing valuable information such as amplitude, phase, damping, and phase velocity.

However, at high frequencies, the wavelength is comparable to the size of the pores, and the solution behaves as piecewise plane waves, one inside each sub-domain. Figure 5 shows the real part of the pressure computed using (8). In this scenario, the plane wave assumption required by Prony's method is invalid, and we need to estimate the effective velocities differently. For that, we propose a simple method based on

counting the number of wavelengths of the solution over the whole domain. To do so, we sample the solution $u_h$ as in the Prony's method, and we compute the phase difference of the pressure from one receiver to the next. In this way, we approximate the number of wavelengths located between one receiver and the next. By repeating the process throughout all receivers, we obtain an estimate of the total number of wavelengths ($\lambda_t$) traveled by the wave. Then, we compute the velocity as:

$$v_{eff} = \frac{L\nu}{\lambda_t}, \qquad (12)$$

where $\nu$ is the frequency in Hertz and $L$ is the distance between the first and last receiver. Experiment 1 in Section 4 compares Prony's method with the one based on counting the number of wavelengths.

## 3.3 High Frequencies Analysis

Approximation of high-frequency solutions by a FEM involves a prohibitive computational cost due to the need of considering several degrees of freedom per wavelength, which translates into large systems of equations. However, our application only requires to estimate the traveltime of the first arriving wave. To do so, we introduce the eikonal equation [47], which we solve with a FMM.
### 3.3.1 Eikonal Equation

We assume the solution of the wave equation (4) takes the following general form:

$$p(t, \boldsymbol{x}) = A(\boldsymbol{x}) e^{-i\omega(T(\boldsymbol{x})+t)}, \quad (13)$$

where A is the amplitude, $\omega$ is the angular frequency and T is a time function (the eikonal), which describes surfaces of constant phase (wavefronts) when T is constant.

Plugging the ansatz (13) into the wave equation (4), separating real and imaginary parts and using the high-frequency assumption ($\omega \to \infty$) in the real part equation, we obtain the so called eikonal equation:

$$|\nabla T(\boldsymbol{x})|^2 = \frac{1}{c(\boldsymbol{x})^2} \quad \text{in } \Omega, \quad (14)$$

where T is the traveltime evaluated at each point $\boldsymbol{x}$ of the domain $\Omega$.

### 3.3.2 Fast Marching Method (FMM)

A FMM is a finite difference scheme that solves boundary value problems with the unknown traveltime satisfying the eikonal equation:

$$\begin{cases} |\nabla T(\boldsymbol{x})|^2 = \frac{1}{c(\boldsymbol{x})^2} \quad \text{in } \Omega, \\ T(\boldsymbol{x}) = f \quad \text{at least on a part of } \partial\Omega, \end{cases} \quad (15)$$

where f is typically 0.

We consider a front given by the initial condition, which is systematically updated in the propagation direction (one grid point at a time), in a downwind fashion from known upwind values. This entropy-satisfying strategy preserves stability in the presence of a wavefront discontinuity, as explained in [56]. The key ingredient of the method is to carefully select the order of traveltime evaluation based on the concept of narrow band: a set of points –that reproduce the shape of the wavefront– with pending traveltime evaluation. Since the wavefront propagates only forward, it crosses each point of the grid only once. The narrow band separates the points with known traveltimes from those with unknown values of T.

Traveltimes in the 3D narrow band for are updated by:

$$\left[ \max(D_{ijk}^{-x}T, -D_{ijk}^{+x}T, 0)^2 + \max(D_{ijk}^{-y}T, -D_{ijk}^{+y}T, 0)^2 + \max(D_{ijk}^{-z}T, -D_{ijk}^{+z}T, 0)^2 \right]^{1/2} = \frac{1}{c_{ijk}}, \quad (16)$$

where $c_{ijk}$ is the velocity at grid point $(i, j, k)$ and $D_{ijk}$ are first order finite difference operators given by:

$$\begin{align} D_{ijk}^{-x} T &= \frac{T_{i,j,k} - T_{i-1,j,k}}{\delta x}, & D_{ijk}^{+x} T &= \frac{T_{i+1,j,k} - T_{i,j,k}}{\delta x} \nonumber \\ D_{ijk}^{-y} T &= \frac{T_{i,j,k} - T_{i,j-1,k}}{\delta y}, & D_{ijk}^{+y} T &= \frac{T_{i,j+1,k} - T_{i,j,k}}{\delta y} \nonumber \\ D_{ijk}^{-z} T &= \frac{T_{i,j,k} - T_{i,j,k-1}}{\delta z}, & D_{ijk}^{+z} T &= \frac{T_{i,j,k+1} - T_{i,j,k}}{\delta z} \tag{17} \end{align}$$

and $\delta_x$, $\delta_y$ and $\delta_z$ are the grid spacings in x, y, and z directions, respectively. We select a grid conforming with the spatial resolution given by the tomography:

$$\delta(\cdot) = \frac{N(\cdot)}{L(\cdot)}, \quad (18)$$

where subscript $(\cdot)$ denote the spatial component, and $N(\cdot)$ and $L(\cdot)$ denote the number of voxels obtained in the tomography, and the dimension of the rock, respectively.

To initialize the narrow band, we compute analytically the traveltime in one end of the rock as:

$$T_{i,j,k} = \frac{c_{i,j,k}}{2 \delta_z}. \quad (19)$$

for $i \in 1, 2, \dots, N_x$, $j \in 1, 2, \dots, N_y$ and $k = 1$.

The FMM algorithm reads:

1. Generate the grid conforming to the voxels.

2. Initialize the narrow band using Equation (19).

3. Do until all determining all traveltimes:

(a) Find the point with minimum traveltime in the narrow band.

(b) Update the narrow band.

(c) Update traveltimes at the points recently added to the narrow band and their neighbours using equation (16).

(d) Compute traveltimes on the bottom part of the domain using Equation (19) with $k = N_z$.

Figure 6 illustrates the sequence 3(a) - 3(b) in 2D. The first order scheme described above is unconditionally stable [55]. For more details on FMM, we refer the reader to [51].

In 1D, a FMM reduces to computing the traveltime from one point of the grid to the next one, using the spacing between the two points and the known velocity at the arrival point.

## 4 Numerical Results

This section contains multiple numerical results. Experiments 1-4 consider 1D periodic and non-periodic
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Example of a narrow band evolution in a two dimensional domain: green squares stand for points where the traveltime is known; red empty circles are points with unknown travel time; black filled circles are the points in the narrow band. (a) A blue triangle identifies the minimum traveltime in the narrow band. (b) The minimum and its neighbors of the narrowband (until correcting all arising discontinuities) are tagged as known traveltime points, let advancing the narrow band.</figcaption></figure>

formations, and assess the performance of the methods for different values of the involved parameters (e.g., number of rock repetitions and number of wavelengths needed in the domain). Experiments 5-8 show the influence of the principal physical quantities (e.g., density, porosity, the size of a pore, and different pore distribution) in the effective velocity. Experiment 9 describes the effective velocity profile computed for four non-periodic formations. Experiment 10 analyzes the lateral extension of the domain for a 2D rock. Experiment 11 exhibits a 2D VTI rock. It also shows the match between the Time-Average formula and FMM in the high-frequency regime. Experiment 12 employs HTI rocks and shows the discrepancy between the Time-Average formula and the FMM. Then, experiment 13 further validates the FMM. Finally, experiment 14 applies the proposed methods to two real rocks.

## 4.1 Samples and Modeling Considerations

We consider samples composed of two materials. We identify solid (S) and fluid (F) pixels with white and black color, respectively (see Table 1).

We generate particular periodic samples by repeating a specific sequence of pixels several times. To refer to a periodical 1D synthetic rock, we introduce the following nomenclature: we first indicate that it is a one-dimensional (1D) periodic (P) rock with the abbreviation 1DP. Next, we write the number of total pixels that compose the sample and finally, inside brackets, the repeated sequence of pixels.

<table><caption>**Table 1** Nomenclature associated to each one-dimensional pixel.</caption><thead><tr><th>Graphical representation</th><th>Material state</th><th>Material properties</th></tr></thead><tbody><tr><td>□</td><td>Solid (S)</td><td>v<sub>solid</sub> &rho;<sub>solid</sub></td></tr><tr><td>■</td><td>Fluid (F)</td><td>v<sub>fluid</sub> &rho;<sub>fluid</sub></td></tr></tbody></table>

Table 2 shows two periodic samples.

<table><caption>**Table 2** Nomenclature and porosity associated with two one-dimensional periodic formations.</caption><thead><tr><th>Graphical representation</th><th>Nomenclature</th><th>Porosity</th></tr></thead><tbody><tr><td>□□□□□</td><td>1DP8(7S1F)</td><td>12.5%</td></tr><tr><td>□□□□□</td><td>1DP8(3S1F)</td><td>25%</td></tr></tbody></table>

We additionally consider four large realistic non-periodic formations composed each of 1203 pixels. Table 3 summarizes the nomenclature and the porosity of each non-periodic formation.

<table><caption>**Table 3** Nomenclature and porosity of non-periodic one-dimensional formations.</caption><thead><tr><th>Nomenclature</th><th>Porosity</th><th>Nomenclature</th><th>Porosity</th></tr></thead><tbody><tr><td>1D Rock 1</td><td>11.31%</td><td>1D Rock 3</td><td>21.36%</td></tr><tr><td>1D Rock 2</td><td>22.61%</td><td>1D Rock 4</td><td>25.85%</td></tr></tbody></table>

Figure 7 describes the one-dimensional model of the formation and identifies the necessary geometric model parameters.

<figure><img src="image_8.png" /><figcaption>**Fig. 7** 1D Model problem. The transmitter (in red) is defined by its central point $T_x$ and size $T_s$. The blue crosses with coordinates $R_{x1}$ and $R_{xn}$ indicate the first and last receiver, respectively.</figcaption></figure>

We denote the domain of the sample by $\Omega_0 = (a, b)$. To truncate the computational domain, we introduce a PML on the outer part of the domain. We automatically select its size $|a - a'| = |b - b'|$ as the minimum between 10% of the domain of the sample and three times the
<table><caption>**Table 5** Nomenclature and porosity associated with 3D periodic formations.</caption><thead><tr><th>Graphical representation</th><th>Nomenclature</th><th>Porosity</th></tr></thead><tbody><tr><td></td><td>3DVTI 12x12x12(1F2S)</td><td>33.3%</td></tr><tr><td></td><td>3DHTI<sub>x</sub> 12x12x12(1F3S)</td><td>25%</td></tr><tr><td></td><td>3DHTI<sub>y</sub> 12x12x12(1F3S)</td><td>25%</td></tr></tbody></table>

<table><caption>**Table 6** Material number, volume fraction, wave velocity and density of the 3D Rock 1 constituents.</caption><thead><tr><th>Material number (i)</th><th>Volume fraction (φ<sub>i</sub>)</th><th>c<sub>i</sub> (m/s)</th><th>ρ<sub>i</sub> (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>1</td><td>0.800185</td><td>6041.98</td><td>2648.0</td></tr><tr><td>2</td><td>0.034559</td><td>2424.87</td><td>1373.0</td></tr><tr><td>3</td><td>0.071186</td><td>1513.27</td><td>1000.0</td></tr><tr><td>4</td><td>0.046283</td><td>2138.54</td><td>1298.4</td></tr><tr><td>5</td><td>0.047787</td><td>1493.02</td><td>1149.2</td></tr></tbody></table>

<table><caption>**Table 7** Material number, volume fraction, wave velocity and density of the 3D Rock 3 constituents.</caption><thead><tr><th>Material number (i)</th><th>Volume fraction (φ<sub>i</sub>)</th><th>c<sub>i</sub> (m/s)</th><th>ρ<sub>i</sub> (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>1</td><td>0.619703</td><td>6645.11</td><td>2710.0</td></tr><tr><td>2</td><td>0.354040</td><td>5161.60</td><td>2197.0</td></tr><tr><td>3</td><td>0.026257</td><td>1513.27</td><td>1000.0</td></tr></tbody></table>

## 4.2 Experiment 1: Prony vs Count Wavelength

In this experiment, we extend the domain at both sides of the original sample to guarantee $10 \lambda_{solid}$ at each side.

Figure 10 shows the theoretical effective velocities computed by equations (2) and (3), and compares the profiles of compressional wave effective velocities

<figure><img src="image_9.png" /><figcaption>**Fig. 10** Effective compressional wave velocities computed with Prony's method (red line) and counting wavelengths (green line).</figcaption></figure>

obtained by the application of Prony's method (red line) and counting the number of wavelengths (green line) for two periodic rocks. We observe that Prony's method converges only in the low-frequency regime, which is the region where the solution is given by a sum of plane waves. The count wavelengths method properly estimates velocities at low frequencies, matching the results obtained by Prony's method. It also performs adequately in the high-frequency regime and converges to the theoretical solution $v_{eff}^{TA}$ given by the Wyllie time-average.
<table><caption>**Table 8** Summary of simulation parameters used in each experiment.</caption><thead><tr><th>Experiment number</th><th>Figure or table number</th><th>Formation nomenclature</th><th>c<sub>solid</sub> (m/s)</th><th>ρ<sub>solid</sub> (kg/m<sup>3</sup>)</th><th>c<sub>fluid</sub> (m/s)</th><th>ρ<sub>fluid</sub> (kg/m<sup>3</sup>)</th><th>Precomputed matrix</th><th>Fitting meshes</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>Figure 10a</td><td>1DP8(3S1F)</td><td rowspan="2">4500</td><td rowspan="2">2800</td><td rowspan="2">800</td><td rowspan="2">1000</td><td rowspan="2">No</td><td rowspan="2">No</td></tr><tr><td>Figure 10b</td><td>1DP8(7S1F)</td></tr><tr><td rowspan="2">2</td><td>Table 9</td><td>1D Rock 1-4</td><td rowspan="2">4500</td><td>1000</td><td rowspan="2">800</td><td rowspan="2">1000</td><td>*</td><td>*</td></tr><tr><td>Figure 11</td><td>1D Rock 4</td><td>*</td><td>Yes</td><td>Yes</td></tr><tr><td rowspan="2">3</td><td>Figure 12a</td><td>1DP8(3S1F)</td><td rowspan="2">4500</td><td rowspan="2">2800</td><td rowspan="2">800</td><td rowspan="2">1000</td><td rowspan="2">No</td><td rowspan="2">Yes</td></tr><tr><td>Figure 12b</td><td>1DP8(7S1F)</td></tr><tr><td rowspan="2">4</td><td>Figure 13a</td><td>1D Rock 1</td><td rowspan="2">4500</td><td rowspan="2">1000</td><td rowspan="2">800</td><td rowspan="2">1000</td><td rowspan="2">Yes</td><td rowspan="2">Yes</td></tr><tr><td>Figure 13b</td><td>1D Rock 2</td></tr><tr><td rowspan="2">5</td><td>Figure 14a</td><td>1DP8(3S1F)</td><td rowspan="2">4500</td><td>*</td><td rowspan="2">800</td><td>*</td><td>Yes</td><td>Yes</td></tr><tr><td>Figure 14b</td><td>1DP8(7S1F)</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">6</td><td>Figure 15a</td><td>1DP8(1S1F)</td><td rowspan="3">4500</td><td>2800</td><td rowspan="3">800</td><td rowspan="3">1000</td><td rowspan="3">Yes</td><td rowspan="3">Yes</td></tr><tr><td>Figure 15b</td><td>1DP8(3S1F)</td><td>1000</td></tr><tr><td></td><td>1DP8(7S1F)</td><td></td></tr><tr><td rowspan="5">7</td><td>Figure 16a</td><td>1DP8(1S1F)</td><td rowspan="5">4500</td><td>2800</td><td rowspan="5">800</td><td rowspan="5">1000</td><td rowspan="5">Yes</td><td rowspan="5">Yes</td></tr><tr><td>Figure 16b</td><td>1DP80(1S1F)</td><td>1000</td></tr><tr><td></td><td>1DP800(1S1F)</td><td></td></tr><tr><td>Figure 17a</td><td>1DP8(3S1F)</td><td>2800</td></tr><tr><td>Figure 17b</td><td>1DP8(7S1F)</td><td></td></tr><tr><td rowspan="2">8</td><td>Figure 18a</td><td>*</td><td rowspan="2">4500</td><td rowspan="2">2800</td><td rowspan="2">800</td><td rowspan="2">1000</td><td rowspan="2">Yes</td><td rowspan="2">Yes</td></tr><tr><td>Figure 18b</td><td></td><td></td><td></td></tr><tr><td>9</td><td>Figure 19</td><td>1D Rock 1-4</td><td>4500</td><td>2800</td><td>800</td><td>1000</td><td>Yes</td><td>Yes</td></tr><tr><td>10</td><td>Figure 20</td><td>2DVTI 12x12 (1F2S)</td><td>4500</td><td>1000</td><td>800</td><td>1000</td><td>Yes</td><td>Yes</td></tr><tr><td rowspan="5">11</td><td>Figure 22a</td><td>1DP12(1S1F)</td><td rowspan="5">4500</td><td>2800</td><td rowspan="5">800</td><td rowspan="5">1000</td><td rowspan="5">Yes</td><td rowspan="5">Yes</td></tr><tr><td>Figure 22b</td><td>2DVTI (1F2S)</td><td>1000</td></tr><tr><td></td><td>2D VTI (1F2S)</td><td></td></tr><tr><td>Table 11</td><td>2D VTI (1F3S)</td><td></td></tr><tr><td></td><td>2D VTI (1F7S)</td><td></td></tr><tr><td rowspan="5">12</td><td>Figure 23a</td><td>2DHTI (1F2S)</td><td rowspan="5">4500</td><td>2800</td><td rowspan="5">800</td><td rowspan="5">1000</td><td rowspan="5">Yes</td><td rowspan="5">Yes</td></tr><tr><td>Figure 23b</td><td>2D HTI (1F2S)</td><td>1000</td></tr><tr><td></td><td>2D HTI (1F3S)</td><td></td></tr><tr><td>Table 12</td><td>2D HTI (1F7S)</td><td></td></tr><tr><td></td><td>3D HTI<sub>x</sub> (1F3S)</td><td></td></tr><tr><td rowspan="2">13</td><td>Figure 24b</td><td>Labyrinth</td><td>4500</td><td>—</td><td>800</td><td>1000</td><td>—</td><td>—</td></tr><tr><td>Figure 25</td><td>3D Rock 1</td><td>*</td><td>*</td><td>*</td><td>*</td><td>Yes</td><td>Yes</td></tr><tr><td>14</td><td>Figure 26</td><td>3D Rock 2</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Symbol * indicates that the corresponding field is a variable explained in its own experiment subsection.

Symbol — indicates that such field is not applicable in the considered experiment.
## 4.3 Experiment 2: Precomputed-Matrix and Fitting vs Non-Fitting Meshes

In this experiment, we assess the accuracy of non-fitting meshes and precomputed matrix technique. We compute the compressional wave effective velocity using Prony's method for the non-periodic formations displayed in Table 3, excited at a frequency of 10 Hz. We extend the domain at both sides of the original sample virtually to guarantee $1 \lambda_{Backus}$ at each side.

Table 9 shows the effective velocity (in m/s) computed over the non-periodic formations defined in Table 3 using the three different methods listed below.

- Method 1 (M1) employs a FEM with a traditional fitting mesh technique.

- Method 2 (M2) uses a FEM with a non-fitting mesh with only one element in each repeated rock and without a precomputed matrix.

- Method 3 (M3) considers a FEM with a non-fitting mesh and a precomputed matrix.

This table also exhibits the Backus effective velocity $v_{eff}^B$ given by equation (2).

<table><caption>**Table 9** Comparison between the theoretical compressional wave effective velocity $v_{eff}^B$ and the computed one (using Prony) from the FEM solution over the non-periodic formations using three different methods: M1 employs fitting meshes; M2 uses a non-fitting mesh without precomputed matrix; and M3 considers a non-fitting mesh with precomputed matrix. The values for the velocity are in m/s.</caption><thead><tr><th>Formation</th><th>v<sub>eff</sub><sup>B</sup></th><th>M1</th><th>M2</th><th>M3</th></tr></thead><tbody><tr><td>1D Rock 1</td><td>2129.870</td><td>2129.866</td><td>2129.870</td><td>2129.870</td></tr><tr><td>1D Rock 2</td><td>1598.210</td><td>1598.222</td><td>1598.220</td><td>1598.220</td></tr><tr><td>1D Rock 3</td><td>1638.169</td><td>1638.205</td><td>1638.199</td><td>1638.199</td></tr><tr><td>1D Rock 4</td><td>1506.608</td><td>1506.513</td><td>1506.547</td><td>1506.547</td></tr></tbody></table>

Table 9 reveals that the computation of the effective velocity through any of the aforementioned methods provides an excellent approximation to the Backus average velocity $v_{eff}^B$ (see equation 2). Moreover, the results obtained with methods two and three coincide. This validates the implementation of the precomputed matrix technique. The solution associated with a non-fitting mesh differs slightly from the one computed with a fitting mesh. To assess these discrepancies, Figure 11 shows the velocity computed with non-fitting meshes against the number of subcells for different values of the fast material density, in the case of 1D Rock 4. For a large number of subcells per element, the error increases as function of the ratio between densities.

<figure><img src="image_10.png" /><figcaption>**Fig. 11** Influence of the number of subcells per element in the estimated effective compressional wave velocity solution. Non-periodic 1D Rock 4 sample excited at 10 Hz for multiple densities.</figcaption></figure>

## 4.4 Experiment 3: Rock-repetitions: One Side vs Two Sides

We consider two periodic rocks. Figure 12 shows the effective velocity (obtained counting wavelengths) in the original rock, as well as the effect of repeating the rock at one or both sides of the original rock.

At low frequencies, the effective velocity is different from the Backus velocity $v_{eff}^B$ when considering only the original rock. But if we repeat it at one or both sides, then the two velocities match. A one-side repetition is computationally more efficient. At high frequencies, all simulations converge to $v_{eff}^{TA}$.

The results of this experiment indicate we should repeat the rock at one side and guarantee a minimum number of wavelengths in the domain.

## 4.5 Experiment 4: Rock-repetitions: Number of Wavelengths

This experiment analyzes the number of times we should repeat the rock. For that, we use the samples 1D Rock 1 and 1D Rock 2.

Figure 13 studies the effect of the number of wavelengths contained in the domain on the effective velocity. We consider three cases: 10, 3, and 1 wavelengths of the solid. Aside from slight differences in the solution, all of them provide good approximations of the quantity of interest.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12** Effective compressional wave velocities computed counting wavelengths on different domain sizes.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Effective compressional wave velocities computed counting wavelengths. The repetition of the rock at only one side guarantees a domain size of $N\lambda_{solid}$.</figcaption></figure>

## 4.6 Experiment 5: Influence of Density

At low frequencies, an increase in the density of the fast material diminishes the effective velocities (see Figure 11). Figure 14 investigates the effect of different density ratios $\rho_{solid}/\rho_{fluid}$ on the compressional wave velocity for two periodic rocks. The one-side repetition of the rock guarantees a domain size of $1 \lambda_{solid}$.

We observe a decrease of the effective compressional wave velocity with the increase of the ratio $\rho_{solid}/\rho_{fluid}$ at low frequencies. The frequency that separates the Backus average zone and the transition zone

takes smaller values when increasing this ratio. At high frequencies, the value of the compressional wave velocity is independent of the densities. These two observations agree with the theoretical limits discussed in Section 2.

## 4.7 Experiment 6: Influence of Porosity

We consider the samples 1DP8(1S1F), 1DP8(3S1F) and 1DP8(7S1F) with porosities equal to 50%, 25%, and 12.5%, respectively. We repeat the rock at one-side to guarantee a domain size of $1 \lambda_{solid}$. Figure
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Effective compressional wave velocities computed counting wavelengths with different density ratios.</figcaption></figure>

15 shows a decrease on the effective compressional wave velocity when the porosity diminishes. The three curves display a qualitatively similar profile –in terms of oscillations and frequency limits– in the transition zone. The frequency that separates the Backus average zone from the transition zone increases when augmenting the porosity.

## 4.8 Experiment 7: Influence of the Size of the Pore

In this experiment, we repeat the rock on one-side to guarantee $1 \lambda_{solid}$. Figure 16 illustrates the behavior

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15** Effective compressional wave velocities computed counting wavelengths in the rock 1DP8(3S1F) –green line–, 1DP8(7S1F) –red line– and 1DP8(1S1F) –dotted gray line–.</figcaption></figure>

of the effective velocity in three synthetic rocks: 1DP8(1S1F), 1DP8(1S1F) and 1DP8(1S1F), all of them of porosity 50% and size of the pores ($d_{fluid}$) equal to 125, 12.5, and 1.25 millimeters, respectively.

Figure 17a considers the synthetic rocks 1DP8(3S1F), 1DP80(3S1F) and 1DP800(3S1F) with porosity 25%. Figure 17b employs the rocks 1DP8(7S1F), 1DP80(7S1F) and 1DP800(7S1F) with porosity 12.5%.

In all case of Figures 16 and 17 we observe a shift of the curves on frequency as the size of the pore decreases. We conclude the shift is independent of the porosity and density. These results agree with the theoretical 1D
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16** Effective compressional wave velocities computed counting wavelengths in samples with different sizes of the pore ($d_{fluid}$).</figcaption></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 17** Effective compressional wave velocities computed counting wavelengths in samples with different sizes of the pore ($d_{fluid}$).</figcaption></figure>

results (see Figure 1): the Backus average zone appears when the ratio between the wavelength and the size of the pore is greater than one. When this ratio is much smaller than one, we are in the time average zone.

## 4.9 Experiment 8: Influence of the Distribution of the Sizes of the Pore

We consider rocks with a specific size of the pore distribution by joining the rocks 1DP400(150T50F) with 1DP400(3T1F) –porosity equal to 25%–, and

1DP400(350T50F) with 1DP400(7T1F) –porosity equal to 12.5%–.

Figure 18 illustrates how the transition zone between low and high frequencies is wide if there are pores with different sizes in the same rock. In this case, the frequency that indicates the change from the Backus average zone to the transition zone is related to the largest size of the pore, whereas the frequency from the transition zone to the time average zone is associated with the smallest size of the pore.
<figure><img src="image_23.png" /><figcaption>**Fig. 18** Effective compressional wave velocities computed counting wavelengths in samples with different size of the pore ($d_{fluid}$).</figcaption></figure>

## 4.10 Experiment 9: 1D Non-periodical Formations

Figure 19 displays the effective velocity profiles over the rocks displayed in Table 3. The profiles converge to the theoretical bounds at low and high frequencies.

To assess FMM's efficiency, we compare its effective velocities with the Wyllie time average $v_{eff}^{TA}$. Table 10 shows that the two methods deliver equivalent results in 1D.

<figure><img src="image_24.png" /><figcaption>**Fig. 19** Effective compressional wave velocities computed counting wavelengths for the non-periodic rocks.</figcaption></figure>

<table><caption>**Table 10** Effective velocities obtained for non-periodical rocks in the high-frequency regime, with Time-Average and FMM</caption><thead><tr><th>Formation</th><th>Porosity</th><th><i>v</i><sup>TA</sup><sub>eff</sub> (m/s)</th><th>FMM (m/s)</th></tr></thead><tbody><tr><td>1D Rock 1</td><td>11.31%</td><td>2954.9672</td><td>2954.9673</td></tr><tr><td>1D Rock 2</td><td>22.61%</td><td>2199.7156</td><td>2199.7156</td></tr><tr><td>1D Rock 3</td><td>21.36%</td><td>2263.5234</td><td>2265.7145</td></tr><tr><td>1D Rock 4</td><td>25.85%</td><td>2049.5007</td><td>2051.2966</td></tr></tbody></table>

## 4.11 Experiment 10: Lateral Extension

This experiment analyzes the importance of the domain size in the perpendicular direction to the main P-wave propagation direction (controlled by the parameter *M*). For this, we test the sample 2DVTI 12x12(1F2S). We
postprocess the solution over one vertical line located at the center of the rock.

Figure 20 reveals that at low frequencies, the domain should be repeated along its lateral dimension (perpendicular direction to the main P-wave propagation direction). In the remaining experiments, we will select M to guarantee at least one $\lambda_{solid}$.

<figure><img src="image_25.png" /><figcaption>Fig. 20 Effective velocities for the sample 2DVTI 12x12 (1F2S) for different values of M.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 21 2D model problem with periodic boundary conditions (BCs) in the direction perpendicular to the source (red line). The original rock is mirrored according to the axis of symmetry (blue line).</figcaption></figure>

It is possible to repeat infinite times the rock in the perpendicular directions to the source. To do so, it is sufficient to mirror the original rock once in 2D (or four times in 3D) and impose periodic boundary conditions (see Figure 21). While this is a viable venue, we did not

explore it in this work. Notice, however, that a similar alternative does not exist in the propagation direction due to the presence of the source.

## 4.12 Experiment 11: VTI Rocks

This experiment compares the solution for the VTI Rock 2DVTI 12x12 (1F2S) defined in Table 4 with the solution computed for the Rock 1DP12(1F2S). We post-process the solution over one vertical line located at the center of the rock.

<figure><img src="image_27.png" /><figcaption>Fig. 22 Comparative of the effective compressional wave velocity between a 2D 12x12 VTI (1F2S) rock and its 1D equivalent.</figcaption></figure>
Figure 22 shows similar results for 1D and 2D. However, in the case of a large contrast between the two densities, the 2D solution exhibits additional oscillations in the transition zone before reaching the high-frequency velocity limit.

Table 11 shows a match between the velocity computed with the FMM and with the Wyllie time-average formula (3).

<table><caption>**Table 11** Effective P-wave velocities computed with Time-Average and FMM for 2D and 3D VTI rocks</caption><thead><tr><th>Nomenclature</th><th>porosity</th><th>v<sup>TA</sup><sub>eff</sub> (m/s)</th><th>FMM (m/s)</th></tr></thead><tbody><tr><td>2D VTI (1F2S)*</td><td>33.3%</td><td>1770.49</td><td>1768.92</td></tr><tr><td>2D VTI (1F3S)*</td><td>25%</td><td>2086.96</td><td>2084.77</td></tr><tr><td>2D VTI (1F7S)*</td><td>12.5%</td><td>2541.18</td><td>2538.43</td></tr><tr><td>3D VTI (1F3S)**</td><td>25%</td><td>2086.96</td><td>2086.96</td></tr></tbody></table>

* Size of 1200x1200 pixels

** Size of 100x100x100 voxels

## 4.13 Experiment 12: HTI Rocks

Figure 23 shows the effective P-wave velocity for the periodic rock 2DHTI 12x12(1F2T) defined in Table 4. In this case, the choice of the analysis line is crucial for the accurate computation of the effective velocity.

The effective compressional wave velocity computed by FMM matches with the value of the fast material, regardless the porosity (see Table 12). On the other hand, the velocity computed by the Wyllie time-average formula (3) is much lower. This result confirms that the use of the Wyllie time-average formula is invalid outside of VTI symmetries, as physically expected.

<table><caption>**Table 12** Effective P-wave velocities computed with Time-Average and FMM. $c_{solid} = 4500$ m/s</caption><thead><tr><th>Nomenclature</th><th>porosity</th><th>v<sup>TA</sup><sub>eff</sub> (m/s)</th><th>FMM (m/s)</th></tr></thead><tbody><tr><td>2D HTI (1F2S)*</td><td>33.3%</td><td>1770.49</td><td>4500.00</td></tr><tr><td>2D HTI (1F3S)*</td><td>25%</td><td>2086.96</td><td>4500.00</td></tr><tr><td>2D HTI (1F7S)*</td><td>12.5%</td><td>2541.18</td><td>4500.00</td></tr><tr><td>3D HTI<sub>x</sub> (1F3S)**</td><td>25%</td><td>2086.96</td><td>4499.99</td></tr><tr><td>3D HTI<sub>y</sub> (1F3S)**</td><td>25%</td><td>2086.96</td><td>4499.99</td></tr></tbody></table>

* Size of 1200x1200 pixels

** Size of 100x100x100 voxels

## 4.14 Experiment 13: The *Labyrinth*

In this experiment, we consider the *labyrinth* depicted in Figure 9. At high frequencies, the first arriving wave

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 23** Effective compressional wave velocity for 2D HTI 12x12 (1F2S) rock. The black solid curve corresponds to the effective P-wave velocity measured along the fluid (L<sub>2,x</sub>), while the dashed red curve is recorded along the solid (L<sub>1,x</sub>).</figcaption></figure>

circumvents the fluid by following the solid red line in Figure 24a. Its effective velocity is approximately $c_{solid}/2$. Figure 24b displays the effective velocity computed with the FMM at the bottom of the domain as a function of space. These results confirm that the FMM is trustworthy and provides physically consistent velocities at all arriving points in a non-trivial domain.
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>**Fig. 24** Sketch of the labyrinth domain and the computed effective velocity at the bottom of the core by the FMM.</figcaption></figure>

## 4.15 Experiment 14: Real rocks

Figures 25 and 26 show the effective P-wave velocity average analyzed along five lines in a cross-section of 3D Rock1 (see Table 6), and 3D Rock 2, respectively.

Table 13 collects the values computed with FMM for the full 3D rocks as well as some of their cross-sections. As expected, the effective compressional wave velocity computed in the high-frequency regime for the 3D rock is in between the one computed for the cross-section and the sound velocity of the fastest material.

<figure><img src="image_32.png" /><figcaption>**Fig. 25** Average of effective P-wave velocity for a cross-section of 3D Rock 1.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Fig. 26** Average of effective P-wave velocity for a cross-section of 3D Rock 2.</figcaption></figure>

<table><caption>**Table 13** Effective P-wave velocities computed by FMM and faster constituent sound velocity</caption><thead><tr><td>Formation</td><td>v<sup>FMM-2D</sup><sub>eff</sub></td><td>v<sup>FMM</sup><sub>eff</sub> (m/s)</td><td>c<sub>fast</sub> (m/s)</td></tr></thead><tbody><tr><td>3D Rock 1</td><td>5946.94</td><td>6031.15</td><td>6041.98</td></tr><tr><td>3D Rock 2</td><td>6415.495</td><td>6607.81</td><td>6645.11</td></tr></tbody></table>

# 5 Conclusions

We propose a set of numerical methods to estimate the effective compressional wave velocity of heterogeneous rocks. At low frequencies, the PML pollutes the values of the computed effective compressional wave ve-
locity. Thus, we extend the domain through rock repre-
tition in order to guarantee at least two wavelengths
in the propagation direction. To make the resulting
problem computationally tractable, we employ non-
fitting meshes, which introduce some controlled numer-
ical errors. When carefully implemented, the resulting
effective velocities exhibit acceptable accuracies. The
Prony's post-processing method provides accurate val-
ues of the effective compressional wave velocity in the
Backus average zone, but the method diverges in the
transition zone, where the behavior of the wave is no
longer given as a sum of plane waves. On the other
hand, the value of the effective compressional wave
velocity estimated by counting the numbers of wave-
lengths converges in the entire frequency spectrum.

At high frequencies, the computation of the effective
velocity by FMM matches with the theoretical limit
expected in VTI rocks. We show the validity of this
numerical method in higher spatial dimensions.

Numerical results, in accordance with the existing
theory, indicate the following: (a) an increase in
the density of the fast material decreases the value
of the effective compressional wave velocity at low
frequencies. On the other hand, effective velocities at
high frequencies become unaltered by these density
changes, as predicted by the theory; (b) an increase
in the porosity produces a decrease in the value of the
effective compressional wave velocity. This decrease is
more pronounced at high frequencies; (c) the presence
of different sizes of the pore produces a widening of
the transition zone since the most significant size of
the pore is related to the frequency regime in which
the Backus average is valid, while the smallest size of
the pore characterizes the frequency for which the time
average holds; (d) the Wyllie time-average formula fails
outside of VTI rocks.

Possible extensions of this work include the use of
parametric solvers like the one presented in [41] for the
low frequency spectrum to analyze simultaneously a
group of rocks.

**Acknowledgements** This research was supported by Rep-
sol, the European Union's Horizon 2020 research and in-
novation programme under the Marie Skłodowska-Curie
grant agreement No 777778 (MATHROCKS), the European
POCTEFA 2014-2020 Project PIXIL (EFA362/19) by the
European Regional Development Fund (ERDF) through the
Interreg V-A Spain-France-Andorra programme, the Project
of the Spanish Ministry of Economy and Competitiveness
with reference MTM-2016-76329-R (AEI/FEDER, EU), the
BCAM "Severo Ochoa" accreditation of excellence (SEV-
2017-0718), and the Basque Government through the BERC
2018-2021 program, the two Elkartek projects ArgIA (KK-
2019-00068) and MATHEO (KK-2019-00085), the grant "Ar-
tificial Intelligence in BCAM number EXP. 2019/00432", and

the Consolidated Research Group MATHMODE (IT1294-19)
given by the Department of Education.

References

1. Abdulle, A., E, W., Engquist, B., Vanden-Eijnden, E.: The heterogeneous multiscale method. Acta Numerica **21**, 1–87 (2012). DOI 10.1017/S0962492912000025

2. Aydin, A.: Upgraded ISRM suggested method for determining sound velocity by ultrasonic pulse transmission technique. Rock Mechanics and Rock Engineering **47**(1), 255–259 (2014). DOI 10.1007/s00603-013-0454-z. URL https://doi.org/10.1007/s00603-013-0454-z

3. Backus, G.: Long-wave elastic anisotropy produced by horizontal layering. Journal of Geophysical Research **67**(11), 4427–4440 (1962). DOI 10.1029/JZ067i011p04427. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JZ067i011p04427

4. Berenger, J.: A perfectly matched layer for the absorption of electromagnetic waves. Journal of computational physics **114**(2), 185–200 (1994)

5. Bergmann, P.G.: The wave equation in a medium with a variable index of refraction. The Journal of the Acoustical Society of America **17**(4), 329–333 (1946). DOI 10.1121/1.1916333. URL https://doi.org/10.1121/1.1916333

6. Berryman, J.G.: Long-wave elastic anisotropy in transversely isotropic media. GEOPHYSICS **44**(5), 896–917 (1979). DOI 10.1190/1.1440984. URL https://doi.org/10.1190/1.1440984

7. Berryman, J.G.: Origin of Gassmann's equations. GEOPHYSICS **64**(5), 1627–1629 (1999). DOI 10.1190/1.1444667. URL https://doi.org/10.1190/1.1444667

8. Biot, M.A.: Theory of propagation of elastic waves in a fluidsaturated porous solid. i. lowfrequency range. The Journal of the Acoustical Society of America **28**(2), 168–178 (1956). DOI 10.1121/1.1908239. URL https://doi.org/10.1121/1.1908239

9. Bonnasse-Gahot, M., Calandra, H., Diaz, J., Lanteri, S.: Hybridizable discontinuous Galerkin method for the 2D frequency-domain elastic wave equations. Geophysical Journal International **213**(1), 637–659 (2017). DOI 10.1093/gji/ggx533. URL https://doi.org/10.1093/gji/ggx533

10. Calo, V.M., Efendiev, Y., Galvis, J., Li, G.: Randomized oversampling for generalized multi-scale finite element methods. Multiscale Modeling & Simulation **14**(1), 482–501 (2016). DOI
10.1137/140988826. URL https://doi.org/10.1137/140988826

11. Cance, P., Capdeville, Y.: Validity of the acoustic approximation for elastic waves in heterogeneous media. GEOPHYSICS 80(4), T161–T173 (2015). DOI 10.1190/geo2014-0397.1. URL https://doi.org/10.1190/geo2014-0397.1

12. Carcione, J.: Wave fields in real media: Wave propagation in anisotropic, anelastic, porous and electromagnetic media, vol. 38. Elsevier (2007)

13. Chapman, M.: Frequency-dependent anisotropy due to meso-scale fractures in the presence of equant porosity. Geophysical Prospecting 51(5), 369–379 (2003). DOI 10.1046/j.1365-2478.2003.00384.x. URL https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-2478.2003.00384.x

14. Chaumont-Frelet, T., Nicaise, S., Pardo, D.: Finite element approximation of electromagnetic fields using nonfitting meshes for geophysics. SIAM Journal on Numerical Analysis 56(4), 2288–2321 (2018). DOI 10.1137/16M1105566. URL https://doi.org/10.1137/16M1105566

15. Chaumont-Frelet, T., Valentin, F.: A multiscale hybrid-mixed method for the helmholtz equation in heterogeneous domains. SIAM Journal on Numerical Analysis 58(2), 1029–1067 (2020). DOI 10.1137/19M1255616. URL https://doi.org/10.1137/19M1255616

16. Cheng, C.H.: Crack models for a transversely isotropic medium. Journal of Geophysical Research: Solid Earth 98(B1), 675–684 (1993). DOI 10.1029/92JB02118. URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92JB02118

17. Cioranescu, D., Donato, P.: An Introduction to Homogenization. Oxford lecture series in mathematics and its applications. Oxford University Press (1999)

18. Ciz, R., Shapiro, S.A.: Generalization of Gassmann equations for porous media saturated with a solid material. GEOPHYSICS 72(6), A75–A79 (2007). DOI 10.1190/1.2772400. URL https://doi.org/10.1190/1.2772400

19. Cnudde, V., Boone, M.: High-resolution x-ray computed tomography in geosciences: A review of the current technology and applications. Earth-Science Reviews 123, 1–17 (2013). DOI https://doi.org/10.1016/j.earscirev. 2013.04.003. URL http://www.sciencedirect.com/science/article/pii/S001282521300069X

20. Dvorkin, J., Nur, A.: Time-average equation revisited. GEOPHYSICS 63(2), 460–464 (1998). DOI

10.1190/1.1444347. URL https://doi.org/10.1190/1.1444347

21. Efendiev, Y., Galvis, J., Li, G., Presho, M.: Generalized multiscale finite element methods: Oversampling strategies. International Journal for Multiscale Computational Engineering 12(6), 465–484 (2014). DOI 10.1615/IntJMultCompEng.2014007646

22. Eshelby, J.D., Peierls, R.E.: The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 241(1226), 376–396 (1957). DOI 10.1098/rspa.1957.0133. URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1957.0133

23. Francfort, G.A., Murat, F.: Homogenization and optimal bounds in linear elasticity. Archive for Rational Mechanics and Analysis 94(4), 307–334 (1986). DOI 10.1007/BF00280908. URL https://doi.org/10.1007/BF00280908

24. Gardner, G.H.F., Gardner, L.W., Gregory, A.R.: Formation velocity and density-the diagnostic basics for stratigraphic traps. GEOPHYSICS 39(6), 770–780 (1974). DOI 10.1190/1.1440465. URL https://doi.org/10.1190/1.1440465

25. Gassmann, F.: ber die elastizitt porser medien (on elasticity of porous media). Veirteljahrsschrift der Naturforschenden Gesellschaft in Zrich 96, 1–23 (1951)

26. Geers, M., Kouznetsova, V., Brekelmans, W.: Multi-scale computational homogenization: Trends and challenges. Journal of Computational and Applied Mathematics 234(7), 2175–2182 (2010). DOI https://doi.org/10.1016/j.cam.2009.08.077. URL http://www.sciencedirect.com/science/article/pii/S0377042709005536. Fourth International Conference on Advanced Computational Methods in ENgineering (ACOMEN 2008)

27. Grechka, V., Kachanov, M.: Effective elasticity of rocks with closely spaced and intersecting cracks. GEOPHYSICS 71(3), D85–D91 (2006). DOI 10.1190/1.2197489. URL https://doi.org/10.1190/1.2197489

28. Hashin, Z., Shtrikman, S.: On some variational principles in anisotropic and nonhomogeneous elasticity. Journal of the Mechanics and Physics of Solids 10(4), 335–342 (1962). DOI https://doi.org/10.1016/0022-5096(62)90004-2. URL http://www.sciencedirect.com/science/article/pii/0022509662900042

29. Hashin, Z., Shtrikman, S.: A variational approach to the theory of the elastic behaviour of multi-
phase materials. Journal of the Mechanics and Physics of Solids **11**(2), 127-140 (1963). DOI https://doi.org/10.1016/0022-5096(63)90060-7.
URL http://www.sciencedirect.com/science/article/pii/0022509663900607

30. Helbig, K.: Anisotropy and dispersion in periodically layered media. Geophysics **49**(4), 364-373 (1984)

31. Hill, R.: The elastic behaviour of a crystalline aggregate. Proceedings of the Physical Society. Section A **65**(5), 349-354 (1952). DOI 10.1088/0370-1298/65/5/307

32. Hovem, J.: Acoustic waves in finely layered media. Geophysics **60**(4), 1217-1221 (1995). DOI 10.1190/1.1443850. URL https://doi.org/10.1190/1.1443850

33. Hudson, J.A.: Wave speeds and attenuation of elastic waves in material containing cracks. Geophysical Journal International **64**(1), 133-150 (1981). DOI 10.1111/j.1365-246X.1981.tb02662.x. URL https://dx.doi.org/10.1111/j.1365-246X.1981.tb02662.x

34. Kuster, G.T., Toksz, M.N.: Velocity and attenuation of seismic waves in twophase media: Part i. theoretical formulations. GEOPHYSICS **39**(5), 587-606 (1974). DOI 10.1190/1.1440450. URL https://doi.org/10.1190/1.1440450

35. Levin, F.K.: Seismic velocities in transversely isotropic media. GEOPHYSICS **44**(5), 918-936 (1979). DOI 10.1190/1.1440985. URL https://doi.org/10.1190/1.1440985

36. Marion, D., Coudin, P.: Fram ray to effective medium theories in stratified media: An experimental study. In: SEG Technical Program Expanded Abstracts 1992, pp. 1341-1343. Society of Exploration Geophysicists (1992)

37. Marion, D., Mukerji, T., Mavko, G.: Scale effects on velocity dispersion: From ray to effective medium theories in stratified media. Geophysics **59**(10), 1613-1619 (1994)

38. Matou, K., Geers, M.G., Kouznetsova, V.G., Gillman, A.: A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. Journal of Computational Physics **330**, 192-220 (2017). DOI https://doi.org/10.1016/j.jcp.2016.10.070. URL http://www.sciencedirect.com/science/article/pii/S0021999116305782

39. Matuszyk, P.J., Torres-Verdín, C., Pardo, D.: Frequency-domain finite-element simulations of 2D sonic wireline borehole measurements acquired in fractured and thinly bedded formations. GEOPHYSICS **78**(4), D193-D207 (2013). DOI 10.1190/geo2012-0397.1. URL https://doi.org/10.1190/geo2012-0397.1

40. Mavko, G., Mukerji, T., Dvorkin, J.: The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media, 2nd edn. Cambridge University Press, Cambridge (2009). URL http://www.worldcat.org/isbn/0521861365

41. Modesto, D., Zlotnik, S., Huerta, A.: Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation. Computer Methods in Applied Mechanics and Engineering **295**, 127-149 (2015). DOI https://doi.org/10.1016/j.cma.2015.03.026. URL http://www.sciencedirect.com/science/article/pii/S004578251500136X

42. Müller, T.M., Caspari, E., Qi, Q., Rubino, J.G., Velis, D., Lopes, S., Lebedev, M., Gurevich, B.: Chapter 3 - acoustics of partially saturated rocks: Theory and experiments. In: J. Ba, Q. Du, J.M. Carcione, H.J. Zhang, T.M. Müller (eds.) Seismic Exploration of Hydrocarbons in Heterogeneous Reservoirs, pp. 45-75. Elsevier (2015). DOI https://doi.org/10.1016/B978-0-12-420151-4.00003-9. URL http://www.sciencedirect.com/science/article/pii/B9780124201514000039

43. Norris, A.: A differential scheme for the effective moduli of composites. Mechanics of Materials **4**(1), 1-16 (1985). DOI https://doi.org/10.1016/0167-6636(85)90002-X.
URL http://www.sciencedirect.com/science/article/pii/016766368590002X

44. Oleinik, O., Shamaev, A., Yosifian, G.: Mathematical Problems in Elasticity and Homogenization. ISSN. Elsevier Science (2009)

45. Osborne, M., Smyth, G.: A modified prony algorithm for exponential function fitting. SIAM Journal on Scientific Computing **16**(1), 119-138 (1995)

46. Prony, C., Gaspard, F., De-Baron, R.: Essai experimental et analytique. Journal de l'Ecole Polytechnique de Paris **1**, 24-76 (1795)

47. Rawlinson, N., Hauser, J., M., S.: Seismic ray tracing and wavefront tracking in laterally heterogeneous media. Advances in Geophysics **49**, 203-273 (2008). DOI https://doi.org/10.1016/S0065-2687(07)49003-3

48. Raymer, L.L., Hunt, E.R., Gardner, J.S.: An improved sonic transit time-to-porosity transform. Society of Petrophysicists and Well-Log Analysts pp. 1-12 (1980)

49. Reuss, A.: Berechnung der fliegrenze von mischkristallen auf grund der plastizittsbedingung fr einkristalle. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift fr Angewandte
Mathematik und Mechanik **9**(1), 49–58 (1929).
DOI 10.1002/zamm.19290090104

50. Romkes, A., Oden, J.T.: Adaptive modeling of wave propagation in heterogeneous elastic solids. Computer Methods in Applied Mechanics and Engineering **193**(6-8), 539–559 (2004). DOI 10.1016/j.cma.2003.10.014

51. Rouy, E., Tourin, A.: A viscosity solutions approach to shape-from-shading. SIAM J. Num. Anal. **29**, 867–884 (1992). DOI https://doi.org/10.1137/0729053

52. Saleh, A.A., Castagna, J.P.: Revisiting the wyllie time average equation in the case of nearspherical pores. GEOPHYSICS **69**(1), 45–55 (2004). DOI 10.1190/1.1649374. URL https://doi.org/10.1190/1.1649374

53. Schoenberg, M., Muir, F.: A calculus for finely layered anisotropic media. GEOPHYSICS **54**(5), 581–589 (1989). DOI 10.1190/1.1442685. URL https://doi.org/10.1190/1.1442685

54. Schot, S.H.: Eighty years of sommerfeld's radiation condition. Historia Mathematica **19**(4), 385 – 401 (1992). DOI https://doi.org/10.1016/0315-0860(92)90004-U. URL http://www.sciencedirect.com/science/article/pii/031508609290004U

55. Sethian, J.A.: Evolution, implementation, and application of level set and fast marching methods for advancing fronts. J. Comp. Phys. **169**, 503–555 (2001). DOI https://doi.org/10.1006/jcph.2000.6657

56. Sethian, J.A., Popovici, A.M.: 3-d traveltime computation using the fast marching method. Geophysics **64**, 516–523 (1999). DOI https://doi.org/10.1190/1.1444558

57. Sommerfeld, A.: Die greensche funktion der schwingungslgleichung. Jahresbericht der Deutschen Mathematiker-Vereinigung **21**, 309 – 352 (1912). URL http://eudml.org/doc/145344

58. Stovas, A., Arntsen, B.: Low frequency waves in finely layered media. In: 65th EAGE Conference & Exhibition (2003)

59. Stovas, A., Ursin, B.: Equivalent time-average and effective medium for periodic layers. Geophysical Prospecting **55**(6), 871–882 (2007). DOI 10.1111/j.1365-2478.2007.00653.x. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2478.2007.00653.x

60. Thomsen, L.: Elastic anisotropy due to aligned cracks in porous rock1. Geophysical Prospecting **43**(6), 805–829 (1995). DOI 10.1111/j.1365-2478.1995.tb00282.x. URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2478.1995.tb00282.x

61. Voigt, W.: Ueber die beziehung zwischen den beiden elasticittsconstanten isotroper krper. Annalen der Physik **274**(12), 573–587 (1889). DOI 10.1002/andp.18892741206

62. Wyllie, M.R.J., Gregory, A.R., Gardner, L.W.: Elastic wave velocities in heterogeneous and porous media. GEOPHYSICS **21**(1), 41–70 (1956). DOI 10.1190/1.1438217. URL https://doi.org/10.1190/1.1438217