Adaptive simulation of wave propagation problems including dislocation sources and random media

Abstract

An adaptive Tikhonov regularization is integrated with an h-adaptive grid-based scheme for simulation of elastodynamic problems, involving seismic sources with discontinuous solutions and random media. The Tikhonov method is adapted by a newly-proposed detector based on the MINMOD limiters and the grids are adapted by the multiresolution analysis (MRA) via interpolation wavelets. Hence, both small and large magnitude physical waves are preserved by the adaptive estimations on non-uniform grids. Due to developing of non-dissipative spurious oscillations, numerical stability is guaranteed by the Tikhonov regularization acting as a post-processor on irregular grids. To preserve waves of small magnitudes, an adaptive regularization is utilized: using of smaller amount of smoothing for small magnitude waves. This adaptive smoothing guarantees also solution stability without over smoothing phenomenon in stochastic media. Proper distinguishing between noise and small physical waves are challenging due to existence of spurious oscillations in numerical simulations. This identification is performed in this study by the MINMOD limiter based algorithm. Finally, efficiency of the proposed concept is verified by: 1) three benchmarks of one-dimensional (1-D) wave propagation problems; 2) P-SV point sources and rupturing line-source including a bounded fault zone with stochastic material properties.

Keywords

adaptive wavelet / adaptive smoothing / discontinuous solutions / stochastic media / spurious oscillations / Tikhonov regularization / minmod limiter

Introduction

The accurate forward numerical simulations of earthquake sources are important for studying the physics of earthquakes, strong ground motions and near fault/field effects. Forward simulations can also provide a constraint for the solutions of inverse problems for the identification of earthquake sources. This is important, since, the inverse problems are naturally ill-posed with several admissible solutions. The direct simulations of the sources, however, are challenging due to developing of singularities and discontinuities around sources. Day et al. [¹] and Dalguer and Day [²] showed that the discontinuities can affect considerably the accuracy of numerical simulations, regardless of the type or the order of numerical discretization. To handle the discontinuous solutions in wave propagation problems, several numerical techniques have been developed, such as: Finite difference (FD) [¹,³–⁶], Finite element (FE) [⁷–¹¹], Spectral element (SE) [¹²–¹⁴], Boundary integral (BI) [¹⁵–¹⁸], Finite volume (FV), and Discontinuous Galerkin (DG) methods [¹⁹–²⁴]. Some of these schemes can inherently handle discontinuities (like the DG and modified FV methods), while others should be modified/improved for proper modeling. In this regard, the aims of this study are:

1) Integration of the Tikhonov-based adaptive-regularization scheme with h-adaptive grid-based numerical methods to simulate: discontinuities, multi-scale responses, random waves in stochastic media, radiated waves from earthquake sources.

2) Developing a new MIMMOD-based solution-indicator for the adaptation of the regularization method; this indicator can work satisfactory on both non-uniform grids and noisy solutions (e.g., those from stochastic-like domains).

The regularization step can be performed by a post-processing stage to filter out and stabilize numerical solutions. The main drawback of the regularization, however, is that fine-scale waves (emitted from earthquake sources or developed in stochastic media) can also be filtered out. To cure this, one effective approach is the adaptive filtering in the spatial domain: adaptive estimations. Generally, artificial (spurious) oscillations develop in numerical simulations of wave equations; distinguishing between these oscillations and real ones is necessary for proper simulations. For this differentiating, two approaches are followed in this study:

1) The wavelet-based graded tree, proposed by Harten [²⁵].

2) The algorithm based on the MINMOD slope limiters, suggested in this study.

The graded tree can detect and decompose effectively a computational domain into different spatial sub-domains in a spatio-resolution representation. This is promising for adaptive simulation and filtering. This approach has successfully been used for adaptive solutions of numerical methods with the Total Variation Diminishing (TVD) feature [²⁵–²⁷]. Most of commonly used numerical methods, however, are not TVD and this leads to developing of spurious oscillations in numerical solutions. For this reason, the graded trees may not be suitable for non-TVD schemes; as the tress is sensitive for spurious oscillations (it is shown in this study). This sensitivity is nearly negligible for the MINMOD limiters; as, they can effectively ignore noise-like fluctuations. There are several definitions of MINMOD limiters with different amount of dissipations. MINMOD limiters with more surrounding points lead to more dissipative results. This dissipation difference of different definitions can be used as a criterion to detect both fine-scale and coarse-scale physical variations from noise-like fluctuations. This is because noises do not have coherency with surrounding data. The limiters can also work directly on non-uniform data, and this makes them effective for detections on grid-based adaptive solutions.

Smoothing on non-uniform grids, however, is an ill-posed problem. For this reason, the famous Tikhonov regularization concept is used to guarantee reliability and stability of estimations directly on the non-uniform grids. For the regularization, a specific case of the Tikhonov regularization, the cubic smoothing spline with variable smoothing parameter is used (to perform an adaptive filtering) [²⁸–³⁰]. For smoothing splines, stability, convergence rate of function estimations and corresponding derivatives [³⁰,³¹], boundary effects [³²] and fast algorithms [³³,³⁴] had been studied. Due to under [³²] or over [³³] smoothing effects, proper ranges of the regularization can be obtained by the trial-error method or the L-curve. This curve exhibits trade-off between estimation errors and smoothness [³⁵].

The Tikhonov-based smoothing methods yield the smoothest result for a distinct error in an estimation [³⁶]. This feature is illustrated in this study by a benchmark problem containing discontinuities. This problem is also simulated by some other commonly used approaches developed for discontinuity handling in stress wave propagation problems. Some of these methods are: the discontinuous time Galerkin [³⁷], Taylor-Galerkin [³⁸], and the-generalized dissipative time integration schemes [³⁹]. The results reveal importance of the regularization and the optimum feature of the Tikhonov method. To study the effects of the adaptive filtering, another 1-D wave propagation problem is presented involving both large and small waves. Performances of the slope limiters and the wavelet-based graded tree are then examined; investigations are performed for data with and without noise effects.

To improve accuracy of numerical simulations, the adaptive estimation can be integrated by multiresolution-based grid adaptations; this extra h-adaptation leads to both method and grid adaptations. Simulations of partial differential equations (PDEs) by the multiresolution analysis (MRA) motivate many researches. MRA-based grid adaptations for PDEs with sharp transition zones improve accuracy of fine responses [⁴⁰–⁵⁴]. In this approach, grid concentrations are in accordance with solution gradients. In most studies, the wavelet theory was only used as a tool to detect high-gradient zones. Thereafter, the spatio-temporal discretization is performed by other commonly used schemes (e.g., the collocation, FD or FV method). This leads to simple and straight forward numerical algorithms. This approach is mainly used for the simulation of parabolic or elliptic systems. This is because, these problems have inherent dissipations [⁵⁵] improving numerical stabilities.

In simulation of wave equations (hyperbolic PDEs), spurious oscillations develop and propagate in computational domains due to the numerical truncation error [⁵⁶] and the lack of the inherent dissipation [⁵⁵]. In wave propagation problems, both numerical dissipation and dispersion phenomena exist [⁵⁶]. Propagating spurious oscillations lead to improper adapted grids (even on smooth regions) and numerical instabilities, especially around discontinuities (this is shown in this study in simulation of earthquake sources).

In conclusion, two effective approaches are utilized to control the Gibb’s phenomenon around discontinuities and spurious oscillations while preserving fine-scale responses, as:

1) Adaptive regularization of noisy solutions [⁵⁵,⁵⁷] on non-uniform grids by the MIMMOD limiters as a post-processing stage. Adaptive smoothing is important to preserve small magnitude responses; such as: high frequency waves radiated from sources or scattered waves developed in stochastic media.

2) Concentrating extra grid points around high gradient solutions to capture fine responses by the MRA. This adaptation is important in seismic source modeling, due to: geometry details of faults, length scale effects: faults can extend to hundreds of kilometers while rupturing would be of the order of meter.

This so-called h-adaptive approach has widely been integrated with different numerical methods; for example, in: High resolution schemes [⁴⁰,⁴⁷,⁵⁸–⁶⁰], Finite element methods [⁶¹–⁶⁴], Mesh-free methods [⁶⁵–⁶⁹], Isogeometric analysis [⁷⁰–⁷³].

Due to source rupturing, there exists bounded zones around sources with crushed materials and soften properties. Variations of properties have stochastic nature leading to wave propagation problems in stochastic media. Such propagation, however, introduces another challenging due to developing of spurious oscillations and numerical instabilities [⁷⁴–⁷⁶]. Magnitude of these oscillations is in accordance with the perturbation size of physical properties. Most of the commonly used numerical methods become unstable on the stochastic media; the methods can handle effectively stochastic media for small perturbations (less than 10% variation in the standard deviation). Our investigations showed that small amount of regularization can stabilize numerical solvers on stochastic-like media [³⁴]. On the other hand, large amount of smoothing can yield over-smoothing and even artificial amplification of propagating waves. For stochastic media, performance of the regularization on numerical simulations can be checked by comparing the numerical scattering attenuation, {{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}} with predicted ones obtained by the single scattering theory [⁷⁷]. This is performed here by a 1-D benchmark.

Wave propagation in stochastic media can mobilize multiscale responses at different scales. In this study, the adaptive regularization approach is used for preserving such solutions. Other proposed approaches considering the scale effect are:

1) The multiresolution-based projection of wave equations in wavelet spaces [⁷⁸,⁷⁹].

2) Different multiscale methods simulating different scales simultaneously [⁸⁰–⁸⁵].

3) Hierarchical enrichment strategies for resolving the fine-scale features [⁸⁶–⁹¹].

Finally, it should be noted that the wavelet-based approach presented in this study is a non-projection method; since operators are not projected into wavelet spaces. The wavelet based projection approach has been used for stress wave propagation problems in several works [⁷⁴–⁷⁶,⁹²–⁹⁵]. However, it suffers from proper projection of nonlinear operators and handling of boundary conditions [⁴³,⁹⁶].

This paper is composed of eight parts. Section 2 devotes to 2-D P-SV wave propagations in media containing earthquake sources (defined by the dislocation theory) and stochastic material properties. Here, construction of stochastic media is reviewed as well as evaluating of the attenuation factor and the predicted one (from the single scattering theory). Section 3 explains the main concept of the multiresolution-based grid adaptation by the interpolating wavelets. Section 4 presents issues related to: ill-posed problems, the Tikhonov regularization, smoothing splines and a benchmark problem. Section 5 explains two different solution detectors: the wavelet-based graded tree and the one based on the MINMOD limiters. Section 6 reviews the spatio-temporal discretization and corresponding algorithm for stress wave problems. Section 7 includes some numerical examples containing three 1-D benchmarks and two 2-D earthquake sources. In the 1-D benchmarks, three examples are studied, as: 1) Performance comparison of the regularization method with other numerical schemes developed for simulation of mechanical waves; 2) Performance of the two detectors in problems with and without noise effects; 3) Some 1-D wave propagation problems in stochastic media to compare their numerical attenuations with predicted ones. The 2-D examples are: a point source and a rupturing earthquake source; the point source problem is simulated twice, with and without a fault zone described by a stochastic medium located around the source. The conclusion is presented at the end of the paper.

Wave equations

The P-SV wave equation and earthquake source modeling

For simulation of earthquake sources, two general approaches have been developed: the kinematic and the dynamic descriptions. In the kinematic description, a predefined spatio-temporal slip evolution is assumed for sources. In the dynamic modeling, rupture evolutions are determined by physical phenomena (e.g., based on applied stress fields and friction laws). In this study, the kinematic approach is used and the sources are simulated by the dislocation theory.

In the following, governing equations of P-SV waves are presented in infinite media with a point source. For the simulation of the infinite media, the concept of absorbing boundaries is used (by considering four surrounding absorbing boundaries). The boundaries are considered explicitly in the wave equations by the approach introduced in Ref. [⁹⁷]. The wave equation is modified by adding a damping term Q\left(x,z\right)\dot{u}\left(x,z,t\right) , where Q\left(x,z\right)\in C^{\mathrm{2}} is an attenuation function. This function is zero in the computational domain and increases gradually while approaching toward the artificial boundaries. This is for guarantee gradual diminishing of outward propagating waves. The modified P-SV wave equations with the absorbing terms are:

\left(\begin{array}{c}\hfill {{\displaystyle \partial }}_x\left(\lambda +\mathrm{2}\mu \right){{\displaystyle \partial }}_x\hfill \\ \hfill \mathrm{0}\hfill \end{array}\begin{array}{c}\hfill \mathrm{0}\hfill \\ \hfill {{\displaystyle \partial }}_z\left(\lambda +\mathrm{2}\mu \right){{\displaystyle \partial }}_z\hfill \end{array}\right)\mathit{u}

+\left(\begin{array}{c}\hfill \mathrm{0}\hfill \\ \hfill {{\displaystyle \partial }}_z\lambda {{\displaystyle \partial }}_x+{{\displaystyle \partial }}_x\mu {{\displaystyle \partial }}_z\hfill \end{array}\begin{array}{c}\hfill {{\displaystyle \partial }}_x\lambda {{\displaystyle \partial }}_z+{{\displaystyle \partial }}_z\mu {{\displaystyle \partial }}_x\hfill \\ \hfill \mathrm{0}\hfill \end{array}\right)\mathit{u}

(1)+\left(\begin{array}{c}\hfill {{\displaystyle \partial }}_z\mu {{\displaystyle \partial }}_z\hfill \\ \hfill \mathrm{0}\hfill \end{array}\begin{array}{c}\hfill \mathrm{0}\hfill \\ \hfill {{\displaystyle \partial }}_x\mu {{\displaystyle \partial }}_x\hfill \end{array}\right)\mathit{u}=\rho {{\displaystyle \partial }}_t^{\mathrm{2}}\mathit{u}+\rho Q\left(x,z\right){{\displaystyle \partial }}_t\mathit{u}+\rho \mathit{f},

where {{\displaystyle \partial }}_Y:={\displaystyle \frac{\partial }{\partial Y}} and {\partial }_Y^2:={\displaystyle \frac{{\partial }^2}{\partial Y^2}} for Y\in \{t,x,z\}; \mathit{u}=\left(\begin{array}{l}{{\displaystyle u}}_x\hfill \\ {{\displaystyle u}}_z\hfill \end{array}\right); \mathit{f}=\left(\begin{array}{l}{f}_x\hfill \\ f_z\hfill \end{array}\right); u_x and u_z denote displacements in the x and z directions, respectively; f_x and f_z denote the components of body forces in the x and z directions, respectively; and \lambda and \mu present the Lame coefficients. For a computational domain \Omega \in \left[\mathrm{0},L\right]\times \left[\mathrm{0},L\right] , the attenuation function Q can be defined as:

Q\left(x,z\right)={{\displaystyle a}}_x\left(e^{b_x{\left(x/L\right)}^{\mathrm{2}}}+e^{b_x\left(\mathrm{1}-{\left(x/L\right)}^{\mathrm{2}}\right)}\right)

(2)+{{\displaystyle a}}_z\left(e^{b_z{\left(z/L\right)}^{\mathrm{2}}}+e^{b_z\left(\mathrm{1}-{\left(z/L\right)}^{\mathrm{2}}\right)}\right),

where a_x and a_z denote attenuation coefficients in the x and z directions, respectively; coefficients b_x and b_z control effective width of absorbing boundaries in the vicinity of such boundaries.

Dislocations in sources (due to discontinuity of displacements on both sides of the sources) are implemented by an equivalent system of double-couple forces. In this context, the equivalent vector of the double-couple body force \mathit{f}\left(t\right):={\left[{{\displaystyle f}}_x,{{\displaystyle f}}_y,{{\displaystyle f}}_z\right]}^{\mathrm{T}} for 3-D sources could be expressed as Ref. [⁹⁸]:

(3)\mathit{f}\left(t\right)=-\mathit{M}\nabla \delta \left(\mathit{r}-{{\displaystyle \mathit{r}}}_{\mathrm{s}}\right)D\left(t\right),

where r is the location vector; r_s shows the source location; D(t) denotes displacement history of a particle on a fault. Function \dot{D}\left(t\right) corresponds to the far-field source time function in a way that {\displaystyle \int \dot{D}\left(t\right)}\mathrm{d}t=\mathrm{1} [⁹⁸]; the function \delta \left(x\right) is the Dirac delta; \nabla is the gradient operator; and M is the symmetric moment-tensor with elements M_ij where i=j=x,y,z. The moment-tensor is functions of the strike (\varphi), dip (\zeta) and rake (\eta ) angles defined by the fault geometry [⁹⁸]. Six independent components of the moment tensor M are:

(4)\begin{array}{l}{{\displaystyle M}}_{xx}=-{{\displaystyle M}}_{\mathrm{0}}\left(\sin \zeta \cos \eta \sin \mathrm{2}\varphi +\sin \mathrm{2}\zeta \sin \eta {\sin }^{\mathrm{2}}\varphi \right),\hfill \\ {{\displaystyle M}}_{yy}={{\displaystyle M}}_{\mathrm{0}}\left(\sin \zeta \cos \eta \sin \mathrm{2}\varphi -\sin \mathrm{2}\zeta \sin \eta {\cos }^{\mathrm{2}}\varphi \right),\hfill \\ {{\displaystyle M}}_{zz}={{\displaystyle M}}_{\mathrm{0}}\left(\sin \mathrm{2}\zeta \sin \eta \right),\hfill \\ {{\displaystyle M}}_{xy}={{\displaystyle M}}_{yx}={{\displaystyle M}}_{\mathrm{0}}\left(\sin \zeta \cos \eta \cos \mathrm{2}\varphi +\sin \mathrm{2}\zeta \sin \eta \sin \varphi \cos \varphi \right),\hfill \\ {{\displaystyle M}}_{xz}={{\displaystyle M}}_{zx}={{\displaystyle -M}}_{\mathrm{0}}\left(\cos \zeta \cos \eta \cos \varphi +\cos \mathrm{2}\zeta \sin \eta \sin \varphi \right),\hfill \\ {{\displaystyle M}}_{yz}={{\displaystyle M}}_{zy}={{\displaystyle -M}}_{\mathrm{0}}\left(\cos \zeta \cos \eta \sin \varphi -\cos \mathrm{2}\zeta \sin \eta \cos \varphi \right),\hfill \end{array}

where {{\displaystyle M}}_{\mathrm{0}}={{\displaystyle \mu }}_{\mathrm{s}}AD is the scalar seismic moment, in which {{\displaystyle \mu }}_{\mathrm{s}} , A, and D are the shear modulus in the source, the fault area, and an average displacement of the fault, respectively.

A rupturing could be simulated by superimposing of several pointwise sources. They are located along the rupturing direction in fault zones; each of which can have its own source time history. Hence, effects of dynamic source process could be captured in both radiation patterns and frequency contents of transient waves.

Due to rupturing of sources, materials around sources are crushed; they have softened mechanical properties with stochastic nature. This means, we face with wave propagation problems in stochastic media. For this reason, this phenomenon is explained briefly in the next sub-section.

Wave propagation in stochastic media

Construction of stochastic media

Let us assume a stochastic parameter \alpha \left(x,z\right) defined as \alpha \left(z,x\right)={{\displaystyle \alpha }}_{\mathrm{0}}\left[\mathrm{1}+\beta \left(x,z\right)\right] where {{\displaystyle \alpha }}_{\mathrm{0}} denotes the background value and \beta \left(x,z\right) presents the relative perturbation with zero mean value over the computational domain, \mathrm{\Omega }\in X\times Z . In the wavenumber domain, autocorrelation of the function \beta \left(x,z\right), denoted by \mathrm{\Phi }\left(x,z\right), has the following relationship with the perturbation function [⁷⁷]:

(5)\mathcal{F}\left[\mathrm{\Phi }\left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]={\displaystyle \frac{\mathrm{1}}{XZ}}{\left|\mathcal{F}\left[\beta \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]\right|}^{\mathrm{2}},

where \mathcal{F}[·] denotes the (2-D) Fourier transform.

For 2-D stochastic media, different autocorrelation functions have been suggested; one of them is the double-exponential function, defined as:

(6a)\mathrm{\Phi }\left(x,y\right)={\xi }^{\mathrm{2}}\exp \left(-\left|x\right|/b\right)\exp \left(-\left|y\right|/a\right),

with the Fourier transform:

(6b)\left|\mathcal{F}\left[\mathrm{\Phi }\left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]\right|={\xi }^{\mathrm{2}}{\displaystyle \frac{\mathrm{4}ab}{\left[\mathrm{1}+{\left(a{{\displaystyle k}}_x\right)}^{\mathrm{2}}\right]\left[\mathrm{1}+{\left(b{{\displaystyle k}}_z\right)}^{\mathrm{2}}\right]}},

where a and b are correlation lengths perpendicular and parallel to the propagation direction of the incident wave, respectively; and {\xi }^{\mathrm{2}} represents variance of fluctuations in the stochastic media. In brief, the following steps can be used to construct a stochastic medium:

1) Select a predefined autocorrelation function for the domain, \mathcal{F}\left[\mathrm{\Phi }\left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right] , and a variance of fluctuations in the stochastic medium (see Eq. (6)).

2) Compute \mathcal{F}\left[\beta \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right] from Eq. (5), as:

(7)\mathcal{F}\left[\beta \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]=\left|\mathcal{F}\left[\beta \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]\right|e^{i\times \phi \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)}=\sqrt{XZ}\sqrt{\mathcal{F}\left[\mathrm{\Phi }\left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)\right]}{e}^{i\times \phi \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right)},

where \varphi \left({{\displaystyle k}}_x,{{\displaystyle k}}_z\right) is random phases with uniform distributions over \left[-\mathrm{\pi },\mathrm{\pi }\right].

3) Compute the (2-D) inverse Fourier transform to obtain the perturbation function \beta \left(x,z\right) .

Generation of 1-D stochastic media is straightforward by using the 1-D Fourier transform and 1-D autocorrelation functions in the abovementioned algorithm.

For numerical simulation of the P-SV wave propagation problems in stochastic media, stochastic features are considered for \rho \left(x,z\right), \mu \left(x,y\right), and \lambda \left(x,y\right) as:

\rho \left(x,z\right)={{\displaystyle \rho }}_{\mathrm{0}}\left[\mathrm{1}+\beta \left(x,z\right)\right],

\mu \left(x,z\right)={{\displaystyle \mu }}_{\mathrm{0}}\left[\mathrm{1}+\beta \left(x,z\right)\right],

(8)\lambda \left(x,z\right)={{\displaystyle \lambda }}_{\mathrm{0}}\left[\mathrm{1}+\beta \left(x,z\right)\right].

Performance of numerical methods on stochastic media

In wave propagation problems on stochastic media, one common approach to investigate robustness of different numerical methods is to compare their attenuation functions {{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}} (due to scattering phenomenon in heterogeneous media) with predicted ones (obtained theoretically with some simplified assumptions).

In numerical simulations, the attenuation function {{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}} due to propagation through a heterogeneous layer with thickness z₀ can be estimated as Ref. [⁷⁷]:

(9){{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}}\left(\omega \right)=-{\displaystyle \frac{\mathrm{2}{{\displaystyle \alpha }}_{\mathrm{0}}}{\omega {{\displaystyle z}}_{\mathrm{0}}}}\ln \left[A\left(\omega \right)/{{\displaystyle A}}_{\mathrm{0}}\left(\omega \right)\right],

where \omega denotes the angular frequency; z₀ is thickness of the stochastic medium; {{\displaystyle A}}_{\mathrm{0}}\left(\omega \right) and A\left(\omega \right) denote amplitude spectrum of waves measured just before and after the heterogeneous layer, respectively (this means, {{\displaystyle A}}_{\mathrm{0}}\left(\omega \right) at z=0 is for incident waves and A\left(\omega \right) at z= z₀ belongs to transmitted waves through the heterogeneous layer); and {{\displaystyle \alpha }}_{\mathrm{0}} is the wave velocity ({{\displaystyle \alpha }}_{\mathrm{0}}=\sqrt{\lambda /\rho }).

Let us assume a 1-D stochastic medium produced by the exponential autocorrelation function, which is defined as: \mathrm{\Phi }\left(x\right)={\xi }^{\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\left|x\right|/b\right) (by assuming a\to \infty in Eq. (6a)). For this stochastic medium, for compression waves, the attenuation function, {{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}} resulted from the single scattering theory is [⁷⁷]:

(10){{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}}\left(k\right)={\xi }^{\mathrm{2}}{\left(\mathrm{1}+\overline{k}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}kb}{\mathrm{1}+\mathrm{4}{b}^{\mathrm{2}}{k}^{\mathrm{2}}}},

where k is the wave number; \overline{k} shows proportionality factor between the speed and density variations; i.e.: {\displaystyle \frac{\delta \alpha \left(x\right)}{\alpha \left(x\right)}}={\displaystyle \frac{\mathrm{1}}{\overline{k}}}{\displaystyle \frac{\delta \rho \left(x\right)}{\rho \left(x\right)}} (where \alpha \left(x\right)=\sqrt{\lambda \left(x\right)/\rho \left(x\right)} ); b is the correlation length; and {\xi }^{\mathrm{2}} denotes variance of fluctuations in the stochastic medium.

By comparing results from Eq. (9) for a numerical method with the predicted one from Eq. (10), one can estimate performance of the numerical method.

In this study to stabilize numerical methods, the regularization concept is used for simulation of permanent displacements around sources and wave propagations in stochastic media. For wave propagations in stochastic media, regularization effects on numerical solutions are investigated numerically by a benchmark problem in Sec.€7. There, a 1-D wave propagation problem and corresponding numerical scattering attenuation ({{\displaystyle Q}}_{\mathrm{s}}^{-\mathrm{1}}) are compared with those resulted from the single scattering theory [⁷⁷].

Finally it should be mentioned that different numerical approaches were developed to cure such artificial oscillations (e.g., around discontinuities or in stochastic media). For a general review, see Appendix A.

Multiresolution-based grid adaptations

Wavelets are known as a mathematical microscope: different local features of data can simultaneously be detected and studied in different resolutions; for details, see Appendix B. In this theory, the scaled-translated wavelet functions \psi \left({\mathrm{2}}^{j}x-n\right) (with scale 2^j at spatial location {\displaystyle \frac{\mathrm{2}n+\mathrm{1}}{{\mathrm{2}}^{j+\mathrm{1}}}} ) measure local variations in data. Intensity of such variations is then reflected in the wavelet coefficients d_j,n. These coefficients can be considered as a criterion for grid adaptations as they measure local variations. Wavelet coefficients of sufficiently large values and corresponding spatial points are preserved while remaining grid points belonging to remaining wavelet functions are omitted. It should be noted that in the non-redundant MRA (based on the pyramid algorithm [⁹⁹]), each wavelet (detail) or scale coefficient is uniquely connected to a specific point of the underlying grid, and hence, the adapted grid would also be unique.

For the D-D wavelet family (see Appendix B), the adaptation criterion is magnitude of wavelet (detail) coefficients. The grid points belonging to detail functions are set to zero (i.e., removed) if corresponding wavelet (detail) coefficients are smaller than predefined thresholds {{\displaystyle ϵ}}_j ({{\displaystyle \in }}_j can be different for each resolution level j), this means:

(11)

where {{\displaystyle \widehat{d}}}_{j,n} denotes thresholded wavelet coefficients. With this thresholding technique, a 1-D uniform grid can be adapted in the multiresolution framework. In Ref. [⁵⁴], it was also shown that, such thresholding has a bounded truncation error which approaches to zero as {{\displaystyle ϵ}}_j\to \mathrm{0} .

The 1-D transform and reduction technique can easily be extendable to 2-D grid points by the dimension-by-dimension transform; wherein, the 1-D transform is performed independently for each direction; for more details see Ref. [⁵⁴].

For resolution of PDEs, the previous procedure of a scalar PDE can be modified to reflect behaviors of all equations. This means, the resultant adapted grid is simply the superposition of all adapted grids; in this case, each adapted grid corresponds to a variable of PDEs.

The 1-D and 2-D grid adaptation algorithms are sensitive for spurious oscillations, developed commonly in numerical solutions of wave equations. Due to the lack of inherent dissipation of wave equations [⁵²], they propagate throughout computational domains and cause improper grid adaptation (and even the numerical instability). These non-dissipative spurious oscillations develop even in problems having smooth solutions. For more discussions and a 2-D example of such erroneous grid adaptation, please see Ref. [⁵²]. In Section 5, effects of noise and spurious oscillations are studied in more detail by some numerical examples.

Tikhonov regularization methods

In the previous section, the MRA-based grid adaptation is explained for adaptive simulations. Such simulations, however, can be sensitive for grid irregularities. It is easy to show that sampling of continuous functions can be viewed as the first kind Fredholm equation. This means, there is a close relationship between the spatial sampling theory (especially the non-uniform one) and ill-posed problems; for more details please see Appendix C. To cure this grid sensitivity, in this study, the Tikhonov method is utilized.

Ill-posed problems and the Tikhonov regularization method

A problem is called ill-posed if it violates one of the following conditions: the solution exists; the solution is unique and stable. The solution is said to be stable, if any perturbation of coefficients, parameters, initial or boundary conditions causes a little change in answers [³⁶].

Ill-posed problems have certain solutions; however, they should be obtained with especial treatments: regularization schemes. In these methods, in general, extra information is used to estimate real answers. By the regularization methods, ill-posed problems are replaced with nearly well-posed ones. In the famous Tikhonov regularization method, the regulated answer is obtained by minimizing a functional obtained by a linear combination of a residual norm \rho {\left(f\right)}^{\mathrm{2}}={\Vert f-y\Vert }_2^2 and additional information \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}, i.e.:

(12)\rho {\left(f\right)}^{\mathrm{2}}+\overline{\lambda }\times \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}={{\displaystyle \min }}_y,

where f and y denote the unknown regularized solution and raw (or noisy) data, respectively; \overline{\lambda } is a smoothing or regularization parameter belonging to range \overline{\lambda }\in \left[\mathrm{0},\infty \right). This parameter controls regularization effects [³⁵,³⁶]. One of the common constraints is smoothness, defined as \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}={\displaystyle \int {\left|{\mathrm{d}}^{m}f\left(x\right)/\mathrm{d}{x}^m\right|}^{\mathrm{2}}}\mathrm{d}x. This assumption for \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}} leads to the famous smoothing spline method of order 2m−1 [³⁰]. For this case, fast algorithms are developed [⁵²,¹⁰⁰,¹⁰¹].

To confine variation range of the smoothing parameter \overline{\lambda }, one can define a new parameter p as: \overline{\lambda }=\left(\mathrm{1}-p\right)/p; in this case, the variation range becomes p\in \left[\mathrm{0},\mathrm{1}\right]. For p=1 and p=0, the interpolation and the least-square based straight line fit are obtained, respectively. For p<1, the interpolating property vanishes while the smoothing feature increases.

For adaptive estimations, a variable smoothing parameter can be used. For this purpose, Eq. (12) can be rewritten as:

(13)\rho {\left(f\right)}^{\mathrm{2}}+\overline{\lambda }\left(x\right)\times \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}={{\displaystyle \min }}_y,

where \overline{\lambda }\left(x\right) denotes the variable smoothing parameter. For the implementation algorithm of Eq. (12) or Eq. (3), see Refs. [³⁴,¹⁰⁰,¹⁰¹].

It can be shown that estimation errors in the Tikhonov-based regularization are comparable with those obtained by the wavelet Galerkin scheme [⁵⁴]. In fact, effects of smoothing splines are similar to those of wavelet-based soft thresholding methods [³⁰].

Proper selection of the smoothing parameter (by computational methods) is not generally reliable. The L-curve approach (a graphical representation) is suggested to investigate regularization (smoothing) effects [³⁵]. In this method, residual norms in estimations are presented against regularization (semi) norms in Log-Log scales; horizontal and vertical axes correspond to the residual and regularization norms (or semi-norms), respectively. With this graph, the trade-off between these two quantities could apparently be investigated, especially the over-smoothing and under-smoothing phenomena. The L-curve method reveals significant role of the Tikhonov regularization scheme in comparison with other regularization methods. The Tikhonov regularization method leads to the smoothest possible result for a certain estimation error [³⁵,³⁶]. This helps investigate performance of Tikhonov-based smoothing and other commonly used methods developed for elastodynamic problems described in Section 7.1 by a 1-D benchmark problem.

For 2-D regularizations, at each grid point, regularized solutions are separately obtained in each direction; i.e., {{\displaystyle \left({{\displaystyle f}}_i\right)}}_z:z\in \left\{x,y\right\} . The average value {{\displaystyle f}}_i=\left\{{{\displaystyle \left({{\displaystyle f}}_i\right)}}_x+{{\displaystyle \left({{\displaystyle f}}_i\right)}}_y\right\}/\mathrm{2} is then used as an estimation of the 2-D regularized solution. This is because, adapted grids do not have a tensor product configuration; hence, the well-known thin plate spline based smoothing with constraint \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}={\displaystyle \iint \left\{{{\displaystyle \left({\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial x^{\mathrm{2}}}}\right)}}^{\mathrm{2}}+\mathrm{2}{{\displaystyle \left({\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial x\partial y}}\right)}}^{\mathrm{2}}+{{\displaystyle \left({\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial y^{\mathrm{2}}}}\right)}}^{\mathrm{2}}\right\}}\mathrm{d}x\mathrm{d}y cannot be used for 2-D adapted grids.

Choosing of the smoothing parameter p (or \overline{\mathit{\lambda }} )

As mentioned in the Introduction section, computational methods developed for assessment of the smoothing parameter p(\overline{\mathit{\lambda }}) may lead to over or under estimated results. For this reason, the trial and error method will be useful. Our study shows that even marginal values of p (p\approx \mathrm{0.99} or p\approx \mathrm{0.98}) can considerably improve numerical stability of wave propagation problems in stochastic-like media [³⁴]. Also for controlling of spurious oscillations (developed around discontinuous or high-gradient solutions) lower p values are recommended, as: 0.8<p<0.95 [⁵¹,⁵²].

A benchmark problem

In the following, estimation errors of a discontinuous test function is investigated; let us assume the test function to be: y\left(x\right)=\sin \left(\mathrm{2}\mathrm{\pi }x\right)-H\left(x-\mathrm{0.5}\right)+\mathrm{0.5}. In simulations, adaptations of both grids and the smoothing parameter are performed; grids are adapted by the D-D wavelet with M=2. For smoothing, the cubic smoothing spline (the Tikhonov regularization where m=2) with constant (p) and variable (p (x)) smoothing parameters are considered. In simulations, assumed parameters are: \in ={\mathrm{10}}^{-\mathrm{3}} (the threshold value); {{\displaystyle J}}_{\max }=\mathrm{10} (the finest resolution level); and {{\displaystyle J}}_{\min }=\mathrm{4} (the coarsest resolution level). In Fig. 1, the results are illustrated; In Figs. 1(a), 1(b), and 1(c), the solid and dashed lines are estimated functions f€(x) and the parameter p(x), respectively. For each result, corresponding estimation error (i.e., \left|f\left(x\right)-y\left(x\right)\right| ) is illustrated in the same column, at the last row.

The results offer that: 1) Adapted grid points are properly concentrated around the discontinuity at different levels of resolution; 2) Tikhonov-based estimations are stable on irregular grids; 3) By the adaptive estimation (with the value range \mathrm{0.9}\le {{\displaystyle p}}_i\le \mathrm{0.99} ), error magnitudes can be controlled both in smooth regions and around discontinuities. In the adaptive estimation, errors in the smooth and the discontinuous zones are near to those obtained with constant parameters {{\displaystyle p}}_i=\mathrm{0.9} and {{\displaystyle p}}_i=\mathrm{0.99}, respectively; 4) Significant errors develop around the discontinuity, due to the Gibb’s phenomenon. This phenomenon and the numerical dissipation can be controlled by grid adaptations; 5) Due to the Runge phenomenon, errors increase near data edges.

For the regularization adaptation, in the next section, two different detectors are presented for distinguishing real solutions and noisy ones in the spatial domain.

Detection of non-smooth zones and parameter adaptations

As mentioned before, grid-based adaptive simulations can resolve more fine resolutions and therefore improve solution accuracy. On the other hand, the grid adaptation can introduce a new source of truncation error and this leads to more sensitivity in grid-adaptations. To control this drawback, the Tikhonov regularization is suggested in this study. The Tikhonov scheme, however, can lead to over-smoothing phenomenon on smooth zones. The method-adaptation enables us to prevent over-smoothing in numerical simulations.

Generally, for proper adaptations, physical solutions should be identified properly even in presence of noise. In the following, two detecting approaches are studied; they are: the wavelet-based graded tree; a detector based on the MINMOD limiter. Performances of them are numerically investigated for data with and without noise.

Wavelet-based graded trees [²⁵]

For a wavelet-based graded tree construction, initially, it is essential to consider distribution of adapted points in different resolution levels by the MRA: the spatio-resolution representation. For a distinct spatial region, when an element (a cell of the space-frequency domain) is included in the tree, all elements of coarser spatial resolutions are also included in that fine resolution. For clarity, such tree construction is explained in more detail in Fig. 2. There, Fig. 2(a) shows an adapted grid (in the below of the illustration), distribution of the adapted grid in different resolution levels and corresponding graded tree. In this spatio-resolution representation, solid points are MRA-based adapted points, solid-thick gray lines are leaves, and thin inclined black lines (solid or dashed ones) represent the tree. In Fig. 2, dashed horizontal or vertical gray lines denote boundaries of the spatio-resolution representation. Figure 2(b) illustrates resolution relocation of the adapted points based on the graded tree concept. This illustration shows that how a spatial domain can be decomposed without overlapping according to tree leaves. This kind of domain decomposition can be used for both grid and method adaptations.

Adaptation of the smoothing parameter p based on the graded tree

One possible algorithm for parameter adaptations can be expressed as:

1)‚Construct a wavelet-based graded tree and decompose the domain based on it.

2)‚For points belonging to the coarsest resolution, level J_min, smoothing parameter p is assumed to have a background value, i.e.: {{\displaystyle p}}_i={{\displaystyle p}}_{\max }={{\displaystyle p}}_{\mathrm{0}}, where {{\displaystyle p}}_i:=p\left({{\displaystyle x}}_i\right) (e.g., p₀=0.99 to minimize smoothing effects on smooth solutions while to preserve the numerical stability).

3)‚For points belonging to the finest resolution, level J_max−1, we set: {{\displaystyle p}}_i={{\displaystyle p}}_{\min } (such zones belong to high gradient or discontinuous solutions which need more regularizations).

4)‚For remaining resolution levels j, where J_min<j<J_max−1, by a linear interpolation, p_i values can be evaluated in a way that {{\displaystyle p}}_i\in \left[{{\displaystyle p}}_{\min },{{\displaystyle p}}_{\max }\right].

For possible choices of p_min and p_max , one can consult Section 4.2 (also see Section 5.2).

The generalized MINMOD limiter and parameter adaptations

For the estimation of the first derivative, different definitions of the generalized MINMOD limiters can be considered, as:

(14a){{\displaystyle D}}_i^{\mathrm{I}}y:=MINMOD\left(\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{1}}}},{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_{i-\mathrm{1}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_i}}\right),

{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y:=MINMOD(\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{2}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{2}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{1}}}},{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_{i-\mathrm{1}}}},

(14b){\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{2}}-{{\displaystyle y}}_{i-\mathrm{2}}}{{{\displaystyle x}}_{i+\mathrm{2}}-{{\displaystyle x}}_{i-\mathrm{2}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_i}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{2}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{2}}-{{\displaystyle x}}_i}},

{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y:=MINMOD(\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{3}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{3}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{2}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{2}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_i-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_i-{{\displaystyle x}}_{i-\mathrm{1}}}},

{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_{i-\mathrm{1}}}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_{i-\mathrm{1}}}},{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{2}}-{{\displaystyle y}}_{i-\mathrm{2}}}{{{\displaystyle x}}_{i+\mathrm{2}}-{{\displaystyle x}}_{i-\mathrm{2}}}},{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{3}}-{{\displaystyle y}}_{i-\mathrm{3}}}{{{\displaystyle x}}_{i+\mathrm{3}}-{{\displaystyle x}}_{i-\mathrm{3}}}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{1}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{1}}-{{\displaystyle x}}_i}},

(14c)\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{2}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{2}}-{{\displaystyle x}}_i}},\theta \hspace{0.17em}{\displaystyle \frac{{{\displaystyle y}}_{i+\mathrm{3}}-{{\displaystyle y}}_i}{{{\displaystyle x}}_{i+\mathrm{3}}-{{\displaystyle x}}_i}},

where,

(15)

and D^j\left({{\displaystyle y}}_i\right):={{\displaystyle D}}_i^{j}y\mathrm{f}\mathrm{o}\mathrm{r}j\in \left\{\mathrm{I},\mathrm{I}\mathrm{I},\mathrm{I}\mathrm{I}\mathrm{I}\right\}; the parameter \theta belongs to the range: \mathrm{1}\le \theta \le \mathrm{2}: smaller values of lead to more dissipative results. Since the aim is to detect real non-smooth regions, small values of can be used to ignore zones belonging to spurious oscillations. The MINMOD limiters with larger stencils (more surrounding data \left\{{{\displaystyle y}}_i\right\}) lead to more dissipative results. Hence, functions {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y and {{\displaystyle D}}_i^{\mathrm{I}}y lead to the most and the least dissipative results, respectively. This difference of dissipations between {{\displaystyle D}}_i^{\mathrm{I}}y and {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y(or {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y) may be used to distinguish real variations from noise-like ones. This is because, it can be assumed that noise has stochastic feature without any coherency with surrounding data. For a zone, when both definitions {{\displaystyle D}}_i^{\mathrm{I}}y and {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}yhave values larger than a threshold, then a physical variation can exist in that zone.

In the following, a detecting algorithm is suggested. Let us assume a function y(x) containing both real variations and noise. The detecting algorithm can be summarized as:

1)‚ Detecting of real non-smooth zones:

Evaluate the first derivative dy/dx by the MINMOD limiter {{\displaystyle D}}_i^{\mathrm{I}}y;

Compute d y/dx by the {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y (or {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y);

Detect zones like {{\displaystyle \mathrm{\Omega }}}_k where both {{\displaystyle |D}}_i^{\mathrm{I}}y| and \left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y\right| (or \left|{{\displaystyle \mathrm{D}}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y\right|) are larger than zero; such domains are likely to have physical variations.

2)‚ Investigation about features of variations: distinguishing physical variations from noise in the zones.

For this purpose, a criterion is needed; here this is based on the {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y (or \left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y\right|) values; a zone {{\displaystyle {{\displaystyle y}}_i\in \mathrm{\Omega }}}_k belongs to the physical data if (similar idea can be found for example in Refs. [¹⁰²,¹⁰³] for threshold estimations):

The first criterion:

(16a)\left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}\left(y\right)\right|\ge \mu -\beta \sigma ,\mathrm{f}\mathrm{o}\mathrm{r}{{\displaystyle y}}_i\in {{\displaystyle \mathrm{\Omega }}}_k,

or the second criterion:

(16b)\left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}\left(y\right)\right|\ge {{\displaystyle \alpha }}_{\mathrm{d}}{{\displaystyle \left(y^{\prime }\right)}}_{\max }^{\mathrm{I}\mathrm{I}},\mathrm{f}\mathrm{o}\mathrm{r}{{\displaystyle y}}_i\in {{\displaystyle \mathrm{\Omega }}}_k,

where \mu and {\sigma }^{\mathrm{2}} denote the average and variance values of \left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y\right|, respectively; \beta is a positive constant in a way that its value may be in the range \beta \in \left[\mathrm{2},\mathrm{4}\right]; {{\displaystyle \alpha }}_{\mathrm{d}} is a positive predefined constant; and {{\displaystyle \left(y^{\prime }\right)}}_{\max }^{\mathrm{I}\mathrm{I}}={{\displaystyle \max }}_i\left(\left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}\left(y\right)\right|\right). Regarding Eqs. (16a) and (16b), it is simple to show that these equations can be related to each other if:

(17){{\displaystyle \alpha }}_{\mathrm{d}}={\displaystyle \frac{\mu -\beta \sigma }{{{\displaystyle \left(y^{\prime }\right)}}_{\max }^{\mathrm{I}\mathrm{I}}}}.

The coefficient {{\displaystyle \alpha }}_{\mathrm{d}} can also be estimated directly by the trial-and-error method.

Once the real non-smooth zones are detected, the smoothing (regularization) parameter p(x) can be evaluated as:

1) ‚Set a background smoothing parameter p=p_max=p₀ for the computational domain.

2)‚ Consider p={{\displaystyle p}}_{\mathrm{0}}-\Delta {{\displaystyle p}}_{\mathrm{N}\mathrm{S}} for real smooth but high gradient zones {{\displaystyle \mathrm{\Omega }}}_k.

3)‚ For discontinuous or high gradient solutions, more regularization may be required; for such zones, gradient of solutions can be used as a detector. For this purpose, evaluate absolute normalized values of {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y (or {{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}\mathrm{I}}y) in {{\displaystyle \mathrm{\Omega }}}_k as: {{\displaystyle \left({{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y\right)}}_N=\left|{{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y\right|/\left|{{\displaystyle \left(y^{\prime }\right)}}_{\max }^{\mathrm{I}\mathrm{I}}\right|.

4)‚ Update {{\displaystyle p}}_i:=p\left({{\displaystyle x}}_i\right) values over {{\displaystyle \mathrm{\Omega }}}_k as: {{\displaystyle p}}_i={{\displaystyle p}}_{\mathrm{0}}-\mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}} -{{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}}{{\displaystyle \left({{\displaystyle D}}_i^{\mathrm{I}\mathrm{I}}y\right)}}_N, where {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}} is a real positive constant.

Based on some numerical results, in general, suitable constant values of p for numerical solutions on adaptive grids may be in the range p\in \left[\mathrm{0.84},\mathrm{0.93}\right] [⁵¹,⁵²]; where smaller values of p are for discontinuous or high gradient solutions. It was shown that even small amount of regularization can have considerable effects on the solution stability [³⁴]. In this regard, a proper value of p₀ would be p₀=0.98; to preserve accuracy of smooth-high gradient solutions \mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}} may belong to the range \mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}}\in\left[\mathrm{0.02},\mathrm{0.04}\right]; and {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}} may belong to the range {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}}\in \left[\mathrm{0.03},\mathrm{0.12}\right].

Performance of the MINMOD limiter and graded tree based algorithms

Noise-like oscillations can degrade performance of different detectors; hence, effects of physical variations and noise in solutions are studied in this subsection.

To study noise effects, two examples are considered with and without noise. In Fig. 3(a), the test function y(x) is illustrated. It contains both smooth-continuous and discontinuous variations, where the continuous one has small amplitude. Such variations can occur around earthquake sources; the discontinuous solution represents permanent displacements around the sources and the smooth variation shows radiated waves. Figures 3(b)– 3(c) and Figs. 3(d)–3(e) correspond to results from the MINMOD limiter and graded tree based approaches, respectively. In Fig. 3(d), non-smooth regions, detected by the MINMOD-based algorithm are shown by gray-filled boxes. For the domain detection and the parameter adaptation, assumed parameters are: (1) For the MINMOD-based approach: p₀=0.98, \alpha =\mathrm{0.0075}, \mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}}=\mathrm{0.02} and {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}}=\mathrm{0.04}; this means \mathrm{0.92}\le {{\displaystyle p}}_i\le \mathrm{0.98} ; (2) For the wavelet-based graded tree: p₀=p_max=0.98, p_min=0.92, {{\displaystyle \in }}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}} (the wavelet threshold), J_max=10 and J_min=4. The results offer that: both the MINMOD limiter and graded tree based approaches can properly detect large and small variations; for smaller variations, larger smoothing parameters are obtained to preserve them during smoothing.

Effects of noise are studied in Fig. 4; there, the variation of the random noise is presented in Fig. 4(a) which is added to the function y(x) (presented in Fig. 3). Regarding the two detectors, parameter values are the same as the previous example (Fig. 3). Performance of the graded tree and MINMOD limiter approaches are illustrated in Figs. 4(b)–4(c) and Figs. 4(d)–4(e), respectively. The results confirm that:

1) The wavelet-based graded tree is sensitive to noise, even by using a large threshold value ({{\displaystyle \in }}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{2}}).

2) The MINMOD based approach can properly handle noise effects.

Performance of different MINMOD limiters (presented in Eq. (14)) in noisy data are also studied by considering the function {{\displaystyle y}}_{\mathrm{1}}\left(x\right)=y\left(x\right)+noise . Results are illustrated in Fig. 5, where the FD operator is the second order central finite difference and the MINMOD limiters are D^{\mathrm{I}}y , D^{\mathrm{I}\mathrm{I}}y, and D^{\mathrm{I}\mathrm{I}\mathrm{I}}y. It is evident that the smoothest limiter D^{\mathrm{I}\mathrm{I}\mathrm{I}}y leads to the smoothest result. As mentioned, this difference in smoothness can be useful to distinguish physical oscillations from noise-like ones.

Algorithm of wavelet-based adaptive methods

The aim of this study is to implement the adaptation concept to capture efficiently both high-gradient and smooth solutions. This concept is employed for both method and grids adaptations. Computing grids and the regularization parameter will be adapted by the wavelet theory and the MINMOD limiter-based algorithm, respectively.

For a system of PDEs, at time step t={{\displaystyle t}}_{n-\mathrm{1}} and spatial location x_i, it is assume that the numerical solution is: {{\displaystyle \mathit{U}}}_i^{n-\mathrm{1}}:=\mathit{U}\left({{\displaystyle \mathit{x}}}_i,{{\displaystyle t}}_{n-\mathrm{1}}\right). The procedure for PDE solution in adaptive framework can be expressed as:

1) Adapt grid points by the D-D wavelet with parameter M=2 based on the previous time step solution \left\{{{\displaystyle \mathit{U}}}_i^{n-\mathrm{1}}\right\}. For points without {{\displaystyle \mathit{U}}}_i^{n-\mathrm{1}} values, they can be estimated by either a local interpolation method (e.g., the Lagrange or spline based scheme) or by smoothing splines.

2) Evaluating spatial derivatives on the adapted grid; this can be done by the Fornberg [¹⁰⁴,¹⁰⁵] fast and recursive scheme. For finite data, to remedy one-sided derivation effects, end padding strategy will be used. For second derivative, the antisymmetric padding [⁵¹,⁵²] and for first derivatives the symmetric extension will be utilized. So, extra non-physical fluctuations, deduced by one sided derivatives, are reduced around data edges. First and second derivatives will have second order accuracies, and so a high-order adaptive numerical scheme is achieved [⁴³,¹⁰⁶]. In antisymmetric and symmetric extensions, respectively the second and first derivatives are equal to zero at the end boundary points. Therefore, in 2-D simulations, the front of the adapted points moves smoothly and gradually.

3) To have semi-discrete representation of PDEs and then using of “method of lines” method, PDEs are first discretized in spatial domain. Standard time-stepping methods can be used to solve resulted ODEs at the time t={{\displaystyle t}}_n (such as the Runge-Kutta 4th order scheme).

4) Adaptive denoising of spurious oscillations directly on non-uniform grids by the Tikhonov scheme. The smoothing parameters (set \left\{p\left({{\displaystyle x}}_i\right)\right\} or \left\{\overline{\lambda }\left({{\displaystyle x}}_i\right)\right\}) are spatially adapted by the algorithm using the MINMOD limiters. Also, for the regularization (smoothing), the constraint is assumed to be \mathrm{\Omega }{\left(f\right)}^{\mathrm{2}}={\displaystyle \int \overline{\lambda }\left(x\right)}{\left|{\mathrm{d}}^{\mathrm{2}}f\left(x\right)/\mathrm{d}{x}^{\mathrm{2}}\right|}^{\mathrm{2}}\mathrm{d}x.

5) Repeating the steps from the beginning.

Fast algorithms are developed for Tikhonov regularization method [³³]. The wavelet transforms have also very fast algorithms. Therefore, both smoothing and adapting procedures are fast and effective. To have cost-effective solvers, solutions are not smoothed or adapted at each time step. It is performed after some time steps (e.g., 10–20 steps) based on the velocity of moving fronts. To capture properly propagation of moving fronts, it is recommended to add some extra new points in different resolution levels around propagating fronts (e.g., one or two points at each side in each resolution level).

Numerical examples

One dimensional benchmark problems

In the following, three 1-D problems are studied:

1) The first example is to compare performance of the Tikhonov-based smoothing with other commonly used approaches developed for numerical simulation of mechanical wave equations. The performance is measured by the error-smoothness representation: the L-curve.

2) In the second problem, performance of the MINMOD-based detector is studied for adaptive smoothing. Here, both localized high-gradient waves and smooth ones with small magnitudes are considered. The aim is proper capturing of both small and large magnitude waves, simultaneously.

3) In the third example, effects of the Tikhonov-based smoothing are investigated on wave propagation problems in stochastic media.

The first benchmark: performance of different approaches for discontinuity simulations

In this Example, the main concern is to study effects of the numerical dissipation and dispersion for different approaches proposed for mechanical wave equations. To investigate only the smoothing effects, all simulations are done on uniform grids/meshes.

The problem is wave propagation in an elastic rod of unit length, containing a discontinuity in initial conditions. It is assumed that: 1) the wave velocity in rod is unit, i.e.: C=\sqrt{E/\rho }=\mathrm{1} , where E and \rho denote the Hooke coefficient and the density, respectively; 2) the rod is only subjected to an imposed initial displacement {{\displaystyle u}}_{\mathrm{0}}\left(x\right) with initial zero velocity ({{\displaystyle \dot{u}}}_{\mathrm{0}}\left(x\right)=\mathrm{0}). Let {{\displaystyle u}}_{\mathrm{0}}\left(x\right) to be:

(18)

This function has a discontinuity at x=0.5, where x\in \left[\mathrm{0},\mathrm{1}\right]. The pattern of this displacement is similar to permanent displacements formed around earthquake sources.

Considered conventional schemes are:

1)The FE method with linear spatial shape functions using the Runge-Kutta 4th order time integration method,

2)The discontinuous time Galerkin scheme (DTG) [³⁷]; in the spatio-temporal domains, the shape functions are also linear,

3)The common (second order) Taylor-Galerkin method (CTG) and higher order one (third order) (HOTG) [³⁸]; for CTG case, corresponding parameters are: \gamma =\mathrm{0.5} and \alpha =\mathrm{0.5}, and for HOTG, two different sets are considered, as: and ,

4)The FE-based discretization having linear shape functions with the α-generalized dissipative time integration scheme (α-GDTI); the assumed values of parameters are: \gamma \in \left\{\mathrm{0.52},\mathrm{0.6},\mathrm{0.7}\right\} and α can have different values ({{\displaystyle \alpha }}_{\min }\le \alpha \le {{\displaystyle \alpha }}_{\max } ) [³⁹],

5)The FE method with linear shape functions; the time-integration is performed by the proposed scheme in Ref. [¹⁰⁷] (FE-(GB-TI)).

In each scheme, the considered parameters are those suggested in corresponding references.

Results of these methods are compared with those obtained from the cubic smoothing spline (SS) in the L-Curve diagram at time t=0.2 (see Fig. 6). There, smoothness and estimation errors are measured by {{\displaystyle \Vert \Omega \Vert }}_{\mathrm{2}}^{\mathrm{2}}={\displaystyle \int {\left|{\mathrm{d}}^{\mathrm{2}}f\left(x\right)/{\mathrm{d}x}^{\mathrm{2}}\right|}^{\mathrm{2}}}\mathrm{d}x and {{\displaystyle \Vert \rho \Vert }}_{\mathrm{2}}^{\mathrm{2}} , respectively. It should be mentioned that, in the FE formulations, to minimize the Gibb’s phenomenon, in each element linear shape functions are used. Number of uniform elements is 255; hence, more than ten elements are used in each wavelength. To minimize oscillations in numerical solutions, it is also assumed \mathrm{\Delta }t\cdot {{\displaystyle L}}_{\mathrm{e}}/C to satisfy the CFL condition; where L_e is the element length and C is the wave velocity (here the time step is assumed to be \mathrm{\Delta }t=\mathrm{0.8}\left\{\left({{\displaystyle L}}_{\mathrm{e}}=\mathrm{1}/\mathrm{255}\right)/\left(C=\mathrm{1}\right)\right\}=\mathrm{0.003}) [¹⁰⁸].

It is evident that for a distinct error in estimation, the Tikhonov regularization (smoothing spline) curve leads to the smoothest possible result: solutions of the other methods are above the SS curve. The numerical results are also compared with the exact solution at t = 0.2 in Fig. 7; there, the light-gray and black lines are the exact and numerical solutions, respectively. The results offer that:

1) ‚Solutions are considerably improved even for small amount of smoothing in the Tikhonov method. This is apparent in top row of Fig. 7, by comparing results of the “SS” method with smoothing parameters p=1 (the FD method without any smoothing effect) and p=0.99 (with small amount of smoothing).

2) ‚Methods “SS”, “TDG”, and “ α-GDTI” may lead to reliable results for discontinuity simulations in elastodynamic problems.

The second benchmark: smoothing effects on waves with different magnitudes

Again, a wave propagation problem in a finite bar is studied. Wave velocity is assumed to be C=1. The bar is subjected to an imposed initial displacement {{\displaystyle U}}_{\mathrm{0}}\left(x\right)= \mathrm{E}\mathrm{x}\mathrm{p}\left(-\mathrm{1000}{\left(x-\mathrm{0.8}\right)}^{\mathrm{2}}\right)+\mathrm{0.1}\mathrm{E}\mathrm{x}\mathrm{p}\left(-\mathrm{200}{\left(x-\mathrm{1.2}\right)}^{\mathrm{2}}\right); hence, waves of both large and small magnitudes are included in solutions. To investigate smoothing effects, both constant and adaptive smoothing schemes are considered. In the latter case, p values are determined by the MINMOD-based algorithm; for such adaptation, the parameters are: p₀=0.98, {{\displaystyle \alpha }}_{\mathrm{d}}=\mathrm{0.0075} , \mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}}=\mathrm{0.01} and {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}}=\mathrm{0.07} . By these parameters, p values are in the range: {{\displaystyle p}}_{\min }\le p\le {{\displaystyle p}}_{\max }, where p_min=0.9 and p_max=0.98. For constant p values, we set: p=p_min=0.9. Numerical results are presented in Fig. 8 at t = 0.4. Figures 8(a)–8(b) and 8(c)–8(d) are obtained from adaptive and non-adaptive regularizations, respectively. The results confirm that:

1)‚The adaptive smoothing approach leads to smaller errors for waves of small and large amplitudes.

2)‚Spurious oscillations cannot affect the MINMOD-limiter based adaptation algorithm (Fig. 8(b)).

The third benchmark problem: wave propagation in a 1-D stochastic medium

A wave propagation problem in an infinite rod is considered. The bar consists of both homogeneous and heterogonous materials. Simulated length of the rod belongs to the range \mathrm{0}\le x\le \mathrm{2}; for ranges and the bar is homogeneous and for \mathrm{0.5}\le x\le \mathrm{1.5} it is heterogonous (this means z₀=1 in Section 2.2.2). In the homogeneous parts, we assume: {{\displaystyle \rho }}_{\mathrm{0}}=\mathrm{1} (density) and {{\displaystyle E}}_{\mathrm{0}}=\mathrm{1} (the modulus of elasticity) as background properties. For generation of the stochastic medium, the 1-D exponential autocorrelation function is used (defined as \mathrm{\Phi }\left(x\right)={\xi }^{\mathrm{2}}\mathrm{E}\mathrm{x}\mathrm{p}\left(-\left|x\right|/b\right) ). At x=0 and x=2, absorbing boundary conditions are considered. It is also assumed, there is a time-dependent boundary excitation in the hom-ogenous part, as: u\left(t\right)={{\displaystyle \mathbf{1}}}_{\left[\mathrm{0},T\right]}\left[\sin \left({\displaystyle \frac{\mathrm{2}\mathrm{\pi }}{T}}t\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sin \left({\displaystyle \frac{\mathrm{4}\mathrm{\pi }}{T}}t\right)\right]; where {{\displaystyle \mathbf{1}}}_{\left[\mathrm{0},T\right]} is the unit box function, with unit value over the interval \left[\mathrm{0},T\right] , where T denotes the period of the excitation in time. Two types of wave propagation problems are studied:

1)‚ Long period waves propagating in the stochastic medium with large correlation length, b (i.e., gradual variation of properties). For this case, assumed parameters are: b=0.30, \xi =\mathrm{0.20} (deviation of fluctuations in the autocorrelation function), T=0.30 (period of the excitation) and p=\left\{\mathrm{0.90},\mathrm{0.98},\mathrm{1.0}\right\} (the smoothing parameter).

2) ‚Waves of small periods in the random medium with fast perturbation in properties (i.e., domains with small correlation length, b). Values of the parameters are: b=0.05, \xi =\mathrm{0.30}, T=0.10, and p=\left\{\mathrm{0.90},\mathrm{0.98},\mathrm{1.0}\right\}.

Effects of the regularization are presented in Fig. 9 where the parameter controls the results. In this figure, Figs. 9(a)–9(c) and 9(d)–9(f) correspond to the problems with the long and small period waves, respectively. The dashed and solid lines denote the scattering attenuation functions resulted from the single scattering theory and numerical simulations, respectively. The attenuation factors for numerical results are calculated by Eq. (9). The results offer that:

1)‚ The smoothing procedure does not change the attenuation function considerably in the case of the long period incident waves (Figs. 9(a)–9(c)),

2)‚ For incident waves of small wave numbers, the attenuation factor is sensitive especially for large values of wave numbers; that is, for large kb values. For this range, the attenuation functions decrease significantly (Figs. 9(d)–9(f)),

3)‚ A considerable smoothing may also lead to artificial amplifications (e.g., see Figs. 9(c) and 9(f)).

In conclusion, authors’ experiences offer that:

1)‚ The regularized numerical method has a long-term stability even for large deviation values of fluctuations in material properties (\xi), most of numerical schemes become unstable for such media;

2)‚ A marginal amount of the regularization can improve considerably the numerical stability (with p values near to one) and preserve numerical stability [³⁴].

In this regard, using of small amount of the regularization (i.e., p\to \mathrm{1} ) is suggested for simulations on stochastic media.

Some 2-D numerical examples with earthquake sources

The following examples are two P-SV source problems selected for simulation of both discontinues and transient solutions. For this reason, a P-SV point source and a P-SV rupturing line source are studied. For the point source example, simulations are done without and with a fault zone. This zone contains materials with stochastic nature and lower average wave velocities than surrounding domain; this is due to crushing and softening of the materials.

In numerical simulations, the main assumptions are:

1)‚ Applying the D-D interpolating wavelet of order 3 (M=2).

2) ‚Decomposing 2-D grids in three levels, where the finest and the coarsest resolution levels are J_max=8 and J_min=5, respectively. This means, for spatial domain \mathrm{\Omega }\in \left[\mathrm{0},L\right]\times \left[\mathrm{0},L\right] , in each dimension, the spatial sampling step is \mathrm{d}x=L/{\mathrm{2}}^{\mathrm{8}} at the finest resolution.

3)‚ Simulation of earthquake sources by the dislocation theory.

4)‚ Using of the MINMOD-based detector.

5)‚ Using of the Tikhonov regularization with the adaptive smoothing parameter.

Example 1: A P-SV point source. The schematic shape of computational domain and corresponding four absorbing boundaries are illustrated in Fig. 10. In this figure, domains {{\displaystyle \mathrm{\Omega }}}_{\mathrm{1}} and {{\displaystyle \mathrm{\Omega }}}_2 denote the back-ground and fault zone (heterogeneous) domains, respectively; “S” shows location of the point source; and {{\displaystyle r}}_i represents receivers. Orientation of the source can be described as: a thrust fault (\eta =\mathrm{\pi }/\mathrm{2} ) with a dip angle (\zeta) placed parallel to the E-W direction (\varphi =\mathrm{\pi }/\mathrm{2}). For this case, components of the moment tensor M (in Eq. (4)) become:

{{\displaystyle M}}_{xx}=-{{\displaystyle M}}_{zz}=-{{\displaystyle M}}_{\mathrm{0}}\sin \mathrm{2}\zeta ,

{{\displaystyle M}}_{xz}={{\displaystyle M}}_{zx}=-{{\displaystyle M}}_{\mathrm{0}}\cos \mathrm{2}\zeta ,

(19){{\displaystyle M}}_{yy}={{\displaystyle M}}_{xy}={{\displaystyle M}}_{yz}=\mathrm{0},

where, x and z correspond to east and vertical directions, respectively. Regarding Eqs. (3) and (4), the equivalent body forces are:

(20)\begin{array}{l}{{\displaystyle f}}_x=-\left\{{{\displaystyle M}}_{xx}{{\displaystyle \partial }}_x\delta \left(r-{{\displaystyle r}}_{\mathrm{s}}\right)+{{\displaystyle M}}_{xz}{{\displaystyle \partial }}_z\delta \left(r-{{\displaystyle r}}_{\mathrm{s}}\right)\right\}D\left(t\right),\hfill \\ {{\displaystyle f}}_y=\mathrm{0},\hfill \\ {{\displaystyle f}}_z=-\left\{{{\displaystyle M}}_{xz}{{\displaystyle \partial }}_x\delta \left(r-{{\displaystyle r}}_{\mathrm{s}}\right)+{{\displaystyle M}}_{zz}{{\displaystyle \partial }}_z\delta \left(r-{{\displaystyle r}}_{\mathrm{s}}\right)\right\}D\left(t\right).\hfill \end{array}

If there is not any structural variation in the y direction, the fault source can be implemented in a 2-D domain as a 90 dip-slip fault activated on the vertical (x-z) plane.

The source displacement function D(t) is assumed to have a ramp excitation model where starting time of fault slip is t_s = 0, the rise-time (slipping duration) is t_r = 0.1 and the maximum amplitude of dislocation is {{\displaystyle D}}_{\infty } (occurred at t=t_r+t_s). It is assumed that the final dislocation ({{\displaystyle D}}_{\infty }) and the average dislocation are the same in the entire fault plane (D=\overline{D} ). In numerical simulations, it is assumed that: {{\displaystyle M}}_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}} ; \Omega \in \left[\mathrm{0},L\right]\times \left[\mathrm{0},L\right] (the computational domain); L=10km, dt=0.0003175 (the time step); and \in =\mathrm{0.25}\times {\mathrm{10}}^{-\mathrm{8}} (the wavelet threshold). For the background domain (W₁) we have: {{\displaystyle \rho }}_{\mathrm{1}}=\mathrm{2200}\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}, {{\displaystyle \left({{\displaystyle V}}_{\mathrm{s}}\right)}}_{\mathrm{1}}=\mathrm{1800}\mathrm{m}/\mathrm{s}, {{\displaystyle \left({{\displaystyle V}}_{\mathrm{p}}\right)}}_{\mathrm{1}}=\mathrm{3150}\mathrm{m}/\mathrm{s}; for the fault zone (W₂), they are: {{\displaystyle \rho }}_{\mathrm{2}}=\mathrm{1830}\left(\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}\right), {{\displaystyle \left({{\displaystyle V}}_{\mathrm{s}}\right)}}_{\mathrm{2}}=\mathrm{1500}\mathrm{m}/\mathrm{s}, {{\displaystyle \left({{\displaystyle V}}_{\mathrm{p}}\right)}}_{\mathrm{2}}=\mathrm{2630}\mathrm{m}/\mathrm{s}. For the adaptive smoothing, we assume: p₀=0.98, {{\displaystyle \alpha }}_{\mathrm{d}}=\mathrm{0.0075} , \mathrm{\Delta }{{\displaystyle p}}_{\mathrm{N}\mathrm{S}}=\mathrm{0.02} and {{\displaystyle \alpha }}_{\mathrm{N}\mathrm{S}}=\mathrm{0.075} . With these parameters, p values are in the range: \mathrm{0.885}\le p\left({{\displaystyle x}}_i\right)\le \mathrm{0.98}. All initial conditions are assumed to be zero. For the absorbing boundary conditions we have: a_x=a_z=30 and b_x=b_z=−120.

In the following, two different simulations are done without and with the fault zone.

The domain without the fault zone: the snapshots of solutions u_x, u_z and corresponding adapted grids are shown in Fig. 11. They are obtained at 0.076, 2.108, 3.203, and 5.0762 s. These snapshots are selected to investigate:

1) ‚Effects of the energy-release in the source at initial time steps and its effects on the grid adaptation (see Figs. 11(a)–11(b));

2)‚ Forming of both transient and permanent displacements;

3)‚ Developing of both P and SV propagating fronts. In Figs. 11(d)–11(e)), these waves are completely separated from each other due to the speed difference;

4)‚ Studying performance of the absorbing boundary conditions (Figs. 11(j)–11(k)).

The results offer that permanent displacements induced by the dislocation equivalent forces have a four-lobed pattern around the source zone (Fig. 11); it diminishes with the distance with the relationship 1/r. The transient propagating waves are also attenuated with the distance with the rate \mathrm{1}/\sqrt{r} [¹⁰⁹]. Hence, in far-fields, only transient displacements are noticeable in wave-fields. Permanent displacements lead to permanent adapted grid points while propagating fronts yield moving adapted points.

Solutions obtained by the proposed scheme with the constant and adaptive smoothing parameter (\overline{\lambda }) are compared in Fig. 12 with the FD method at the receivers r₁ and r₄. The FD method is based on the five point central formulation without any smoothing. All of the solutions are obtained with the same time steps (dt=0.0003175). It is clear that:

1)‚ Numerical dissipation and dispersion exist obviously in the FD method.

2)‚ The dissipation effects are considerable around singular solutions (by comparing solutions in the receiver r₁ and r₄).

3) The dispersion phenomenon affects all of the computational domain.

4) The adaptive smoothing leads to smaller dissipation and more permanent displacement near the point source.

Numerical solutions of the adaptive solver and the FD method (without filtering) around the point source are compared in Fig. 13. This comparison is important since around the source, singular solutions develop. It is evident that in the FD method, propagating spurious oscillations lead to the numerical instability. Developing of such oscillations in the vicinity of a discontinuous solution is known as the Gibb’s phenomenon. These non-dissipative spurious oscillations propagate throughout the computational domain due to advection feature of the wave equation.

The domain with the fault zone: The fault zone with a stochastic feature and thickness z₀=1.2km is generated with autocorrelation function defined by the double exponential description, Eq. (6) [⁷⁷]. The horizontal and vertical correlation lengths in the fault zone are assumed to be a=40 and b=60, respectively. The standard deviation of fluctuations, \xi , is assumed to be 20% for the density (p) and the Lame parameters (\lambda and \mu). The numerical results are presented in Fig. 14. This figure offers that:

1) Some waves, denoted by T are trapped in the fault zone.

2) Multi-reflected waves exist in the zone.

3) The multi-reflected waves of different frequencies transmit and leakage continuously to the back-ground domain.

4) Due to the stochastic nature of materials in the fault zone, stochastic waves develop there.

5) In the fault zone, the trapped waves propagate with smaller velocities due to the crushed and softening feature of the fault zone materials.

6) By the adaptive smoothing, it is possible to preserve high-frequency small amplitude waves; due to these waves extra adapted points are concentrated (see Figs. 11 and 14).

Example 2: A rupturing line source. The rupturing phenomenon is simulated along a line source by superposition of several point sources along the source. In the numerical simulation, the main assumptions are: the rupture velocity is 0.9 times the shear wave velocity; the spatial length of rupturing is {{\displaystyle \mathrm{\Delta }}}_{\mathrm{s}}=\mathrm{40}\times \left(\mathrm{10000}/{\mathrm{2}}^{\mathrm{8}}\right)\mathrm{m}; this associates with 40 dislocation points located in the finest resolution; rupturing direction is parallel to the axis (horizontal direction); initiating point of rupturing is at the same location of the point source “S” in the previous example; each fault segment has the same particle displacement time histories (the ramp function), released with different time delays along the line source; each point source has the same released energy; this means for each segment, the energy is: {{\displaystyle M}}_{\mathrm{0}}={\displaystyle \frac{{{\displaystyle M}}_{\mathrm{0}}^{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}{\mathrm{40}}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}/\mathrm{40} . Other parameters are the same as the previous point source example.

Snapshots of u_x and u_z solutions and corresponding adapted grids are illustrated in Fig. 15. In these illustrations, the snapshots are taken at 0.4730, 2.108, 3.203, and 5.0762 s. The snapshots are chosen to study:

1) Developing of rupturing (by comparing Figs. 15 (a)–15(b), 15(d)–15(e), 15(g)–15(h) and 15(j)–15(k)).

2) Forming of both SV and P waves during source rupturing (snapshots Figs. 15 (b)–15(c)),

3) Directivity effects (Figs. 15 (d)–15(e)),

4) Performance of absorbing boundaries and existence of permanent displacements (snapshots (j,k)).

The results offer that: 1) The adapted grid points are concentrated properly around the source line; 2) A discontinues permanent displacement (with a horizontal extension) is developed along the faulting direction in the u_x solutions; 3) For the u_z solutions, the permanent displacements are only concentrated at the ends of the rupturing line; 4) As the total energy is evenly distributed along the source line and released with some time delays, the transient waves radiated from the source have smaller amplitudes and lower frequencies than those of the point source.

Conclusions and future studies

Conclusions

A grid-based adaptive scheme reinforced by an adaptive regularization is proposed for simulation of stress wave propagation problems, such as: Seismic sources with permanent displacements and emitted waves; wave propagations through stochastic media.

The main steps of simulations are: The grids are adapted via the MRA; the Tikhonov-based regularization is used as a post-processor for controlling the Gibb’s phenomenon and spurious oscillations; to preserve physical waves with small and large magnitudes, an adaptive estimation is proposed by using an adaptive smoothing directly on non-uniform grids. A new detector is proposed based on the MINMOD limiters for identifications real and noisy solutions directly on non-uniform grids. This is especially important for simulations of propagating waves in stochastic media. The numerical results offer that:

1) The main novelty is preserving of both large and small physical waves in homogenous and heterogeneous (stochastic) media by the MINMOD-based detector used for adaptive regularizations performed directly on non-uniform grids.

2) To distinguish physical small waves from numerical (artificial) oscillations, two approaches are examined as detectors: the wavelet-based graded trees, the MINMOD limiter-based algorithm (proposed in this work). The former is sensitive to noise, while the latter can properly detect fine features even in noisy solutions.

3) The wavelet graded tree is suitable for numerical methods having the TVD feature (such as: essentially non-oscillatory (ENO), weighted ENO (WENO), central or central-upwind high resolution methods). For schemes without the TVD property, other approaches should be considered.

4) Regarding stochastic media, long-waves are not so sensitive for smoothing (in the post-processing stage), while short-waves are. Small amount of the regularization can guarantee solution stability (Fig. 13) and preserve predicted scattering attenuations in stochastic media (see Fig. 6).

5) A proper adaptive regularization can capture more precisely real solutions, like those around singular solutions: seismic sources (e.g., permanent displacement around sources; see Fig. 12).

6) In the one of the one-dimensional benchmarks, the proposed method is compare with other commonly used methods developed fundamentally for elastodynamic problems. The Tikhonov-based spatial filtering leads to the smoothest results (see Fig. 7). There, some schemes are based on the weak-form formulations, while in the proposed scheme the strong form is directly resolved. In the weak formulations, solutions have C⁰ continuity, while in the proposed approach, it is C².

7) The proposed method is simple and straightforward as all calculations are done in the physical spaces (not wavelet spaces).

Finally, it should be mentioned that disadvantages of this regularization-based smoothing is proper selection of the regularization (smoothing) parameter. In fact, the trial and error method would be useful; in general, the recommended range would be \mathrm{0.85}\le P\le \mathrm{0.98} .

Future studies

This study can be improved by considering: thermal coupling effects in sources [⁵⁹,¹¹⁰], coupled waves including magneto-electro-thermal behavior (important for earthquake early warning) [¹¹¹–¹¹³], coupled waves in porous media [¹¹⁴].

References

[1] Day S M, Dalguer L A, Lapusta N, Liu Y. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture. Journal of Geophysical Research. Solid Earth, 2005, 110(B12): B12307

[2] Dalguer L A, Day S M. Comparison of fault representation methods in finite difference simulations of dynamic rupture. Bulletin of the Seismological Society of America, 2006, 96(5): 1764–1778

[3] Dalguer L A, Day S M. Staggered-grid split-node method for spontaneous rupture simulation. Journal of Geophysical Research. Solid Earth, 2007, 112: B02302

[4] Day A, Steven M. Three-dimensional finite difference simulation of fault dynamics: Rectangular faults with fixed rupture velocity. Bulletin of the Seismological Society of America, 1982, 72: 705–727

[5] Madariaga R, Olsen K, Archuleta R. Modeling dynamic rupture in a 3D earthquake fault model. Bulletin of the Seismological Society of America, 1998, 88: 1182–1197

[6] Moczo P, Kristek J, Galis M, Pazak P, Balazovjech M. Finite-difference and finite-element modeling of seismic wave propagation and earthquake motion. Acta Physica Slovaca, 2007, 57(2): 177–406

[7] Duan B, Day S M. Inelastic strain distribution and seismic radiation from rupture of a fault kink. Journal of Geophysical Research. Solid Earth, 2008, 113(B12): B12311

[8] Galis M, Moczo P, Kristek J A. 3-D hybrid finite-difference—finite-element viscoelastic modelling of seismic wave motion. Geophysical Journal International, 2008, 175(1): 153–184

[9] Barall M. A grid-doubling finite-element technique for calculating dynamic three-dimensional spontaneous rupture on an earthquake fault. Geophysical Journal International, 2009, 178(2): 845–859

[10] Ely G P, Day S M, Minster J B. A support-operator method for 3-D rupture dynamics. Geophysical Journal International, 2009, 177(3): 1140–1150

[11] Aagaard B T, Heaton T H, Hall J F. Dynamic earthquake ruptures in the presence of lithostatic normal stresses: Implications for friction models and heat production. Bulletin of the Seismological Society of America, 2001, 91(6): 1756–1796

[12] Kaneko Y, Lapusta N, Ampuero J P. Spectral element modeling of spontaneous earthquake rupture on rate and state faults: Effect of velocity-strengthening friction at shallow depths. Journal of Geophysical Research. Solid Earth, 2008, 113: B09317

[13] Kaneko Y, Ampuero J P, Lapusta N. Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics. Journal of Geophysical Research. Solid Earth, 2011, 116(B10): B10313

[14] Galvez P, Ampuero J P, Dalguer L A, Somala S N, Nissen-Meyer T. Dynamic earthquake rupture modelled with an unstructured 3-D spectral element method applied to the 2011 M9 Tohoku earthquake. Geophysical Journal International, 2014, 198(2): 1222–1240

[15] Tada T, Madariaga R. Dynamic modelling of the flat 2-D crack by a semi-analytic BIEM scheme. International Journal for Numerical Methods in Engineering, 2001, 50(1): 227–251

[16] Lapusta N, Rice J R, Ben-Zion Y, Zheng G. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction. Journal of Geophysical Research. Solid Earth, 2000, 105(B10): 23765–23789

[17] Andrews D J. Dynamic plane-strain shear rupture with a slip-weakening friction law calculated by a boundary integral method. Bulletin of the Seismological Society of America, 1985, 75: 1–21

[18] Das S. A numerical method for determination of source time functions for general three-dimensional rupture propagation. Geophysical Journal International, 1980, 62(3): 591–604

[19] Benjemaa M, Glinsky-Olivier N, Cruz-Atienza V M, Virieux J, Piperno S. Dynamic non-planar crack rupture by a finite volume method. Geophysical Journal International, 2007, 171(1): 271–285

[20] Benjemaa M, Glinsky-Olivier N, Cruz-Atienza V M, Virieux J. 3-D dynamic rupture simulations by a finite volume method. Geophysical Journal International, 2009, 178(1): 541–560

[21] Dumbser M, Käser M, De La Puente J. Arbitrary high-order finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D. Geophysical Journal International, 2007, 171(2): 665–694

[22] Dumbser M, Käser M. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes — II. The three-dimensional isotropic case. Geophysical Journal International, 2006, 167(1): 319–336

[23] Titarev V A, Toro E F. ADER: Arbitrary high order Godunov approach. Journal of Scientific Computing, 2002, 17(1/4): 609–618

[24] Käser M, Dumbser M. An Arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes — I. The two-dimensional isotropic case with external source terms. Geophysical Journal International, 2006, 166(2): 855–877

[25] Harten A. Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Communications on Pure and Applied Mathematics, 1995, 48(12): 1305–1342

[26] Cohen A, Kaber S M, Muller S, Postel M. Fully Adaptive multiresolution finite volume schemes for conservation laws. Mathematics of Computation, 2001, 72(241): 183–226

[27] Müller S. Adaptive Multiscale Schemes for Conservation Laws. Berlin: Springer, 2003

[28] Reinsch C H. Smoothing by spline functions. Numerische Mathematik, 1967, 10(3): 177–183

[29] Reinsch C H. Smoothing by spline functions. II. Numerische Mathematik, 1971, 16(5): 451–454

[30] Unser M. Splines: A perfect fit for signal and image processing. IEEE Signal Processing Magazine, 1999, 16(6): 22–38

[31] Ragozin D L. Error bounds for derivative estimates based on spline smoothing of exact or noisy data. Journal of Approximation Theory, 1983, 37(4): 335–355

[32] Loader C. Smoothing: Local Regression Techniques. In: Gentle J E, Härdle W K, Mori Y, eds. Handbook of Computational Statistics. Berlin Heidelberg: Springer, 2012, 571–96

[33] Hutchinson M F, de Hoog F R. Smoothing noisy data with spline functions. Numerische Mathematik, 1985, 47(1): 99–106

[34] Yousefi H, Ghorashi S S, Rabczuk T. Directly simulation of second order hyperbolic systems in second order form via the regularization concept. Communications in Computational Physics, 2016, 20(1): 86–135

[35] Hansen P C. Rank-deficient and discrete ill-posed problems: Numerical aspects of linear inversion. Society for Industrial and Applied Mathematics, 1998

[36] Petrov Y P, Sizikov V S. Well-Posed, Ill-Posed, and Intermediate Problems with Applications. Leiden: Koninklijke Brill NV, 2005

[37] Li X D, Wiberg N E. Structural dynamic analysis by a time-discontinuous Galerkin finite element method. International Journal for Numerical Methods in Engineering, 1996, 39(12): 2131–2152

[38] Youn S K, Park S H. A new direct higher-order Taylor-Galerkin finite element method. Computers & Structures, 1995, 56(4): 651–656

[39] Hilber H M, Hughes T J R, Taylor R L. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics, 1977, 5(3): 283–292

[40] Alves M A, Cruz P, Mendes A, Magalhães F D, Pinho F T, Oliveira P J. Adaptive multiresolution approach for solution of hyperbolic PDEs. Computer Methods in Applied Mechanics and Engineering, 2002, 191(36): 3909–3928

[41] Cruz P, Mendes A, Magalhães F D. Using wavelets for solving PDEs: An adaptive collocation method. Chemical Engineering Science, 2001, 56(10): 3305–3309

[42] Cruz P, Mendes A, Magalhães F D. Wavelet-based adaptive grid method for the resolution of nonlinear PDEs. AIChE Journal. American Institute of Chemical Engineers, 2002, 48(4): 774–785

[43] Jameson L, Miyama T. Wavelet analysis and ocean modeling: A dynamically adaptive numerical method “WOFD-AHO”. Monthly Weather Review, 2000, 128(5): 1536–1549

[44] Alam J M, Kevlahan N K R, Vasilyev O V. Simultaneous space–time adaptive wavelet solution of nonlinear parabolic differential equations. Journal of Computational Physics, 2006, 214(2): 829–857

[45] Bertoluzza S, Castro L. Adaptive Wavelet Collocation for Elasticity: First Results. Pavia, 2002

[46] Griebel M, Koster F. Adaptive wavelet solvers for the unsteady incompressible Navier-Stokes equations. In: Malek J, Nečas J, Rokyta M, eds. Advances in Mathematical Fluid Mechanics. Berlin: Springer, 2000, 67–118

[47] Santos J C, Cruz P, Alves M A, Oliveira P J, Magalhães F D, Mendes A. Adaptive multiresolution approach for two- dimensional PDEs. Computer Methods in Applied Mechanics and Engineering, 2004, 193(3–5): 405–425

[48] Vasilyev O V, Kevlahan N K R. An adaptive multilevel wavelet collocation method for elliptic problems. Journal of Computational Physics, 2005, 206(2): 412–431

[49] Ma X, Zabaras N. An adaptive hierarchical sparse grid collocation algorithm for the solution of stochastic differential equations. Journal of Computational Physics, 2009, 228(8): 3084–3113

[50] Mehra M, Kevlahan N K R. An adaptive wavelet collocation method for the solution of partial differential equations on the sphere. Journal of Computational Physics, 2008, 227(11): 5610–5632

[51] Yousefi H, Noorzad A, Farjoodi J. Simulating 2D waves propagation in elastic solid media using wavelet based adaptive method. Journal of Scientific Computing, 2010, 42(3): 404–425

[52] Yousefi H, Noorzad A, Farjoodi J. Multiresolution based adaptive schemes for second order hyperbolic PDEs in elastodynamic problems. Applied Mathematical Modelling, 2013, 37(12–13): 7095–7127

[53] Bürger R, Ruiz-Baier R, Schneider K. Adaptive multiresolution methods for the simulation of waves in excitable media. Journal of Scientific Computing, 2010, 43(2): 261–290

[54] Holmström M. Solving hyperbolic PDEs using interpolating wavelets. SIAM Journal on Scientific Computing, 1999, 21(2): 405–420

[55] Gottlieb D, Hesthaven J S. Spectral methods for hyperbolic problems. Journal of Computational and Applied Mathmatics, 2001, 128(1–2): 83–131

[56] Cebeci T, Shao J P, Kafyeke F, Laurendeau E. Computational Fluid Dynamics for Engineers: From Panel to Navier-Stokes Methods with Computer Programs. Long Beach, California: Springer, 2005

[57] Hesthaven J S, Gottlieb S, Gottlieb D. Spectral Methods for Time-Dependent Problems. Cambridge: Cambridge University Press, 2007

[58] Yousefi H, Tadmor E, Rabczuk T. High resolution wavelet based central schemes for modeling nonlinear propagating fronts. Engineering Analysis with Boundary Element, 2019, 103: 172–195

[59] Yousefi H, Taghavi Kani A, Mahmoudzadeh Kani I. Multiscale RBF-based central high resolution schemes for simulation of generalized thermoelasticity problems. Frontiers of Structural and Civil Engineering, 2019, 13(2): 429–455

[60] Yousefi H, Taghavi Kani A, Mahmoudzadeh Kani I. Response of a spherical cavity in a fully-coupled thermo-poro-elastodynamic medium by cell-adaptive second-order central high resolution schemes. Underground Space, 2018, 3(3): 206–217

[61] Nguyen-Xuan H, Nguyen-Hoang S, Rabczuk T, Hackl K. A polytree-based adaptive approach to limit analysis of cracked structures. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 1006–1039

[62] Budarapu P R, Gracie R, Bordas S P A, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[63] Nguyen-Xuan H, Liu G R, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 252–273

[64] Badnava H, Msekh M A, Etemadi E, Rabczuk T. An h-adaptive thermo-mechanical phase field model for fracture. Finite Elements in Analysis and Design, 2018, 138: 31–47

[65] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[66] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2008, 197(6-8): 641–658

[67] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[68] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[69] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81: 48–71

[70] Anitescu C, Hossain M N, Rabczuk T. Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshes. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 638–662

[71] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1157–1178

[72] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger K U, Bazilevs Y, Rabczuk T. Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47–48): 3410–3424

[73] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S P A, Rabczuk T. Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21–22): 1892–1908

[74] Hong T K, Kennett B L N. Scattering attenuation of 2D elastic waves: theory and numerical modeling using a wavelet-based method. Bulletin of the Seismological Society of America, 2003, 93(2): 922–938

[75] Hong T K, Kennett B L N. Modelling of seismic waves in heterogeneous media using a wavelet-based method: Application to fault and subduction zones. Geophysical Journal International, 2003, 154(2): 483–498

[76] Hong T K, Kennett B L N. Scattering of elastic waves in media with a random distribution of fluid-filled cavities: Theory and numerical modelling. Geophysical Journal International, 2004, 159(3): 961–977

[77] Roth M, Korn M. Single scattering theory versus numerical modelling in 2-D random media. Geophysical Journal International, 1993, 112(1): 124–140

[78] Persson P O, Runborg O. Simulation of a waveguide filter using wavelet-based numerical homogenization. Journal of Computational Physics, 2001, 166(2): 361–382

[79] Engquist B, Runborg O. Wavelet-based numerical homogenization with applications. Multiscale and Multiresolution Methods, 2002, 20: 97–148

[80] Weinan E, Engquist B, Huang Z. Heterogeneous multiscale method: A general methodology for multiscale modeling. Physical Review. B, 2003, 67(9): 092101

[81] Weinan E, Engquist B, Li X, Ren W, Vanden-Eijnden E. Heterogeneous multiscale methods: A review. Communications in Computational Physics, 2007, 2: 367–450

[82] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[83] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[84] Talebi H, Silani M, Bordas S P, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[85] Engquist B, Holst H, Runborg O. Multiscale methods for wave propagation in heterogeneous media over long time. Numerical Analysis of Multiscale Computations, 2012, 82: 167–86

[86] Rabczuk T, Areias P, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[87] Amiri F, Anitescu C, Arroyo M, Bordas S P A, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[88] Nanthakumar S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[89] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23-24): 1391–1411

[90] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[91] Ghorashi S S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[92] Hong T K, Kennett B L N. A wavelet-based method for simulation of two-dimensional elastic wave propagation. Geophysical Journal International, 2002, 150(3): 610–638

[93] Hong T K, Kennett B L N. On a wavelet-based method for the numerical simulation of wave propagation. Journal of Computational Physics, 2002, 183(2): 577–622

[94] Wu Y, McMechan G A. Wave extrapolation in the spatial wavelet domain with application to poststack reverse-time migration. Geophysics, 1998, 63(2): 589–600

[95] Operto S, Virieux J, Hustedt B, Malfanti F. Adaptive wavelet-based finite-difference modelling of SH-wave propagation. Geophysical Journal International, 2002, 148(3): 476–498

[96] Beylkin G, Keiser J M. An adaptive pseudo-wavelet approach for solving nonlinear partial differential equations. Wavelet Analysis and Its Applications, 1997, 6: 137–97

[97] Sochacki J, Kubichek R, George J, Fletcher W R, Smithson S. Absorbing boundary conditions and surface waves. Geophysics, 1987, 52(1): 60–71

[98] Lay T, Wallace T C. Modern Global Seismology. San Diego, California: Academic Press, 1995

[99] Mallat S G. A Wavelet Tour of Signal Processing. San Diego: Academic Press, 1998

[100] Eilers P H C. A perfect smoother. Analytical Chemistry, 2003, 75(14): 3631–3636

[101] Stickel J J. Data smoothing and numerical differentiation by a regularization method. Computers & Chemical Engineering, 2010, 34(4): 467–475

[102] Hadjileontiadis L J, Panas S M. Separation of discontinuous adventitious sounds from vesicular sounds using a wavelet-based filter. IEEE Transactions on Biomedical Engineering, 1997, 44(12): 1269–1281

[103] Hadjileontiadis L J, Liatsos C N, Mavrogiannis C C, Rokkas T A, Panas S M. Enhancement of bowel sounds by wavelet-based filtering. IEEE Transactions on Biomedical Engineering, 2000, 47(7): 876–886

[104] Fornberg B. Classroom note: Calculation of weights in finite difference formulas. SIAM Review, 1998, 40(3): 685–691

[105] Fornberg B. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation, 1988, 51(184): 699–706

[106] Jameson L. A wavelet-optimized, very high order adaptive grid and order numerical method. SIAM Journal on Scientific Computing, 1998, 19(6): 1980–2013

[107] Noh G, Bathe K J. An explicit time integration scheme for the analysis of wave propagations. Computers & Structures, 2013, 129: 178–193

[108] Hulbert G M, Chung J. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. Computer Methods in Applied Mechanics and Engineering, 1996, 137(2): 175–188

[109] Vidale J, Helmberger D V, Clayton R W. Finite-difference seismograms for SH waves. Bulletin of the Seismological Society of America, 1985, 75: 1765–1782

[110] Zhuang X, Huang R, Liang C, Rabczuk T. A coupled thermo-hydro-mechanical model of jointed hard rock for compressed air energy storage. Mathematical Problems in Engineering, 2014, 2014: 1–11

[111] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[112] Ghasemi H, Park H S, Rabczuk T. A multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[113] Das B. Problems and Solutions in Thermoelasticity and Magneto-thermoelasticity. Cham: Springer International Publishing, 2017

[114] Donoho D L. Interpolating Wavelet Transforms. Stanford: Stanford University, 1992

[115] Harten A, Engquist B, Osher S, Chakravarthy S R. Uniformly high order accurate essentially non-oscillatory schemes, III. Journal of Computational Physics, 1997, 131(1): 3–47

[116] Katerinaris J A. ENO and WENO schemes. In: Dougalis V A, Kampanis N A, Ekaterinaris J A, eds. Effective Computational Methods for Wave Propagation. New York: Chapman and Hall/CRC, 2008, 521–92

[117] Gu Y, Wei G W. Conjugate filter approach for shock capturing. Communications in Numerical Methods in Engineering, 2003, 19(2): 99–110

[118] Wei G W, Gu Y. Conjugate filter approach for solving Burgers’ equation. Journal of Computational and Applied Mathematics, 2002, 149(2): 439–456

[119] Fatkullin I, Hesthaven J S. Adaptive high-order finite-difference method for nonlinear wave problems. Journal of Scientific Computing, 2001, 16(1): 47–67

[120] Leonard B P. Locally modified QUICK scheme for highly convective 2-D and 3-D flows. In: Taylor C, Morgan K, eds. Numerical Methods in Laminar and Turbulent Flow. Swansea: Pineridge Press, 1987, 35–47

[121] Leonard B P. Simple high-accuracy resolution program for convective modelling of discontinuities. International Journal for Numerical Methods in Fluids, 1988, 8(10): 1291–1318

[122] Leonard B P. The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection. Computer Methods in Applied Mechanics and Engineering, 1991, 88(1): 17–74

[123] Leonard B P. Bounded higher-order upwind multidimensional finite-volume convection-diffusion algorithms. In: Minkowycz W J, Sparrow E M, eds. Advances in Numerical Heat Transfer. Taylor and Francis, 1996, 1–57

[124] Darwish M S, Moukalled F H. Normalized variable and space formulation methodology for high-resolution schemes. Numerical Heat Transfer Part B-Fundamentals, 1994, 26(1): 79–96

[125] Bürger R, Kozakevicius A. Adaptive multiresolution WENO schemes for multi-species kinematic flow models. Journal of Computational Physics, 2007, 224(2): 1190–1222

[126] Wang L L. Foundations of Stress Waves. Amsterdam: Elsevier, 2007

[127] Day S M. Effect of a shallow weak zone on fault rupture: Numerical simulation of scale-model experiments. Bulletin of the Seismological Society of America, 2002, 92(8): 3022–3041

[128] Galis M, Moczo P, Kristek J, Kristekova M. An adaptive smoothing algorithm in the TSN modelling of rupture propagation with the linear slip-weakening friction law. Geophysical Journal International, 2010, 180(1): 418–432

[129] Ampuero J P. A physical and numerical study of earthquake nucleation. Dessertation for the Doctoral Degree. Paris: University of Paris, 2002 (In French)

[130] Festa G, Vilotte J P. Influence of the rupture initiation on the intersonic transition: Crack-like versus pulse-like modes. Geophysical Research Letters, 2006, 33(15): L15320

[131] Festa G, Vilotte J P. The Newmark scheme as velocity-stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics. Geophysical Journal International, 2005, 161(3): 789–812

[132] Trangenstein J A. Numerical solution of hyperbolic partial differential equations. Cambridge: Cambridge University Press, 2009

[133] LeVeque R J. Numerical methods for conservation laws. 2nd ed. Basel: Birkhäuser Verlag, 1992

[134] Cockburn B, Karniadakis G, Shu C W. The development of discontinuous galerkin methods. In: Cockburn B, Karniadakis G, Shu C-W, eds. Discontinuous Galerkin Methods. Heidelberg: Springer Berlin, 2000, 3–50

[135] Banks J W, Henshaw W D. Upwind schemes for the wave equation in second-order form. Journal of Computational Physics, 2012, 231(17): 5854–5889

[136] Hughes T J R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Englewood Cliffs: Prentice Hall, 1987

[137] Shafiei M, Khaji N. Simulation of two-dimensional elastodynamic problems using a new adaptive physics-based method. Meccanica, 2014, 49(6): 1353–1366

[138] Yousefi H, Noorzad A, Farjoodi J, Vahidi M. Multiresolution-based adaptive simulation of wave equation. Applied Mathematics & Information Sciences, 2012, 6: 47S–58S

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature