Next Article in Journal
Energy Mechanism and Acoustic Emission Characteristics in Rock-Backfill Composite Structure Specimens under Multi-Level Cyclic Loads: Cement-Tailings Ratio Effect
Next Article in Special Issue
What Extra Information Can Be Provided by Multi-Component Seismic Data: A Case Study of 2D3C Prospecting of a Copper–Molybdenum Mine in Inner Mongolia, China
Previous Article in Journal
Inhibiting Mechanism of High pH on Molybdenite Flotation. An Experimental and DFT Study
Previous Article in Special Issue
Real-Time Ambient Seismic Noise Tomography of the Hillside Iron Oxide–Copper–Gold Deposit
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Full-Waveform Modeling of Complex Media Seismic Waves for Irregular Topography and Its Application in Metal Ore Exploration

1
Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
2
University of Chinese Academy of Sciences, Beijing 100049, China
*
Author to whom correspondence should be addressed.
Minerals 2024, 14(7), 664; https://doi.org/10.3390/min14070664
Submission received: 27 May 2024 / Revised: 16 June 2024 / Accepted: 17 June 2024 / Published: 27 June 2024
(This article belongs to the Special Issue Seismics in Mineral Exploration)

Abstract

:
Seismic exploration has caught widespread attention in metal ore exploration due to its higher resolution. However, the presence of topography and complex underground structures in metal ore exploration complicates seismic records. Therefore, it is essential to apply a numerical simulation method suitable for metal ore exploration to study the propagation law of seismic waves in shallow and ore-forming zones, providing reliable theoretical support for multi-component seismic techniques. In particular, the presence of topography generates strong-amplitude surface waves, scattered waves, and converted waves, which consistently distort seismic records and affect the imaging accuracy of the metallogenic belts. Additionally, the propagation of seismic waves is also affected by the anisotropy and viscoelasticity of the underground medium. This paper proposes an elastic wave finite-difference numerical simulation method suitable for irregularly topographical and complex medium conditions, named the comprehensive parameter correction method, which implements a free-surface boundary condition based on the concept of medium averaging. It is algorithmically simple and implies no additional computational costs. Meanwhile, the results obtained by this method are highly consistent with those of the spectral element method, demonstrating its accuracy. By presenting several numerical simulation cases and illustrating the impact of topography and medium conditions on seismic records, this paper demonstrates the necessity of considering irregularly topographical and complex medium conditions in metal ore exploration. In conclusion, the numerical simulation method we propose provides a solid theoretical foundation for the application of seismic exploration methods in metal ore exploration.

1. Introduction

With the development of industrialization, global demand for metal mineral resources is rising. Traditional metal ore exploration techniques, such as gravitational, magnetic, and electrical exploration, have advantages in shallow prospecting (<500 m), but are limited in deeper exploration areas (0.5 km–1.5 km) [1,2]. Seismic exploration, with its great depth and high resolution, has attracted significant attention in metal ore exploration. However, complex underground structures and topography can distort seismic records. In this context, it is necessary to implement numerical simulations of seismic waves to aid in understanding wave propagation phenomena and geophysical characteristics of ore deposits. Consequently, in recent years, many scholars have conducted seismic numerical simulations on mineral resource exploration [3,4,5,6,7,8].
The spectral element method (SEM) [9,10] and the finite difference method (FDM) [11,12] are two commonly used numerical simulation methods in seismic exploration. SEM faces challenges such as high computational cost and difficulty in designing high-quality unstructured grids. Compared with the SEM, the FDM offers the advantages of higher computational efficiency, lower memory requirements, and easier implementation. Among the FDMs, the staggered-grid FDM has been widely used in seismic exploration [12,13,14]. Therefore, in this study, we will employ the staggered-grid FDM to conduct numerical simulations.
In metal mineral exploration, complex topographical conditions are often encountered, producing a variety of complex seismic wave propagation phenomena such as scattered waves, site effects, and mode conversions between surface and body waves [15]. These phenomena distort seismic records and subsequently affect the imaging accuracy of the metallogenic belt. Therefore, to more accurately simulate the propagation of seismic waves in ore-forming zones, it is essential to consider the influence of complex topography.
In the presence of irregular topography, implementing a free-surface boundary condition (FSBC) is crucial for accurate numerical simulation. Within the framework of the FDM, the previous FSBC was represented by applying a difference approximation of displacement derivatives [11,16,17,18,19,20]. With the development of the staggered-grid FDM, the implementation of an FSBC shifts to the velocity–stress staggered grid. When dealing with complex topography, two common strategies in the FDM are staircase approximation [21,22,23] and specific coordinate transformation techniques [24,25,26,27]. The latter implements an FSBC in a curvilinear coordinate system instead of the traditional Cartesian coordinate system. In contrast, the staircase approximation is easier to implement because it discretizes the irregular topography in a conventional Cartesian coordinate system, making it suitable for any topography. The staircase discretization methods include the improved vacuum formulation [28,29,30], generalized stress-image method [14,31,32], and parameter-modification method [22,33,34]. Compared with the other two methods, the parameter-modification method has higher simulation accuracy and greater convenience, showing promising prospects particularly in near-surface numerical simulation [23,35].
In addition to the impact of complex topography, complex subsurface medium conditions also pose significant challenges in metal ore seismic exploration. Traditional metal ore seismic exploration often assumes isotropy. However, anisotropy and viscoelasticity are more common in geological formations. Considering anisotropy and viscoelasticity allows for a more accurate description of the physical phenomena of seismic waves. Anisotropy causes seismic waves to propagate at different speeds in different directions [36,37]. For viscoelastic media, they cause the attenuation of seismic waves, which is specifically characterized by the decay of their amplitude and phase [38]. To obtain more accurate numerical results, it is necessary to consider the influence of complex medium properties during numerical simulation.
In this paper, we address the challenges in metal ore exploration, including irregular topography and complex media, by proposing a comprehensive parameter correction method (CPCM). The CPCM aims to analyze seismic wave propagation in shallow layers and ore-forming zones, offering theoretical support for subsequent seismic imaging and inversion. Through staircase approximation and the correction of constitutive relations and density at the free surface, the CPCM efficiently deals with the numerical simulation of irregular topography. For scenarios involving anisotropic, viscoelastic, or anisotropic–viscoelastic media, we derived the corresponding mathematical expressions of an FSBC using the CPCM. Subsequently, we will present numerical experiments to demonstrate the feasibility of the CPCM, and compare its simulation results with those of the SEM to validate its accuracy. Additionally, we will show the impact of topography and complex media on the seismic wavefield, confirming the importance of considering complex topography and media in metal ore seismic exploration. Finally, we extended our investigation to the Half Mile Lake deposit model [4], testing the performance of the CPCM in addressing practical problems.

2. Method of Free-Surface Parameter Correction for Different Media

In this section, we will derive the FSBC for various media within the staggered-grid FD scheme. The derivation process primarily focuses on elastic medium conditions, with subsequent derivations for anisotropic, viscoelastic, and anisotropic–viscoelastic media, being essentially similar to elastic media. Therefore, we will only provide corresponding notes and final derivation results.

2.1. Elastic Scenario

First, we modified the constitutive relations for the free surface in elastic media. Here, we only considered a 2D case, but a 3D case can be derived similarly. The constitutive relations for isotropic elastic media in the 2D case are the following:
σ x x σ z z σ x z = λ + 2 μ λ 0 λ λ + 2 μ 0 0 0 2 μ ε x x ε z z ε x z ,
where (σxx, σxz, σzz) represent the stress tensor, and (εxx, εxz, εzz) denote the strain tensor. λ and μ are the Lamé parameters. Then, following the averaged medium approach proposed by Moczo et al. (2002) [39], we assumed two completely elastic isotropic half-spaces separated at z = 0. Based on the traction and displacement continuity conditions at the discontinuous interface z = 0, stress tensor and strain tensor components were divided into two categories.
The first category includes continuous stress and strain tensor components, as the following:
σ C = σ z z , σ x z T , ε C = ε x x .
The second category includes discontinuous stress and strain tensor components, as the following:
σ D = σ x x , ε D = ε z z , ε x z T ,
where the subscript “C” denotes the continuous components, and “D” denotes the discontinuous components. Substituting Equations (2) and (3) into Equation (1), we then obtained the following simplified constitutive relationship:
σ D σ C = K L L T M ε C ε D ,
where the superscript “T” denotes the transpose of the matrix. K, L, and M are matrices of elastic parameters, with their detailed components listed as the following:
K = λ + 2 μ , L = λ 0 , M = λ + 2 μ 0 0 2 μ .
Then, rewriting Equation (4) and relocating the discontinuous variables of stress and strain tensor components to the left side, we derived the following equation:
σ D ε D = L M 1 K L M 1 L T M 1 M 1 L T σ C ε C .
We denote the parameters in the upper medium with the superscript “u” and those in the lower medium with the superscript “d”. Applying the constitutive relationship from Equation (6) to the upper and lower half-spaces, we obtained the following formula:
σ D i ε D i = L i ( M i ) 1 K i L i ( M i ) 1 ( L i ) T ( M i ) 1 ( M i ) 1 ( L i ) T σ C i ε C i , i u , d
where
K i = λ i + 2 μ i , L i = λ i 0 , M i = λ i + 2 μ i 0 0 2 μ i , i u , d .
For the continuous stress and strain tensor components, they are continuous at the interface z = 0. Accordingly, the average values of continuous stress and strain tensor components at the interface are
σ ¯ C ε ¯ C = σ C u ε C u = σ C d ε C d .
Next, we denote Q as the discontinuous variable across the interface, i.e., QuQd. We define the averaging function as
Q ¯ = 1 2 Q u + Q d ,
where the superscript “–” denotes the average value of the medium parameters at the interface. So, according to Equation (10), for the discontinuous stress and strain, the average at the interface can be expressed as
σ ¯ D ε ¯ D = 1 2 σ D u ε D u + σ D d ε D d = 1 2 L u ( M u ) 1 K u L u ( M u ) 1 ( L u ) T ( M u ) 1 ( M u ) 1 ( L u ) T + L d ( M d ) 1 K d L d ( M d ) 1 ( L d ) T ( M d ) 1 ( M d ) 1 ( L d ) T σ ¯ C ε ¯ C = ( L M 1 ) ¯ ( K L M 1 L T ) ¯ ( M 1 ) ¯ ( M 1 L T ) ¯ σ ¯ C ε ¯ C .
Then, we rewrote Equation (11) following the format of Equation (4), resulting in the following constitutive relationship for the averaged medium at the interface:
σ ¯ D σ ¯ C = K * L * L * T M * ε ¯ C ε ¯ D ,
where
K * = ( K L M 1 L T ) ¯ + ( L M 1 ) ¯ ( M 1 ) ¯ 1 ( M 1 L T ) ¯ , L * = ( L M 1 ) ¯ ( M 1 ) ¯ 1 , M * = ( M 1 ) ¯ 1 .
Substituting Equation (8) into Equation (13), we obtained the averaged elastic parameter matrix as
K * = 2 μ d λ d + μ d λ d + 2 μ d + 2 μ u λ u + μ u λ + + 2 μ + + χ , L * = λ d λ u + 2 μ u + λ u λ d + 2 μ d λ d + 2 μ d + λ u + 2 μ u 0 , M * = 2 λ d + 2 μ d λ u + 2 μ u λ d + 2 μ d + λ u + 2 μ u 0 0 4 μ d μ u μ d + μ u ,
where
χ = λ d λ u + 2 μ u + λ u λ d + 2 μ d 2 2 λ d + 2 μ d λ u + 2 μ u λ d + 2 μ d + λ u + 2 μ u .
Thus, we can express the constitutive relationship for the averaged medium as
σ ¯ x x σ ¯ z z σ ¯ x z = 2 μ d λ d + μ d λ d + 2 μ d + 2 μ u λ u + μ u λ + + 2 μ + + χ λ d λ u + 2 μ u + λ u λ d + 2 μ d λ d + 2 μ d + λ u + 2 μ u 0 λ d λ u + 2 μ u + λ u λ d + 2 μ d λ d + 2 μ d + λ u + 2 μ u 2 λ d + 2 μ d λ u + 2 μ u λ d + 2 μ d + λ u + 2 μ u 0 0 0 4 μ d μ u μ d + μ u ε ¯ x x ε ¯ z z ε ¯ x z .
Next, considering the upper medium as air and the lower medium as the Earth, the interface (z = 0) transforms into a free-surface boundary. Setting μu = 0 and substituting it into Equations (14)–(16), we then obtained the following corrected interface constitutive relationship:
σ ¯ x x σ ¯ z z σ ¯ x z = 2 μ d λ d + μ d λ d + 2 μ d + χ λ d λ u + λ u λ d + 2 μ d λ d + 2 μ d + λ u 0 λ d λ u + λ u λ d + 2 μ d λ d + 2 μ d + λ u 2 λ u λ d + 2 μ d λ d + 2 μ d + λ u 0 0 0 0 ε ¯ x x ε ¯ z z ε ¯ x z ,
where
χ = λ u λ d + 2 μ d + λ d λ u 2 2 λ u λ d + 2 μ d λ u + λ d + 2 μ d .
If λu = 0, Equation (18) becomes meaningless. Therefore, based on the vacuum approximation theory and mathematical limit concept, we let λu approach an infinitesimal value but not equal to zero. This approximation implies that the elastic modulus of the atmosphere above the surface is much smaller than that of the Earth below, in terms of the physical properties at the air–Earth interface. Therefore, by substituting λu → 0 into Equations (17) and (18), we obtained the following:
σ ¯ x x σ ¯ z z σ ¯ x z = 2 μ d λ d + μ d λ d + 2 μ d 0 0 0 0 0 0 0 0 ε ¯ x x ε ¯ z z ε ¯ x z .
Removing all superscripts, we then obtained the following corrected constitutive relationship at the free-surface boundary:
σ x x σ z z σ x z = 2 μ λ + μ λ + 2 μ 0 0 0 0 0 0 0 0 ε x x ε z z ε x z .
For the calculation of density ρ at the interface, we used the arithmetic mean of the densities of the upper and lower media. Since the upper half-space is air, the density is 0 based on the vacuum approximation, and the average density at the interface is ρ/2. Combining this with Equation (20), we obtained the modified constitutive relations for the stress tensor and density at the free surface:
σ x x z = 0 = 2 μ ( λ + μ ) λ + 2 μ ε x x , σ z z z = 0 = 0 , σ x z z = 0 = σ z x z = 0 = 0 , ρ z = 0 = ρ 2 .
In practical exploration environments, the Earth’s surface is not perfectly flat and horizontal. Therefore, after addressing numerical simulation issues related to a horizontal free surface, it became essential to further consider the influence of topography.
In the process of simulating the propagation of subsurface elastic waves using the standard staggered-grid FDM, since the discrete grids were all standard rectangular grids, it was necessary to approximate the irregular topography of the surface through staircase approximation. To address this, we adopted the irregular surface parameter-modified method proposed by Cao & Chen (2018) [21] for isotropic media. In the case of isotropic media, during grid partitioning, as illustrated in Figure 1a, we assigned components of the normal stress tensor and medium parameters (λ, μ, and ρ) to integer grid points (i, j). Shear stress components were allocated to non-integer grid points at the center of the grid cells (i + 1/2, j + 1/2). The x-direction velocity component (vx) and the z-direction velocity component (vz) were assigned to the remaining horizontal and vertical half-grid points (i, j + 1/2) and (i + 1/2, j). For medium parameters not located at integer grid points, we employed the parameter-averaging method proposed by Moczo et al. (1997) [39] for the calculation:
ρ x = ρ i , j + ρ i + 1 , j 2 , ρ z = ρ i , j + ρ i , j + 1 2 , < μ > x z = 1 4 1 μ i , j + 1 μ i + 1 , j + 1 μ i , j + 1 + 1 μ i + 1 , j + 1 1 .
We set up three types of grid cells, labeled as H, VL, and VR, to represent different grid types under three distinct topographical conditions. In the H cell, it is a horizontal boundary grid cell (Figure 1a). The VL cell represents a vertical boundary grid cell with air on the left side (Figure 1b), while the VR cell represents a vertical boundary grid cell with air on the right side (Figure 1c). The constitutive relationships at the free surface for these cells need to be reconstructed through coordinate axis rotations. For the VL grid cell, the free surface is equivalent to rotating the H grid cell counterclockwise by 90°, and for the VR grid cell, it is equivalent to rotating the H grid cell clockwise by 90°. The modified constitutive relationships and density at the free surface for these three types of grids are the following:
H : σ x x = 2 μ ( λ + μ ) λ + 2 μ ε x x σ z z = 0 σ x z = 2 μ ε x z ρ x = ρ x 2 ρ z = ρ z , VL : σ x x = 0 σ z z = 2 μ ( λ + μ ) λ + 2 μ ε z z σ x z = 2 μ ε x z ρ x = ρ x ρ z = ρ z 2 , VR : σ x x = 0 σ z z = 2 μ ( λ + μ ) λ + 2 μ ε z z σ x z = 0 ρ x = 0 ρ z = ρ z 2 .
Subsequently, in the grid cells within an irregular topography, as illustrated in Figure 1d, we categorized the topography transition cells into inner corner grid cells (IL, IR) and outer corner grid cells (OL, OR), each subject to specific treatment. Both required density averaging. For the constitutive relationships at these grid cells, we simplified the inner corner grid cells by treating them as internal grid points without specific adjustments. Regarding the outer corner grid cells, we set all normal stresses to zero, a method similar to the approach used by Hayashi et al. (2001) [40].
To implement the above method in the staggered-grid FDM program, it was necessary to redefine the medium parameters. We needed to transform the parameter treatment of special grid cells under irregular topography into corrections for these newly defined medium parameters. The FDM format for the staggered grid with the redefined medium parameters is the following:
v x t = b x * ( σ x x + σ x z ) , v z t = b z * ( σ z x + σ z z ) ,
and
σ x x t = E 1 x x v x x + E 2 x x v z z , σ z z t = E 1 z z v x x + E 2 z z v z z , σ x z t = < μ * > x z v x z + v z x .
The redefined medium parameters include the modified reciprocal densities ( b x * and b z * ) and parameters ( E 1 x x , E 2 x x , E 1 z z , E 2 x x , < μ * > x z ) associated with the medium’s Lamé parameters, specifically represented as the following:
b x = 2 b i , j b i + 1 , j b i , j + b i + 1 , j , b z = 2 b i , j b i , j + 1 b i , j + b i , j + 1 , b x * = b x , b z * = b z , E 1 x x = E 2 z z = λ + 2 μ , E 2 x x = E 1 z z = λ , < μ * > x z = < μ > x z ,
where bx and bz are the harmonic mean values of the reciprocal density b in the x- and z- directions for the upper and lower media. For grid points on the complex surface condition, we implemented this by modifying these newly defined medium parameters. For example, the aforementioned VR boundary was achieved through the following parameter correction:
b x * = b x , b z * = 2 b z , E 2 z z = 2 μ ( λ + μ ) λ + 2 μ , E 1 x x = E 2 x x = E 1 z z = 0 , < μ * > x z = < μ > x z .
There was no need to deliberately set bx and <μ>xz to zero because the adjacent points of this grid cell had parameters with a value of zero. According to the harmonic mean, automatic zeroing could be achieved. Therefore, there was no need to distinguish the relative positional relationship between the grid cells in the air and on the ground during the simulation calculation, simplifying the computational process. Similarly, the treatment at the transition cells of irregular topography was achieved by modifying the medium parameters. For example, external corner points (OL or OR) were implemented through the following parameter correction:
b x * = b x , b z * = 2 b z , E 2 z z = E 1 x x = E 2 x x = E 1 z z = 0 , < μ * > x z = < μ > x z 2 .
Therefore, we achieved a unified computation for the complex free surface and the internal region of the model. This method exhibits its good applicability to different Poisson’s ratios of the medium and is suitable for both frequency and time domains [22]. It can be directly embedded into existing FDM programs without increasing computational costs or memory requirements. Subsequently, this method for handling irregular topography can be further extended to anisotropic media, viscoelastic media, and anisotropic–viscoelastic media based on the same derivation principles.

2.2. Anisotropic Scenario

For anisotropic media, vertically transverse isotropic (VTI) media [41] and horizontally transverse isotropic (HTI) media [42,43] are the most commonly used types. Note that if we did not specify the medium conditions, they were assumed to be isotropic or elastic by default. For example, anisotropic media were considered to be anisotropic elastic media unless otherwise stated. Taking VTI media as an example (with a similar derivation for HTI media), their constitutive relationship is the following:
σ x x σ z z σ x z = c 11 c 13 0 c 13 c 33 0 0 0 c 44 ε x x ε z z ε x z .
Adopting the same derivation process as isotropic elastic, we obtained the modified constitutive relationship for VTI media at the free surface as the following:
σ x x σ z z σ x z = c 11 c 33 c 13 2 2 c 33 0 0 0 0 0 0 0 0 ε x x ε z z ε z x .
For parameter correction in irregular topography, we employed the same approach as outlined in Section 2.1. The constitutive relations and densities at the free surface for the rotated H, VL, and VR grid cells are the following:
H : σ x x = c 11 c 33 c 13 2 2 c 33 ε x x σ z z = 0 σ x z = 2 < c 44 > x z ε x z ρ x = ρ x 2 ρ z = ρ z , VL : σ x x = 0 σ z z = c 11 c 33 c 13 2 2 c 33 ε z z σ x z = 2 < c 44 > x z ε x z ρ x = ρ x ρ z = ρ z 2 , VR : σ x x = 0 σ z z = c 11 c 33 c 13 2 2 c 33 ε z z σ x z = 0 ρ x = 0 ρ z = ρ z 2 .
It is important to note here that we used harmonic averaging for the anisotropic elastic coefficients, which can be expressed as
< c 44 > x z = 1 4 1 c 44 i , j + 1 c 44 i + 1 , j + 1 c 44 i , j + 1 + 1 c 44 i + 1 , j + 1 1 .
The handling of the inner and outer corner grid cells was similar to the approach described in Section 2.1 and will not be further elaborated here.

2.3. Viscoelastic Scenario

For viscoelastic media, the derivation of free-surface parameter correction was slightly different. Firstly, based on the research by Carcione et al. (1988) [44], the constitutive relationship for viscoelastic media is presented as the following:
σ x x σ z z σ x z = λ + 2 μ λ 0 λ + μ 2 μ 0 λ λ + 2 μ 0 λ + μ 2 μ 0 0 0 2 μ 0 0 2 μ ε x x ε z z ε x z l = 1 L 1 e 1 l ( 1 ) l = 1 L 2 e 11 l ( 2 ) l = 1 L 2 e 13 l ( 2 ) ,
and the memory variable equations which govern seismic wave energy attenuation are
e 1 l ( 1 ) t = φ 1 l ( ε x x + ε z z ) 1 τ σ l ( 1 ) e 1 l ( 1 ) , e 11 l ( 2 ) t = 0.5 φ 2 l ε x x ε z z 1 τ σ l ( 2 ) e 11 l ( 2 ) , e 13 l ( 2 ) t = φ 2 l ε x z 1 τ σ l ( 2 ) e 13 l ( 2 ) ,
with
φ ν l = 1 τ σ l ν 1 τ ε l ( ν ) τ σ l ( ν ) l = 1 L ν τ ε l ( ν ) τ σ l ( ν ) 1 ,   ν = 1   or   2 ,
where ν = 1 or 2 represents dilatational or shear mode, respectively; Lν is the relaxation mechanism parameter; τ ε l ( ν ) and τ σ l ( ν ) are the strain and stress relaxation times, respectively; and e l ( ν ) , e 11 l ( ν ) , and e 13 l ( ν ) are viscoelastic memory variables. Therefore, in the derivation process of free surface parameter correction for viscoelastic media, there was a slight modification compared to the previous derivation. We set the viscoelastic memory variable matrix as
W 1 = l = 1 L 1 e 1 l ( 1 ) , W 2 = l = 1 L 2 e 11 l ( 2 ) l = 1 L 2 e 13 l ( 2 ) .
Hence, the constitutive relationship for Equation (25) can be expressed as
σ D σ C = K L S 1 S 2 L T M S 3 S 4 ε C ε D W 1 W 2 ,
Matrices K, L, and M are the same as in Equation (5), and matrix S can be represented as
S 1 = λ + μ , S 2 = 2 μ 0 , S 3 = λ + μ 0 , S 4 = 2 μ 0 0 2 μ .
The subsequent derivation process was similar, and eventually, we obtained the constitutive relationship correction for the free surface in viscoelastic media:
σ x x σ z z σ x z = 2 μ ( λ + μ ) λ + 2 μ 0 0 μ ( λ + μ ) λ + 2 μ 2 μ ( λ + μ ) λ + 2 μ 0 0 0 0 0 0 0 0 0 0 0 0 0 ε x x ε z z ε x z l = 1 L 1 e 1 l ( 1 ) l = 1 L 2 e 11 l ( 2 ) l = 1 L 2 e 13 l ( 2 ) .
For viscoelastic parameter correction in irregular terrain, we adopted the same approach as outlined in Section 2.1. The constitutive relations and densities at the free surface for the rotated H, VL, and VR grid cells are the following:
H : σ x x = 2 γ ε x x + γ l = 1 L 1 e 1 l ( 1 ) + 2 γ l = 1 L 2 e 11 l ( 2 ) σ z z = 0 σ x z = 2 μ ε x z + 2 μ l = 1 L 2 e 13 l ( 2 ) ρ x = ρ x 2 ρ z = ρ z , VL : σ x x = 0 σ z z = 2 γ ε z z + γ l = 1 L 1 e 1 l ( 1 ) 2 γ l = 1 L 2 e 11 l ( 2 ) σ x z = 2 μ ε x z + 2 μ l = 1 L 2 e 13 l ( 2 ) ρ x = ρ x ρ z = ρ z 2 , VR : σ x x = 0 σ z z = 2 γ ε z z + γ l = 1 L 1 e 1 l ( 1 ) 2 γ l = 1 L 2 e 11 l ( 2 ) σ x z = 0 ρ x = 0 ρ z = ρ z 2 .
Here, we defined μ (λ + μ)/(λ + 2μ) as γ, and the treatment of inner and outer corner grid cells was similar to the previous method in Section 2.1.

2.4. Anisotropic–Viscoelastic Scenario

For the correction derivation of the FSBC in anisotropic–viscoelastic media, the process was similar to that of viscoelastic media. First, based on the constitutive model and wave equation proposed by Carcione (1995) [45] for anisotropic–viscoelastic media, the constitutive relationship for anisotropic–viscoelastic media is given by
σ x x σ z z σ x z = c 11 c 13 0 ξ 2 c 44 0 c 13 c 33 0 ξ 2 c 44 0 0 0 c 44 0 0 c 44 ε x x ε z z ε x z l = 1 L 1 e 1 l ( 1 ) l = 1 L 2 e 11 l ( 2 ) l = 1 L 2 e 13 l ( 2 ) ,
where
ξ = c 11 + c 33 2 c 44 .
Finally, we obtained the corrected constitutive relationship for the free surface in anisotropic viscoelastic media:
σ x x σ z z σ z x = c 11 c 13 c 13 2 2 c 33 0 0 ξ 2 1 c 13 c 33 c 44 2 1 + c 13 c 33 0 0 0 0 0 0 0 0 0 0 0 0 0 ε x x ε z z ε z x l = 1 L 1 e 1 l ( 1 ) l = 1 L 2 e 11 l ( 2 ) l = 1 L 2 e 13 l ( 2 ) .
In irregular topography, we similarly derived the constitutive relation corrections for the H, VL, and VR points at the free surface as the following:
H : σ x x = p ε x x + q l = 1 L 1 e 1 l ( 1 ) + s l = 1 L 2 e 11 l ( 2 ) σ z z = 0 σ x z = c 44 ε x z + c 44 l = 1 L 2 e 13 l ( 2 ) ρ x = ρ x 2 ρ z = ρ z , VL : σ x x = 0 σ z z = p ε z z + q l = 1 L 1 e 1 l ( 1 ) s l = 1 L 2 e 11 l ( 2 ) σ x z = c 44 ε x z + c 44 l = 1 L 2 e 13 l ( 2 ) ρ x = ρ x ρ z = ρ z 2 , VR : σ x x = 0 σ z z = p ε z z + q l = 1 L 1 e 1 l ( 1 ) s l = 1 L 2 e 11 l ( 2 ) σ x z = 0 ρ x = 0 ρ z = ρ z 2 ,
where
p = c 11 c 13 c 13 2 2 c 33 , q = ξ 2 1 c 13 c 33 , s = c 44 1 + c 13 c 33 .
The treatment of inner and outer corner points was similar to that discussed earlier, so we will not elaborate further here.
Thus, we successfully derived the method for optimizing surface media parameters under irregularly topographical and different medium conditions. This method integrates both topographical and medium considerations, which we refer to as the comprehensive parameter correction method. Using this method, it becomes possible to comprehensively simulate seismic wave propagation under various geological conditions, providing a more accurate and reliable tool for complex media metal ore exploration in irregular topography. In the following sections, we will further apply the CPCM to numerical simulations of metal ore based on some case studies.

3. Numerical Experiments of a Simple Model

In this section, we will explore the effects of topographical and medium conditions on numerical simulations of metal ore exploration by applying the CPCM to two simple metal ore models. Simultaneously, we will compare the simulation results with those obtained from the SEM to demonstrate the accuracy of this approach.

3.1. The Influence of Topography

We designed a simple metal ore model, denoted as model A, with a flat surface. The dimensions of the model are 1500 m × 1000 m. Within the depth range of 600 m to 700 m underground, there exists a large rectangular metal ore deposit (depicted by the red anomalous region in Figure 2a). The velocity and density parameters for the metal ore and their wall rock were independently set based on the parameter ranges summarized by Malehmir et al. (2012) [46]. The specific velocity and density parameters are shown in Table 1.
Subsequently, to further explore the influence of irregular topography, we applied a sinusoidal undulating surface to model A, forming model B, as illustrated in Figure 2b. The velocity and density parameters remained the same as those of model A. Through this model adjustment, we could more clearly assess the impact of topography on numerical simulations of metal ore exploration.
The source was positioned at a distance of 250 m on the model surface (indicated by the red pentagram in Figure 2a), serving as a single vertical force source. The source’s time function adopted a Ricker wavelet with a dominant frequency of 60 Hz, and the grid spacing was set to 1 m. The time sampling interval was 6 × 10–5 s, and the total simulation time was 1 s. In addition to the free surface, the other three boundaries of the model used Convolutional Perfectly Matched Layer (CPML) [47] to absorb spurious reflected waves. We first set the surface of model A separately as either an absorbing boundary or a free-surface boundary for numerical simulation. A 6th-order staggered-grid FDM was employed for the simulations. Figure 3a,b show the synthetic shot gathers (vx component) under these two different surface conditions. In the shot gathers of both conditions, a very clear direct P wave can be observed, along with P–P, P–S, and S–S waves generated by reflections from the ore body, which we marked in the figures. Here, the P–P wave represents the reflected P wave generated by the direct P wave in the ore body, the P–S wave indicates the converted S wave generated by the direct P wave in the ore body, and the S–S wave is the reflected S wave by the direct S wave. These three types of reflected waves are prominent, while the S–P wave (P wave reflected from the ore body by the direct S wave) is less observable due to its small amplitude and overlapping with other waveforms. Additionally, reflections from the right boundary (RRB) of the rectangular metal ore model were observed in both scenarios due to the incident angle. However, compared to the absorbing boundary, as shown in Figure 3b, the free surface introduced surface waves and multiple reflections (MRs) from the surface and ore body, which may have potentially interfered with subsequent data interpretation. Similarly, comparable waveform recordings were observed in the shot gathers of the vz component, as depicted in Figure 4a,b. In the corresponding wavefield snapshots, more pronounced wavefield differences were also observed (Figure 5a,b).
To further explore the influence of topography on numerical simulations of metal ore, we applied model B for numerical simulation. The source was placed at the same location as in model A, and the receivers were positioned along the surface, with all simulation parameters remaining consistent. Synthetic shot gathers of the horizontal and vertical particle velocity components (vx and vz) under irregular topography are shown in Figure 3c and Figure 4c, respectively. In these figures, we observed corresponding P–P, P–S, and S–S reflected waves, noting significant differences from the records observed with a horizontal surface. These differences arose from the complex surface reflections and the data acquisition on irregular topography. Figure 5c shows wavefield snapshots of the metal ore model B, indicating that the wavefield became more complex under complex topographical conditions, further distorting the seismic records. Therefore, in numerical simulations of seismic exploration for metal ore, careful consideration of topographical conditions is crucial for enhancing the reliability of the simulation results and the accuracy of geological interpretation.

3.2. The Influence of Complex Medium Conditions

In this section, we will discuss the influence of complex medium conditions on numerical simulations of metal ore exploration. For this purpose, we adjusted the full space medium conditions in model B into three different types, anisotropic media, viscoelastic media, and anisotropic–viscoelastic media, with detailed parameters listed in Table 2. The parameters QKappa and Qμ represent the bulk quality factor and the shear quality factor, respectively. Through simulations, we obtain synthetic shot gathers under these three different medium conditions, as shown in Figure 6 and Figure 7. We have marked the waveform characteristics with black arrows, where qP represents the propagation of longitudinal wave under anisotropic conditions, and qS denotes the propagation of shear wave in an anisotropic medium. Therefore, qP–P, qP–S, qS–P, and qS–S represent the reflection waves of the metal ore under anisotropic conditions, respectively. Comparing these records with those under elastic medium conditions in Figure 3c and 4c, we notice a significant amplitude attenuation and energy weakening of metal ore reflection waves in viscoelastic media. Introducing anisotropic media (including anisotropic–viscoelastic media) leads to notable differences in recorded reflection waves compared to elastic media, potentially causing errors in subsequent data processing and interpretation. Figure 8 depicts wavefield snapshots under the three media conditions, where noticeable energy attenuation was observed in viscoelastic media compared to elastic media, while anisotropic media directly altered wavefield propagation compared to isotropic media. Therefore, different medium conditions have a significant impact on seismic wave propagation, potentially leading to inaccuracies in subsequent seismic data processing and interpretation. So, considering complex medium factors is necessary for accurately reflecting the underground wavefield characteristics in seismic exploration for metal ores.

3.3. Accuracy Analysis

To verify the accuracy of the CPCM in numerical simulations, we used the simulation results of model B for single trace analysis. We extracted the single trace waveforms received at the surface at offsets of 250 m, 650 m, and 1050 m, labeled as R1, R2, and R3, respectively (Figure 2b). We used the SPECFEM2D software package [9,48,49,50,51] as a reference for the CPCM. The SPECFEM2D software package employs the SEM for the numerical simulation of elastic waves. Because the SEM naturally satisfies the FSBC, utilizing this software package enables precise numerical simulation results at the Earth’s surface. Comparing our results with those from SPECFEM2D could effectively validate the accuracy and reliability of the CPCM in simulations involving irregular topography.
Figure 9 illustrates the waveform records from receivers R1 to R3 under the corresponding four different medium conditions (elastic, viscoelastic, anisotropic, and anisotropic–viscoelastic). By meticulously comparing the waveforms obtained from the CPCM with those from the SPECFEM2D method, we observed a remarkable consistency between the two. We calculated the L2 misfit norm for the single trace records obtained by the two methods and recorded the error results in Table 3. It is noteworthy that the L2 misfit norms of the waveforms recorded by each receiver under both methods were less than 2%. This significant result thoroughly demonstrates the high accuracy of the CPCM in numerical simulations under irregularly topographical and complex medium conditions.

4. Numerical Simulation for the Half Mile Lake Deposit

In this section, we will primarily investigate the application effectiveness of the CPCM in the numerical simulations of complex metal ore models. During the model design process, we referred to the model in Bellefleur et al.’s 2012 paper [4] and generated a new complex model with topography. We then applied anisotropic viscoelastic medium conditions to this model and analyzed the theoretical shot gather data.

4.1. Geological Background and 2D Model

In this study, the reference model we used was the Half Mile Lake deposit model proposed by Bellefleur et al. (2012) [4]. This ore deposit is located in New Brunswick, Canada, and is a volcanic-hosted massive sulfide deposit, primarily composed of a pyrrhotite breccia matrix with laterally continuous layers of pyrite and pyrrhotite [52]. To better illustrate the simulation results for the metal ore area, we selected a portion of the model for analysis. The model spans a length of 3 km and a depth of 1.6 km, as shown in Figure 10. In this model, the strata hosting the ore deposit are inclined approximately 45°, exhibiting three distinct mineralized zones. In the process of modeling surrounding lithologies, we disregarded lithological units with limited thickness or lateral extent, such as minor veins. The model only includes rock types with sufficiently continuous properties to be determined. Based on the summary of lithological units and their physical properties in the article by Bellefleur et al. [4], we labeled the different lithologies in Figure 10 and listed the velocity–density parameters corresponding to each numbered lithology in Table 4.
During the numerical simulation, a single-point vertical force was located at (x = 2500 m, z = 72 m), and the source time function was a Ricker wavelet with a dominant frequency of 60 Hz. The grid size was set to 1 m, and the simulation time was 1.5 seconds. Initially, we set the surface as an absorbing boundary and conducted finite-difference (FD) modeling under elastic conditions. Subsequently, considering real topographical and complex medium factors, as the original text did not define the anisotropy and viscoelasticity of the subsurface medium, we assumed the Thomsen parameters of anisotropy to be ε = 0.2 and δ* = –0.2, with attenuation quality factors set to QKappa = 80 and Qμ = 40. By analyzing the differences in numerical simulations between the two scenarios, we demonstrated the feasibility of our proposed method in practical complex model situations.

4.2. FD Modeling Results

Figure 11 shows the synthetic shot gathers with the surface set as an absorbing boundary. Without the interference of surface waves and complex media, the reflection waves from the ore body are quite clear. Under actual topographical conditions and anisotropic–viscoelastic media conditions, we obtained the synthetic shot gather records shown in Figure 12. From the figure, reflections from the ore body can be observed, as indicated by the white arrows. However, due to the prominent Rayleigh waves, the reflections generated by the underground ore body were further overshadowed, making the reflected wave energy less apparent compared to that obtained under absorbing boundary conditions. Additionally, in numerical simulation results under actual topographical conditions, staircase diffractions were observed, which also attenuated some of the reflected wave energy, as indicated by the red arrows in Figure 12. Staircase diffractions are caused by the staircase discretization algorithm, which does not exist in reality. Therefore, in the numerical simulations, we used smaller grid spacings to discretize or adopt appropriate methods to suppress staircase diffractions. This helped minimize the impact of staircase diffractions on seismic records. Further discussion on this topic will be provided in the Discussion section. In conclusion, the CPCM demonstrates effective application in numerical simulations of complex models, providing more realistic synthetic shot gathers and offering better guidance for subsequent data processing and interpretation.

5. Discussion

In the process of numerical simulation of metal ore models with complex terrain conditions, we adopted the staircase discretization approach to handle terrain issues. However, when we used a stepped discretization strategy to represent free surfaces with large grid spacing, staircase diffractions were inevitably generated, disrupting the wavefield and recordings, as indicated by the red arrows in Figure 3c. When conducting numerical simulation tests on the Half Mile Lake deposit, we found that this staircase diffraction phenomenon was obvious, potentially obscuring some of the layer reflections. This is a common problem in the FD modeling of complex topography. Since these staircase diffractions were algorithmically generated, we reduced their occurrence significantly by densifying the grid, i.e., reducing the grid spacing, albeit at the cost of increased computational burden. Another approach is to employ a superposition concept. For instance, Zhou et al. (2023) [23] proposed a probabilistic superposition-based discrete model design method to optimize the staircase diffraction issue in simulating irregularly topographical surface waves, which could be applied to numerical simulations of metal ore exploration in the future. In the future, we will conduct further research on suppressing staircase diffractions.

6. Conclusions

In studies of numerical simulations for the seismic exploration of metal ores, the complexity of topography and subsurface media often poses challenges to accurately simulating wavefields. To address this issue, this paper proposes a comprehensive parameter correction method (CPCM) based on an averaged medium approach and the limit theorem. The CPCM effectively simulates wave propagation characteristics in complex surface and medium conditions. Regarding its simulation accuracy, the CPCM is comparable to the spectral element method. It achieves simulation of irregular topography under different medium conditions by adjusting the constitutive relationships and density of the surface. Another major advantage of the CPCM is its high compatibility with existing FD programs, requiring no significant increase in computational resources to achieve high-precision simulation.
Through numerical experiments, we found significant differences in the simulation process of seismic exploration for metal ores, depending on whether irregularly topographical and complex medium conditions were considered. These differences can lead to substantial errors during the seismic data inversion and interpretation stage. Therefore, this paper emphasizes the importance of considering complex surface and medium conditions in the exploration process of metal ores. In conclusion, the CPCM demonstrates broad prospects for application in the field of metal ore exploration. It can assist in improving the accuracy of locating metal ore deposits and provides robust technical support for research and practice in related fields.

Author Contributions

Conceptualization, W.S. and S.H.; methodology, W.S. and X.Z.; software, W.S. and X.Z.; validation, W.S.; formal analysis, W.S.; investigation, W.S.; resources, S.H.; writing—original draft preparation, W.S.; writing—review and editing, W.S., S.H., and X.Z.; supervision, S.H.; funding acquisition, S.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China (2021YFA0716901, 2022YFB3904601) and the National Natural Science Foundation of China (Grant No. 42174160).

Data Availability Statement

The data presented in this study are available on request from the corresponding author due to internal policy.

Acknowledgments

The authors appreciate the feedback and suggestions from all reviewers. In addition, the authors thank Wang Hao for his valuable feedback on the content, which helped to improve the completeness of this article.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Liu, G.D.; Hao, T.Y. Searching of hidden mineral deposits by geophysical methods. Chin. J. Geophys. 1995, 38, 850–854. (In Chinese) [Google Scholar]
  2. Witherly, K. The evolution of minerals exploration over 60 years and the imperative to explore undercover. Lead. Edge 2012, 31, 292–295. [Google Scholar] [CrossRef]
  3. Bohlen, T.; Müller, C.; Milkereit, B. 5. Elastic Seismic-Wave Scattering from Massive Sulfide Orebodies: On the Role of Composition and Shape. In Hardrock Seismic Exploration; Society of Exploration Geophysicists: Houston, TX, USA, 2003; Volume 70–89. [Google Scholar]
  4. Bellefleur, G.; Malehmir, A.; Müller, C. Elastic finite-difference modeling of volcanic-hosted massive sulfide deposits: A case study from Half Mile Lake, New Brunswick, Canada. Geophysics 2012, 77, Wc25–Wc36. [Google Scholar] [CrossRef]
  5. Dehghannejad, M.; Malehmir, A.; Juhlin, C.; Skyttä, P. 3D constraints and finite-difference modeling of massive sulfide deposits: The Kristineberg seismic lines revisited, northern Sweden. Geophysics 2012, 77, Wc69–Wc79. [Google Scholar] [CrossRef]
  6. Malinowski, M.; Schetselaar, E.; White, D.J. 3D seismic imaging of volcanogenic massive sulfide deposits in the Flin Flon mining camp, Canada: Part 2-Forward modeling. Geophysics 2012, 77, Wc81–Wc93. [Google Scholar] [CrossRef]
  7. Pang, Y.; Yan, L.; Liu, Y.; Tang, L.; Zhu, R.; Liu, G. Seismic Wave Finite-Difference forward Modeling for Orogenic Gold Deposits. Minerals 2022, 12, 1465. [Google Scholar] [CrossRef]
  8. Malehmir, A.; Urosevic, M.; Bellefleur, G.; Juhlin, C.; Milkereit, B. Seismic methods in mineral exploration and mine planning—Introduction. Geophysics 2012, 77, Wc1–Wc2. [Google Scholar] [CrossRef]
  9. Komatitsch, D.; Tromp, J. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophys. J. Int. 1999, 139, 806–822. [Google Scholar] [CrossRef]
  10. Komatitsch, D.; Vilotte, J.-P.; Vai, R.; Castillo-Covarrubias, J.M.; Sánchez-Sesma, F.J. The spectral element method for elastic wave equations—Application to 2-D and 3-D seismic problems. Int. J. Numer. Methods Eng. 1999, 45, 1139–1164. [Google Scholar] [CrossRef]
  11. Kelly, K.R.; Ward, R.W.; Treitel, S.; Alford, R.M. Synthetic Seismograms: A Finite-Difference Approach. Geophysics 1976, 41, 2–27. [Google Scholar] [CrossRef]
  12. Virieux, J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics 1986, 51, 889–901. [Google Scholar] [CrossRef]
  13. Virieux, J. SH-wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics 1984, 49, 1933–1942. [Google Scholar] [CrossRef]
  14. Levander, A.R. Fourth-order finite-difference P-SV seismograms. Geophysics 1988, 53, 1425–1436. [Google Scholar] [CrossRef]
  15. Tarrass, I.; Giraud, L.; Thore, P. New curvilinear scheme for elastic wave propagation in presence of curved topography. Geophys. Prospect. 2011, 59, 889–906. [Google Scholar] [CrossRef]
  16. Alterman, Z.; Karal, F.C., Jr. Propagation of elastic waves in layered media by finite difference methods. Bull. Seismol. Soc. Am. 1968, 58, 367–398. [Google Scholar]
  17. Alterman, Z.S.; Rotenberg, A. Seismic waves in a quarter plane. Bull. Seismol. Soc. Am. 1969, 59, 347–368. [Google Scholar] [CrossRef]
  18. Ilan, A.; Ungar, A.; Alterman, Z. An Improved Representation of Boundary Conditions in Finite Difference Schemes for Seismological Problems. Geophys. J. Int. 1975, 43, 727–745. [Google Scholar] [CrossRef]
  19. Ilan, A.; Loewenthal, D. Instability of Finite Difference Schemes due to Boundary Conditions in Elastic Media*. Geophys. Prospect. 1976, 24, 431–453. [Google Scholar] [CrossRef]
  20. Vidale, J.E.; Clayton, R.W. A Stable Free-Surface Boundary-Condition for Two-Dimensional Elastic Finite-Difference Wave Simulation. Geophysics 1986, 51, 2247–2249. [Google Scholar] [CrossRef]
  21. Cao, J.; Chen, J.B. A parameter-modified method for implementing surface topography in elastic-wave finite-difference modeling. Geophysics 2018, 83, T313–T332. [Google Scholar] [CrossRef]
  22. Cao, J.; Chen, J.B.; Dai, M.X. An adaptive free-surface expression for three-dimensional finite-difference frequency-domain modelling of elastic wave. Geophys. Prospect. 2018, 66, 707–725. [Google Scholar] [CrossRef]
  23. Zhou, X.H.; Huo, S.D.; Wang, H.; Dong, S.L.; Liang, Y.; Cao, J. Model parameter design for modeling surface topography in VTI elastic finite-difference near-surface simulations. Geophysics 2023, 88, C33–C52. [Google Scholar] [CrossRef]
  24. Zhang, W.; Zhang, Z.; Chen, X. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids. Geophys. J. Int. 2012, 190, 358–378. [Google Scholar] [CrossRef]
  25. Zhang, W.; Zhuang, Y.; Chung, E.T. A new spectral finite volume method for elastic wave modelling on unstructured meshes. Geophys. J. Int. 2016, 206, 292–307. [Google Scholar] [CrossRef]
  26. Tessmer, E.; Kosloff, D.; Behle, A. Elastic wave propagation simulation in the presence of surface topography. Geophys. J. Int. 1992, 108, 621–632. [Google Scholar] [CrossRef]
  27. Hestholm, S.; Ruud, B. 3D free-boundary conditions for coordinate-transform finite-difference seismic modelling. Geophys. Prospect. 2002, 50, 463–474. [Google Scholar] [CrossRef]
  28. Zeng, C.; Xia, J.; Miller, R.D.; Tsoflias, G.P. An improved vacuum formulation for 2D finite-difference modeling of Rayleigh waves including surface topography and internal discontinuities. Geophysics 2012, 77, T1–T9. [Google Scholar] [CrossRef]
  29. Zahradník, J.Í.; Moczo, P.; Hron, F.E. Testing four elastic finite-difference schemes for behavior at discontinuities. Bull. Seismol. Soc. Am. 1993, 83, 107–129. [Google Scholar]
  30. Bohlen, T.; Saenger, E.H. Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves. Geophysics 2006, 71, T109–T115. [Google Scholar] [CrossRef]
  31. Robertsson, J.O.A. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography. Geophysics 1996, 61, 1921–1934. [Google Scholar] [CrossRef]
  32. Kristek, J.; Moczo, P.; Archuleta, R.J. Efficient Methods to Simulate Planar Free Surface in the 3D 4th-Order Staggered-Grid Finite-Difference Schemes. Stud. Geophys. Et Geod. 2002, 46, 355–381. [Google Scholar] [CrossRef]
  33. Mittet, R. Free-surface boundary conditions for elastic staggered-grid modeling schemes. Geophysics 2002, 67, 1616–1623. [Google Scholar] [CrossRef]
  34. Xu, Y.X.; Xia, J.H.; Miller, R.D. Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach. Geophysics 2007, 72, Sm147–Sm153. [Google Scholar] [CrossRef]
  35. Dong, S.L.; Zhou, X.H.; Chen, J.B. Finite-difference modeling with topography using 3D viscoelastic parameter-modified free-surface condition. Geophysics 2023, 88, T211–T226. [Google Scholar] [CrossRef]
  36. Fryer, G.J.; Frazer, L.N. Seismic waves in stratified anisotropic media. Geophys. J. Int. 1984, 78, 691–710. [Google Scholar] [CrossRef]
  37. Petrov, I.B.; Golubev, V.I.; Petrukhin, V.Y.; Nikitin, I.S. Simulation of Seismic Waves in Anisotropic Media. Dokl. Math. 2021, 103, 146–150. [Google Scholar] [CrossRef]
  38. Robertsson, J.O.; Blanch, J.O.; Symes, W.W. Viscoelastic finite-difference modeling. Geophysics 1994, 59, 1444–1456. [Google Scholar] [CrossRef]
  39. Moczo, P.; Kristek, J.; Vavryčuk, V.C.; Archuleta, R.J.; Halada, L. 3D Heterogeneous Staggered-Grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities. Bull. Seismol. Soc. Am. 2002, 92, 3042–3066. [Google Scholar] [CrossRef]
  40. Hayashi, K.; Burns, D.R.; Toksöz, M.N. Discontinuous-Grid Finite-Difference Seismic Modeling Including Surface Topography. Bull. Seismol. Soc. Am. 2001, 91, 1750–1764. [Google Scholar] [CrossRef]
  41. Postma, G.W. Wave Propagation in a Stratified Medium. Geophysics 1955, 20, 780–806. [Google Scholar] [CrossRef]
  42. Crampin, S. Seismic-wave propagation through a cracked solid: Polarization as a possible dilatancy diagnostic. Geophys. J. Int. 1978, 53, 467–496. [Google Scholar] [CrossRef]
  43. Crampin, S. Effective anisotropic elastic constants for wave propagation through cracked solids. Geophys. J. Int. 1984, 76, 135–145. [Google Scholar] [CrossRef]
  44. Carcione, J.M.; Kosloff, D.; Kosloff, R. Wave propagation simulation in a linear viscoelastic medium. Geophys. J. Int. 1988, 95, 597–611. [Google Scholar] [CrossRef]
  45. Carcione, J.M. Constitutive Model and Wave-Equations for Linear, Viscoelastic, Anisotropic Media. Geophysics 1995, 60, 537–548. [Google Scholar] [CrossRef]
  46. Malehmir, A.; Durrheim, R.; Bellefleur, G.; Urosevic, M.; Juhlin, C.; White, D.J.; Milkereit, B.; Campbell, G. Seismic methods in mineral exploration and mine planning: A general overview of past and present case histories and a look into the future. Geophysics 2012, 77, Wc173–Wc190. [Google Scholar] [CrossRef]
  47. Komatitsch, D.; Martin, R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics 2007, 72, SM155–SM167. [Google Scholar] [CrossRef]
  48. Komatitsch, D.; Vilotte, J.-P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bull. Seismol. Soc. Am. 1998, 88, 368–392. [Google Scholar] [CrossRef]
  49. Tromp, J.; Komatitsch, D.; Liu, Q.Y. Spectral-element and adjoint methods in seismology. Commun. Comput. Phys. 2008, 3, 1–32. [Google Scholar]
  50. Xie, Z.; Komatitsch, D.; Martin, R.; Matzen, R. Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML. Geophys. J. Int. 2014, 198, 1714–1747. [Google Scholar] [CrossRef]
  51. Komatitsch, D.; Xie, Z.; Bozdağ, E.; Sales de Andrade, E.; Peter, D.; Liu, Q.; Tromp, J. Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion. Geophys. J. Int. 2016, 206, 1467–1478. [Google Scholar] [CrossRef]
  52. Adair, R.N. Stratigraphy, structure, and geochemistry of the Halfmile Lake massive-sulfide deposit, New Brunswick. Explor. Min. Geol. 1992, 1, 151–166. [Google Scholar]
Figure 1. (ac) show the basic fundamental grid cells in 2D topography for staircase discretization. (a) H grid cell; (b) VL grid cell; (c) VR grid cell; (d) Topography implementation method. The black solid arrows indicate the conventional coordinate system, while the red dashed arrows represent the coordinate system after rotation. Black solid-line squares represent the positions of staggered half-grid points, and the purple solid line represents the free surface. IL and IR denote inner corner grid cells; OL and OR denote outer corner grid cells.
Figure 1. (ac) show the basic fundamental grid cells in 2D topography for staircase discretization. (a) H grid cell; (b) VL grid cell; (c) VR grid cell; (d) Topography implementation method. The black solid arrows indicate the conventional coordinate system, while the red dashed arrows represent the coordinate system after rotation. Black solid-line squares represent the positions of staggered half-grid points, and the purple solid line represents the free surface. IL and IR denote inner corner grid cells; OL and OR denote outer corner grid cells.
Minerals 14 00664 g001
Figure 2. Simple metal ore models. (a) Model A. (b) Model B. The red anomalous region represents the metal ore deposit. The red pentagram on the surface indicates the position of the seismic source. The green inverted triangles represent the locations of the surface receivers (R1, R2, and R3) respectively.
Figure 2. Simple metal ore models. (a) Model A. (b) Model B. The red anomalous region represents the metal ore deposit. The red pentagram on the surface indicates the position of the seismic source. The green inverted triangles represent the locations of the surface receivers (R1, R2, and R3) respectively.
Minerals 14 00664 g002
Figure 3. Synthetic shot gathers of the horizontal particle velocity components (vx). (a) Surface with absorbing boundary. (b) Surface with horizontal free-surface boundary. (c) Surface with irregular topography. The red arrows indicate staircase diffraction.
Figure 3. Synthetic shot gathers of the horizontal particle velocity components (vx). (a) Surface with absorbing boundary. (b) Surface with horizontal free-surface boundary. (c) Surface with irregular topography. The red arrows indicate staircase diffraction.
Minerals 14 00664 g003
Figure 4. Synthetic shot gathers of the vertical particle velocity components (vz). (a) Surface with absorbing boundary. (b) Surface with horizontal free-surface boundary. (c) Surface with irregular topography.
Figure 4. Synthetic shot gathers of the vertical particle velocity components (vz). (a) Surface with absorbing boundary. (b) Surface with horizontal free-surface boundary. (c) Surface with irregular topography.
Minerals 14 00664 g004
Figure 5. Snapshots of wavefields at different time points (t = 0.2 s, 0.267 s, and 0.334 s). (a) Snapshot with surface as an absorbing boundary. (b) Snapshot with surface as a horizontal free-surface boundary. (c) Snapshot with surface featuring irregular topography.
Figure 5. Snapshots of wavefields at different time points (t = 0.2 s, 0.267 s, and 0.334 s). (a) Snapshot with surface as an absorbing boundary. (b) Snapshot with surface as a horizontal free-surface boundary. (c) Snapshot with surface featuring irregular topography.
Minerals 14 00664 g005
Figure 6. Synthetic shot gathers of the horizontal particle velocity components (vx). (a) Viscoelastic media. (b) Anisotropic media. (c) Anisotropic–viscoelastic media.
Figure 6. Synthetic shot gathers of the horizontal particle velocity components (vx). (a) Viscoelastic media. (b) Anisotropic media. (c) Anisotropic–viscoelastic media.
Minerals 14 00664 g006
Figure 7. Synthetic shot gathers of the vertical particle velocity components (vz). (a) Viscoelastic media. (b) Anisotropic media. (c) Anisotropic–viscoelastic media.
Figure 7. Synthetic shot gathers of the vertical particle velocity components (vz). (a) Viscoelastic media. (b) Anisotropic media. (c) Anisotropic–viscoelastic media.
Minerals 14 00664 g007
Figure 8. Snapshots of wavefields at different time points (t = 0.2 s, 0.267 s, and 0.334 s). (a) Snapshot with viscoelastic media. (b) Snapshot with anisotropic media. (c) Snapshot with anisotropic–viscoelastic.
Figure 8. Snapshots of wavefields at different time points (t = 0.2 s, 0.267 s, and 0.334 s). (a) Snapshot with viscoelastic media. (b) Snapshot with anisotropic media. (c) Snapshot with anisotropic–viscoelastic.
Minerals 14 00664 g008
Figure 9. Vertical and horizontal components of particle velocity at R1, R2, and R3. (a) Elastic media. (b) Viscoelastic media. (c) Anisotropic media. (d) Anisotropic–viscoelastic media. The left half of the figure shows the horizontal component of the particles (vx), while the right half shows the vertical component (vz). The figure compares the numerical simulation results obtained using the CPCM with those obtained using the spectral element method applied by the SPECFEM2D version 8.1.0 software package.
Figure 9. Vertical and horizontal components of particle velocity at R1, R2, and R3. (a) Elastic media. (b) Viscoelastic media. (c) Anisotropic media. (d) Anisotropic–viscoelastic media. The left half of the figure shows the horizontal component of the particles (vx), while the right half shows the vertical component (vz). The figure compares the numerical simulation results obtained using the CPCM with those obtained using the spectral element method applied by the SPECFEM2D version 8.1.0 software package.
Minerals 14 00664 g009
Figure 10. Composite 2D geological section of the Half Mile Lake deposit model.
Figure 10. Composite 2D geological section of the Half Mile Lake deposit model.
Minerals 14 00664 g010
Figure 11. Synthetic shot gathers of the Half Mile Lake deposit model with the surface as an absorbing boundary. (a) Horizontal particle velocity component (vx). (b) Vertical particle velocity component (vz). “OR” indicates ore body reflection.
Figure 11. Synthetic shot gathers of the Half Mile Lake deposit model with the surface as an absorbing boundary. (a) Horizontal particle velocity component (vx). (b) Vertical particle velocity component (vz). “OR” indicates ore body reflection.
Minerals 14 00664 g011
Figure 12. Synthetic shot gathers of the Half Mile Lake deposit model with real topographical and anisotropic–viscoelastic media conditions. (a) Horizontal particle velocity component (vx). (b) Vertical particle velocity component (vz).
Figure 12. Synthetic shot gathers of the Half Mile Lake deposit model with real topographical and anisotropic–viscoelastic media conditions. (a) Horizontal particle velocity component (vx). (b) Vertical particle velocity component (vz).
Minerals 14 00664 g012
Table 1. Velocity and density parameters of the simple metal ore models.
Table 1. Velocity and density parameters of the simple metal ore models.
LithologyVp (m/s)Vs (m/s)ρ (kg/m3)
Metal ore680045004700
Wall rock560033002820
Table 2. Parameters of the media under different conditions. “ε” and “δ*” represent Thomsen parameters.
Table 2. Parameters of the media under different conditions. “ε” and “δ*” represent Thomsen parameters.
MediumThomsen ParametersQuality Factors
εδ*QKappaQμ
Viscoelastic4020
Anisotropic0.2–0.2
Anisotropic–viscoelastic0.2–0.24020
Table 3. The L2 misfit norms for different receivers under four different medium conditions in Figure 9. The four medium conditions were elastic, viscoelastic, anisotropic, and anisotropic–viscoelastic.
Table 3. The L2 misfit norms for different receivers under four different medium conditions in Figure 9. The four medium conditions were elastic, viscoelastic, anisotropic, and anisotropic–viscoelastic.
MediumReceiverL2 Misfit Norms
vx (%)vz (%)
ElasticityR11.190.094
R21.2791.148
R31.4082.081
ViscoelasticityR10.6430.151
R20.9970.508
R32.0210.893
AnisotropyR11.0420.085
R21.4981.065
R31.0091.393
Anisotropic–viscoelasticityR10.9020.186
R21.0570.545
R33.741.651
Table 4. Lithological units and parameters of velocity and density for Half Mile Lake deposit. These data are referenced from the study by Bellefleur et al. (2012) [4].
Table 4. Lithological units and parameters of velocity and density for Half Mile Lake deposit. These data are referenced from the study by Bellefleur et al. (2012) [4].
LithologyVp (m/s)Vs (m/s)Density (kg/m3)
Rhyolite580033702580
Mafic volcanics580033302690
Sediments582033802680
Q–F porphyry583033502460
Felsic volcanics580033702580
Stringer zone560032502820
Sulfide569038403420
Red stripe583033502460
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Su, W.; Huo, S.; Zhou, X. Full-Waveform Modeling of Complex Media Seismic Waves for Irregular Topography and Its Application in Metal Ore Exploration. Minerals 2024, 14, 664. https://doi.org/10.3390/min14070664

AMA Style

Su W, Huo S, Zhou X. Full-Waveform Modeling of Complex Media Seismic Waves for Irregular Topography and Its Application in Metal Ore Exploration. Minerals. 2024; 14(7):664. https://doi.org/10.3390/min14070664

Chicago/Turabian Style

Su, Wenchao, Shoudong Huo, and Xuhui Zhou. 2024. "Full-Waveform Modeling of Complex Media Seismic Waves for Irregular Topography and Its Application in Metal Ore Exploration" Minerals 14, no. 7: 664. https://doi.org/10.3390/min14070664

APA Style

Su, W., Huo, S., & Zhou, X. (2024). Full-Waveform Modeling of Complex Media Seismic Waves for Irregular Topography and Its Application in Metal Ore Exploration. Minerals, 14(7), 664. https://doi.org/10.3390/min14070664

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop