A numerical extension of White's theory of P-wave attenuation to non-isothermal poroelastic media.

Mesoscopic P-wave attenuation in layered, partially saturated thermo-poroelastic media is analyzed by combining the theories of Biot poroelasticity and Lord-Shulman thermoelasticity (BLS). The attenuation is quantified by estimating the quality factor Q. The mesoscopic attenuation effect, commonly referred to as wave-induced fluid flow (WIFF), is the process that converts fast compressional and shear waves into slow diffusive Biot waves at mesoscopic heterogeneities larger than the pore scale, but much smaller than the dominant wavelengths. This effect was first modeled in White's isothermal theory by quantifying the seismic response of a periodic sequence of planar porous layers that are alternately saturated with gas or water. This work presents a numerical extension of White's theory for the non-isothermal case in this type of sequence. For this purpose, an initial-boundary-value problem (IBVP) for the BLS wave propagation equations is solved using the finite element method, where the particle velocity field is recorded at uniformly distributed receivers. The quality factor is estimated using spectral-ratio and frequency-shift methods. The Q-estimates show that thermal effects influence the attenuation of the P-wave and the velocity dispersion compared to the isothermal case.

P- and S-wave simulation using a Cole-Cole model to incorporate thermoelastic attenuation and dispersion.

In thermoelastic wave attenuation, such as that caused by heterogeneities much smaller than the wavelength, e.g., Savage theory of spherical pores, the shape of the relaxation peak differs from that of the Zener (or standard linear solid) mechanical model. In these effective homogeneous media, the anelastic behavior is better represented by a stress-strain relation based on fractional derivatives; particularly, P- and S-wave dispersion and attenuation is well described by a Cole-Cole equation. We propose a time-domain algorithm for wave propagation based on the Grünwald-Letnikov numerical derivative and the Fourier pseudospectral method to compute the spatial derivatives. As an example, we consider Savage theory and verify the algorithm by comparison with the analytical solution in homogeneous media based on the frequency-domain Green function. Moreover, we illustrate the modeling performance with wave propagation in a two half-space medium where one section is lossless and the other is a Cole-Cole medium. This apparently simple example, which does not have an analytical solution, shows the complexity of the wavefield that characterizes a single flat interface.

Finite-element numerical simulations of seismic attenuation in finely layered rocks.

P-wave conversion to slow diffusion (Biot) modes at mesoscopic (small-scale) inhomogeneities in porous media is believed to be the most important attenuation mechanisms at seismic frequencies. This study considers three periodic thin layers saturated with gas, oil, and water, respectively, a realistic scenario in hydrocarbon reservoirs, and perform finite-element numerical simulations to obtain the wave velocities and quality factors along the direction perpendicular to layering. The results are validated by comparison to the Norris-Cavallini analytical solution, constituting a cross-check for both theory and numerical simulations. The approach is not restricted to partial saturation but also applies to relevant properties in reservoir geophysics, such as porosity and permeability heterogeneities. This paper considers two cases, namely, the same rock skeleton and different fluids, and the same fluid and different dry-rock properties. Unlike the two-layer case (two fluids), the results show two relaxation peaks and the agreement between numerical and analytical solutions is excellent.

Waves at a fluid-solid interface: Explicit versus implicit formulation of boundary conditions using a discontinuous Galerkin method.

An accurate solution of the wave equation at a fluid-solid interface requires a correct implementation of the boundary condition. Boundary conditions at fluid-solid interface require continuity of the normal component of particle velocity and traction, whereas the tangential components vanish. A main challenge is to model interface waves, namely, the Scholte and leaky Rayleigh waves. This study uses a nodal discontinuous Galerkin (dG) finite-element method with the medium discretized using an unstructured uniform triangular meshes. The natural boundary conditions in the dG method are implemented by (1) using an explicit upwind numerical flux and (2) by using an implicit penalty flux and setting the modulus of rigidity of the acoustic medium to zero. The accuracy of these methods is evaluated by comparing the numerical solutions with analytical ones, with source and receiver at and away from the interface. The study shows that the solutions obtained from the explicit and implicit boundary conditions provide the correct results. The stability of the dG scheme is determined by the numerical flux, which also implements the boundary conditions by unifying the numerical solution at shared edges of the elements in an energy stable manner.

A Simulation of a COVID-19 Epidemic Based on a Deterministic SEIR Model.

An epidemic disease caused by a new coronavirus has spread in Northern Italy with a strong contagion rate. We implement an SEIR model to compute the infected population and the number of casualties of this epidemic. The example may ideally regard the situation in the Italian Region of Lombardy, where the epidemic started on February 24, but by no means attempts to perform a rigorous case study in view of the lack of suitable data and the uncertainty of the different parameters, namely, the variation of the degree of home isolation and social distancing as a function of time, the initial number of exposed individuals and infected people, the incubation and infectious periods, and the fatality rate. First, we perform an analysis of the results of the model by varying the parameters and initial conditions (in order for the epidemic to start, there should be at least one exposed or one infectious human). Then, we consider the Lombardy case and calibrate the model with the number of dead individuals to date (May 5, 2020) and constrain the parameters on the basis of values reported in the literature. The peak occurs at day 37 (March 31) approximately, with a reproduction ratio R0 of 3 initially, 1.36 at day 22, and 0.8 after day 35, indicating different degrees of lockdown. The predicted death toll is approximately 15,600 casualties, with 2.7 million infected individuals at the end of the epidemic. The incubation period providing a better fit to the dead individuals is 4.25 days, and the infectious period is 4 days, with a fatality rate of 0.00144/day [values based on the reported (official) number of casualties]. The infection fatality rate (IFR) is 0.57%, and it is 2.37% if twice the reported number of casualties is assumed. However, these rates depend on the initial number of exposed individuals. If approximately nine times more individuals are exposed, there are three times more infected people at the end of the epidemic and IFR = 0.47%. If we relax these constraints and use a wider range of lower and upper bounds for the incubation and infectious periods, we observe that a higher incubation period (13 vs. 4.25 days) gives the same IFR (0.6 vs. 0.57%), but nine times more exposed individuals in the first case. Other choices of the set of parameters also provide a good fit to the data, but some of the results may not be realistic. Therefore, an accurate determination of the fatality rate and characteristics of the epidemic is subject to knowledge of the precise bounds of the parameters. Besides the specific example, the analysis proposed in this work shows how isolation measures, social distancing, and knowledge of the diffusion conditions help us to understand the dynamics of the epidemic. Hence, it is important to quantify the process to verify the effectiveness of the lockdown.

Numerical simulation of wave-induced fluid flow seismic attenuation based on the Cole-Cole model.

The acoustic behavior of porous media can be simulated more realistically using a stress-strain relation based on the Cole-Cole model. In particular, seismic velocity dispersion and attenuation in porous rocks is well described by mesoscopic-loss models. Using the Zener model to simulate wave propagation is a rough approximation, while the Cole-Cole model provides an optimal description of the physics. Here, a time-domain algorithm is proposed based on the Grünwald-Letnikov numerical approximation of the fractional derivative involved in the time-domain representation of the Cole-Cole model, while the spatial derivatives are computed with the Fourier pseudospectral method. The numerical solution is successfully tested against an analytical solution. The methodology is applied to a model of saline aquifer, where carbon dioxide (CO2) is injected. To follow the migration of the gas and detect possible leakages, seismic monitoring surveys should be carried out periodically. To this aim, the sensitivity of the seismic method must be carefully assessed for the specific case. The simulated test considers a possible leakage in the overburden, above the caprock, where the sandstone is partially saturated with gas and brine. The numerical examples illustrate the implementation of the theory.

Synthetic waveforms of axial motion in a borehole with drill string.

This work studies the wave motion in a fluid-filled borehole in the presence of drill string and geological formation. The synthetic waveforms are obtained by a three-dimensional axis-symmetric full-wave numerical simulation in a two-dimensional multi-domain where the medium is uniform with respect to the azimuth. The discretization is performed in cylindrical coordinates. In order to simulate the waves at the origin (axis of the polar radius), a very small radius is used to avoid the singularity. The free-surface and rigid boundary conditions are tested and it is shown that the rigid one constitutes the best approximation. The simulations provide the amplitude distribution and motion diagrams in the borehole vertical cross-sections and at the outer boundary, away from the borehole. Propagation in the presence of hard and soft formations is analysed. The dispersion, the amplitude, and the orbital polarization of the modes excited by a point source acting in the fluid inside a drillstring are considered and examples of comparison with literature results obtained using multi-modal analysis are shown. The proposed approach is more general than the multi-modal analysis, since it allows for arbitrary variations of the properties in the plane of symmetry.

Reflection and transmission coefficients of a single layer in poroelastic media.

Wave propagation in poroelastic media is a subject that finds applications in many fields of research, from geophysics of the solid Earth to material science. In geophysics, seismic methods are based on the reflection and transmission of waves at interfaces or layers. It is a relevant canonical problem, which has not been solved in explicit form, i.e., the wave response of a single layer, involving three dissimilar media, where the properties of the media are described by Biot's theory. The displacement fields are recast in terms of potentials and the boundary conditions at the two interfaces impose continuity of the solid and fluid displacements, normal and shear stresses, and fluid pressure. The existence of critical angles is discussed. The results are verified by taking proper limits-zero and 100% porosity-by comparison to the canonical solutions corresponding to single-phase solid (elastic) media and fluid media, respectively, and the case where the layer thickness is zero, representing an interface separating two poroelastic half-spaces. As examples, it was calculated the reflection and transmission coefficients for plane wave incident at a highly permeable and compliant fluid-saturated porous layer, and the case where the media are saturated with the same fluid.

Elastic surface waves in crystals--part 2: cross-check of two full-wave numerical modeling methods.

We obtain the full-wave solution for the wave propagation at the surface of anisotropic media using two spectral numerical modeling algorithms. The simulations focus on media of cubic and hexagonal symmetries, for which the physics has been reviewed and clarified in a companion paper. Even in the case of homogeneous media, the solution requires the use of numerical methods because the analytical Green's function cannot be obtained in the whole space. The algorithms proposed here allow for a general material variability and the description of arbitrary crystal symmetry at each grid point of the numerical mesh. They are based on high-order spectral approximations of the wave field for computing the spatial derivatives. We test the algorithms by comparison to the analytical solution and obtain the wave field at different faces (stress-free surfaces) of apatite, zinc and copper. Finally, we perform simulations in heterogeneous media, where no analytical solution exists in general, showing that the modeling algorithms can handle large impedance variations at the interface.

Wave simulation in biologic media based on the Kelvin-Voigt fractional-derivative stress-strain relation.

The acoustic behavior of biologic media can be described more realistically using a stress-strain relation based on fractional time derivatives of the strain, since the fractional exponent is an additional fitting parameter. We consider a generalization of the Kelvin-Voigt rheology to the case of rational orders of differentiation, the so-called Kelvin-Voigt fractional-derivative (KVFD) constitutive equation, and introduce a novel modeling method to solve the wave equation by means of the Grünwald-Letnikov approximation and the staggered Fourier pseudospectral method to compute the spatial derivatives. The algorithm can handle complex geometries and general material-property variability. We verify the results by comparison with the analytical solution obtained for wave propagation in homogeneous media. Moreover, we illustrate the use of the algorithm by simulation of wave propagation in normal and cancerous breast tissue.

Elastic surface waves in crystals. Part 1: review of the physics.

We present a review of wave propagation at the surface of anisotropic media (crystal symmetries). The physics for media of cubic and hexagonal symmetries has been extensively studied based on analytical and semi-analytical methods. However, some controversies regarding surfaces waves and the use of different notations for the same modes require a review of the research done and a clarification of the terminology. In a companion paper we obtain the full-wave solution for the wave propagation at the surface of media with arbitrary symmetry (including cubic and hexagonal symmetries) using two spectral numerical modeling algorithms.

Elastic medium equivalent to Fresnel's double-refraction crystal.

In 1821, Fresnel obtained the wave surface of an optically biaxial crystal, assuming that light waves are vibrations of the ether in which longitudinal vibrations (P waves) do not propagate. An anisotropic elastic medium mathematically analogous to Fresnel's crystal exists. The medium has four elastic constants: a P-wave modulus, associated with a spherical P wave surface, and three elastic constants, c(44), c(55), and c(66), associated with the shear waves, which are mathematically equivalent to the three dielectric permittivity constants epsilon(11), epsilon(22), and epsilon(33) as follows: mu(0)epsilon(11)<==>rho/c(44), mu(0)epsilon(22)<==>rho/c(55), mu(0)epsilon(33)<==>rho/c(66), where mu(0) is the magnetic permeability of vacuum and rho is the mass density. These relations also represent the equivalence between the elastic and electromagnetic wave velocities along the principal axes of the medium. A complete mathematical equivalence can be obtained by setting the P-wave modulus equal to zero, but this yields an unstable elastic medium (the hypothetical ether). To obtain stability the P-wave velocity has to be assumed infinite (incompressibility). Another equivalent Fresnel's wave surface corresponds to a medium with anomalous polarization. This medium is physically unstable even for a nonzero P-wave modulus.