Seismic Wave Finite-Difference Forward Modeling for Orogenic Gold Deposits

: The demand for deep prospecting has led to an increase in the enthusiasm for seismic techniques in mineral exploration. Reﬂection seismology applications in the base metal industry have achieved success. For orogenic gold deposits, however, their applicable conditions remain to be investigated. This paper simulated seismic wave propagation based on a ﬁnite-difference algorithm with an accuracy of eighth order in space and second order in time to investigate the factors inﬂuencing the reﬂection seismic exploration results. Then, the paper assessed the algorithm’s feasibility for orogenic gold deposits, taking the giant Zaozigou deposit in central China as an example. The forward modeling showed that the petrophysical properties, dimensions, and dip of targets signiﬁcantly affected the seismic exploration results. In the Zaozigou model, shallowly dipping orebodies were well imaged with precise extension and thickness. Steeply dipping orebodies were recognized but their thickness information was lost. Steeply dipping orebodies at depth were not detectable under a surface conﬁguration. These problems could be effectively solved by increasing the array length and using vertical seismic proﬁling methods. For small orebodies, multiwave and multicomponent seismic techniques offered more valuable information in terms of mineral exploration. In conclusion, it was possible to locate orogenic gold deposits using the reﬂection seismology method.


Introduction
Most current mineral deposit explorations and mining all over the world are carried out in shallow areas within several hundred meters under the surface. There have been significant declines in the discovery rate and quality of newfound deposits over the past few decades [1,2]. It has become increasingly difficult to discover giant mineral deposits in the traditional search spaces (depths less than 500 m) due to the depletion of shallow deposits [3,4]. The progressive but persistent move to deeper areas of mineral deposit exploration puts forward new requests for exploration methods and techniques [5,6]. Owing to their high resolution at greater depths, seismic methods are receiving increasing attention [3,7,8]. Seismic exploration techniques have been used in the hydrocarbon resource industry for almost 100 years and have a very mature technical system [9,10]. Among these methods, the reflection seismology method has become one of the most important tools in the hydrocarbon sector and has played a dominant role in the delineation of oil and gas deposits [11][12][13][14][15]. Additionally, a wide range of case studies on the applications of reflection seismology in the base metal and volcanogenic massive sulfide deposits have achieved success with varying degrees, such as in the Dapai deposit, China [7,[16][17][18], the Bathurst Mining Camp, Canada [19], the Laisvall deposit, Sweden [20], and the Archean A two-dimensional model of an ideal homogeneous isotropic linear elastic material was considered in this paper [46]. Derived from Newton's second law, the differential equations of motion describe the relationship between the wavefield displacement and the stress state of an elastic material per unit volume in the small deformation approximation [47]: where ρ is the density of the material; u x and u z are the horizontal and vertical particle displacements, respectively; τ xx , τ zz , τ xz , and τ zx are the components of the stress tensor, and τ zx is equal to τ xz ; and f x and f z are the body force components in different directions per unit mass. According to the generalized Hooke's law [48], the linear relationship between stress components and strain components can be described though constitutive equations: where λ and µ are the Lame parameters, which characterize the elastic properties of the medium, and e x , e z , and γ xz are strain components defined by geometric equations as: The separation of P-waves and S-waves can be realized though the divergence and curl operators of the displacement vector.

Numerical Implementation
The finite-difference time-domain method is often used to give numerical solutions to continuous partial derivative problems [49]. An eighth-order accurate FD operator was adopted to approximate the first derivative in the space domain, and a second-order accuracy approximation was used for the second derivative in the time domain [50]: where ∆x denotes the special mesh spacing and ∆t is the time step; u n i,j represents the displacement at position (i∆x, j∆z) at time n∆t; and a k is the element of the finite-difference coefficient array a = 4 5 , − 1 5 , 4 105 , − 1 280 . A cascade of the absorbing boundary conditions and an exponential-damping attenuation layer was used to remove the effects of artificial boundaries [51]. Considering the calculation requirement brought by the FD algorithm, we applied a parallel computational architecture with a GPU to speed up our elastic wave simulation [52][53][54][55][56]. The forward modeling code was developed based on the work of Weiss and Shragge [57] on Madagascar, an open-source software package for multidimensional seismic data analysis.

Forward Modeling for Factors
Petrophysical properties, dimensions and dip are three primary factors that significantly affect the seismic response from ore deposits in different ways [42,58]. In order to facilitate a more effective seismic survey design and interpretation, this section presents the effects of these factors, mainly taking the conventional P wave exploration method as an example.
Primary attention was given to the backpropagating wavefield in this study, because these events might only be received at surface or in shallow boreholes and could contribute to exploration [43]. The models were simplified by eliminating the interference from the free surface [42,43].

Petrophysical Property
Seismic wave reflection occurs when a wave impedance difference exists [59]. The reflection coefficient (R), calculated from the impedance contrast [60], was used as a direct parameter here. Sulfide veins consisting of pyrite, arsenopyrite chalcopyrite, galena, sphalerite, and stibnite are the primary exploration targets of orogenic gold deposits and yield a wide range of wave velocities and densities, as summarized by Salisbury et al. [58]. In this set of simulations, the orebody was simplified into a sulfide vein with a P-wave velocity of 5.3 km/s and a density of 4.1 g/cm 3 . Gabbro, anorthosite, and quartz-mica schist were used as the host rocks in models A-C. The elastic wave velocities and densities for these models are shown in Table 1. The models were discretized with a grid spacing of 4 m. The source wavelets used to generate the synthetic seismic data were Ricker signals with a Minerals 2022, 12, 1465 4 of 17 dominant frequency of 60 Hz and were sampled every 0.4 ms. The FD modeling results showed that the reflections in model A were hard to recognize due to the low reflection coefficient (Figure 1a,b). Model B and model C showed distinct records of reflected waves with a stepwise increasing reflection coefficient (Figure 1c-f). yield a wide range of wave velocities and densities, as summarized by Salisbury et al. [58]. In this set of simulations, the orebody was simplified into a sulfide vein with a P-wave velocity of 5.3 km/s and a density of 4.1 g/cm 3 . Gabbro, anorthosite, and quartz-mica schist were used as the host rocks in models A-C. The elastic wave velocities and densities for these models are shown in Table 1. The models were discretized with a grid spacing of 4 m. The source wavelets used to generate the synthetic seismic data were Ricker signals with a dominant frequency of 60 Hz and were sampled every 0.4 ms. The FD modeling results showed that the reflections in model A were hard to recognize due to the low reflection coefficient (Figure 1a,b). Model B and model C showed distinct records of reflected waves with a stepwise increasing reflection coefficient (Figure 1c-f).   Table 1.

Dimensions
Dimensions, which refer to the size and extent of a target, can influence the imaging result significantly. Because metamorphism and deformation can cause orebodies to become contorted and discontinuous and because even economic ore deposits can be small relative to seismic wavelengths [35,58,61], the effects of dimensions needed to be further  Table 1.

Dimensions
Dimensions, which refer to the size and extent of a target, can influence the imaging result significantly. Because metamorphism and deformation can cause orebodies to become contorted and discontinuous and because even economic ore deposits can be small relative to seismic wavelengths [35,58,61], the effects of dimensions needed to be further clarified for orogenic gold deposits. For exploration experiences in sedimentary basins, it is generally believed that the seismic reflection method has inherent limitations of resolution in both the vertical and horizontal directions. We considered the limitations of resolution by mainly drawing on those widely used in the hydrocarbon and base metal industries. The minimum thickness t min and diameter d min that can be resolved by reflection are described by the quarter-wavelength criterion and the first Fresnel zone, respectively, as where λ presents the average wavelength in the backgrounds, and z is the depth of the target [62]. Two models ( Figure 2) were designed to show the effects of dimensions. The elastic parameters used in these models are listed in Table 2. P-wave sources with a primary frequency of 60 Hz were placed at the surface with an average P-wave wavelength of 78 m. The models were discretized through 5 by 5 m cells, and the time step was 0.5 ms. Receivers were also located at the surface with a space of 5 m.   Table 2.
To verify the effects of the target dimensions shown in synthetic seismograms, the reverse time migration (RTM) method was introduced to image the orebodies based on traditional compressional wave seismic exploration. A total of 101 sets of single-shot records were used with a shot interval of 50 m. Figure 3a shows that the imaging result of model D clearly delineated the extension range and thickness of Orebody Ib, whereas Orebody Ia was imaged as a point. In the image of model E (Figure 3b), the additional Orebody II was presented like a curve and Orebody Ia failed to be imaged.  Table 2. Two Type I orebodies were involved in model D (Figure 2a-d). The smaller one (Orebody Ia) at about a 1300 m depth was 84 m wide and 25 m thick. The larger one (Orebody Ib) was located at a depth of 1800 m and was 3160 m in length and 68 m in thickness ( Figure 2a). The thickness of the smaller one was slightly higher than t min but the length was lower than d min , while the larger one achieved the resolution requirements in both directions. The synthetic seismogram ( Figure 2b) shows that Orebody Ia generated one simple event. Reflections occurred twice when the incident wave propagated through Orebody Ib. The signal from the top interface was clearly distinguished from that reflected from the bottom interface, roughly characterizing the shape of the reflector. The snapshot at 0.4 s ( Figure 2c) shows that Orebody Ia performed as a new point source, which caused the energy to spread in all directions, forming a circular wavefront. The events corresponding to Orebody Ib ( Figure 2d) had their strongest amplitudes in the down-dip direction. At both ends, new circle wavefronts were clearly observed. These events propagated backward to the surface and were tangential to the reflected events on the synthetic record ( Figure 2b), which helped to determine the extension of reflector.
Model E (Figure 2e-h) was composed of model D and an additional Type II orebody (Orebody II) placed closely under the smaller Type I orebody, which was 1762 m long and 34 m thick ( Figure 2e). The thickness of Orebody II was slightly higher than t min and the length was much higher than d min . The Type II orebody generated a stronger reflection than the Type I orebodies because of its higher reflection coefficient. The weak energy corresponding to Orebody Ia was dimly observed at 0.34 s (white arrow in Figure 2g). At 0.51 s (Figure 2h), nevertheless, this upgoing energy was interfered with by that reflected from the underlying Orebody II and could hardly be discerned. The snapshot displays a combination of prominent reflection and weak diffraction from Orebody II. The waves propagating backward formed two groups of wavefronts, which intersected with each other because of the bend of the orebody. As was expected from the snapshots, the seismogram recorded at the surface ( Figure 2f) showed an absence of signals relating to Orebody Ia. Signals from the underlying Orebody II took on a cross shape. The events corresponding to the top and the bottom of the orebody were so close that it was very difficult to distinguish them.
To verify the effects of the target dimensions shown in synthetic seismograms, the reverse time migration (RTM) method was introduced to image the orebodies based on traditional compressional wave seismic exploration. A total of 101 sets of single-shot records were used with a shot interval of 50 m. Figure 3a shows that the imaging result of model D clearly delineated the extension range and thickness of Orebody Ib, whereas Orebody Ia was imaged as a point. In the image of model E (Figure 3b), the additional Orebody II was presented like a curve and Orebody Ia failed to be imaged.

Dip
Two models containing a shallowly dipping (25 • ) and a steeply dipping (70 • ) orebody were designed in this part to reveal the effects of dip (Figure 4a,e). The corresponding elastic parameters are displayed in Table 3. The spatial and temporal sampling intervals were 5 m and 0.5 ms. The source wavelets injected at the surface were Ricker signals whose center frequency was 70 Hz.

Dip
Two models containing a shallowly dipping (25°) and a steeply dipping (70°) orebody were designed in this part to reveal the effects of dip (Figure 4a,e). The corresponding elastic parameters are displayed in Table 3. The spatial and temporal sampling intervals were 5 m and 0.5 ms. The source wavelets injected at the surface were Ricker signals whose center frequency was 70 Hz.    Table 3.
Direct waves, reflected waves, and diffracted waves were observed in the synthetic data (Figure 4b,f) for both models. Little difference between the two seismograms was found except for some nuance in the curvature corresponding to reflections. The dip of the orebodies could not be directly identified through the records. Snapshots at 0.325 s (Figure 4c,g) showed that the two models generated apparently different wavefronts. The reflected energy was mainly concentrated in the down-dip direction, that is, the shallowly dipping orebody generated near-horizontal reflected wavefronts, while the steeply dip-  Table 3. Direct waves, reflected waves, and diffracted waves were observed in the synthetic data (Figure 4b,f) for both models. Little difference between the two seismograms was found except for some nuance in the curvature corresponding to reflections. The dip of the orebodies could not be directly identified through the records. Snapshots at 0.325 s (Figure 4c,g) showed that the two models generated apparently different wavefronts. The reflected energy was mainly concentrated in the down-dip direction, that is, the shallowly dipping orebody generated near-horizontal reflected wavefronts, while the steeply dipping orebody induced near-horizontal wavefronts.
The RTM imaging approach was also applied to estimate the exploration effect. A total of 61 P-wave sources were distributed at the surface to generate single-shot records with the source interval being 50 m. A sharp and clear imaging result was obtained for the shallowly dipping orebody. By contrast, the steeply dipping orebody was unclear on the image, especially the lower part. Based on the propagation direction of reflected waves, we designed another model (model H). The model contained the identical orebody used in model G, but a longer receiver array of 4.5 km was used to generate multiple shot records for the RTM image. The imaging result of this steeply dipping orebody in model H was relatively clear ( Figure 5).

Geology and 2D Model
The giant Zaozigou orogenic Au-Sb deposit in central China is located in the northwesternmost part of the West Qinling Orogen ( Figure 6) [35,63,64]. It hosts 118 t of Au, averaging 3.42 g/t and > 0.12 Mt Sb, averaging 0.99%, representing one of the largest active gold producers in the West Qinling Orogen [65,66]. The Zaozigou deposit is hosted in the metasedimentary rocks of the Triassic Gulangdi Formation slate, and Triassic dacite dikes are emplaced into the unit [67,68]. Two styles of orebodies have been identified at Zaozigou based on their geometries (Figures 6 and 7). Orebodies with the 'Au' prefix (e.g., Au1) are the steeply dipping (60°-90°) orebodies. Orebodies with the 'M' prefix (e.g., M9) are the shallowly dipping orebodies, most of which dip between 19° and 40°. The largest Au1 orebody is known to be 800 m in strike length and extends in the down-dip direction for at least 1600 m. The shallowly dipping orebodies are commonly approximately 100 m in strike length and extend approximately 200 m vertically. The ore minerals of the Zaozigou deposit are primarily stibnite, pyrite, and arsenopyrite, followed by chalcopyrite, galena, sphalerite, and native gold. The gangue minerals mainly consist of quartz, sericite, calcite, feldspar, and dolomite with minor epidote, apatite, titanite, zircon, and monazite [65].
After a dozen years of miniyng, the resources in the shallow area have been nearly exhausted. There is an urgent need to search for new resources at depth. In addition, the

Geology and 2D Model
The giant Zaozigou orogenic Au-Sb deposit in central China is located in the northwesternmost part of the West Qinling Orogen ( Figure 6) [35,63,64]. It hosts 118 t of Au, averaging 3.42 g/t and >0.12 Mt Sb, averaging 0.99%, representing one of the largest active gold producers in the West Qinling Orogen [65,66]. The Zaozigou deposit is hosted in the metasedimentary rocks of the Triassic Gulangdi Formation slate, and Triassic dacite dikes are emplaced into the unit [67,68]. Two styles of orebodies have been identified at Zaozigou based on their geometries (Figures 6 and 7). Orebodies with the 'Au' prefix (e.g., Au1) are the steeply dipping (60 • -90 • ) orebodies. Orebodies with the 'M' prefix (e.g., M9) are the shallowly dipping orebodies, most of which dip between 19 • and 40 • . The largest Au1 orebody is known to be 800 m in strike length and extends in the down-dip direction for at least 1600 m. The shallowly dipping orebodies are commonly approximately 100 m in strike length and extend approximately 200 m vertically. The ore minerals of the Zaozigou deposit are primarily stibnite, pyrite, and arsenopyrite, followed by chalcopyrite, galena, sphalerite, and native gold. The gangue minerals mainly consist of quartz, sericite, calcite, feldspar, and dolomite with minor epidote, apatite, titanite, zircon, and monazite [65]. 022, 12, x FOR PEER REVIEW 10 of 18 Figure 6. Geological sketch map of the Zaozigou Au-Sb deposit (modified after [66]).
The 2D model was a 1.3 km deep by 1.1 km long profile (Figure 7), in which the geologic interfaces were constrained by field observations and internal reports of the Zaozigou mine. A stibnite vein was chosen to represent all the orebodies for modeling, given that stibnite is the primary ore mineral of the Zaozigou deposit [65,[69][70][71]. Petrophysical properties were assigned to each individual lithology present in the model. Densities were measured by the mines. The P-wave velocities of slate and dacite were obtained from Salisbury et al. [58] and that of stibnite was obtained from Liu et al. [72]. Then, the S-wave velocities were estimated from the corresponding P-wave velocities using a Vp/Vs value of √3. Detailed elastic parameters are shown in Table 4.
After the construction of the model, we populated this section with the P-wave velocity, S-wave velocity, and density per meter. A P-wave source centered at 130 Hz, which was obtained by bedrock-coupled dynamite or vibroseis [45] and sampled every 0.1 ms, was used to generate synthetic seismic data. The source was located at X = 1 km and Z = 0 km (see the pentagram in Figure 7). An absorbing boundary condition was applied. The elastic displacement components were recorded at the surface with 1101 receivers placed per meter.   Figure 6 at the Zaozigou deposit (modified after [65,66]).

Finite-Difference Modeling Resultes
The synthetic seismograms (Figure 8a,b) showed significant reflected signals, especially multiple sets of reflected P-waves and converted S-waves, which were generated when waves went through the steeply dipping orebodies. Reflected P-waves from the top After a dozen years of miniyng, the resources in the shallow area have been nearly exhausted. There is an urgent need to search for new resources at depth. In addition, the steeply dipping orebodies in Zaozigou trend northeast, while the shallowly dipping orebodies, which were discovered during mining, typically trend northwest (Figure 6). At the planning stage, the exploration lines were designed to be perpendicular to the strike direction of the steeply dipping orebodies ( Figure 6). Consequently, the drillholes were arranged along the northwest direction, leading to poor control of the shallowly dipping orebodies. It is of great significance for the estimation of reserves and mining planning of the Zaozigou deposit to determine the extension range of these orebodies. However, a large number of metal materials introduced during production restrict the traditional electric and magnetic methods. In view of the situations mentioned above, the seismic exploration method may be a good choice for the Zaozigou deposit. In the next section, we simulated the application of the seismic method in Zaozigou with a surface configuration by finite-difference modeling and then analyzed the expected results.
The 2D model was a 1.3 km deep by 1.1 km long profile (Figure 7), in which the geologic interfaces were constrained by field observations and internal reports of the Zaozigou mine. A stibnite vein was chosen to represent all the orebodies for modeling, given that stibnite is the primary ore mineral of the Zaozigou deposit [65,[69][70][71]. Petrophysical properties were assigned to each individual lithology present in the model. Densities were measured by the mines. The P-wave velocities of slate and dacite were obtained from Salisbury et al. [58] and that of stibnite was obtained from Liu et al. [72]. Then, the S-wave velocities were estimated from the corresponding P-wave velocities using a Vp/Vs value of √ 3. Detailed elastic parameters are shown in Table 4. After the construction of the model, we populated this section with the P-wave velocity, S-wave velocity, and density per meter. A P-wave source centered at 130 Hz, which was obtained by bedrock-coupled dynamite or vibroseis [45] and sampled every 0.1 ms, was used to generate synthetic seismic data. The source was located at X = 1 km and Z = 0 km (see the pentagram in Figure 7). An absorbing boundary condition was applied. The elastic displacement components were recorded at the surface with 1101 receivers placed per meter.

Finite-Difference Modeling Resultes
The synthetic seismograms (Figure 8a,b) showed significant reflected signals, especially multiple sets of reflected P-waves and converted S-waves, which were generated when waves went through the steeply dipping orebodies. Reflected P-waves from the top and bottom of the M6 orebody were identified on the vertical component record (Figure 8a). Energy diffracted from the endpoints of the M6 and Au16 orebodies was also observed on this record. These diffracted events were tangential to reflected events. They partially overlapped with the events reflected from the subface of the M6 orebody and appeared to be branches of the main reflection events. The S-waves were more distinct on the horizontal component record (Figure 8b). Diffractions from the M4 and M7 orebodies were observed. Additionally, the diffractions from the endpoints of the M6 orebody clearly showed two parallel events. The prefix "R" represents reflected events, and "D" represents diffracted events. The following "P" and "S" indicate P-waves and S-waves. For example, RP-M6 represents the reflected P-waves from M6 orebody. See rock properties in Table 4. A video on the vertical component of the simulated wavefield containing both P-and S-waves is available online at doi: www.mdpi.com/xxx/s1. The P-wave reverse-time migration method was used to predict the exploration results of the Zaozigou deposit. A total of 131 sources were successively placed on the surface with an interval of 10 m. For each shot, the record length was 0.4 s. Most geological interfaces in our model were observed in the final imaging result (Figure 9). Among the steeply dipping orebodies, the Au1 (top 700 m), Au9, Au121, and Au122 orebodies were well imaged. Only the top and bottom were visible for the remaining steeply dipping orebodies (Au16 and Au144). The portion of the Au1 orebody below a 700 m depth failed to be imaged. The shallowly dipping orebodies were imaged with considerable detail in terms of their extension, thickness, and shape. The prefix "R" represents reflected events, and "D" represents diffracted events. The following "P" and "S" indicate P-waves and S-waves. For example, RP-M6 represents the reflected P-waves from M6 orebody. See rock properties in Table 4. A video on the vertical component of the simulated wavefield containing both P-and S-waves is available online at Supplementary Materials.
On the snapshot of the compressional wave wavefield at 0.066 s (Figure 8c), reflections from the top and bottom boundaries of the M6 orebody, as well as those from the left and right boundaries of the Au1 orebody, were the most conspicuous events. The amplitudes of the reflected waves from the M6 orebody had their maximum in the down-dip direction. The endpoints of the M6 and Au16 acted as new sources, which generated diffracted waves propagating throughout the whole medium. The snapshot of the shear wave wavefield at 0.121 s (Figure 8d) showed that a large number of converted S-waves were generated. Several parallel events were induced by diffraction from the endpoints of the M6 orebody. It is worth noting that upgoing waves documenting orebodies at depth were generated through multiple reflections and diffractions in this model.
The P-wave reverse-time migration method was used to predict the exploration results of the Zaozigou deposit. A total of 131 sources were successively placed on the surface with an interval of 10 m. For each shot, the record length was 0.4 s. Most geological interfaces in our model were observed in the final imaging result (Figure 9). Among the steeply dipping orebodies, the Au1 (top 700 m), Au9, Au121, and Au122 orebodies were well imaged. Only the top and bottom were visible for the remaining steeply dipping orebodies (Au16 and Au144). The portion of the Au1 orebody below a 700 m depth failed to be imaged. The shallowly dipping orebodies were imaged with considerable detail in terms of their extension, thickness, and shape. Minerals 2022, 12, x FOR PEER REVIEW 13 of 18 Figure 9. RTM imaging result of the Zaozigou model. Orebodies with yellow margin were successfully imaged in our experiment, while orebodies with white margin failed to be imaged.

Effects of Factors
The effects of the petrophysical properties on the reflected wavefields were investigated by comparison of the finite-difference simulations for three simplified models (models A to C) with various lithologies of wall rocks (Figure 1). The comparison among these models showed that the reflected energy increased with the reflection coefficient R. This implied that the R coefficient played an essential role in the reflection seismic exploration of the orogenic gold deposits by controlling the intensity of the reflected energy, which was similar to the effects in base metal and volcanogenic massive sulfide deposits [21]. Specifically, Salisbury et al. [62] proposed a minimum coefficient of R = 0.06 after allowance has been made for noise. Apparently, this rule also fit our models, that is, when R was higher than 0.06, significant reflections were generated at the contact (Figure 1c-f). Thus, it is strongly recommended to obtain density and velocity values under realistic geological conditions and assess whether the exploration target meets the requirement before conducting seismic surveys.
The forward modeling and imaging experiments indicated that the dimensions of the targets significantly affected the recognition precision of the seismic waves (Figures 2 and  3). If the exploration target meets the resolution requirements in both the vertical and horizontal directions, such as for Orebody Ib in model D, it will be easily identified by seismic surveys, with a sharp top and bottom (Figure 2b and Figure 3a). When the target is much higher than the resolution in a certain direction and slightly higher than the resolution requirement in the other direction, such as for Orebody II in model E, it may be recognized as a surface (Figure 3b). Moreover, if target is close to the resolution in one direction but lower than the resolution in the other direction, such as for Orebody Ia in model D, it will influence the wave propagation as a diffractor [21] and be imaged as a point [62]. However, if this target is underlain by a strong reflector, it may not be discovered on seismograms and imaging results (Figure 2f and Figure 3b).
The dip of the interface was the third factor concerned in this paper. The comparison between models F and G showed that the wave propagation direction mainly depended on the orebody dip ( Figure 4). Shallowly dipping interfaces, such as the orebody in model F, generally led to convergent reflected energy (Figure 4c). Thus, a short range of receivers

Effects of Factors
The effects of the petrophysical properties on the reflected wavefields were investigated by comparison of the finite-difference simulations for three simplified models (models A to C) with various lithologies of wall rocks (Figure 1). The comparison among these models showed that the reflected energy increased with the reflection coefficient R. This implied that the R coefficient played an essential role in the reflection seismic exploration of the orogenic gold deposits by controlling the intensity of the reflected energy, which was similar to the effects in base metal and volcanogenic massive sulfide deposits [21]. Specifically, Salisbury et al. [62] proposed a minimum coefficient of R = 0.06 after allowance has been made for noise. Apparently, this rule also fit our models, that is, when R was higher than 0.06, significant reflections were generated at the contact (Figure 1c-f). Thus, it is strongly recommended to obtain density and velocity values under realistic geological conditions and assess whether the exploration target meets the requirement before conducting seismic surveys.
The forward modeling and imaging experiments indicated that the dimensions of the targets significantly affected the recognition precision of the seismic waves (Figures 2 and 3). If the exploration target meets the resolution requirements in both the vertical and horizontal directions, such as for Orebody Ib in model D, it will be easily identified by seismic surveys, with a sharp top and bottom (Figures 2b and 3a). When the target is much higher than the resolution in a certain direction and slightly higher than the resolution requirement in the other direction, such as for Orebody II in model E, it may be recognized as a surface (Figure 3b). Moreover, if target is close to the resolution in one direction but lower than the resolution in the other direction, such as for Orebody Ia in model D, it will influence the wave propagation as a diffractor [21] and be imaged as a point [62]. However, if this target is underlain by a strong reflector, it may not be discovered on seismograms and imaging results (Figures 2f and 3b).
The dip of the interface was the third factor concerned in this paper. The comparison between models F and G showed that the wave propagation direction mainly depended on the orebody dip ( Figure 4). Shallowly dipping interfaces, such as the orebody in model F, generally led to convergent reflected energy (Figure 4c). Thus, a short range of receivers alone was sufficient for the exploration target, whereas steeply dipping interfaces, such as the target in model G, generally gave rise to diverging wavefronts (Figure 4g), inducing relatively weak reflected signals above the target (black arrows in Figure 4f,g). This means that a longer source-detector distance was needed to obtain adequate information about the target. The difference in the propagation direction was difficult to identify on seismic records but was clearly exhibited in the final imaging results. A shallowly dipping orebody (e.g., an orebody with a dip of 25 • ) could be clearly imaged on the seismic profile (Figure 4d), whereas a steeply dipping orebody (e.g., an orebody with a dip of 70 • ) could not obtain an imaging result with sharp borders or even be invisible on the resulting image (Figure 4h). The significantly improved imaging result of the steeply dipping orebody in model H indicated that increasing the length of the receiver array could address the issue induced by large dips. In practice, clear imaging results may require smaller dips given that the limited signal-to-noise ratio of field data and interference from other units were not considered in our simple models.

Implications for the Reflection Seismology Method on Orogenic Gold Deposits
Lithology units (metasedimentary rocks, granitic rocks, and sulfides) such as in the Zaozigou model are common in orogenic gold deposits [35,73]. Our calculation exhibited that the reflection coefficients among these units were higher than 0.06, fulfilling the requirement for significant reflection mentioned above. As a result, it should be possible to use reflection seismic techniques for orogenic gold exploration and delineation from the perspective of petrophysical properties. Our FD modeling showed that multiple reflections and diffractions formed a complex wavefield recorded on the surface (Figure 8a,b). These signals presented a wealth of information on subsurface structures, offering the opportunity of detecting orebodies through a proper data acquisition and processing method, which was further verified by the following RTM imaging (Figure 9). In the Zaozigou model, most of the shallowly dipping orebodies were well imaged with a precise extension and thickness ( Figure 9). The reflections from the M4 and M7 orebodies were not obvious on the synthetic record, but the diffracted energy caused by the endpoints of these orebodies was clearly observed (Figure 8b). The endpoints of the shallowly dipping orebodies impacted the wavefield as new sources inducing diffracted waves, which propagated to the surface. This energy was recorded individually or formed a tangential relationship with the reflected energy on the seismograms, contributing to the determination of the spatial extent of the shallowly dipping orebodies. Therefore, for shallowly dipping orebodies, which are difficult to determine by boreholes in orogenic gold deposits, the diffracted events corresponding to their endpoints are good indicators for their spatial extent.
Most of the steeply dipping orebodies of the Zaozigou deposit were imaged but lost important thickness information due to the low lateral resolution of the seismic waves ( Figure 9). The Au16 and Au144 orebodies were not clearly imaged because the propagation of efficient signals to the surface was impeded by the shallowly dipping orebodies above them. The imaging of the Au1 orebody below 700 m failed because its large dip led to little useful information being recorded on the surface within the limited receiver array set used in this model. From the preceding Zaozigou model, it is clear that steeply dipping orebodies in orogenic gold deposits can be detected using the reflection seismology method. However, it is difficult to identify their thickness. In addition, steeply dipping orebodies at depth may not be visualized. Moreover, the diverging and relatively weak reflected signals of steeply dipping orebodies will be significantly disturbed by strong reflectors above them. This will lead to imaging failure for small steeply dipping orebodies (e.g., the Au16 and Au144 orebodies at Zaozigou in Figure 9).
For steeply dipping orebodies at depth, which were undetectable in the Zaozigou model, a long source-to-receiver distance that allows more information to be recorded at the surface may address this issue. Additionally, vertical seismic profiling (VSP) techniques [74,75] can provide much information on steeply dipping orebodies at depth through downgoing waves recorded by receivers vertically or near-vertically arranged in a borehole. Furthermore, with a higher lateral resolution [76], they enable the obtaining of the thickness information of steeply dipping orebodies. In conclusion, VSP may be more promising for the exploration of steeply dipping orebodies.
Our modeling showed that the distance between the two reflected S-waves corresponding to the top and bottom (or left and right) of orebodies was longer than that of the P-waves, leading to a better recognition ability of S-waves (Figure 8a,b). Meanwhile, S-waves have a higher resolution relative to their corresponding P-waves, owing to their slower velocities, that is, their shorter wavelengths, according to relations (6)- (7). This implies that S-waves have an advantage in identifying the shape of small-sized objects. Furthermore, based on our observation that most S-wave events have significant amplitudes on the horizontal component record (Figure 8b), we hold the opinion that multiwave and multicomponent seismic techniques [10,77,78] can offer more valuable information about small-scale orebodies in mineral exploration.

Conclusions
The FD forward modeling experiments showed that the petrophysical properties, dimensions, and dip of targets affect the results of reflection seismic exploration. Petrophysical properties determine the intensity of reflections and hence have an essential impact on the detectability of orebodies. The dimensions of targets influence the recognition precision of the morphology of an orebody. The dip of orebodies affects the amount of information recorded at the surface through controlling the propagation direction of reflected waves and further influences the imaging results.
From the results of the Zaozigou model presented, we conclude that it was possible to locate and image orogenic gold deposits using the reflection seismology method. Shallowly dipping orebodies could be well imaged with precise extension and thickness. Steeply dipping orebodies could be recognized but lost important thickness information. Steeply dipping orebodies at depth were not detectable under a surface configuration. This problem could be effectively solved by increasing the array length or using VSP seismic methods. Additionally, VSP seismic methods enable the obtaining of the thickness informing of steeply dipping orebodies. For small orebodies, multiwave and multicomponent seismic techniques can offer more valuable information in mineral exploration.