Pre‐stack depth imaging techniques for the delineation of the Carosue Dam gold deposit, Western Australia

In this study, we explore the latest generation of seismic imaging algorithms (migration) that have been successful in the oil and gas exploration industry and apply them in the more challenging hard rock environment. Seismic migration is a crucial processing step required to build an accurate image of the Earth's subsurface. The seismic method applied in a hard rock environment has a specific set of challenges: complex geological settings often comprised of steeply dipping interfaces and heterogeneities (faults, fracture zones, thrusts, etc.); spatially variable zones of alteration, low intrinsic signal‐to‐noise ratio; complex near‐surface conditions (the weathered overburden has a very high contrast in seismic properties with base formations). Here, we present how a dense source–receiver 3D grid in combination with the latest pre‐stack depth imaging techniques can image the geology with a remarkable level of detail and to depths as shallow as the top of fresh rock. We also propose a comprehensive velocity model building strategy applied specifically to the Carosue Dam deposit, Western Australia, which is mainly characterized by volcaniclastic and volcanic rocks within the Carosue Basin disrupted by a complex system of gold‐bearing faults. To arrive at an optimized velocity model for pre‐stack depth imaging, we combine drillhole data with tomographic refinement. Excellent correlation between resultant seismic images and refraction tomography, vertical seismic profiling and drillholes can be attributed to the integrated velocity model building workflow used for depth migration.


INTRODUCTION
The golden era of geophysical discoveries of shallow, rich mineral deposits has well passed. Exploration targets have now moved to a much greater depth window. In general, depths below 500 m are difficult to investigate with precision by any other geophysical method except the seismic 1950s when the first high-resolution seismic study was conducted by Berson (1957), numerous and valuable seismic studies followed, demonstrating the high penetration, resolution and precision of reflection seismic in detecting complex geological underground forms at all depths. The history, development and application of modern seismic methods for the exploration of mineral resources across all continents have been well summarized by Malehmir et al. (2012). While the potential of seismic was clearly demonstrated through several early field studies (Durrheim & Maccelari, 1991;Milkereit et al., 1992;Pretorius et al., 1989), the mining industry was reluctant to embrace it on a larger scale, certainly incomparable to the utilization of potential field and electromagnetic methods. Until 1994, seismic for mineral exploration was exclusively conducted along a profile (2D). The first breakthrough came with the areal (3D) seismic design recorded in South Africa (De Wet & Hall, 1994). Subsequently, Eaton et al. (1997) used 3D seismic for the delineation of a deep ore deposit in Canada, while Pretorius et al. (1997) conducted the first 3D survey for the purpose of underground mine planning and development. These early works demonstrated that the complexity of the hard rock environment is best assessed and delineated by 3D reflection seismic. Instead of flourishing, a significant downturn in the application of 3D seismic followed until 2011-2012. Since then, a slow but steady rise in the use of 3D reflection seismic for exploration in brown fields, to extend the existing mine life and discover new resources at a greater depth, has occurred (Urosevic, 2013;Urosevic et al., 2017).
Extensive structural complexity typically encountered in a hard rock environment requires high-precision measurements. In the case of seismic measurements that translate into requirements for high spatial and temporal resolution as recognized in early seismic works (Stolz et al., 2004;Urosevic et al., 2005). The high spatial density of sensors is of key importance for imaging structurally controlled mineralization . Unfortunately, a high-resolution, highfidelity seismic approach requires a large number of sensors to be simultaneously deployed over a large territory which makes such investigations very expensive. It is possible, based on reciprocity, to partially compensate for a small sensor count by introducing additional source positions or repeating specific shooting patterns. Indeed, such an approach characterized hard rock 3D seismic surveys in the previous decade (Urosevic et al., 2017;Ziramov et al., 2015Ziramov et al., , 2016. One to two thousand receivers were typically deployed for 3D surveys at the time using moving, overlapping source-receiver templates. While this approach made 3D seismic more accessible and affordable to the mineral industry, it suffered from several drawbacks. Insufficient offset range is probably the most important issue that in some limited cases could have affected the imaging of steep dips at depth, perhaps outside of the mining window but still of interest to structural geologists. The potential for high sensor density designs for imaging the complex geometry of salt domes was hinted at in 1990 in the oil industry, while the first (primitive) nodal systems referred to distributed cable recordings dating back to the 1970s (Dean et al., 2018). The real 'explosion' in the production and utilization of modern nodal systems is recorded only in the middle of the last decade in the oil sector. The benefit of a high receiver count for imaging complex structures in a hard rock environment at all depths was quickly realized in the minerals industry and deployed even for 2D surveys (Naghizadeh et al., 2019). Similarly, modern vibroseis and electro-kinetic vibrating sources provide a broad-band seismic signal that is required for the utilization of novel inversion and imaging techniques (Naghizadeh et al., 2019;Pertuz et al., 2022;Singh et al., 2022). A combination of a large number of sensors (high data density) and broad-band seismic sources enables the successful implementation of novel processing (signal enhancements, interpolations, etc.), imaging techniques such as reverse time migration (RTM) and inversions such as full-waveform inversion (FWI). This led to early investigations into the potential of FWI in hard rock environments (Adamczyk et al., 2014(Adamczyk et al., , 2015, while investigations into imaging techniques were documented by several authors (Brodic et al., 2021;Broüning et al., 2020;Ding & Malehmir, 2021;Hloušek et al., 2015;Singh et al., 2022).
The successful application of the seismic reflection method is to the first degree controlled by the contrast and the distribution of elastic properties of the geological environment (Bohlen et al., 2003). However, mineral composition and more elusive rock characteristics such as alteration and rock texture play an important role in recorded reflectivity (Schetselaar et al., 2019). Our ability to delineate reflectivity also depends on its geometry and volume, as well as the presence of fractures and joints, grain alignment and crystal symmetry, etc. Therefore, reflectivity in hard rock is controlled by many factors and it is not readily predictable even in well explored mine sites. This possibly peculiar reflectivity pattern of hard rocks is further exacerbated by high structural complexities. In terms of gold exploration in Western Australia, both the reflectivity patterns and excessively complex underground structures are commonplace (Stolz et al., 2004). The necessity to utilize high-precision imaging techniques, such as Kirchhoff pre-stack depth migration (PSDM) to delineate complex undergraduate structures, was recognized early in the application of seismic for gold exploration across Yilgarn craton of Western Australia (Urosevic et al., 2007;Urosevic & Evans, 2007a). Consequently, the application of seismic for gold exploration in Western Australia must be carefully staged. Hence, seismic data acquisition is typically proceeded by a seismic feasibility study looking into the seismic rock properties. Positive results will grant project continuation. After the acquisition, data processing and analysis should adopt the best approach for imaging complex geological structures at all scales with the maximum possible precision.
Such an approach was carried out at the Carosue Dam site. Extensive rock property measurements were conducted prior to 3D data acquisition. As the primary objective of the survey was to refine the existing geological model needed for the discovery of additional gold resources within the Carosue Dam Operations, maximum precision seismic imaging was required. That included the delineation of numerous smallsized, steeply dipping brittle-ductile faults that facilitate the emplacement of gold mineralization at this site.
To accomplish such a task, we evaluated several approaches to velocity model building, composed of refraction and iterative reflection tomography, to construct high-precision seismic images using concurrently pre-stack Kirchhoff depth migration and reverse time migration. We describe a complete depth imaging workflow rooted in the standard oil and gas industry to depth imaging tailored to a hard rock environment. This workflow is applied to a high-density, high-quality, large (50 km 2 ) 3D seismic survey shot over the complex gold mineralization target in WA (Carouse Dam). The resultant images correlated very well with the existing downhole logs, vertical seismic profiling (VSP) and drillhole information that provide validation of the results. High-quality seismic images obtained in this study gave rise to high-quality geological analysis and interpretation that are hoped to define new drilling targets soon.

Geological background
The northwest-trending Carosue Basin is located in the Kurnalpi Terrane of Swager (1997) and the Murrin Murrin domain of Cassidy et al. (2006) (Cassidy et al., 2006). The Menangina domain mostly comprises mafic volcanic and intrusive rocks, ultramafic intrusions, felsic volcanic rocks and minor fine-grain sedimentary rocks. East of the Carosue Basin is the older Murrin Murrin domain (Cassidy et al., 2006) that comprises mafic and ultramafic volcanic and subvolcanic rocks with minor epiclastic sedimentary rocks. The contact between the Carosue Basin and the Murrin Murrin domain is obscured locally by sandstone, siltstone and polymictic conglomerate of the younger Yilgangi Basin, which is intruded by ca. 2.66 Ga syenite and monzodiorite (Nelson, 1996). The eastern part of the Carosue Basin and isolated remnants of the Yilgangi Basin in the Carosue Dam area are strongly deformed with tight folding and stretching fabrics.
The general architecture of the two largest gold deposits in the Carosue Basin, Karari and Whirling Dervish, is steep northeast-dipping tabular to lenticular beds that are moderately foliated with localized steeply northward-plunging transposition folds and associated shearing. The stratigraphic sequence ( Figure 1) from west to east in the area covered by the 3D seismic survey is (i) mostly northeast-dipping Menangina mafic-ultramafic succession, (ii) a mostly northeast-dipping, quartz-bearing turbidite sequence with local disruption by folding and thrust faulting (iii) the informally termed central mineralized zoned (CMZ; Witt et al., 2009) comprising mostly volcaniclastic sandstones, minor andesite and minor conglomerate, (iv) a variably thick (50-250 m) sequence of andesite, intruded by monzonite, syenite and lamprophyre (the Eastern Intrusive Complex or EIC; Witt et al., 2009) that forms the hanging wall to most of the gold deposits in the Carosue Basin, (v) intercalated volcaniclastic sandstone, conglomerate and minor andesite and (vi) a thick sequence of fine-to medium-grained epiclastic sediments. The margin between the Carosue Basin and the Murrin Murrin rocks to the east is obscured by epiclastic sandstones and polymictic conglomerates of the Yilgangi Basin that are strongly foliated and folded. This zone of high strain is considered to be the likely trace of the regionally significant Keith Kilkenny Shear Zone.
Within the Karari and Whirling Dervish gold deposits, the EIC succession has undergone high shear strain and hydrothermal alteration. Across this structural domain, the reversal in younging directions (east-facing in the CMZ; overturned west-facing stratigraphy east of the EIC) demonstrates the deposits are sited on a district-scale, faulted-out fold closure with a shallow north plunge. Shearing in the CMZ is generally concentrated on the contacts between the mechanically weaker andesitic tuffaceous conglomerates or lamprophyre dykes and more competent sandstone or monzonite units. Gold mineralization is generally within these sheared domains and in the heavily fractured and pervasively altered zones in the competent units (see Witt et al., 2009, for further details on the gold mineralization). Numerous northto-northeast trending, steeply dipping brittle-ductile faults dissect the Carosue Basin and offset mineralization. The Osman Fault has dextrally offset the Karari and Whirling Dervish deposits by around 1 km. Similar chloritic faults have developed as late reactivation of earlier-formed mineralized shear zones. The Resurrection fault and Young Fella fault, which in places includes cataclasite fault gouge, are such examples.

Rock properties: seismic feasibility study
The seismic exploration strategy for Carosue Dam utilized a 3-phased approach to determine whether the local geology was amenable to modern seismic reflection imaging. The first phase completed in 2017 involved the characterization of the rock properties and predicting seismic reflectivity via measurements of P-wave velocity (Vp) and specific gravity on the core samples. The second stage of seismic activity involved the acquisition of 2D and borehole seismic, which demonstrated excellent reflectivity from key lithologies and structures within the area of interest. This is illustrated in F I G U R E 2 Rock property measurements in a mineralized zone associated with the Whirling Dervish orebody (borehole WDX047), with the highlighted area showing strong P-impedance changes from wireline and VSP measurements Figure 2, by comparing wireline and VSP measurements against sampled mineralized zones. WDX047 is the deepest borehole acquired in the Whirling Dervish area, intercepting mineralization at 1.1-1.6 km of depths. Synthetic traces computed from P-impedance from wireline logging show high reflectivity in areas of gold-bearing monzonite dykes. Following the successful 2D and wireline programme, the HiSeis team proceeded with the third stage of the seismic investigation, whereby a 50-km 2 3D seismic survey was acquired over an area covering the Karari and Whirling Dervish deposits.

3D seismic data acquisition
Acquisition of 3D seismic utilized two vibroseis trucks, each having a mass of 27,200 kg, in flip flop shooting order. Data acquisition geometry utilized a rolling swath design using INOVA's Hawk wireless seismic acquisition system. Recording parameters are listed in Table 1.
The acquisition phase took place in early 2019 with the final survey design illustrated in Figure 3, showing the distribution of azimuths and offsets of the survey. The active shouting template highlights a maximum inline of 2.2 km and a maximum crossline offset of 1.2 km. It is a common occurrence to acquire hard rock seismic surveys through operating mine sites, with access constraints around key infrastructure, including the mining operations, pits, waste dumps and tailings storage facilities. To reduce acquisition footprint, offset regularization was applied prior to data imaging.
Seismic data were processed shortly after acquisition, and imaging was carried out using a Kirchhoff pre-stack time migration (PSTM). This initial volume showed encouraging results, which provided new insights into the geology of the area. In this study, we have revisited imaging of the data with the aim of trialling pre-stack depth migration algorithms, in order to improve the imaging of small complex structures that are known to exist in the target zone but were not fully resolved with time imaging. Our expectation was that depth imaging will not only increase the imaging accuracy but will provide better resolution and hence enable less ambiguous geologic interpretation. The re-processing of the Carosue Dam data was conducted through the following steps: seismic data processing, pre-stack depth imaging and velocity model building.

Seismic data processing
The reflection seismic method faces many challenges in a hard rock geological environment. Often high noise levels due to mining activity, a complex overburden, overthrust zones, fracture zones and steeply dipping stratigraphy require us to constantly advance our approach to seismic processing and imaging. The Carosue Dam environment is known to exhibit all of these challenges. In order to address these, a tailored processing sequence was designed and is shown in Table 2. The seismic processing can be grouped into the following categories: • Signal-to-noise (S/N) improvements • Wavelet shaping and conditioning (deconvolution) • Near-surface static corrections • Corrections for positioning (imaging) The higher data density characteristic of hard rock survey designs provided proper wavefield sampling in the spatial domain, which supported the effective application of sourcegenerated noise removal techniques. Refraction statics were calculated using refraction tomography-based inversion. The process involves measuring the arrival times of refracted seis-mic energy. These travel times are then inverted for seismic velocities at each layer and their corresponding depths. This is because the angle of the refracted event directly depends on the ratio of velocities between the two layers (Snell's law). To improve the accuracy of the velocity-depth model, several iterations of tomographic updates were utilized. A refraction and residual static that corresponds to smooth topography (floating datum) has been applied to compensate for smallscale velocity anomalies in the near surface that imaging could not correct. This is one of the main differences relative to legacy processing, where a full static solution that shifts the data to a final flat datum has been applied prior to imaging. Shot gathers before and after data processing are illustrated in Figure 4. It can be observed that the conditioning steps have increased signal-to-noise ratio (SNR) and vertical resolution of the seismic records.

Seismic data imaging
Seismic imaging (migration) is the final step applied to the seismic data that moves reflection events to their true spatial location. Since migration is essentially an inverse process to seismic modelling, its natural computation region is depth.  In practice, however, migration can be applied in either time or depth, before or after CDP stacking. These processes are known as pre-stack and post-stack migrations, respectively. Time migration produces an image in two-way travel time, utilizing a simplified extrapolation of the wavefield. Because of that, it is a more robust process when it comes to imperfections or smaller errors in the velocity model. Depth imaging directly outputs images in depth and accommodates for lateral velocity changes as it incorporates ray bending through a heterogeneous and anisotropic medium, mainly attributed to a more accurate wavefield extrapolation and imaging principle. In this study, we focused on pre-stack 3D depth migration techniques that in recent years have been successfully used in oil and gas exploration. The challenge in a hard rock environment is how to appropriately apply them in a complex geological environment that produces discontinuous reflections with low SNR, such as that experienced at Carosue Dam.
If the basic migration assumptions are met, the quality and accuracy of the migrated image largely depend upon the type of migration. The following are common basic migration assumptions (Dondurur, 2018): • All reflection events are a result of primary reflections, and the data are devoid of noise components and S waves. • Velocities at every sample location of the seismic event are known. • There are no out-of-plane reflections (applicable for 2D migration).
Seismic migration is a reflection imaging process, and we must make sure that the input data only contain reflections. During the seismic data conditioning stage, all non-reflected waves such as surface, air-wave and diving waves are attenuated.
The second assumption is immensely more difficult to accomplish. Vertical velocity distribution can be obtained from a reliable sonic log or velocities from vertical seismic profiling (VSP). However, borehole velocities provide only a one-dimensional function and do not necessarily match with the migration velocity, which undergoes 3D wave propagation. The only way to assess a 3D velocity model is through the surface seismic method.
Imaging structural complexities in the near surface (top 1 km) of the Carosue Basin are of critical importance for exploration and, consequently, time imaging has to be replaced by depth imaging techniques to accommodate for the lateral velocity variations of the subsurface (Buzlukov & Landa, 2013). Lateral velocity variations are often associated with complex geology and the depth migration algorithm has the ability to accommodate such features, including steep dips (Yilmaz, 2001).

METHODS
Kirchhoff's integral migration is a ray-based migration method widely used in both time and depth domains due to its effectiveness and robustness. In the depth domain, it is based on solutions of the wave equation, assuming a highfrequency approximation. Thus, seismic waves approximate rays and ray paths, with the assumption that the scale of the structure is greater than the seismic wavelength (Etgen et al., 2009;Jones, 2003). However, ray-based migration methods rely on a gently varying smooth velocity field for calculating travel times, which make these less accurate than wavefield extrapolation-based methods such as reverse time migration (RTM).
The Kirchhoff depth migration (KDM) process is computed in two stages: 1. Travel time table calculation, using dynamic ray tracing or a solution to the Eikonal equation. 2. Collection of the associated data samples within the migration aperture and their summation.
Our application of the KDM algorithm utilized ray tracing in a 60 × 60 m grid. Offset bin increment was 50 m, with an anti-aliasing filter applied. The maximum migration aperture was 6 km, limited by an 80˚migration angle.
Due to the efficiency of its implementation and computation, we use KDM to iteratively update the velocity model and RTM for the final imaging solution.
RTM is the most effective way to image complex structures and handle strong velocity variations. It uses the wave equation to model complete wavefronts. In general, the RTM method is based on two key steps: 1. Forward propagation of the source wavefield and backward propagation of the receiver (recorded) wavefield. 2. Construction of the image by applying the imaging condition.
RTM is able to solve most imaging challenges but is characterized by expensive computational costs and significant computer memory demands. An increase in maximum frequency causes the run times to increase exponentially. The concept of wavefield extrapolation migration was first proposed 50 years ago (Claerbout, 1971), and it was only in recent years, with the expansion of digital technology, that the RTM method became a widely used imaging method in oil and gas seismic exploration. With the advent of commercial cloud computing, these algorithms are now accessible on a broader scale.
RTM is a two-way wave extrapolation method and employs both the down-going and up-going wavefields. The imaging condition is obtained by cross-correlating the two wave fields at each time step (Claerbout, 1971). However, it is expected that due to differences in the amplitudes of up-going and down-going wavefields, the correlation can introduce longperiod artefacts, especially in areas of strong velocity contrast F I G U R E 5 RTM artefact removal. Due to strong velocity contrast in the near surface, the reflectors in RTM stacked volume are curving beyond 90˚(a). After the application of post-migration mute the artefact disappears (b), the reflectors extend to the TOFR, without curving. Straight and steeply dipping reflectors are verified by a monzonite wireframe from borehole measurements. The mute function is illustrated on a single migrated shot record (c), the iso-velocity of 5700 m/s that approximates TOFR, and a monzonite body is overlain on RTM stack with artefacts removed (d) (Robein, 2010). We observe significant RTM noise in the shallow, near the top of fresh rock (TOFR). This was in turn mitigated through the application of a Laplacian filter and trace muting in the pre-summation process ( Figure 5).
The algorithm is applied on processed shot gathers, producing a migrated volume for each shot. To test computation speed and the parameters for RTM, we migrated 2% of all shot records in the southern region of the survey. Figure 6 shows an assessment of the maximum frequency to be used for the final migration run. The final parameters were a maximum frequency of 60 Hz and an aperture of 6 km to match the KDM imaging distance. It was found that sorting the RTM data into common image point (CIP) gathers comes at a great computational cost as well. For this reason, we have run RTM on the full data set only once using the final velocity model and final choice of parameters. Final depth migrations were computed using a depth interval of 6 m and a CDP binning grid of 7.5×7.5 m. Both algorithms migrate data from a smooth surface directly to the final flat datum. This excludes the need to use a static shift and an assumed replacement velocity to shift data to a final flat datum, reducing ambiguity caused by near-surface velocities.

Handling the near surface
The success of pre-stack depth migration methods largely depends on the accuracy of the velocity model. One of the critical issues the seismic method faces in a hard rock environment is the effect of the low-velocity zone associated with weathering and its consequences on imaging quality. In many hard rock seismic projects, the velocity change from overburden to the fresh rock is drastic and varies over such small distances that even depth imaging cannot fully resolve this problem. In this study, the near-surface velocity model was built using the approach introduced by Singh et al. (2019). The method involves applying a floating datum refraction static to first-arrival times, which are then used in refraction tomography to compute a velocity model for imaging.
A non-linear, first-arrival tomographic travel time inversion (Zhang & Toksöz, 1998) was run until convergence to a unique solution was achieved. Upon each iteration, ray tracing was used to compute synthetic travel times and minimize the difference between the measured times and synthetic first break picks. The process resulted in a detailed velocity model of the near surface (Figure 7) that was used to compute F I G U R E 6 Maximum frequency test for RTM. Highlighted (blue) is the position of the section. Top left is a plan view of the stacked volume. A total of 2% of shot records have been migrated using 40, 50 and 60 Hz maximum frequency. Highlighted in the square is the Karari deposit zone with an expected east-dipping reflector. The reflector is best resolved with a maximum frequency of 60 Hz F I G U R E 7 Velocity model from refraction tomography, confirming NW-SE trend of stratigraphic sequences illustrated in geology map. Velocity-depth slice at 300 m below surface. Low-velocity areas (blue to green) relate to mining infrastructure (pits, waste dumps or tailings dams) refraction tomo-statics, as well as for the interpretation of the TOFR.
As refraction tomography provided velocities down to the TOFR, the initial velocity model has been extended using a pre-stack time migration (PSTM) velocity model converted to interval velocity and smoothed in a 600 × 600 × 600 m (X, Y, Z) grid. A large smoother was used to avoid the introduction of velocity errors in the initial velocity model for migration.

Velocity model building
Velocity updates using pre-stack reflection tomography have been widely used in the oil and gas seismic industry. It requires large enough source-receiver (S-R) offsets, good reflectivity and CIP gathers that have been migrated with an appropriate initial velocity model. The initial velocity model needs to be accurate in the near surface, because low trace fold and small S-R offsets mean the reflection tomography will not be able to update the velocity model.
In Figure 8 we are revisiting the deepest borehole WDX047, where we compare refraction topography using raw first arrivals, to the initial velocity model, final velocity model and the sonic log. We notice that due to different angles of wave propagation, the values of velocities show a substantial difference. Velocity overlays show that values vary from 4800 to 6100 m/s. All models converge at approximately 200 m of depth and a velocity value of 5700 m/s. This iso-velocity is extracted from refraction tomography and used for the interpretation of the TOFR (Figure 16). Down to 200 m depth our initial velocity model was built using the Singh et al. (2019) approach, which is consequently a smooth version of raw refraction tomography, while the deeper velocity is guided by interval velocities from PSTM with a maximum value set to 6100 m/s, based on the maximum velocity of sonic logs. This smoothed initial velocity model will be updated using several iterations of KDM followed by grid-based reflection tomography.
The pre-stack reflection tomography method uses multiple iterations of residual moveout (RMO) analysis on a sparse (60×60 m) grid and depth-migrated CIP gathers to update the initial model. The grid-based tomography method has been described by Woodward et al. (2008). The process is repeated until the velocity error between input and output model is minimized (Tanis et al., 2006). To QC velocity updates a γ attribute was computed. It is a unitless measure of the RMO depth error, represented as a fraction of the reference depth of the primary event on a CIP.
where Z is the depth of the event, Z 0 is the depth at the referenced offset, γ is the gamma factor, and X is the offset.  Figure 9 shows a γ attribute of the full survey at a depth slice of 2.5 km and sample of CIP gathers for the first iteration of velocity model building. Pre-conditioning of gathers involved dip filtering and coherency filters to remove remaining coherent noise and improve S/N ratio, respectively. A block diagram of the process is illustrated in Figure 10. Iterations were focused on correcting velocity discrepancies observed in shallow zones and then progressively worked deeper to update the entire velocity model. Figure 11 shows a comparison of the final KDM stack and final velocity model. We can notice the correlation between reflectors and velocity anomalies. A comparison of CIP gathers shows that the final velocity model better flattens reflectors than the initial model ( Figure 12).

RESULTS
One of the main aims of our analysis was to evaluate the effectiveness and accuracy of depth imaging solutions in a hard rock environment. For this purpose, we validated depth imaging results against the legacy pre-stack time migration (PSTM), refraction tomography, drillholes, modelled wireframes, magnetics and vertical seismic profiling (VSP)-migrated section.
In Figures 13 and 14, we compare the initial PSTM, Kirchhoff depth migration (KDM) and reverse time migration (RTM). The PSTM volume was processed shortly after acquisition of the seismic, independently of this study. Both depth-imaged sections have enhanced steep dips, better resolved the near surface and contained fewer migration artefacts when compared with the legacy PSTM. In comparison to KDM, the RTM results have better-defined steep dips, particularly in the near surface. The higher wavenumber of KDM is beneficial in resolving discreet contacts, while RTM has clearer reflector terminations and faults, such as the Young Fella fault highlighted in Figure 14. The monzonite wireframes modelled from drilling show a good correlation with KDM and RTM sections and confirm the steeper dips that depth imaging resolves well. This correlation is a testament F I G U R E 1 1 Final velocity model depth slice at 1.5 km (a) and W-E section view (b), with corresponding final KDM depth slice (c) and section view (d). The correlation between anomalies in the velocity model and KDM stack migrated with the final velocity model is highlighted with a red arrow to an effective high-density seismic design and an appropriate depth imaging workflow being applied.

DATA INTERPRETATION AND DISCUSSION
To validate the depth imaging results, we compared them against other seismic outputs. Refraction tomography uses first-arrival measurements to compute a detailed near-surface velocity model. The correlation between steep dips observed in the depth imaging volumes and a depth slice from the velocity model shows an impressive correlation ( Figure 15). The near surface is the hardest zone to image, and the match between these two seismic outputs is accurate to within a couple of metres.
In Figure 16, we compare a blended instantaneous phase and energy attribute from the reverse time migration (RTM) volume at a depth slice 778 m below surface with the Carosue Dam geological map as well as a top of fresh rock iso-surface extracted from the refraction tomography with airborne mag-netics. The four geological domains -Footwall Sedimentary Sequence, Central Mineralized Zone, Eastern Intrusive Complex and an Andesite Sequence -are easily recognized in the seismic data sets.
An overlay of drillhole measurements and the depth imaging results is shown in Figure 17a. The 'shape' of the drillhole represents acoustic impedance (AI) while colour corresponds to lithology. A reflection on the seismic image is expected where a change in AI is present. We can observe that reflectors broadly correlate well to AI changes. In some instances, changes in AI are related to the presence of mineralized structures. The high reflectivity domain in the reverse time migration volume correlates with mineralized monzonite intrusions.
Next, we compare our results with vertical seismic profiling (VSP) data acquired during the preliminary de-risking stage of the project in borehole WDX047. Another good correlation between surface and VSP migrated sections was achieved, although the resolution of the VSP is much higher (Figure 17b). This is expected because of the one-way travel path through the attenuative overburden and altogether shorter F I G U R E 1 2 KDM gathers with the corresponding velocity model. CIP gathers migrated with the initial velocity model (a), using the final velocity model (b). Gathers are converted to time using migration velocities and scaled to depth using a final velocity model for comparison purposes. Initial velocity model with superimposed KDM stack and highlighted gather location, W-E section view (c). Final velocity model with final KDM stack (d) propagation distance in VSP, which means less absorption of the higher frequencies. Overall, the match with drillhole data gives us confidence in the velocity model used for imaging as well as in the interpretation of this complex volume.
Several structures have been easily mapped in the RTM volume, with the most significant including the Osman Fault, the Resurrection fault and the reactivated Young Fella fault. By imaging these as well as other geological features, the 3D seismic survey was found to successfully verify, update and extend the current geological model, map key lithologies and generate new drill targets. These can then be used to accelerate subsurface characterization to inform decisions throughout the life of mine.
Finally, we compare the run time of fast-track pre-stack time migration (PSTM), Kirchhoff depth migration (KDM) and RTM in reference to the imaging uplift. We use cloud computing technology to efficiently run processing and imaging flows on a multi-node system. The advantage of cloud computing is in its scalability, as we can use a large number of nodes for most processes. For example, using 100 nodes for 1 h results in the same computing cost as using 1 node for 100 h, except the latter option is not nearly as time-efficient.
Before computing a full-depth imaging process, it is important to estimate node hours required for the job, which will reflect the run cost. PSTM and KDM node hours are estimated by running a single offset plane and multiplying the run time with the total number of offset planes. RTM run time is assessed from a single shot, or for better accuracy, a selection of shots belonging to the same swath. In Figure 6, we illustrated the RTM imaging outcome using only 2% of available shot records, with different maximum frequencies imaged. The estimated full run is extrapolated to be 50 times longer.
Our final PSTM, KDM and RTM runs, using 50 nodes, each node consisting of 30 cores, 256 Gb Random Access Memory (RAM), took 50, 70 and 100 h, respectively. This makes RTM ∼50% more expensive to run, compared with KDM. Even so, considering the uplift in imaging steep reflectors and faults, we are encouraged to regularly run RTM on our 3D data.
The velocity model building process is what makes depth imaging a significantly longer to run technique compared F I G U R E 1 3 Comparison of imaging results in the north region of the survey. Inline A-A' section has no vertical exaggeration; depth slice is at 1.5 km from surface. Legacy PSTM is scaled to depth using smoothed velocity from depth imaging. KDM and RTM show improved imaging of steep dips particularly in near surface and better-resolved reflectors in depth slice. PSTM suffers from migration artefacts (migration swings) and unreliable near surface. Highlighted in white arrows are east-dipping reflectors indicating a good correlation with monzonite models measured from boreholes with time imaging. It included 10 iterations of refraction tomography and four iterations of sparse KDM each followed by reflection tomography. It took us two months to achieve the final velocity model for depth imaging, twice as long as the time processing and fast-track PSTM. However, the imaging uplift from both KDM and RTM migrations, with respect to precision of reflector positioning in 3D space is of essence in the mining industry, therefore justifying the extended imaging timeline.

CONCLUSIONS
A tailored approach to the application of depth imaging techniques applied to seismic data recorded over a complex geo-logical environment in a mineral exploration context is proposed and implemented in this study. Pre-stack depth imaging (Kirchhoff and reverse time migration (RTM)) achieved a high degree of correlation between depth migration solutions compared with drillhole data, magnetics, refraction tomography and VSP. This further highlights the importance and role of depth imaging for the delineation of complex hard rock environments. The proposed workflow for velocity model building was of key importance for the successful application of depth imaging techniques and ultimately in developing an efficient methodology for hard rock data imaging.
This study demonstrated that both PSDM algorithms when coupled with appropriate velocity model building showed significant uplift when compared with the initial time imaging, especially in the critical top 1 km and as shallow as the TOFR. F I G U R E 1 4 Comparison of imaging results across Whirling Dervish and Karari mine. Crossline B-B' section has no vertical exaggeration; depth slice is at 500 m from the surface. Legacy PSTM is scaled to depth using smoothed velocity from depth imaging. KDM and RTM show excellent correlation with existing ore and monzonite surfaces from drilling. Highlighted in white arrows are faults visible on RTM volume in the Karari basin F I G U R E 1 5 Comparison of KDM (a) and RTM (b) stacked volume in the area of the Karari deposit. Depth slice at 300 m below the surface corresponds to refraction tomography. Both volumes show a good correlation against velocity anomalies on refraction tomography. White arrows on RTM volume highlight the contact of steep reflectors with higher velocities on the depth slice F I G U R E 1 6 Observation of depth slice and comparison with other data sets. Red line demarcates high reflectivity mafic and possible domain bounding Keith Kilkenny thrust. Overall, a good correlation is observed between different data sets with respect to the separation of lithologies and major structures. Rock description in the geological map is given in Figure 1 A specific velocity model building technique, that involves drillhole data integration and tomographic refinement proved effective for hard rock seismic imaging producing accurate images even in the near surface. Excellent correlation between resultant seismic images, refraction tomography, VSP and drillholes can be attributed to the appropriate velocity model used for depth migration. Although the RTM method is computationally more expensive than Kirchhoff's depth migration, it is preferred when detailed subsurface imaging is required.
This study sets a foundation for the future application of depth imaging approaches for the exploration of mineral resources. While we still work on improving the performance of our imaging approach applied to a hard rock environment, it is becoming clear that depth imaging will play a crucial role in the mineral sector due to its accuracy and resolution. Specific challenges associated with each new hard rock environment mean that some site-specific optimization may still be needed for improving the performance of depth imaging. Hence, at this stage, each depth imaging project should be considered as a unique study to extract maximum value from the seismic data set.
The proposed depth imaging methods led to improved ore identification which enabled the interment of an improved geological model that will guide drilling for future mine expansion.

F I G U R E 1 7
Overlay of deviated WDEX047 drillhole with RTM W-E section (a) and embedded with a VSP section (b). The shape of the drillhole is the acoustic impedance (AI) and colour corresponds to lithology. Reflections on the seismic image are expected where we have a change in acoustic impedance. We can observe that reflectors broadly correlate to AI changes and in particular with monzonite intrusions (red surface). In addition, a good correlation exists between RTM and VSP section, where VSP exhibits a higher resolution

A C K N O W L E D G E M E N T S
We thank Northern Star Limited for permission to publish the work and for their support throughout the project. We acknowledge the support from Hisies team, especially from Ian James, Alexey Artemov, Reece Cunnold and Sana Zulic during various stages of the project, since the seismic feasibililty study began in 2017.

D A T A AVA I L A B I L I T Y S T A T E M E N T
Research data are not shared.