Advanced BEM-based methodologies to identify and simulate wave fields in complex geostructures

To enhance the applicability of BEM for geomechanical modeling numerically optimized BEM models, hybrid FEM-BEM models, and parallel computation of seismic Full Waveform Inversion (FWI) in GPU are implemented. Inverse modeling of seismic wave propagation in inhomogeneous and heterogeneous half-plane is implemented in Boundary Element Method (BEM) using Particle Swarm Optimization (PSO). The Boundary Integral Equations (BIE) based on the fundamental solutions for homogeneous elastic isotropic continuum are modified by introducing mesh-dependent variables. The variables are optimized to obtain the site-specific impedance functions. The PSO-optimized BEM models have significantly improved the efficiency of BEM for seismic wave propagation in arbitrarily inhomogeneous and heterogeneous media. Similarly, a hybrid BEM-FEM approach is developed to evaluate the seismic response of a complex poroelastic soil region containing underground structures. The far-field semi-infinite geological region is modeled via BEM, while the near-field finite geological region is modeled via FEM. The BEM region is integrated into the global FEM system using an equivalent macro-finite-element. The model describes the entire wave path from the seismic source to the local site in a single hybrid model. Additionally, the computational efficiency of time domain FWI algorithm is enhanced by parallel computation in CPU and GPU

Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion

We introduce a technique to compute exact anelastic sensitivity kernels in the time domain using parsimonious disk storage. The method is based on a reordering of the time loop of time-domain forward/adjoint wave propagation solvers combined with the use of a memory buffer. It avoids instabilities that occur when time-reversing dissipative wave propagation simulations. The total number of required time steps is unchanged compared to usual acoustic or elastic approaches. The cost is reduced by a factor of 4/3 compared to the case in which anelasticity is partially accounted for by accommodating the effects of physical dispersion. We validate our technique by performing a test in which we compare the $K_\alpha$ sensitivity kernel to the exact kernel obtained by saving the entire forward calculation. This benchmark confirms that our approach is also exact. We illustrate the importance of including full attenuation in the calculation of sensitivity kernels by showing significant differences with physical-dispersion-only kernels

Rapid Seismic Waveform Modeling and Inversion with Universal Neural Operators

Seismic waveform modeling is a powerful tool for determining earth structure models and unraveling earthquake rupture processes, but it is computationally expensive. We introduce a scheme to vastly accelerate these calculations with a recently developed machine learning paradigm called the neural operator. Once trained, these models can simulate a full wavefield for arbitrary velocity models at negligible cost. We use a U-shaped neural operator to learn a general solution operator to the 2D elastic wave equation from an ensemble of numerical simulations performed with random velocity models and source locations. We show that full waveform modeling with neural operators is nearly two orders of magnitude faster than conventional numerical methods, and more importantly, the trained model enables accurate simulation for arbitrary velocity models, source locations, and mesh discretization, even when distinctly different from the training dataset. The method also enables efficient full-waveform inversion with automatic differentiation

GPU-accelerated discontinuous Galerkin methods on hybrid meshes: applications in seismic imaging

Seismic imaging is a geophysical technique assisting in the understanding of subsurface structure on a regional and global scale. With the development of computer technology, computationally intensive seismic algorithms have begun to gain attention in both academia and industry. These algorithms typically produce high-quality subsurface images or models, but require intensive computations for solving wave equations. Achieving high-fidelity wave simulations is challenging: first, numerical wave solutions may suffer from dispersion and dissipation errors in long-distance propagations; second, the efficiency of wave simulators is crucial for many seismic applications. High-order methods have advantages of decreasing numerical errors efficiently and hence are ideal for wave modelings in seismic problems. Various high order wave solvers have been studied for seismic imaging. One of the most popular solvers is the finite difference time domain (FDTD) methods. The strengths of finite difference methods are the computational efficiency and ease of implementation, but the drawback of FDTD is the lack of geometric flexibility. It has been shown that standard finite difference methods suffer from first order numerical errors at sharp media interfaces. In contrast to finite difference methods, discontinuous Galerkin (DG) methods, a class of high-order numerical methods built on unstructured meshes, enjoy geometric flexibility and smaller interface errors. Additionally, DG methods are highly parallelizable and have explicit semi-discrete form, which makes DG suitable for large-scale wave simulations. In this dissertation, the discontinuous Galerkin methods on hybrid meshes are developed and applied to two seismic algorithms---reverse time migration (RTM) and full waveform inversion (FWI). This thesis describes in depth the steps taken to develop a forward DG solver for the framework that efficiently exploits the element specific structure of hexahedral, tetrahedral, prismatic and pyramidal elements. In particular, we describe how to exploit the tensor-product property of hexahedral elements, and propose the use of hex-dominant meshes to speed up the computation. The computational efficiency is further realized through a combination of graphics processing unit (GPU) acceleration and multi-rate time stepping. As DG methods are highly parallelizable, we build the DG solver on multiple GPUs with element-specific kernels. Implementation details of memory loading, workload assignment and latency hiding are discussed in the thesis. In addition, we employ a multi-rate time stepping scheme which allows different elements to take different time steps. This thesis applies DG schemes to RTM and FWI to highlight the strengths of the DG methods. For DG-RTM, we adopt the boundary value saving strategy to avoid data movement on GPUs and utilize the memory load in the temporal updating procedure to produce images of higher qualities without a significant extra cost. For DG-FWI, a derivation of the DG-specific adjoint-state method is presented for the fully discretized DG system. Finally, sharp media interfaces are inverted by specifying perturbations of element faces, edges and vertices

High-performance tsunami modelling with modern GPU technology

PhD ThesisEarthquake-induced tsunamis commonly propagate in the deep ocean as long waves and develop into sharp-fronted surges moving rapidly coastward, which may be effectively simulated by hydrodynamic models solving the nonlinear shallow water equations (SWEs). Tsunamis can cause substantial economic and human losses, which could be mitigated through early warning systems given efficient and accurate modelling. Most existing tsunami models require long simulation times for real-world applications. This thesis presents a graphics processing unit (GPU) accelerated finite volume hydrodynamic model using the compute unified device architecture (CUDA) for computationally efficient tsunami simulations. Compared with a standard PC, the model is able to reduce run-time by a factor of > 40. The validated model is used to reproduce the 2011 Japan tsunami. Two source models were tested, one based on tsunami waveform inversion and another using deep-ocean tsunameters. Vertical sea surface displacement is computed by the Okada model, assuming instantaneous sea-floor deformation. Both source models can reproduce the wave propagation at offshore and nearshore gauges, but the tsunameter-based model better simulates the first wave amplitude. Effects of grid resolutions between 450-3600 m, slope limiters, and numerical accuracy are also investigated for the simulation of the 2011 Japan tsunami. Grid resolutions of 1-2 km perform well with a proper limiter; the Sweby limiter is optimal for coarser resolutions, recovers wave peaks better than minmod, and is more numerically stable than Superbee. One hour of tsunami propagation can be predicted in 50 times on a regular low-cost PC-hosted GPU, compared to a single CPU. For 450 m resolution on a larger-memory server-hosted GPU, performance increased by ~70 times. Finally, two adaptive mesh refinement (AMR) techniques including simplified dynamic adaptive grids on CPU and a static adaptive grid on GPU are introduced to provide multi-scale simulations. Both can reduce run-time by ~3 times while maintaining acceptable accuracy. The proposed computationally-efficient tsunami model is expected to provide a new practical tool for tsunami modelling for different purposes, including real-time warning, evacuation planning, risk management and city planning

Rapid three-dimensional multiparametric MRI with quantitative transient-state imaging

Novel methods for quantitative, transient-state multiparametric imaging are increasingly being demonstrated for assessment of disease and treatment efficacy. Here, we build on these by assessing the most common Non-Cartesian readout trajectories (2D/3D radials and spirals), demonstrating efficient anti-aliasing with a k-space view-sharing technique, and proposing novel methods for parameter inference with neural networks that incorporate the estimation of proton density. Our results show good agreement with gold standard and phantom references for all readout trajectories at 1.5T and 3T. Parameters inferred with the neural network were within 6.58% difference from the parameters inferred with a high-resolution dictionary. Concordance correlation coefficients were above 0.92 and the normalized root mean squared error ranged between 4.2% - 12.7% with respect to gold-standard phantom references for T1 and T2. In vivo acquisitions demonstrate sub-millimetric isotropic resolution in under five minutes with reconstruction and inference times < 7 minutes. Our 3D quantitative transient-state imaging approach could enable high-resolution multiparametric tissue quantification within clinically acceptable acquisition and reconstruction times.Comment: 43 pages, 12 Figures, 5 Table

A SEISMOLOGIC STUDY OF THE NORTHERN MISSISSIPPI EMBAYMENT

Part 1: Crustal structure in the New Madrid Seismic Zone (NMSZ) is investigated through a detailed study of explosion data obtained from the Embayment Seismic Excitation Experiment. The data show a distinct anisotropy in distance attenuation for both P and S waves in the range from 0 to 200km distance. Waves that propagate northward from the 1,134kg Marked Tree, Arkansas, explosion attenuate quickly with distance until a range of about 100km from the source where high-amplitude, high-phase velocity critical reflections from the boundary between the middle crust and rift pillow structure produce high amplitude waves. Propagation southward from the 2,268kg Mooring, Tennessee blast shows less distance attenuation compared to northward propagation. Reflections from the middle crust-lower crust boundary occur but do not significantly increase in amplitude with distance and travel with slower apparent phase velocity than observed for the northward propagation data set. A smooth velocity model is developed using a stabilized Weichert-Herglotz travel time inversion using first arrival travel times. Then an inversion using the travel time of both direct and middle crustal reflected waves is developed to obtain a 2D inhomogeneous-layered isotropic crustal model. The result reveals that there is a significant southwest dip to the top of the middle crust interface in the vicinity of the NMSZ, consistent with previously inferred changes in the thickness of the rift pillow model. This 2D feature characterizes the local wave propagation along the Reelfoot Rift and demonstrates the need for an improvement of the current Central United States velocity model.Part 2: Obtaining reliable empirical Greens functions (EGFs) from ambient noise by seismic interferometry requires homogenously distributed noise sources. However, it is difficult to attain this condition since ambient noise data usually contains highly correlated signals from earthquakes or other transient sources from human activities. Removing these transient signals is one of the most essential steps in the whole data processing flow to obtain EGFs. We propose to use a denoising method based on the continuous wavelet transform to achieve this goal. The noise level is estimated in the wavelet domain for each scale by determing the 99% confidence level of the empirical probability density function of the noise wavelet coefficients. The correlated signals are then removed by an efficient soft thresholding method. The same denoising algorithm is also applied to remove the noise in the final stacked cross-correlogram. A complete data processing workflow is provided with the overall data processing procedure divided into four stages: (1) single station data preparation, (2) removal of earthquakes and other transient signals in the seismic record, (3) spectrum whitening, cross-correlation and temporal stacking, and (4) remove the noise in the stacked cross-correlogram to deliver the final EGF. The whole process is automated to make it accessible for large datasets. Synthetic data constructed with a recorded earthquake and recorded ambient noise is used to test the denoising method. We then apply the new processing workflow to data recorded by the USArray Transportable Array stations near the New Madrid Seismic Zone where many seismic events and transient signals are observed. We compare the EGFs calculated from our workflow with commonly used time domain normalization method and our results show improved signal-to-noise ratios. The new workflow can deliever reliable EGFs for further studies.Part 3: We incorporate seismic ambient noise data recorded by different temporary and permanent broadband stations around the northern Mississippi Embayment from 1990 to 2018 to develop a crustal shear wave velocity (Vs) model for this area with full waveform ambient noise tomography. Empirical Greens functions at periods between 8 and 40s for all the possible pairs of stations are extracted by using a new seismic ambient noise data processing flow based on the continuous wavelet transform. Synthetic waveforms are then calculated through a heterogeneous Earth model using a GPU-enabled collocated finite-difference code. The cross-correlation time shifts between the synthetic waveforms and the extracted empirical Greens functions are used to construct the velocity updated kernel by using the adjoint method. Starting from the Central United States Velocity Model, the shear wave velocity model is then iteratively updated with the Vs kernel calculated in each iteration. Checkerboard tests show that perturbations in the top 30km of the crust are well recovered but amplitude recovery ability gradually decreases for deeper structure. We find that velocity lows characterize the Reelfoot Rift Graben and Rough Creek Graben separated by a high velocity crust. High velocity anomalies are observed under the Ozark Uplift and Paducah Gravity Lineament. A low velocity area previously interpreted as the Missouri Batholith is observed between them. A massive high velocity body in the southeast Mississippi Embayment is observed and is explained by the faulting as well as partly mafic intrusion. The Ouachita-Appalachian Thrust Front is clearly observed with a thinner crustal layer underneath. The rift pillow is well observed in the final tomography model along the Reelfoot Rift in the lower crust. The final inverted velocity model is consistent with local geological features and can be used for other seismological studies such as earthquake source determination and earthquake hazard assessment