3D Heterogeneous Media Research Articles

An efficient probabilistic workflow for estimating induced earthquake parameters in 3D heterogeneous media

Abstract. We present an efficient probabilistic workflow for the estimation of source parameters of induced seismic events in three-dimensional heterogeneous media. Our workflow exploits a linearized variant of the Hamiltonian Monte Carlo (HMC) algorithm. Compared to traditional Markov chain Monte Carlo (MCMC) algorithms, HMC is highly efficient in sampling high-dimensional model spaces. Through a linearization of the forward problem around the prior mean (i.e., the “best” initial model), this efficiency can be further improved. We show, however, that this linearization leads to a performance in which the output of an HMC chain strongly depends on the quality of the prior, in particular because not all (induced) earthquake model parameters have a linear relationship with the recordings observed at the surface. To mitigate the importance of an accurate prior, we integrate the linearized HMC scheme into a workflow that (i) allows for a weak prior through linearization around various (initial) centroid locations, (ii) is able to converge to the mode containing the model with the (global) minimum misfit by means of an iterative HMC approach, and (iii) uses variance reduction as a criterion to include the output of individual Markov chains in the estimation of the posterior probability. Using a three-dimensional heterogeneous subsurface model of the Groningen gas field, we simulate an induced earthquake to test our workflow. We then demonstrate the virtue of our workflow by estimating the event's centroid (three parameters), moment tensor (six parameters), and the earthquake's origin time. Using the synthetic case, we find that our proposed workflow is able to recover the posterior probability of these source parameters rather well, even when the prior model information is inaccurate, imprecise, or both inaccurate and imprecise.

Fully scalable solver for frequency-domain visco-elastic wave equations in 3D heterogeneous media: A controllability approach

We develop a controllability strategy for the computation of frequency-domain solutions of the 3D visco-elastic wave equation, in the perspective of seismic imaging applications. We generalize the controllability results for such equations beyond the sound-soft scattering (obstacle) problem. We detail the conjugate gradient implementation and show how an inner elliptic problem needs to be solved to compute the Riesz representative of the gradient at each iteration. We select a spectral-element spatial discretization and a fourth-order Runge-Kutta time discretization. We implement the controllability method in the framework of the SEM46 full waveform modeling and inversion software, to inherit for its excellent scalability which relies on an efficient domain decomposition algorithm. We perform a series of numerical experiments to validate the approach and illustrate its scalability up to more than fifteen hundred cores. In this case, with an elapsed time of less than 50 minutes, we solve a problem on a cubic domain containing up to 160 wavelengths in each direction, involving more than 1.7 billion unknowns.

Efficient numerical solution of the volume integro-differential equation for time-varying temperature fields in heterogeneous media

An efficient numerical method for calculation of time-varying temperature fields in an infinite homogeneous medium containing heterogeneous inclusions is presented. The medium is subjected to distributed or localized heat sources. The problem is reduced to the integro-differential equation for the temperature fields in the regions occupied by the inclusions only. For numerical solution, this equation is discretized using Gaussian radial functions for approximation of spatial-dependence and piece-wise constant functions for approximation of time-dependence of the temperature fields. For Gaussian approximating functions, the elements of the matrix of the discretized problem are calculated in explicit analytical forms. The iterative minimal residual method and fast Fourier transform technique are used for solution of the discretized problems at a sequence of time moments. The numerical algorithms for solution of the discretized equations are presented for 1D, 2D, and 3D heterogeneous media. Numerical solutions are compared with exact and approximate solutions in the cases of spherical layered inclusions subjected to distributed and localized heat sources.

Three-dimensional flow of a power-law nanofluid within a cubic domain filled with a heat-generating and 3D-heterogeneous porous medium

Three-dimensional flow of a power-law nanofluid within a cubic domain filled with a heat-generating and 3D-heterogeneous porous medium

Radiation dose calculation in 3D heterogeneous media using artificial neural networks.

External beam radiotherapy (EBRT) treatment planning requires a fast and accurate method of calculating the dose delivered by a clinical treatment plan. However, existing methods of calculating dose distributions have limitations. Monte Carlo (MC) methods are accurate but can take too long to be clinically viable. Deterministic approaches are quicker but can be inaccurate under certain conditions, particularly near heterogeneities and air interfaces. Neural networks trained on MC-derived data have the potential to reproduce dose distributions that agree closely with the MC method while being significantly quicker to deploy. In this work we present a framework for training machine learning models capable of directly calculating the dose delivered to a point in three-dimensional (3D) heterogeneous media given only spatially local information. The framework consists of three parts. First, we describe a novel method of randomly generating 3D heterogeneous geometries using simplex noise. Dose distributions for training were obtained by importing these geometries into a MC simulation. The second and third parts of the framework are precalculated data channels, aligned with the patient computed tomography (CT) image, to be used as input to the model. These data channels are a computationally efficient way of encoding the parameters of an incident radiation beam while also allowing the model to learn from data that would otherwise be outside of its receptive field. We demonstrate the viability of the framework by a training small, fully connected neural network model to reproduce dose distributions from megavoltage photon beams. The trained network displayed excellent agreement with MC dose distributions in randomly generated geometries with an average gamma index (3%/3mm) pass rate of 94.7% and an average error of 1.45% of peak dose. Finally, the network was used to calculate the dose in a patient CT image, on which the network was not trained, producing similarly impressive results. A novel method of generating training data for learned radiation dosimetry models has been introduced, along with preprocessing steps that allow even simple models to reproduce accurate dose distributions for EBRT. More importantly, we have demonstrated that a model trained using the proposed framework can generalize from the training data to predicting the therapeutic dose in realistic media.

Fully Scalable Solver for Frequency-Domain Visco-Elastic Wave Equations in 3D Heterogeneous Media: A Controllability Approach

Fully Scalable Solver for Frequency-Domain Visco-Elastic Wave Equations in 3D Heterogeneous Media: A Controllability Approach

Exploring Physical Links between Fluid Injection and Nearby Earthquakes: The 2012 Mw 4.8 Timpson, Texas, Case Study

ABSTRACTIn this work, we integrate a fluid-flow model of 3D deformable porous media with a dynamic rupture model of earthquakes in 3D heterogeneous geologic medium. The method allows us to go beyond fault failure potential analyses and to examine how big an earthquake can be if part of a fault reaches failure due to fluid injection. We apply the method to the 17 May 2012 Mw 4.8 Timpson, Texas, earthquake as a case study. The simulated perturbations of pore pressure and stress from wastewater injection at the time of the mainshock are high enough (several MPa) to trigger an earthquake. Dynamic rupture modeling could reproduce the major observations from the Mw 4.8 event, including its size, focal mechanism, and aftershock sequence, and thus building a more convincing physical link between fluid injection and the Mw 4.8 earthquake. Furthermore, parameter space studies of dynamic rupture modeling allow us to place some constraints on fault frictional properties and background stresses. For the Timpson case, we find that a dynamic friction coefficient of ∼0.3, a value of ∼0.1 m for the critical slip distance in the slip-weakening friction law, and uniform effective normal stress are associated with the Timpson earthquake fault. By reproducing main features of the aftershock sequence of the mainshock, we also demonstrate that the method has potential to become a predictive tool for fluid injection design in the future.

PyKonal: A Python Package for Solving the Eikonal Equation in Spherical and Cartesian Coordinates Using the Fast Marching Method

AbstractThis article introduces PyKonal: a new open-source Python package for computing travel times and tracing ray paths in 2D or 3D heterogeneous media using the fast marching method for solving the eikonal equation in spherical and Cartesian coordinates. Compiled with the Cython compiler framework, PyKonal offers a Python application program interface (API) with execution speeds comparable to C or Fortran codes. Designed to be accurate, stable, fast, general, extensible, and easy to use, PyKonal offers low- and high-level API functions for full control and convenience, respectively. A scale-independent implementation allows problems to be solved at micro, local, regional, and global scales, and precision can be improved over existing open-source codes by combining different coordinate systems. The resulting code makes state-of-the-art computational capabilities accessible to novice programmers and is efficient enough for modern research problems in seismology.

Fully implicit hybrid two-level domain decomposition algorithms for two-phase flows in porous media on 3D unstructured grids

Simulation of subsurface flows in porous media is difficult due to the nonlinearity of the operators and the high heterogeneity of material coefficients. In this paper, we present a scalable fully implicit solver for incompressible two-phase flows based on overlapping domain decomposition methods. Specifically, an inexact Newton-Krylov algorithm with analytic Jacobian is used to solve the nonlinear systems arising from the discontinuous Galerkin discretization of the governing equations on 3D unstructured grids. The linear Jacobian system is preconditioned by additive Schwarz algorithms, which are naturally suitable for parallel computing. We propose a hybrid two-level version of the additive Schwarz preconditioner consisting of a nested coarse space to improve the robustness and scalability of the classical one-level version. On the coarse level, a smaller linear system arising from the same discretization of the problem on a coarse grid is solved by using GMRES with a one-level preconditioner until a relative tolerance is reached. Numerical experiments are presented to demonstrate the effectiveness and efficiency of the proposed solver for 3D heterogeneous medium problems. We also report the parallel scalability of the proposed algorithms on a supercomputer with up to 8,192 processor cores.

An efficient and high accuracy finite-difference scheme for the acoustic wave equation in 3D heterogeneous media

In this paper we developed a new explicit compact high-order finite difference scheme to solve the 3D acoustic wave equations with spatially variable acoustic velocity. The boundary conditions for the second derivatives of spatial variables have been derived by using the wave equation and the boundary conditions themselves. Theoretical analysis shows that the new scheme has an accuracy order of O(τ2) + O(h4), where τ is the time step and h is the grid size. Combined with Richardson extrapolation or Runge–Kutta method, the new method can be improved to 4th-order accuracy in time, which has been implemented and verified in this paper. Four numerical experiments are conducted to validate the efficiency and accuracy of the new scheme. The stability of the new scheme has been proved by an energy method, which shows that the new scheme is conditionally stable with a Currant–Friedrichs–Lewy (CFL) number which is slightly lower than that of the Padé approximation based method. However, the new scheme is much simpler to implement.

Estimation of seismic source parameters in 3D elastic media using the reciprocity theorem

We have developed a novel method based upon reciprocity principles to simultaneously estimate the location of a seismic event and its source mechanism in 3D heterogeneous media. The method finds double-couple (DC) and non-DC mechanisms of microearthquakes arising from localized induced and natural seismicity. Because the method uses an exhaustive search of the 3D elastic media, it is globally convergent. It does not suffer from local minima realization observed with local optimization methods, including Newton, Gauss-Newton, or gradient-descent algorithms. The computational efficiency of our scheme is derived from the reciprocity principle, in which the number of 3D model realizations corresponds to the number of measurement receivers. The 3D forward modeling is carried out in the damped Fourier domain with a 3D finite-difference frequency-domain fourth- and second-order code developed to simulate elastic waves generated by seismic sources defined by forces and second-order moment density tensors. We evaluate the results of testing this new methodology on synthetic data for the Raft River geothermal field, Idaho, as well as determine its applicability in designing optimal borehole monitoring arrays in a fracking experiment at the Homestake Mine, South Dakota. We also find that the method proposed here can retrieve the moment tensors of the space distributed source with data arising from spatially restricted arrays with limited aperture. The effects of uncertainties on the source parameter estimation are also examined with respect to data noise and model uncertainty.

Adjoint computational methods for 2D inverse design of linear transport equations on unstructured grids

We address the problem of inverse design of linear hyperbolic transport equations in 2D heterogeneous media. We develop numerical algorithms based on gradient-adjoint methodologies on unstructured grids. While the flow equation is compulsorily solved by means of a second order upwind scheme so to guarantee sufficient accuracy, the necessity of using the same order of approximation when solving the sensitivity or adjoint equation is examined. Two test cases, including Doswell frontogenesis, are analysed. We show the convenience of using a low order method for the adjoint resolution, both in terms of accuracy and efficiency. An analytical explanation for this fact is also provided in the sense that, when employing higher order schemes for the adjoint problem, spurious high frequency numerical components slow down the convergence process.

3DMRT: A Computer Package for 3D Model‐Based Seismic Wave Propagation

ABSTRACT The prediction of accurate source‐to‐receiver travel times and wave paths through heterogeneous media is of major interest in global seismology and microseismic communities. Many algorithms have been proposed to address this problem, among which eikonal solvers have the best accuracy but lack computational efficiency. To facilitate the use of eikonal solvers with a high performance and visualization ability, this article presents a free, open‐source, graphical package named 3DMRT. 3DMRT propagates wavefronts through a 3D heterogeneous medium. Starting with a geologic model, the package first constructs a gridded velocity model. The application implements nine eikonal solvers and provides one‐point and multipoint raytracing functionality. In addition, 3DMRT is complemented with an additional command line tool that allows integration into other programs for further applications.

Source characteristics of the March 16, 2014 Mw 6.7 earthquake and its implications for the Mw 8.2 Pisagua mainshock

The March 16, 2014 Mw 6.7 earthquake played an important role in the complex seismic sequence leading to the nucleation of the April 1st, 2014 Mw 8.2 Pisagua earthquake in northern Chile. The Mw 6.7, an upper plate reverse faulting event with nodal planes highly rotated counterclockwise with respect to the strike of the megathrust at its location, is analyzed. For each nodal plane of the regional W-phase centroid moment tensor solution, the spatiotemporal slip distribution is inverted using near-field records. The optimal space and time smoothing constraints are determined objectively based on the Akaike's Bayesian Information Criterion. We analyze seismicity on the region affected by the Mw 6.7 and Mw 8.2 Pisagua mainshock, with accurate hypocenter locations obtained using a 3D heterogeneous medium. We also performed regional moment tensor inversion of events (M≥4.0) occurred before the Mw 8.2. Our results suggest the sub-vertical fault plane dipping southward (dip∼70°) as being the causative fault of the March 16 Mw 6.7 earthquake. The rupture propagated to the northwest, lasted 15s, and yielded a total seismic moment of 1.36×1019Nm (Mw 6.7). Estimated slip is characterized by total rupture length of ∼22km, and suggest shallow slip distributed between ∼ 3km, down to 17km depth. A comparison between coastal tide gauge records and forward tsunami modeling show similarity of the time series, thus supporting the plausibility of our estimated slip models. Our interpretation suggests that the rupture of the Mw 6.7 intraplate could be explained by stress accommodation between two segments along-dip with different mechanical properties at the interplate boundary, and inside the upper plate.

Direct frequency-domain 3D acoustic solver with intermediate data compression benchmarked against time-domain modeling for full-waveform inversion applications

We have developed a fast direct solver for numerical simulation of acoustic waves in 3D heterogeneous media. The Helmholtz equation is approximated by a 27-point finite-difference stencil of second-order accuracy that is optimized to reduce the numerical dispersion. Due to the optimization, dispersion errors less than 1% are achieved with model discretization of only five points per wavelength. Wave-propagation problem requires solving a large system of linear equations with complex sparse symmetric coefficient matrix with seismic shots representing the right-hand side. The first step is the triangular factorization of the coefficient matrix followed by solving systems of linear equations with triangular coefficient matrices as a second step. Having defined the triangular factors, the second step is very cheap, and its linear scaling with respect to the number of shots is the main advantage of direct methods. To reduce memory consumption and computational time at the factorization step, the lower triangular factors are compressed using a low-rank approximation of their nonzero blocks. The compression enables rapid solving of systems of more than [Formula: see text] equations corresponding to realistic geophysical models. Accuracy and performance comparison of our solver with a highly optimized time-domain solver proves that these approaches complement each other — depending on the problem size and computing configuration, either solver may be preferable.

Meshless local B-spline collocation method for heterogeneous heat conduction problems

Several numerical issues still pertain in the modeling of heterogeneous heat conduction, particularly from the viewpoints of material discontinuity and its handling and the presence of heat source. In this paper, a meshless local B-spline collocation method is presented for unsteady heat conduction problems of heterogeneous media. Unknown field variables are approximated by using B-spline basis functions within overlapped compact domains covering the geometry of materials. The present method is a truly meshless approach. The proposed approach is mainly coming with the following advantage that it is straightforward in dealing with discontinuity across the interface of heterogeneous materials. Treatment of discontinuity by using non-crossing interface compact domains allows material discontinuity to be handled geometrically without enforcing additional term/function at the interface. Several heat conduction problems in 2D and 3D heterogeneous media with arbitrary discontinuity shapes are considered. Attention is given for heterogeneous heat conduction problem accompanied by the presence of crack as well. The analysis is then completed by simulating effect of heat generation, in particular which produces high temperature rise inside a heterogeneous structure/component. Simulation results show that the proposed method is a simple and accurate numerical technique for solving unsteady heat conduction problems of heterogeneous media.

Tracking a Foam Front in a 3D, Heterogeneous Porous Medium

Foam is to be used as a blocking agent for confining a pollutant source zone and avoid spreading in an aquifer. To this end, it is necessary to determine where injected foam flows and stays inside a porous medium. This study examines the use of electrical resistivity tomography for this purpose. Foam is injected in a large-scale 3D heterogeneous porous medium (0.84 × 0.84 × 0.84 m). During the injection, electrical resistivity tomography measurements are performed. We show that combining a large number of measurements with inversion techniques allows for the monitoring of a foam front in 3D during the injection process.

A hierarchy of models for simulating experimental results from a 3D heterogeneous porous medium

In this work we examine the dispersion of conservative tracers (bromide and fluorescein) in an experimentally-constructed three-dimensional dual-porosity porous medium. The medium is highly heterogeneous ($\sigma_Y^2=5.7$), and consists of spherical, low-hydraulic-conductivity inclusions embedded in a high-hydraulic-conductivity matrix. The bi-modal medium was saturated with tracers, and then flushed with tracer-free fluid while the effluent breakthrough curves were measured. The focus for this work is to examine a hierarchy of four models (in the absence of adjustable parameters) with decreasing complexity to assess their ability to accurately represent the measured breakthrough curves. The most information-rich model was (1) a direct numerical simulation of the system in which the geometry, boundary and initial conditions, and medium properties were fully independently characterized experimentally with high fidelity. The reduced models included; (2) a simplified numerical model identical to the fully-resolved direct numerical simulation (DNS) model, but using a domain that was one-tenth the size; (3) an upscaled mobile-immobile model that allowed for a time-dependent mass-transfer coefficient; and, (4) an upscaled mobile-immobile model that assumed a space-time constant mass-transfer coefficient. The results illustrated that all four models provided accurate representations of the experimental breakthrough curves as measured by global RMS error. The primary component of error induced in the upscaled models appeared to arise from the neglect of convection within the inclusions. Interestingly, these results suggested that the conventional convection-dispersion equation, when applied in a way that resolves the heterogeneities, yields models with high fidelity without requiring the imposition of a more complex non-Fickian model.

Comparison of enriched meshless finite volume and element free Galerkin methods for the analysis of heterogeneous media

Meshless methods are considered as powerful numerical methods for computational mechanics, but due to high-order continuity of solution space in the standard meshless methods, they may lead to undesired oscillations in derivative fields in heterogeneous media. Towards an efficient meshless computation in heterogeneous media, the enrichment technique is already used in combination with different methods, such as element free Galerkin (EFG) method. In this work, it is demonstrated that the enriched formulation of the meshless finite volume method (FVM) can be considered as a potential alternative. In contrast to the EFG method, the meshless FVM is based on local weak form and utilizes a Petrov–Galerkin procedure. In this paper, a detailed comparison is made between the enriched formulation of these two methods, and then, their capability to the analysis of 1D and 2D heterogeneous media is investigated. It is demonstrated that enriched EFG method is more accurate for the analysis of 1D heterogeneous problems. However, for the analysis of 2D heterogeneous problems, enriched meshless FVM reveals more accuracy in both displacement and stress fields predictions.

Efficient and accurate numerical simulation of acoustic wave propagation in a 2D heterogeneous media

In this paper, a compact fourth-order finite difference scheme is derived to solve the 2D acoustic wave equation in heterogenous media. The Padé approximation is used to obtain fourth-order accuracy in both temporal and spatial dimensions, and the alternating direction implicit (ADI) technique is used to reduce the computational cost. Due to the non-constant wave velocity, the conventional ADI method is hard to implement as the algebraic manipulation cannot be used here. A novel numerical strategy is proposed in this work so that the compact scheme still maintains fourth-order accuracy in time and space. The fourth-order convergence order was firstly proved by theoretical error analysis, then was confirmed by numerical examples. It was shown that the proposed method is conditionally stable with a Courant–Friedrichs–Lewy (CFL) condition that is comparable to other existing finite difference schemes. Several numerical examples were solved to demonstrate the efficiency and accuracy of the new algorithm.