3D Heterogeneous Media Research Articles

An iterative solver for the 3D Helmholtz equation

We develop a frequency–domain iterative solver for numerical simulation of acoustic waves in 3D heterogeneous media. It is based on the application of a unique preconditioner to the Helmholtz equation that ensures convergence for Krylov subspace iteration methods. Effective inversion of the preconditioner involves the Fast Fourier Transform (FFT) and numerical solution of a series of boundary value problems for ordinary differential equations. Matrix-by-vector multiplication for iterative inversion of the preconditioned matrix involves inversion of the preconditioner and pointwise multiplication of grid functions. Our solver has been verified by benchmarking against exact solutions and a time-domain solver.

Residual-driven online multiscale methods for acoustic-wave propagation in 2D heterogeneous media

Common applications, such as geophysical exploration, reservoir characterization, and earthquake quantification, in modeling and inversion aim to apply numerical simulations of elastic- or acoustic-wave propagation to increasingly large and complex models, which can provide more realistic and useful results. However, the computational cost of these simulations increases rapidly, which makes them inapplicable to certain problems. We apply a newly developed multiscale finite-element algorithm, the generalized multiscale finite-element method (GMsFEM), to address this challenge in simulating acoustic-wave propagation in heterogeneous media. The wave equation is solved on a coarse grid using multiscale basis functions that are chosen from the most dominant modes among those computed by solving relevant local problems on a fine-grid representation of the model. These multiscale basis functions are computed once in an off-line stage prior to the simulation of wave propagation. Because these calculations are localized to individual coarse cells, one can improve the accuracy of multiscale methods by revising and updating these basis functions during the simulation. These updated bases are referred to as online basis functions. This is a significant extension of previous applications of similar online basis functions to time-independent problems. We tested our new algorithm and numerical results for acoustic-wave propagation using the acoustic Marmousi model. Long-term developments have a strong potential to enhance inversion algorithms because the basis functions need not be regenerated everywhere. In particular, recomputation of basis functions is required only at regions in which the model is updated. Thus, our method allows faster simulations for repeated calculations, which are needed for inversion purpose.

A fast marching algorithm for the factored eikonal equation

The eikonal equation is instrumental in many applications in several fields ranging from computer vision to geoscience. This equation can be efficiently solved using the iterative Fast Sweeping (FS) methods and the direct Fast Marching (FM) methods. However, when used for a point source, the original eikonal equation is known to yield inaccurate numerical solutions, because of a singularity at the source. In this case, the factored eikonal equation is often preferred, and is known to yield a more accurate numerical solution. One application that requires the solution of the eikonal equation for point sources is travel time tomography. This inverse problem may be formulated using the eikonal equation as a forward problem. While this problem has been solved using FS in the past, the more recent choice for applying it involves FM methods because of the efficiency in which sensitivities can be obtained using them. However, while several FS methods are available for solving the factored equation, the FM method is available only for the original eikonal equation.In this paper we develop a Fast Marching algorithm for the factored eikonal equation, using both first and second order finite-difference schemes. Our algorithm follows the same lines as the original FM algorithm and requires the same computational effort. In addition, we show how to obtain sensitivities using this FM method and apply travel time tomography, formulated as an inverse factored eikonal equation. Numerical results in two and three dimensions show that our algorithm solves the factored eikonal equation efficiently, and demonstrate the achieved accuracy for computing the travel time. We also demonstrate a recovery of a 2D and 3D heterogeneous medium by travel time tomography using the eikonal equation for forward modeling and inversion by Gauss–Newton.

Experimental evaluation of the applicability of phase, amplitude, and combined methods to determine water flux and thermal diffusivity from temperature time series using VFLUX 2

• Recent McCallum et al., 2012 , Luce et al., 2013 methods added to VFLUX. • Recent methods are tested with non-steady flows, and in 3D heterogeneous media. • Spatial or temporal changes in thermal diffusivity output suggests errors in flux. • Amplitude ratio, or combined amplitude ratio and phase shift methods recommended. Vertical fluid exchange between surface water and groundwater can be estimated using diurnal signals from temperature time series methods based on amplitude ratios ( A r ), phase shifts (Δ ϕ ), or combined use of both ( A r Δ ϕ ). The A r , Δ ϕ , and A r Δ ϕ methods are typically applied in conditions where one or more of their underlying assumptions are violated, and the reliability of the various methods in response to non-ideal conditions is unclear. Additionally, A r Δ ϕ methods offer the ability to estimate thermal diffusivity ( κ e ) without assuming any thermal parameters, although the value of such output has not been broadly tested. The A r , Δ ϕ , and A r Δ ϕ methods are tested under non-steady, 1D flows in sand column experiments, and multi-dimensional flows in heterogeneous media in numerical modeling experiments. Results show that, in non-steady flow conditions, estimated κ e values outside of a plausible range for streambed materials (0.028–0.180 m 2 d −1 ) coincide with time periods with erroneous flux estimates. In heterogeneous media, sudden changes of κ e with depth also coincide with erroneous flux estimates. When (known) fluxes are variable in time, poor identification of Δ ϕ leads to poor flux estimates from Δ ϕ and A r Δ ϕ methods. However, when fluxes are steady, or near zero, A r Δ ϕ methods provide the most accurate flux estimates. This comparison of A r , Δ ϕ and A r Δ ϕ methods under non-ideal conditions provides guidance on their use. In this study, A r Δ ϕ methods have been coded into a new version of VFLUX, allowing users easy access to recent advances in heat tracing.

Finite-difference modeling of SH-wave conversions in shallow shear-wave refraction surveying

The shallow shear-wave refraction method works successfully in an area with a series of horizontal layers. Complex near-surface geology, however, may not fit into the assumption of a series of horizontal layers. It is theoretically inevitable that a plane SH-wave undergoes wave-type conversions along an interface in an area of non-horizontal layers. One real example has shown that the shallow SH-wave refraction method provides velocities of a converted wave rather than SH-wave. Moreover, it is impossible to identify the converted wave by refraction data itself. In this paper, we implement numerical simulation for conversion of SH- to P-wave in 3D heterogeneous medium with the finite-difference method. An SH-wave source excitation method that we give in the numerical simulation is testified, which can only generate SH-wave without P-wave. The numerical modeling results demonstrate that the conversion of the SH-wave to other wave-types will occur in an area of non-horizontal layers. All the converted P-wave arrivals are shown reversed polarity like S-wave arrivals in the modeling of reverse of the source and we have clarified the peculiar properties of converted P-waves from the S-wave. Our numerical simulation results confirm that velocities calculated from an SH-wave refraction survey are velocities of converted waves. Therefore, special attention should be paid to this pitfall in the real world.

Ground-Motion Simulation in the Lower Tagus Valley Basin

Throughout history, the Lower Tagus Valley (LTV) region has been shaken by several earthquakes, including some with moderate to large magnitudes and with sources located inside the basin, for example the 1344 (M6.0) and 1909 (M6.0) earthquakes. Previous simulations (Bezzeghoud et al. Natural Hazard 69: 1229–1245, 2011) have revealed strong amplification of the amplitude waves in the region, because of the effect of the low-velocity sediments that fill the basin. The model used in those simulations was updated in this work by including new high-resolution geophysical and geotechnical data available for the area (seismic reflection, aeromagnetic, gravimetric, deep wells, standard penetration tests, and geological data). To contribute to improved assessment of seismic hazard in the LTV, we simulated propagation of seismic waves produced by moderate earthquakes in a 3D heterogeneous medium by using elastic finite-difference wave propagation code. The method, successfully used by Grandin et al. (Geophys J Int 171: 1144–1161, 2007), involves evaluation of the seismic potential of known faults in the area studied and three-dimensional seismic ground motion modelling by use of finite difference methods. On the basis of this methodology, we calculated the ground motion for the April 23, 1909, Benavente (Portugal) earthquake (Mw = 6.0) in dense grid points, then computed the synthetic isoseismic map of the area by use of appropriate relationships between seismic intensity (MMI) and peak ground velocity (PGV). The synthetic results, in contrast with available macroseismic and instrumental data, enable validation of the source models proposed for the area, identification of the sources of historical earthquakes, and could also indicate which areas are more exposed to seismic ground motion.

Effect of geometry on concentration polarization in realistic heterogeneous permselective systems.

This study extends previous analytical solutions of concentration polarization occurring solely in the depleted region, to the more realistic geometry consisting of a three-dimensional (3D) heterogeneous ion-permselective medium connecting two opposite microchambers (i.e., a three-layer system). Under the local electroneutrality approximation, the separation of variable methods is used to derive an analytical solution of the electrodiffusive problem for the two opposing asymmetric microchambers. The assumption of an ideal permselective medium allows for the analytic calculation of the 3D concentration and electric potential distributions as well as a current-voltage relation. It is shown that any asymmetry in the microchamber geometries will result in current rectification. Moreover, it is demonstrated that for non-negligible microchamber resistances, the conductance does not exhibit the expected saturation at low concentrations but instead shows a continuous decrease. The results are intended to facilitate a more direct comparison between theory and experiments, as now the voltage drop is across a realistic 3D and three-layer system.

Finite-difference approximation of wave equation: a study case of the SIMA velocity model

Synthetic seismograms enable to model the theoretical seismic response of the Earth interior due to different structural features and changes in the physical properties of crust and mantle. This approximation provides a best understanding of the real seismic data recorded in field experiments. In this paper, we are showing the development and application of a new scheme based on a multi-order explicit finite-difference algorithm for acoustic waves in a 2D heterogeneous media. The results of the modeling are compared with the seismic data acquired within the SIMA project providing new insight about the internal structure of the subsurface allowing improving the velocity model obtained in previous works.

Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media

Modeling and imaging techniques for geophysics are extremely demanding in terms of computational resources. Seismic data attempt to resolve smaller scales and deeper targets in increasingly more complex geologic settings. Finite elements enable accurate simulation of time-dependent wave propagation in heterogeneous media. They are more costly than finite-difference methods, but this is compensated by their superior accuracy if the finite-element mesh follows the sharp impedance contrasts and by their improved efficiency if the element size scales with wavelength, hence with the local wave velocity. However, 3D complex geologic settings often contain details on a very small scale compared to the dominant wavelength, requiring the mesh to contain elements that are smaller than dictated by the wavelength. Also, limitations of the mesh generation software may produce regions where the elements are much smaller than desired. In both cases, this leads to a reduction of the time step required to solve the wave propagation and significantly increases the computational cost. Local time stepping (LTS) can improve the computational efficiency and speed up the simulation. We evaluated a local formulation of an LTS scheme with second-order accuracy for the discontinuous Galerkin finite-element discretization of the wave equation. We tested the benefits of the scheme by considering a geologic model for a North-Sea-type example.

Direct and inverse problems in electromagnetic sounding of three-dimensional heterogeneous medium

The methods for solving three-dimensional (3D) direct and inverse problems of electromagnetic sounding are considered. It is shown that the method of integral equations is an efficient instrument for mathematical modeling of electromagnetic fields in 3D heterogeneous media. Adaptation of the integral equation technique to the solution of inverse problems is demonstrated.

Scale Dependence of Effective Hydraulic Conductivity Distributions in 3D Heterogeneous Media: A Numerical Study

Upscaling procedures and determination of effective properties are of major importance for the description of flow in heterogeneous porous media. In this context, we study the statistical properties of effective hydraulic conductivity (Keff) distributions and their dependence on the coarsening scale. First, we focus on lognormal stationary isotropic media. Our results suggest that Keff is lognormally distributed independently on the coarsening scale. The scale dependence of the mean and variance of Keff are in agreement with recent analytical derivations obtained using coarse graining filtering techniques. In the second part, we focus on binary media, analysing the dependence of Keff distributions on the coarsening scale and also on the high-K facies volume fraction p. When p is near the percolation threshold pc, the decrease of the normalized variance with the coarsening scale is remarkably (102 times) slower compared to the situation in which p far from pc, but also compared to the cases of lognormal media studied before. This result permits to assess the degree of difficulty that systems with p near pc pose for upscaling procedures. Also we point out in terms of Keff statistics the relative influence of the coarsening scale and of the high-K facies connectivity.

Toward a realistic deformation model of the 2008 magmatic intrusion at Etna from combined DInSAR and GPS observations

The combination of (i) DInSAR data, capable of observing deformation pattern at a spatial resolution unachievable with other sparse geodetic measurements, (ii) continuous GPS data, able to provide temporal constraints on source evolution, and (iii) numerical modeling procedures, appropriate to consider a non-uniform opening distribution of a source embedded in a 3D heterogeneous medium, allowed us to infer a complex and realistic deformation model of the magmatic intrusion that occurred in the northern flank of Etna on 13 May 2008. Numerical modeling of ground deformation data defines a near-vertical dyke intruded for 2.5 km starting from a depth of 1400 m asl right below the summit craters and reaching shallow crust level in the northern flank. From the estimated opening distribution of the propagating magma-filled crack, which reached a maximum value of about 2 m, a volumetric expansion of crustal rocks of about 5.3 × 10 6 m 3 was obtained. Also, we clarified the temporal evolution of the northward magmatic intrusion, which lasted just over 5 h with an initial magma propagation velocity of about 1.2 km/h, and decreased to about 0.24 km/h as the driving pressure lowered due to the effusive activity started at southern vents.

Simulation of seismic wave propagation in 3D heterogeneous media: a parallel computing approach

The use of parallel computing makes it feasible to simulate realistic seismological events, whose reconstruction requires wide domains, high frequencies and the introduc- tion of the dissipation terms. The propagation problem of seismic waves is a key feature of the earthquake dynamics that we are interested in numerically modeling and simu- lating. In particular, in this work we present several preliminary results about the load balancing for the parallel resolution of a simulation of the propagation of seismic waves in a 3D heterogeneous medium. The Finite Element Method is employed for the spatial discretization by using non–structured tetrahedral meshes. The Newmark method is used for the time discretization. With the aim to study a priori the load balancing, we intro- duce two performance indices: closing nodes and node balancing. In particular, the first one estimates the amount of processor data, the latter provides information on work–load distribution. The variation of these indices as functions of the number of processors and of the number of nodes of the grid is then investigated.

Coupled simulation of near-wellbore and reservoir models

It is well known that some near-wellbore physics such as formation damage, acid stimulation or near-well reservoir heterogeneity can have a big impact on well productivity. To evaluate this impact, near-wellbore numerical models have been developed to simulate near-wellbore physics using very fine gridblocks. However, near-well models are usually developed standalone and are decoupled from reservoir models. Using a standalone near-well model, which does not take into account production scenarios, does not allow us to correctly predict well injectivity or productivity. In this paper, we present a technique for the coupled simulation of a near-well flow model and a reservoir model in a simple and consistent way. This approach is applied to multiphase flow in 3D heterogeneous media. In this approach, data between these two models are updated through boundary conditions for the near-well model and numerical productivity indices (PI) for the reservoir model. Examples, including an example of two-phase flow in heterogeneous media, simulation of drilling induced formation damage in multi-layered reservoir and water shut-off treatment for a horizontal well, show that this coupled modelling gives quite satisfactory results.

An efficient Monte Carlo-SubSampling approach to first passage problems

In the realm of the Continuous Time Random Walk (CTRW) and in conjunction with the Monte Carlo (MC) approach, we consider the transport of a chemical or radioactive pollutant in a 3D heterogeneous medium, focusing on the first passage time (FPT), defined as the time required by the walkers representative of the dangerous particles to travel within the medium until crossing a target disk, thus entering another medium which should instead remain clean, such as a water well or an aquifer. The advantage of the MC approach is the possibility of simulating different features of the travel such as different waiting-time probabilities, space dependent jump lengths, absorption and desorption phenomena, Galilei invariant and variant velocities driven by external forces. When the computer time required for collecting a suitable number of target crossings is excessive, we propound to hybridize the MC approach with the recent SubSampling computational procedure, usefully applied in the engineering reliability field to computing very small failure probabilities in short computer time. To tackle the FPT problem we iteratively consider groups of a few thousands of walkers: in each iteration we select a fraction of them closer to the target, ignoring the remaining ones, and then restore the group by creating with the MC technique new walkers even more close to the target. The successive groups have large conditional probabilities which can be estimated one after the other in short time and whose product yields the total probability. By so doing the FPT can be actually computed in much shorter times: we report examples of Gaussian and anomalous distributions in which the reduction with respect to a pure MC computation is of orders of magnitude.

Simulation of acoustic wavefields in heterogeneous media: A robust method for automatic suppression of numerical dispersion

Numerical dispersion limits the application of numerical simulation methods for solving the acoustic wave equation in large-scale computation. The nearly analytic discrete method (NADM) and its improved version for suppressing numerical dispersion were developed recently. This new method is a refinement of two previous methods and further increases the ability of suppressing numerical dispersion for modeling acoustic wave propagation in 2D heterogeneous media, which uses acoustic wave displacement, particle velocity, and their gradients to restructure the acoustic wavefield via the truncated Taylor expansion and the high-order interpolation approximate method. For the method proposed, we investigate its implementation and compare it with the higher-order Lax-Wendroff correction (LWC) scheme, the original nearly analytic discrete method (NADM) and its im-proved version with regard to numerical dispersion, computational costs, and computer storage requirements. The numerical dispersion relations provided by the refined algorithm for 1D and 2D cases are analyzed, as well as the numerical results obtained by this method against the exact solution for the 2D acoustic case. Numerical results show that the refined method gives no visible numerical dispersion for very large spatial grid increments. It can simulate high-frequency wave propagation for a given grid interval and automatically suppress the numerical dispersion when the acoustic wave equation is discretized, when too few samples per wavelength are used, or when models have large velocity contrasts. Unlike the high-order LWC methods, our present method can save substantial computational costs and memory requirements because very large grid increments can be used. The refined method can be used for the simulation of large-scale wavefields.

Segmentally Iterative Ray Tracing in Complex 2D and 3D Heterogeneous Block Models

Abstract We describe a complex geologic model as an aggregate of arbitrarily shaped blocks separated by cubic splines in 2D and triangulated interfaces in 3D. Recently we have introduced a segmentally iterative ray-tracing (SIRT) method based on Fermat’s principle of stationary travel time, which has been documented to be robust and fast for a complex block model with a constant velocity defined in each block. In this work, we extend the constant velocity to a generally continuous distribution with an analytical expression of travel time, and develop SIRT in the redefined velocity distribution. As a three-point perturbation scheme, SIRT requires an explicit analytical travel time between two intersection points expressed as a function of coordinates of the two points. In these situations, we derive a general midpoint perturbation formula, and further a detailed perturbation formula for familiar media with a constant velocity gradient. SIRT is a scheme in which we perturb the intersection points of an initial-guess ray path in sequence by the first-order explicit formulas instead of using traditional iterative methods. A key consideration is the fact that the number of intersection points may be variable during the iteration process. Numerical tests demonstrate that SIRT is effective in implementing kinematic two-point ray tracing in complex 2D and 3D heterogeneous media.

Integrated reservoir, petrophysical, and seismic simulation of CO2 storage in coal beds

In this paper, we present results of an integrated simulation study to assess the effectiveness of seismic imaging for monitoring [Formula: see text] sequestration. We considered two scenarios. In the first, injected [Formula: see text] remained confined within a shallow coal formation. In the second, the sequestered [Formula: see text] gas leaked through a semi-permeable shale layer to an overlying sand unit sealed above. In both scenarios, the [Formula: see text] injection process resulted in enhanced production of coal-bed methane. The reservoir and seismic simulations required the construction of 3D geologic and facies models, the estimations of seismic velocities based on rock physics correlations, and the development of geostatistical dual-porosity reservoir descriptions of the coal and overlying shale and sand units. We ran the 3D reservoir simulations with GEM 2008.10, an equation-of-state compositional reservoir simulator from Computer Modeling Group with the capability to model adsorption of injected [Formula: see text] in coal, desorption of [Formula: see text], and matrix shrinkage-swelling effects. The reservoir description was derived from a 3D geostatistical model of the Big George coal in the Powder River Basin. We used Batzle-Wang, Gassmann, and Mavko-Dvorkin-Mukerji relationships to incorporate the effects of mineral contents and changing reservoir properties on P-wave velocity. We built synthetic seismograms by simulating the propagation of acoustic waves through the 3D heterogeneous media, and used reverse time migration to create 3D seismic images corresponding to the state of the reservoir at the beginning and end of ten years of [Formula: see text] injection. The resulting seismic images clearly identified the regions of [Formula: see text] gas saturation, closely matching the gas-saturation profiles predicted by the reservoir simulator. Further research will be necessary to fine-tune an integrated reservoir simulation, rock physics modeling, and seismic imaging workflow to provide highly advanced imaging, monitoring, and verification capabilities for [Formula: see text] sequestration. This study is a first step toward this goal.

Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy

Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy

Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy

Employing a 3D finite element method, we develop an algorithm to calculate gravity changes due to pressurized sources of any shape in elastic and inelastic heterogeneous media. We consider different source models, such as sphere, spheroid and sill, dilating in elastic media (homogeneous and heterogeneous) and in elasto-plastic media. The models are oriented to reproduce the gravity changes and the surface deformation observed at Campi Flegrei caldera (Italy), during the 1982-1984 unrest episode. The source shape and the characteristics of the medium have great influence on the calculated gravity changes, leading to very different values for the source densities. Indeed, the gravity residual strongly depends upon the shape of the source. Non negligible contributions also come from density and rigidity heterogeneities within the medium. Furthermore, if the caldera is elasto- plastic, the resulting gravity changes exhibit a pattern similar to that provided by a low effective rigidity. Even if the variation of the source volumes is quite similar for most of the models considered, the density inferred for the source ranges from ?400 kg/m3 (super critical water) to ?3300 kg/m3 (higher than trachytic basalts), with drastically different implications for risk assessment.