G R A V I T Y M O D E L L I N G I N T H E W E S T E R N B U S H V E L D C O M P L E X, S O U T H A F R I C A, U S I N G I N T E G R A T E D G E O P H Y S I C A L D A T A Stephen John Coomber A dissertation submitted to the Faculty of Science, University of the Witwatersrand, Johannesburg, in fulfilment of the requirements for the degree of Master of Science. Johannesburg, 2008 DECLARATION I declare that this dissertation is my own, unaided work. It is being submitted for the degree of Master of Science in the University of the Witwatersrand, Johannesburg. It has not been submitted before for any degree or examination in any other University. ________________________________ Stephen John Coomber _______ day of _________________________ 2008 ii ABSTRACT A 10 km x 10 km study area in the western Bushveld Complex, south of the Pilanesberg Complex, was selected for testing the inversion of vertical component gravity (Gz) data to determine the geometry of the Bushveld Complex/Transvaal Supergroup contact. This contact has a density contrast of ~0.350 g.cm-3 making it a suitable target for gravity inversion. The resulting 3D gravity model agrees well with the 3D seismic interpretation, indicating that the depths determined from the seismic data are appropriate. The gravity inversion could be extended laterally to investigate regions without seismic data coverage. This methodology may prove useful where upwellings in the floor of the Bushveld Complex distort seismic data, but can be imaged by gravity inversions. The Gz dataset was created from converted Airborne Gradient Gravity (AGG) data, combined with upward continued ground Gz gravity data, providing extensive coverage. This combined dataset was used in an interactive, iterative 3D gravity inversion methodology used to model the geometry of the Bushveld Complex/Transvaal Supergroup contact and densities of the Bushveld Complex, Transvaal Supergroup and Iron-Rich Ultramafic Pegmatoids (IRUPs). The resulting 3D gravity model provides an acceptable first-pass model of the Bushveld Complex/Transvaal Supergroup contact. In the shallow south-west region of the study area, the steeply dipping contact was determined from borehole intersections. 3D seismic data was the only constraint towards the north-east, where the contact flattens out to a sub-parallel contact, at ~2 000 m depth. In the north-western section, the Bushveld Complex/Transvaal Supergroup contact is fault-bounded by a conjugate set of the Rustenburg Fault, causing the Bushveld to onlap the Transvaal sediments. In the southern region, the contact changes as the conjugate fault dies out, and the Bushveld Complex becomes layered/sub-parallel to Transvaal sediments. This, and other geological features (e.g. faulting, folding, dykes), can be explained in relation to the regional tectonic history, relating to motion along the Thabazimbi-Murchison iii Lineament (TML). Pre-Bushveld emplacement NW-SE far-field stress caused NW trending extensional features in the region (e.g. Rustenburg Fault). Re-orientation of the compressive force to NE-SW, in syn- to post-emplacement, caused compressive features in the region (e.g. open folds with axes trending NW). Ground gravity data (100 m x 100 m station- and line-spacing) were also inverted to obtain a 3D model of the overburden, constrained by borehole data. However, the inversion failed to satisfy the gravity data and borehole data simultaneously, relating to difficulties in modelling the regional gravity field and the gradational nature of the weathered contact. Several rapid variations in overburden thickness were mapped, with particular success in the high frequency ground gravity survey (30 m x 30 m station- and line-spacing) with the identification of a deeply weathered (~10 m deep) channel relating to an mapped fault. iv Dedicated to my Mom and Dad, for all the love, support and friendship v ACKNOWLEDGEMENTS As in any major academic work, this would have not have been possible without the support of many incredible people and organisations. My supervisor, Sue Webb, gave ever-positive advice, insight and support. Her knowledge of the field of geophysics is dumb-founding and awe-inspiring. Plus, she let me row as much as I wanted. Gordon Chunnett (Anglo Platinum) provided incredible financial support and permission to work on the project. Dr. Andy Rompel (Anglo Technical Division) provided access to technical support and allowed me to finalise my Masters during working hours, which all allowed for a considerably less stressed two years of study. Many technical experts gave their advice and above-and-beyond-the-call-of-duty support for the project, especially Kevin Fisher and Nick Parker (Geosoft), Peter Fullagar and Glenn Pears (VPmg), Jennifer Levett (Gocad), Des Fitzgerald and Dominik Argast (Intrepid) and Paul Wouters (Anglo Platinum). Tim and Bronwyn Chalke also provided a clear understanding on the structure of an MSc and many crash-courses in various computer programs. The benefactors of the SEG?s Lucien Lacoste Scholarship supplied a most unexpected and immensely fulfilling grant. Stephanie Scheiber, Shawn Letts and Mark Hamilton provided incredible friendship and inspirational insight and ideas into the world of geophysics. Guys, thanks for all the adventures, all the laughs and Steph, especially, thanks for all the tea. Finally, to the three people who made the last two years count. To Mom and Dad, there are no words for how much I appreciate and love you guys. To Michelle, I?ve got you under my skin? vi TABLE OF CONTENTS DECLARATION................................................................................................... ii ABSTRACT........................................................................................................... iii ACKNOWLEDGEMENTS................................................................................... vi LIST OF FIGURES................................................................................................ xi LIST OF TABLES.................................................................................................. xx CHAPTER 1:? .................................................................... 1?INTRODUCTION 1.1? ............................................................................................. 1?Geology 1.2? ........................................................................................ 3?Geophysics 1.2.1? ....................................................................... 3?Geophysical Methods 1.2.2? ..................................................................... 5?Geophysical Modelling 1.3? .................................................................. 6?Goals of the Dissertation 1.4? ....................................................................... 7?Dissertation Structure CHAPTER 2:? ........................................ 9?METHOD ? GROUND GRAVITY 2.1? ....................................................................................... 9?Introduction 2.1.1? .................................................................................. 9?Gravity Theory 2.1.2? .................................................................... 10?Gravitational Potential 2.1.3? ............................................... 11?Reductions to Ground Gravity Data 2.1.4? ..................................................... 12?Factors Affecting Rock Density 2.1.5? ..................................................................... 13?Gravity Measurement 2.2? ......................................... 14?Scintrex CG-3 Autograv Gravity Meter 2.2.1? ............................................................ 16?In-Field Setup & Operation 2.3? .......................... 17?Differential Global Positioning Systems (DGPS) 2.3.1? ............................................................................................. 18?Theory 2.3.2? ....................................................................................... 20?Instrument 2.3.3? ............................................................................... 20?Quality Control 2.3.4? .......................................................................... 21?Downloading Data 2.3.5? ................................................................................................. 21?Data 2.4? ................................................................................. 21?Survey Design 2.5? ............................................................................... 23?Quality Control 2.6? ............................................................................ 24?Office Processing 2.6.1? .......................................................................... 24?Downloading Data 2.6.2? .............................................................................. 24?Drift Correction 2.6.3? ......................................................................... 25?Latitude Correction 2.6.4? .......................................................................... 26?free-air Correction 2.6.5? ........................................................................ 27?Bouguer Correction 2.6.6? .......................................................................... 27?Terrain Correction 2.6.7? .......................................................................... 28?Gravity Anomalies 2.7? ......................................................... 29?Final Ground Gravity Dataset CHAPTER 3:? ................................................................................ 30? METHOD ? AIRBORNE FULL TENSOR GRADIENT GRAVITY 3.1? ..................................................................................... 30?Introduction vii 3.1.1? ....................................... 30?History of Airborne Gravity Instruments 3.1.2? .................................................................................... 32?FTG Theory 3.2? ...................................................................................... 36?FALCON? 3.3? ........................................................................................ 38?Air-FTG? 3.4? ............................................................... 42?Air-FTG Survey Design? 3.4.1? ................................................................................... 42?Flight Height 3.4.2? ............................................................................. 43?Production Lines 3.4.3? ......................................................................................... 44?Tie Lines 3.4.4? ............................................................................... 45?Flight Direction 3.4.5? .............................................................. 45?Styldrift Flight Parameters 3.5? ...................................................................... 45?Air-FTG Processing? 3.5.1? ...................................................................... 46?High Rate Processing 3.5.2? .......................................................................... 46?Terrain Correction 3.5.3? ................................................................ 46?Spike and Shift Removal 3.5.4? ......................................................................................... 46?Levelling 3.5.5? ................................................................... 47?Full Tensor Processing 3.6? ....................................................... 47?Airborne FTG Gravity Dataset CHAPTER 4:? ................................. 48?METHOD ? GROUND MAGNETICS 4.1? ..................................................................................... 48?Introduction 4.1.1? ............................................................. 48?Local Magnetic Anomalies 4.1.2? ..................................... 49?Magnetic Properties of Rocks & Minerals 4.1.3? .............................................. 50?Induced & Remanent Magnetisation 4.1.4? .......................................... 52?Diurnal Variations & Magnetic Storms 4.1.5? ..................................................................... 54?Aeromagnetic Survey 4.2? ..................................................................................... 56?Instruments 4.2.1? ............. 56?Geometrics G856AX Proton Precession MagnetometerTM 4.2.2? ...................................................................... 57?Garmin? GPS 12 XL 4.3? ................................................................................. 58?Survey Design 4.4? ............................................................................... 59?Quality Control 4.5? ....................................................................................... 60?Processing 4.5.1? .......................................................................... 60?Downloading Data 4.5.2? .............................. 60?Corrections & Interpolation of Field Readings CHAPTER 5:? ....... 62?DATA ENHANCEMENT FOR INTERPRETATION 5.1? ..................................................................................... 62?Introduction 5.1.1? ......................................................................... 62?Geophysical Filters 5.1.2? ........................................................ 64?Processing Software Packages 5.2? ...................................................................... 64?Vertical Continuation 5.2.1? ............... 65?Lowering Noise Levels during Downward Continuation 5.3? ..................................................... 67?Integrating FTG data to G dataz 5.4? ......................................... 68?Merging Airborne and Ground Datasets 5.4.1? .............................................................................. 69?Final G Datasetz 5.5? .................................................................................... 70?Sun-shading 5.6? ....................................................................... 71?Euler Deconvolution 5.6.1? ............................................................................................. 71?Theory 5.6.2? .............................................................................. 72?Structural Index viii 5.6.3? ............................................... 73?Application of Euler Deconvolution 5.6.4? ................................ 74?Euler Deconvolution Applied to Profile Data CHAPTER 6:? ............................................... 76?GEOPHYSICAL DATASETS 6.1? .............................................................................. 76?Geological Map 6.2? .................................................................... 78?Regional Gravity Data 6.3? .................................................................. 79?Ground Gravity Dataset 6.3.1? .................................................................... 81?Overburden Thickness 6.3.2? ............................................................................ 82?Gravity Traverses 6.3.3? ..................................... 83?Gravity Traverses ? Vertical Continuation 6.3.4? ............................................... 86?Upward Continued Ground G Dataz 6.4? ....................................................... 87?Airborne FTG Gravity Dataset 6.4.1? ................................................... 89?Downward Continued AGG Data 6.4.2? ...................................................................... 90?Integrated AGG Data 6.4.3? ................................................... 90?Error in Integrated Airborne Data 6.4.4? .............................................................................. 92?Final G Datasetz 6.5? ..................................................................... 94?Aeromagnetic Dataset 6.5.1? ................................................................................. 96?Cultural Noise 6.5.2? ..................................................... 97?Sun-shaded Aeromagnetic Data 6.5.3? ................ 98?Euler Deconvolution Applied to Aeromagnetic Dataset 6.6? ............................................................. 101?Ground Magnetic Dataset 6.6.1? ....................................................................... 102?Magnetic Traverses 6.6.2? ................................ 103?Magnetic Traverses ? Vertical Continuation 6.7? ......................................................................... 105?Topographic Data 6.8? ................................................................................. 106?Seismic Data 6.9? ............................................................................... 108?Borehole Data 6.10? ....................................................................... 109?Rock Property Data CHAPTER 7:? ........................................ 111?GEOPHYSICAL MODELLING 7.1? ........................................................................................... 111?Theory 7.2? ................................. 113?3-Dimensional Gravity Forward Modelling 7.2.1? ....................................................................... 115?Rectangular Prisms 7.3? ................................................. 117?3-Dimensional Gravity Inversion 7.3.1? ............................................................................ 118?Linear Inversion 7.3.2? .................................................................... 118?Non-Linear Inversion 7.3.3? ........................................... 121?3-Dimensional Inversion Algorithms 7.4? ...................................................... 122?VPmg Inversion Methodology 7.4.1? ....................................................................................... 122?Overview 7.4.2? ................................................................ 123?Model Parameterisation 7.4.3? ...................................................................... 126?Inversion Algorithm 7.4.4? ............................................................... 128?Geometrical Constraints 7.4.5? ....................................................... 131?Physical Property Constraints 7.5? .................................................................................... 132?Conclusion CHAPTER 8:? ...................... 133?APPLIED GEOPHYSICAL MODELLING 8.1? .......................................................... 133?2.5-Dimensional Modelling 8.1.1? ........................................................................ 134?Gravity Modelling 8.1.2? ...................................................................... 137?Magnetic Modelling ix 8.2? ..................................................... 139?3-Dimensional Starting Model 8.2.1? ....................................................................... 144?Forward Modelling 8.3? ................................................ 147?VPmg Inversion: Regional Model 8.3.1? ............................. 149?Inversion 1: Homogeneous Density Inversion 8.3.2? ? ....................................................................................... 152? Inversion 2: Heterogeneous Transvaal Supergroup Density Inversion 8.3.3? ................. 156?Inversion 3: Transvaal Supergroup Contact Inversion 8.3.4? ............ 157?Inversion 4: Heterogeneous Bushveld Density Inversion 8.3.5? .............. 161?Inversions 5-7: Homogeneous IRUP Density Inversion 8.3.6? ........................................................... 166?Summary of the Inversions 8.4? ............................................ 168?VPmg Inversion: Overburden Model 8.4.1? ..................................................... 168?3-Dimensional Starting Model 8.4.2? ... ...................................................................................................... 169? Inversion 1: Heterogeneous Bushveld Complex Density Inversion 8.4.3? ..... 172?Inversion 2: Bushveld Complex Contact Geometry Inversion 8.4.4? .................................................... 176?Inverted Overburden Thickness 8.5? .................................................................................... 177?Conclusion CHAPTER 9:? .............................. 179?INTERPRETATION & DISCUSSION 9.1? ............................................................................ 179?Tectonic Setting 9.1.1? .......... 179?Structural Deformation in the western Bushveld Complex 9.1.2? ..... 182?Geological History of the Thabazimbi-Murchison Lineament 9.1.3? .......................................... 185?Intrusion of the Pilanesberg Complex 9.2? .............................. 186?Interpretation of Regional Modelled Features 9.2.1? ................. 186?Contact of Bushveld Complex/Transvaal Supergroup 9.2.2? ............................................................................................. 189?Folds 9.2.3? ............................................................................................. 193?Faults 9.2.4? ................................................................ 195?Magnetic Interpretation 9.3? ................................... 198?Interpretation of Local Modelled Features 9.3.1? .................................................... 198?Ground Gravity Survey Results 9.3.2? ........................ 199?High Resolution Ground Gravity Survey Results 9.4? ........................................... 201?Gravity Modelling as a Seismic Tool 9.5? ................................................................................... 204?Conclusions CHAPTER 10:? ................................................................................... 205? CONCLUSION AND SUGGESTIONS FOR FURTHER WORK APPENDIX A .................................................................................................... 208 APPENDIX B .................................................................................................... 209 REFERENCES ..................................................................................................210 x LIST OF FIGURES Figure 1.1 ? Simplified geological map of the Bushveld Complex and surrounding Transvaal Supergroup, including the Rustenburg Layered Suite (mafic phase), Pilanesberg Complex alkaline intrusion, the Thabazimbi-Murchison Lineament (TML), and the location of the study area (after Cairncross and Dixon (1995)). . 1? Figure 1.2 ? (right) Regional Bouguer gravity data over the Bushveld Complex, with the northern, western and eastern limbs clearly visible. (left) Zooming in on the western Bushveld Complex, with the outline of the study region (solid box), approximate Transvaal Supergroup contact (dashed line) and approximate outline of the Pilanesberg Complex (dotted line). Data courtesy of the South African Council for Geoscience............................................................................ 4? Figure 2.1 ? Equipotential surfaces caused by different substratum. (a) Uniform substratum causes equal gravitational field vectors across the equipotential surface. (b) An ore body of arbitrary shape and with density greater than the substratum (i.e. ? > ? ) results in a mass excess causing unequal field lines across an up-warped equipotential surface (after Corner, 1993). 2 1 ....................... 11? Figure 2.2 ? Overhead view of the Scintrex Autograv CG-3, showing Control Console, Mechanical Tilt Meters (X and Y axis) and Data Output serial port (Scintrex, 1995)................................................................................................... 15? Figure 2.3 ? Schematic diagram of the orbital paths of the 24-satellite constellation of the GPS system, orbiting the Earth at an altitude of approximately 20 200 km. Note: diagram not to scale (Garmin, 2007)......................................................... 17? Figure 2.4 ? Ground gravity traverses (crosses) and grids (dots) relative to outline of AGG survey area................................................................................................. 22? Figure 3.1 ? a) Schematic co-ordinate system showing the three principle components of the gravity field (G ? blue, G ? green, G ? red), with derivatives in the x, y and z direction of each component creating the nine components of the gravity tensor. b) The components of the FTG gravity field, which describe the rate of change of the gravity vector as one moves in the direction of the vector. c) The FTG matrix, T , showing the components which are equal (i.e. T = T , T = T , and T = T ), hence one of each of these components can be considered redundant (Bell Geospace, 2006). x y z ij xy yx xz zx yz zy ....................................................................... 33? Figure 3.2 ? Total gravity field (G ) and FTG gravity components for a buried synthetic cube (outline of cube presented. Parameters: dimensions 1 000 m x 1 000 m x 1 000 m; depth 25 m; density contrast 1.0 g.cm ). z -3 ............................ 35? Figure 3.3 ? Cessna Grand Caravan, modified to house the FALCON system for airborne surveys (BHP Billiton, 2007). ? .............................................................. 38? Figure 3.4 ? a) Schematic diagram of the Air-FTG instrument, showing the arrangement of the accelerometers within each spinning disk, b) image of the Air-FTG system, showing positioning of CGIs, and c) geometric arrangement of the GGIs (Murphy, 2004). ? ? ..............................................................................40? Figure 3.5 ? Comparing surveys flown at a constant height or draped, as well as the actual height flown for a draped survey in regions of steep topography (Syberg, 1972). .................................................................................................................. 43? xi Figure 3.6 ? The effect of line-spacing on T data, where lines have been removed from an original survey. a) Original survey using 50 m line-spacing, b) simulated survey using 150 m line-spacing, and c) simulated survey using 250 m line-spacing. The target anomaly (length: 350 m) is still resolved in b) but already shows aliasing in c). The negative anomaly (circled) all but disappears in c), (Murphy, 2004). zz ......................................................................................... 44? Figure 4.1 ? South-North schematic profile of a magnetic anomaly over a magnetic dyke with a dipolar field. The dyke is striking E-W, and situated in the southern hemisphere (Roux, 1980). ................................................................................... 49? Figure 4.2 ? The effect of remanence on S-N profile: the induced magnetization vector and remanent magnetization vector combine to produce a resultant vector, commonly different to the present-day field (Roux, 1980)................................. 51? Figure 4.3 ? Dynamos in the upper atmosphere generating magnetic fields which causes the diurnal (USGS, 2005). ....................................................................... 52? Figure 4.4 ? Schematic diurnal variation, showing the vertical component of the magnetic field through a 24 hour period at varying magnetic latitudes (Roux, 1980). .................................................................................................................. 53? Figure 4.5 ? Schematic magnetic storm, showing the vertical component of the magnetic field. Note the erratic nature of the field and the change in scale from Figure 4.4. ........................................................................................................... 54? Figure 4.6 ? JetRanger Bell 206 IIIB, with horizontally mounted magnetometers (Fugro, 2005). ..................................................................................................... 55? Figure 4.7 ? Circuit diagram of a proton precession magnetometer (Telford et al., 1990). .................................................................................................................. 56? Figure 4.8 ? Ground magnetic traverses, relative to the outline of the AGG data...... 58? Figure 4.9 ? Plot of base station data for the two days spent on the ground magnetic survey. Note: surveys were approximately one year apart, accounting for the ~100 nT shift from Day 1 to Day 2. This shift was corrected for in the survey data. ..................................................................................................................... 59? Figure 5.1 ? Schematic example of a spatial 3 x 3 smoothing filter kernel. Note, the 6 x 4 input grid is reduced to a 4 x 2 filtered grid, after Cooper (2004c). ............... 1? Figure 5.2 ? The suture method, with circle (radius: R) traversing the boundary of the ground gravity survey, interpolating between the ground gravity and calculated airborne G data. Note: this figure does not show the final sutured image.z ...... 69? Figure 5.3 ? Euler solutions over a synthetic dyke S-N profile using computer program Euler (Cooper, 2004b). Upper panel: Original data (solid black line) and Reduction-to-Pole (RTP, dashed red line) calculated from original data. Middle panel: Horizontal gradient (solid black line) and vertical gradient (dashed red line) data calculated from original data. Lower panel: Euler solutions for a structural index of 1 (i.e. green dots), with overlying position of synthetic dykes. Note the spread of the solutions, both horizontally and vertically.............................................................................................................. 75? Figure 6.1 ? Geological map (after Walraven (1981)) showing ground gravity survey (grid and traverses), AGG survey and aeromagnetic survey. ............................. 77? xii Figure 6.2 ? Regional Bouguer gravity data over the Bushveld Complex. (right) Northern, western and eastern limbs of the Bushveld Complex are clearly visible as gravity high anomalies. (left) Zooming in on the western limb of the Bushveld Complex, showing the outline of the AGG survey area (solid box), and the approximate position of the Transvaal Supergroup contact (dashed line) and Pilanesberg Complex (dotted line). Data courtesy of the South African Council for Geoscience..................................................................................................... 78? Figure 6.3 ? a) free-air gravity map of the ground survey grid, with 100 m x 100 m grid and 30 m x 30 m grid in the eastern corner. b) Bouguer gravity map of the ground survey grid, with similar grid positions. c) DTM of the ground gravity survey region....................................................................................................... 80? Figure 6.4 ? Thickness of the overburden approximated from the DTM and borehole data. Borehole positions represented by crosses. ............................................... 82? Figure 6.5 ? Bouguer gravity profiles (station spacing: 150 m) for Line 1 (solid black) and Line 2 (solid red), and corresponding elevations for Line 1 (dashed black) and Line 2 (dashed red), running from NE to SW. The gravity anomaly in Line 2 is related to a traversed IRUP. ......................................................................... 83? Figure 6.6 ? Line 1: Profiles of free-air ground gravity data (black), airborne gravity data flown at 80 m (red), and airborne gravity data downward continued 80 m (blue). .................................................................................................................. 84? Figure 6.7 ? Line 2: Profiles of free-air ground gravity data (black), airborne gravity data flown at 80 m (red), and airborne gravity data downward continued 80 m (blue). The anomaly in the downward continued profile in the NE is a result of edge effects from FFTs in the downward continuation. ..................................... 85? Figure 6.8 ? a) Original ground free-air G data, upward continued to b) 20 m, c) 40 m and d) 80 m (units: mgal). z ............................................................................... 86? Figure 6.9 ? Final levelled airborne FTG gravity grid, showing T data and xy- recording positions (showing lines and tie lines). Data courtesy of Anglo Platinum. zz ............................................................................................................. 87? Figure 6.10 ? The individual tensor components of the AGG data (units: E?tv?s).... 88? Figure 6.11 ? Original AGG T data, flown at 80 m above ground level downward continued to b) 70 m, c) 60 m, and d) 40 m (units: E?tv?s). zz .............................. 89? Figure 6.12 ? a) The integrated airborne data (units: arbitrary), and b) the integrated airborne data, divided by 10 000 and with a DC shift applied to fit the upward continued ground data (units: mgal). Note the position of the ground gravity grid is outlined in b). ........................................................................................... 90? Figure 6.13 ? a) Difference between the upward continued ground data and the shifted airborne G data, and b) histogram of the grid provided in a). Note the trend in the difference not resolved by the DC shift using least-squares. z ........... 91? Figure 6.14 ? a) Final merged G dataset, 80 m height (units: mgal), and b) zooming in to the ground gravity grid section, with sun-shading from the NE applied. Note how much smoother the upward continued data is relative to the integrated AGG data. z ........................................................................................................... 93? Figure 6.15 ? (top) Regional aeromagnetic map of the Pilanesberg region, data courtesy of the South African Council for Geoscience. (bottom) Aeromagnetic xiii map of the study region, sun-shaded from NE. Data courtesy of Anglo Platinum. ............................................................................................................. 95? Figure 6.16 ? Sun-shaded colour aeromagnetic maps zooming in on areas overlaying QuickBird orthophotographs. a) Village causing the dappled features in the magnetic data, b) Power lines and a dump causing other man-made noise in the aeromagnetic data................................................................................................ 97? Figure 6.17 ? a) Grey-scale aeromagnetic data, b) the same data, sun-shaded from inclination 45?, azimuth 45?, highlighting features perpendicular to 45?, and c) sun-shaded from inclination 90?, highlighting the edges of all magnetic sources. ............................................................................................................................. 98? Figure 6.18 ? Euler solutions (structural index: 1) overlaying aeromagnetic data (sun- shaded from 90? inclination) with solutions for varying depth ranges: green (5- 15 m), blue (15-25 m), red (25-35 m) and black (35-45 m). Inset: Zoom of the solutions over two prominent dykes, with depths dominated by 5-15 m solutions. ........................................................................................................................... 100? Figure 6.19 ? Ground magnetic profiles locations (black), with associated ground profile data (red), and gravity stations (black crosses), overlaying aeromagnetic data in the study region. .................................................................................... 101? Figure 6.20 ? Magnetic profiles of Line 1 (black) and Line 2 (red), running SW to NE. A 100 nT shift was applied to Line 2, allowing for easier viewing. The significant anomalies in the SW relate to the E-W trending structure.............. 102? Figure 6.21 ? Line 1: Profiles of ground magnetic data (black), airborne magnetic data downward continued 18 m, offset 150 nT (blue), airborne magnetic data flown at 20 m (red), offset 300 nT. ................................................................... 104? Figure 6.22 ? Line 2: Profiles of ground magnetic data (black), airborne magnetic data downward continued 18 m, offset 150 nT (blue), and airborne magnetic data flown at 20 m, offset 300 nT (red). Note: change of y-axis scale from Figure 6.21. ....................................................................................................... 104? Figure 6.23 ? DTM of the study region, determined from aeromagnetic survey (data courtesy of Anglo Platinum). Outline of AGG survey, ground gravity stations and geological boundaries included. ................................................................. 105? Figure 6.24 ? Final G gravity map showing position of seismic NE-SW inline profiles and perpendicular cross profiles. The inline profile in bold is shown (left) with the trace of the Transvaal Supergroup contact (yellow). Data courtesy Anglo Platinum. z ................................................................................................ 107? Figure 6.25 ? a) Converted G free-air gravity map with borehole positions and outline of the ground gravity survey (borehole data courtesy of Anglo Platinum and PTM), b) 3D view of Transvaal Supergroup contact (interpolated from seismic profiles) and boreholes in the study region, with only 27 shallow boreholes intersecting the contact in the south-western region. z ....................... 109? Figure 7.1 ? Comparing the forward and inverse method. A represents the measured anomaly, and A the calculated anomaly. Parameters p , p ,? are the source attributes (e.g. depth, thickness, density), Blakely (1995). 0 1 2 ............................... 112? xiv Figure 7.2 ? Arbitrary 3D body with density ?(x?,y?,z?) and unit vector r ........................................................................................................................... 114? ? pointing from an element of the body to the calculation point P(x,y,z), (Blakely, 1995). Figure 7.3 ? Schematic diagram of a 3D approximation of a body, using rectangular prisms (Blakely, 1995). ..................................................................................... 116? Figure 7.4 ? Cross-section model of a sedimentary basin. Basin is approximated by rectangular blocks, extending infinitely perpendicular to the profile and with one block per gravity station (Bott, 1960). .............................................................. 119? Figure 7.5 ? 3D model, defined by Cordell and Henderson (1968), for iterative inversions. Block thicknesses are relative to a reference surface and the calculated gravity data occur on an observation surface................................... 120? Figure 7.6 ? Comparing geometry variation methods for two sets of similar geological units (lithological boundaries represented by red and blue lines). a) Conventional fixed mesh (black lines), with invariable horizontal cell divisions and cells assigned according to the lithology occupying the greatest percentage of the cell. b) Adaptive mesh with variable horizontal cell divisions (solid horizontal lines) and sub-cells (dotted horizontal lines), after Fullagar and Pears (2007). Note how the adaptive mesh represents the geology more closely. .... 125? Figure 7.7 ? a) Gently dipping geological boundaries with pierce points near the prism centre, boundary B is fixed and boundary A can vary, b) Steeply dipping geological boundary with pierce point near the prism edge, contact A is within range of pierce point P (after Fullagar and Pears (2007))................................. 129? Figure 7.8 ? Radius of influences, with R defined by the depth of the pierce point below the topography and R defined by the distance to nearest pierce point. Geometry changes within the radius of influence are damped during inversion (after Fullagar and Pears (2007)). 1 2 ..................................................................... 131? Figure 8.1 ? 2.5D gravity model of Line 1, with measured data (green line) and calculated data (black line). Inversions produced a density contrast of 0.333 g.cm .-3 ................................................................................................................ 135? Figure 8.2 ? 2.5D gravity model of Line 2, with measured data (green line) and calculated data (black line). Inversions produced a density contrast of 0.334 g.cm . An IRUP was modelled with an inverted density contrast of 0.338 g.cm . -3 - 3 ........................................................................................................................ 136? Figure 8.3 ? Profiles of a) Line 1 across Dyke 1, b) Line 2 across Dyke 1, and c) Line 1 across Dyke 2. ................................................................................................ 138? Figure 8.4 ? Simplified geological map (after Walraven (1981)) of the study area, within the AGG survey area (i.e. the white box). The Transvaal Supergroup contact was simplified to a single contact (dashed white line), separating the Bushveld Complex (green) and Transvaal Supergroup (vertical white lines). All other structure was removed. ............................................................................ 141? Figure 8.5 ? Plan view of the study area, showing the Transvaal Supergroup contact with IRUPs labelled 1-14. ................................................................................. 142? Figure 8.6 ? 3D geological starting model of the study region, showing topographic surface (i.e. surface separating air from the Bushveld Complex and Transvaal Supergroup), IRUPs, and Transvaal Supergroup contact. Dimensions of the xv bounding box are shown and apply to all subsequent images of the regional model. Vertical exaggeration: 3. ..................................................................... re 8.7 ? a) 3D voxet of the starting model, with the transparent Transvaal . 143? Figu ergroup .......... 144? Figu . 5? Figu 46? Figu e 48? Figu ? Figu 1? Figu .153? Figu ? Figu 5? Figu up .157? Figu blue Figu ? Figu 0? Figu l 3? Supergroup contact and cross-sections showing the prisms defining the geophysical units: air (blue), Bushveld Complex (orange), Transvaal Sup (lime) and IRUPs (red). b) 3D voxet model showing the discretized 200 x 200 m prisms of the Transvaal Supergroup (lime) and IRUPs (red). Vertical Exaggeration: 3. ...................................................................................... re 8.8 ? a) Observed data, b) calculated data from the starting model, and c) residual data, taking the difference between the observed and calculated data Note, the change in mgal range......................................................................... 14 re 8.9 ? Histogram of residual data for the starting model, bin size: 0.25 mgal, with normal Gaussian overlay........................................................................... 1 re 8.10 ? Schematic diagram showing the inversion steps carried out in the project to create a final model. As the features undergoing inversion becom shallower, the RMS misfit increases................................................................. 1 re 8.11 ? a) Observed data, b) calculated data after inversion 1, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 150 re 8.12 ? Histogram of residual data after inversion 1, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.9........ 15 re 8.13 ? Plan view of the density of the inverted heterogeneous Transvaal Supergroup (units: g.cm ). Density properties are constant in each vertical prism, extending from the Transvaal Supergroup contact to half-space. -3 ......... re 8.14 ? a) Observed data, b) calculated data after inversion 2, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 154 re 8.15 ? Histogram of residual data after inversion 2, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.12...... 15 re 8.16 ? 3D model showing the minor differences between the original Transvaal Supergroup contact (yellow) and inverted Transvaal Supergro contact (green). Vertical exaggeration: 3. ........................................................ re 8.17 ? Plan view of the inverted densities of the heterogeneous Bushveld Complex (units: g.cm ). The heterogeneous Transvaal Supergroup appears as it falls below the density range (?<2.870 g.cm ), and the IRUPs appear red as they fall above the density range (? = 3.200 g.cm ). -3 -3 -3 ........................................ 158? re 8.18 ? a) Observed data, b) calculated data after inversion 4, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 159 re 8.19 ? Histogram of residual data after inversion 4, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.15...... 16 re 8.20 ? 3D slice of inverted homogeneous IRUP densities, with background Bushveld Complex density and Transvaal Supergroup (units: g.cm ). Vertica exaggeration: 3. -3 ................................................................................................. 16 xvi Figure 8.21 ? a) Observed data, b) calculated data after inversions 5-7, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 164? Figure 8.22 ? Histogram of residual data after inversions 5-7, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.19. ........................................................................................................................... 165? Figure 8.23 ? a) DTM of the ground gravity survey, b) Plan view of gravity stations (black dots) and boreholes (labelled) used to model the starting overburden surface (yellow)................................................................................................. 169? Figure 8.24 ? Plan view of the inverted heterogeneous Bushveld Complex density, compared to regional heterogeneous Bushveld Complex (units: g.cm ). Regional density variation and sloping Transvaal Supergroup contact are accounted for. -3 .................................................................................................... 170? Figure 8.25 ? a) Observed data, b) calculated data after inversion 1, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 171? Figure 8.26 ? Histogram of residual data after inversion 1, bin size: 0.125 mgal, with normal Gaussian overlay................................................................................... 171? Figure 8.27 ? Plan view of the original Bushveld Complex contact (yellow) compared to the inverted surface (red), and showing borehole constraint points on the surface (black dots). The regional density model is shown in the background. Inset: 3D view of the original (transparent yellow) and inverted (red) Bushveld Complex contact. Vertical exaggeration: 10.................................................... 173? Figure 8.28 ? a) Observed data, b) calculated data after inversion 1, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. ....................................................................................... 174? Figure 8.29 ? Histogram of residual data after inversion 2, bin size: 0.125 mgal, with normal Gaussian overlay................................................................................... 175? Figure 8.30 ? Overburden thickness, using the constrained, inverted bedrock contact and DTM. .......................................................................................................... 177? Figure 9.1 ? Regional map of the western Bushveld Complex showing study area (red box) and the fold axes of open folds, with synclines (S1, S2, S3) and anticlines (A1, A2) labelled (after Walraven (1976)). The dashed anticline axis in the north relates to the Crocodile River Fragment (Ferguson, 1973). .................... 180? Figure 9.2 ? The position of the TML (after Cawthorn and Walraven (1998)), extending from the Molopo Farms Complex in Botswana (considered a satellite intrusion of the Bushveld Complex) to Mpumalanga (east of Potgietersrus), cutting across the northern boundaries of the eastern and western limbs of the Bushveld Complex. ........................................................................................... 182? Figure 9.3 ? (left) Synthetic 2D T response across a Bushveld Complex layered with the Transvaal Supergroup (after Antoine (2004)) using Grav2DC (Cooper, 2003) and SignProc (Cooper, 2000), (right) Gridded T map over the study area, showing the response across a similar geological scenario (black line), from AGG data. zz zz ......................................................................................................... 187? xvii Figure 9.4 ? (left) Synthetic 2D T response from a Bushveld Complex onlapping the Transvaal Supergroup (after Antoine (2004)) using Grav2DC (Cooper, 2003) and SignProc (Cooper, 2000), (right) Gridded T map over the study area, showing the response across a similar geological scenario (black lines), from AGG data. zz zz ......................................................................................................... 188? Figure 9.5 ? Geological map (after Walraven (1981)) of the study area, showing the axis of syncline 3 (S3), NE of the study area, and the inferred axis of anticline 3 (A3), SW of the study area................................................................................ 190? Figure 9.6 ? Schematic cross-section (trending from SW to NE) representing possible deformational history of the Bushveld Complex in the study region. a) Intrusion of the Bushveld Complex (layers originally horizontal, Letts (2007)) into the Transvaal Supergroup (dipping ~10?, Cawthorn, 1998) adjacent to the Rustenburg Fault, b) NE-SW far-field stress causes gentle open folds, forming a syncline-anticline pair (i.e. S3-A3), followed by ~10? tilting and erosion (Cawthorn, 1998 and Bumby, 1998), and c) Present tilted fold, overlaying interpreted seismic profile, with Transvaal Supergroup contact drawn in yellow. ........................................................................................................................... 191? Figure 9.7 ? Geological map (after Walraven (1981)) of the study area, showing short wavelength, NE-trending folds (wavelength: ~5 km), in the Transvaal Supergroup, as seen in the Bushveld Complex/Transvaal Supergroup contact. ........................................................................................................................... 193? Figure 9.8 ? a) Plan view of the inverted Transvaal Supergroup heterogeneous density, delineating the Rustenburg Fault conjugate set (white dashed line), and b) Plan view of the inverted Bushveld Complex heterogeneous density, delineating the E-W-trending fault (white dashed line). The Transvaal Supergroup (? ? 2.620 g.cm ) is imaged in blue along the SW edge as it falls below the colour range for the density scale, while IRUPs (? ? 3.091 g.cm ) occur in red as they fall above the colour range. -3 -3 ............................................... 194? Figure 9.9 ? Vertically sun-shaded aeromagnetic data with lines of Euler solutions delineating dykes. Colours represent the depth of dykes (green = 5-15 m, blue = 15-25 m, red = 25-35 m). Pilanesberg dykes trend NW-SE, with a predominant depth of 5-15m, while the E-W trending dyke shows a depth of 25-35m. Note: the depth of dykes varies within the specified range. A positive magnetic anomaly sub-parallel to the UG2 sub-outcrop is represented by the yellow line. ........................................................................................................................... 196? Figure 9.10 ? Overburden thickness (in colour, transparent) based on boreholes, with warm colours for thick overburden and cool colours for thick overburden. Overlays QuickBird orthophotograph and with plotted dextral faults (pink lines), unspecified faults (black lines) and Pilanesberg dykes (blue lines) based on geological mapping from boreholes and field exposure. .................................. 198? Figure 9.11 ? Residual Bouguer gravity grid (colour) from 30 m x 30 m line- and station-spacing survey, overlaying aeromagnetic data (greyscale). The surrounding outline of the 100 m x 100 m line- and station-spacing ground gravity survey (black line) also plotted. Features mapped in geological field survey and from aeromagnetic data: dextral faults (pink lines), IRUPs (red xviii polygons) and delineated dykes (blue). (inset) Profile across the linear gravity low relating to dextral fault (density values from Mar? et al. (2002) and density inversion results). .............................................................................................. 200? Figure 9.12 ? Comparison of a) dipping Bushveld Complex, and b) proposed flattening Bushveld Complex. The dipping Bushveld Complex shows the economic Critical Zone dipping into the Earth, beyond the depths of current shaft technology (i.e. >4 km). However, the flattening Bushveld Complex shows the dip of the Critical Zone decreasing down-dip, making it shallower than previously expected and within minable range. Note abbreviations: Bushveld Complex (BC), Upper Zone (UZ), Main Zone (MZ), Critical Zone (CZ), Lower Zone (LZ), Transvaal Supergroup (TVL).................................... 202? xix LIST OF TABLES Table 1.1 ? Generalized stratigraphy of the western Bushveld Complex and underlying Pretoria Group (after Eales and Cawthorn (1996), Reczko et al. (1997)). Note abbreviations: Bushveld Complex (BC), Transvaal Supergroup (TVL). ................................................................................................................... 2? Table 2.1 ? Densities of rocks in the study area (Mar? et al., 2002; Telford et al., 1990). .................................................................................................................. 12? Table 2.2 ? Scintrex CG-3 Autograv Setup parameters chosen for field work, October 2006 and May 2007............................................................................................. 16? Table 4.1 ? Magnetic susceptibilities of minerals and rocks in the study region (Telford et al., 1990). Note minor variations between rocks of the Bushveld Complex and Transvaal Supergroup, but these are small in comparison to the variation with intrusives, which have a prominent magnetic signature. ............. 50? Table 5.1 ? Structural indices for simple potential field bodies (after Thompson (1982) and Geosoft (2006))................................................................................. 73? Table 5.2 ? Parameters of synthetic dykes used in Euler deconvolution.................... 74? Table 6.1 ? Statistics of ground gravity survey........................................................... 81? Table 6.2 ? Statistics of histogram for the difference between upward continued ground G data and the equivalent section of the shifted airborne G data.z z ....... 91? Table 6.3 ? Parameters for the four initial solution categories, based on depth and with user-defined spatial uncertainty restrictions. .............................................. 99? Table 6.4 ? Starting densities of the regional model (after Mar? et al. (2002) and Ashwal et al. (2005)). ....................................................................................... 110? Table 7.1 ? Comparing the forward modelling run time of VPmg and Grav3D for a model with the same number of data points, using a 1 Ghz Pentium III with 512 MB RAM, after Fullagar (2007b). .................................................................... 126? Table 8.1 ? Modelling parameters for Line 1 and Line 2 gravity profiles................ 134? Table 8.2 ? Range of parameters for 2.5D profiles presented in Figure 8.3. First row shows range, second row (italics) shows value used in presented model......... 139? Table 8.3 ? Statistics for the histogram of the residual data, starting model (units: mgal). ................................................................................................................ 146? Table 8.4 ? Starting densities of the regional model (after Mar? et al. (2002) and Ashwal et al. (2005)), with inverted densities. ................................................. 149? Table 8.5 ? Statistics for the histogram of the residual data after inversion 1 (units: mgal). ................................................................................................................ 151? Table 8.6 ? Statistics for the histogram of the residual data after inversion 2 (units: mgal). ................................................................................................................ 155? Table 8.7 ? Statistics for the histogram of the residual data after inversion 4 (units: mgal). ................................................................................................................ 160? Table 8.8 ? IRUPs inverted densities, including change in density (using the starting density of 3.200 g.cm ), average inverted density and average change in density.-3 ........................................................................................................................... 162? xx xxi Table 8.9 ? Statistics for the histogram of the residual data after inversions 5-7 (units: mgal). ................................................................................................................ 165? Table 8.10 ? Inversion parameters and results for the regional model. Note, the Number of Input Iterations are given, number of iterations successfully run are given in brackets. .............................................................................................. 167? Table 8.11 ? Statistics for the histogram of the residual data after inversion 1 (units: mgal). ................................................................................................................ 172? Table 8.12 ? Statistics for the histogram of the residual data after inversion 2 (units: mgal). ................................................................................................................ 175? Table 9.1 ? Summary of published far-field stresses in the region of the Bushveld Complex, associated with the motion along the TML, including strain ellipsoids and related geological events (after Du Plessis and Walraven (1990)). Note abbreviations: TML (Thabazimbi-Murchison Lineament), BC (Bushveld Complex), RLS (Rustenburg Layered Suite, i.e. Bushveld Complex mafic phase), TVL (Transvaal Supergroup). .............................................................. 183? CHAPTER 1: INTRODUCTION 1.1 Geology The mafic phase of the Bushveld Complex (2 058 ? 0.8 Ma, Buick et al. (2001)) is 7- 9 km thick and covers an area of approximately 65 000 km2, making it the largest known layered igneous complex in the world (Eales and Cawthorn, 1996), (Figure 1.1). It is also host to the world?s largest reserves of Platinum Group Elements (PGEs), chromium and vanadium and a host of other minerals (Viljoen and Reimold, 2002). Seismic Line (Odgers et al., 1993) Figure 1.1 ? Simplified geological map of the Bushveld Complex and surrounding Transvaal Supergroup, including the Rustenburg Layered Suite (mafic phase), Pilanesberg Complex alkaline intrusion, the Thabazimbi-Murchison Lineament (TML), and the location of the study area (after Cairncross and Dixon (1995)). 1 The study region is in the vicinity of the farm Styldrift, south of the Pilanesberg, approximately 25.5?S, 27?E. Geophysical exploration is being carried out in the region as a joint venture between Anglo Platinum, Bafokeng Rasimone Platinum Mine (BRPM) and Platinum Group Metals (PTM). Stratigraphically, the study region coincides with outcrop of the Main Zone of the Rustenburg Layered Suite (Table 1.1). At its thinnest, the Main Zone has been eroded down to a thickness of ~500 m, below which the PGE-bearing Merensky Reef and UG2 chromitite layer are found in the Critical Zone. The contact with the underlying Transvaal Supergroup is a prominent density and velocity contrast. The Transvaal Supergroup outcrops in the south-west of the study area and rapidly dips (~75?) to the north-east down to ~2 000 m depth, after which it flattens out. Similar geometry of the Transvaal Supergroup contact has also been noted in Odgers et al. (1993) and Davidson and Chunnett (1999), indicating that the Bushveld Complex may flatten out towards its centre, to within reach of today?s shaft technology (i.e. 4 km depth). Table 1.1 ? Generalized stratigraphy of the western Bushveld Complex and underlying Pretoria Group (after Eales and Cawthorn (1996), Reczko et al. (1997)). Note abbreviations: Bushveld Complex (BC), Transvaal Supergroup (TVL). Stratigraphic Unit Ave. Thickness (m) Major Rock Types Economic Deposits Upper Zone (BC) 1 700 Gabbro, olivine diorite, anorthosite Vanadium, Iron Main Zone (BC) 3 400 Gabbro, anorthosite, norite, pyroxenite Dimension stone Upper Critical Zone (BC) Anorthosite, norite, pyroxenite, chromitite (Merensky Reef & UG2) PGEs, chromium (minor Au, Ni, Cu) Lower Critical Zone (BC) 1 300 Pyroxenite, harzburgite, dunite Lower Zone (BC) 800 Pyroxenite, harzburgite, dunite Cobalt Marginal Zone (BC) Highly variable Norite Pretoria Group (TVL) 3 000 Shale, quartzite 2 Various other geological features occur in the region, including faulting, folding and dykes. Firstly, the contact with the underlying Transvaal Supergroup is structurally complex, resulting from a conjugate set of the Rustenburg Fault. The Rustenburg Fault itself is a pre-Bushveld Complex feature, trending NW, and showing possible reactivation during the intrusion of the nearby Pilanesberg Complex (~1 300 Ma, Cawthorn (1988)). Syn-Bushveld large-scale folding, with axes trending NW, are also seen in the region. This folding is related to far-field stresses associated with sinistral motion along the Thabazimbi-Murchison Lineament (TML), (Du Plessis and Walraven, 1990; Walraven, 1974; Walraven and Darracott, 1976). Two small-scale features occurring in the region are syn-Bushveld Iron Rich Ultramafic Pegmatoids, i.e. IRUPs, (Reid and Basson, 2002; Scoon and Mitchell, 2004) and two sets of post- Bushveld dykes, trending NW (related to the Pilanesberg intrusion) and E-W (late- stage occurrence). These features are modelled using one or more of the geophysical datasets available (e.g. gravity, magnetic). 1.2 Geophysics 1.2.1 Geophysical Methods The predominant gabbros, anorthosites and norites found in the Main Zone are easily weathered, hence, the study area has ~5-10 m overburden, with little outcrop. This makes geophysics an important tool for determining structure in the area. The study area has a comprehensive suite of overlapping geophysical datasets, including: regional gravity (Figure 1.2), Airborne Gravity Gradiometry (AGG), 3D seismic, aeromagnetic, terrain and borehole datasets. 3 Northern Limb Eastern Limb 100 km Pilanesberg Complex mgal Figure 1.2 ? (right) Regional Bouguer gravity data over the Bushveld Complex, with the northern, western and eastern limbs clearly visible. (left) Zooming in on the western Bushveld Complex, with the outline of the study region (solid box), approximate Transvaal Supergroup contact (dashed line) and approximate outline of the Pilanesberg Complex (dotted line). Data courtesy of the South African Council for Geoscience. The high density mafic rocks of the Bushveld Complex correspond to large amplitude gravity anomalies, while the less dense rocks of the Transvaal Supergroup, alkaline Pilanesberg Complex and cover sequences (e.g. Karoo, Cape and Waterberg Supergroups) appear as gravity lows. Developments in AGG technology have made this method a viable and widely-used option for data collection (Murphy, 2004; Nabighian et al., 2005). Bell Geospace developed the Air-FTG? system, which is capable of successfully measuring the Full Tensor Gradient (FTG) gravity field from an airborne platform (Hammond and Murphy, 2003; Murphy, 2004; Murphy et al., 2006; Zuidweg and Mumaw, 2007). In order to obtain a greater understanding of the underlying geology, Bell Geospace was contracted to perform a 10 km x 10 km airborne FTG survey, which formed the extent of this study region. A large number of ground gravity data points were also collected to constrain near-surface overburden geometry, model the Bushveld 4 Complex/Transvaal Supergroup contact in 2D, and for use in statistically determining the correct mgal-range of airborne Gz data, converted from airborne FTG data. Two 3D seismic surveys have been carried out in the study region. The first covers a region of 6 km x 3 km, while the second covers 10 km x 5 km, with a 5 km x 0.5 km overlap between the two. The AGG survey fully overlaps the seismic surveys. While seismic data are now routinely collected in the Bushveld Complex for mine planning, the deeper portions of these data are rarely interpreted. One enigmatic question in the Bushveld Complex has been the geometry of the contact between the Bushveld Complex and the underlying Transvaal. A research seismic line (Figure 1.1) suggests that this contact flattens out towards the interior (Odgers et al., 1993). The implications of a flattened geometry are two-fold: the mafic rocks of the Bushveld Complex might be connected, as proposed by Cawthorn et al. (1998), and the mafic rocks could be lifted closer to surface in regions of updoming (Campbell, 1990). A constraint on the depth to prominent contacts, such as the Bushveld Complex/Transvaal Supergroup contact, is required for the seismic interpretation to correlate these reflectors. If boreholes do not extend to the depth of the contact, an alternative method must be found to estimate this information. This project aims to demonstrate that gravity modelling can provide these constraints. 1.2.2 Geophysical Modelling Improvements to 3D potential field inversion algorithms render the method increasingly practical. Grav3D (Li and Oldenburg, 1998) is an industry-standard inversion program, allowing rock property inversion in a mesh of rigid cubes. Each cube is assigned a density which varies during the inversion. Rock units may then be interpreted from the inverted densities. VPmg offers an applied approach to the inversion method, allowing the user to build a 3D geological starting model in Gocad. The model is divided into vertical prisms, with geological boundaries forming variable inter-prism geological units. Geometry and property inversions may 5 then be run on the starting model, constrained by a priori knowledge such as borehole intersections and measured density values (Fullagar et al., 2000; Fullagar and Pears, 2007; Fullagar et al., 2004). At the time of undertaking these inversions, FTG inversions were not fully developed in VPmg. Hence, the FTG gravity data was converted to vertical component Gz data (e.g. Hinks et al. (2004), Lane (2004b)), for use in Gz gravity inversions. Additional geophysical datasets (especially, airborne and ground magnetics) were used to model geological features such as faults, dykes and Iron Rich Ultramafic Pegmatoids (IRUPs). These features present a problem to mine planning and operations and need to be accurately mapped and modelled in order to maximise safety and minimise loss-of-ground when mining. 1.3 Goals of the Dissertation Gravity surveys have numerous advantages for modelling contacts with large density contrasts: the method is cheaper and takes considerably less time, compared to seismic surveys. The advantage of an airborne survey is that it is more rapid and cost effective than collecting ground data. The project proposes to test converted airborne Gz data in inverting for suitable contacts, beyond the depth of borehole data, to assist with constraining these contacts in seismic interpretations. This will also test the suitability of running gravity inversions in adjacent regions, as a first-pass modelling process, for the structural geometry of the Bushveld Complex/Transvaal Supergroup contact. AGG data is integrated to free-air Gz data, using upward continued ground gravity data to statistically correlate the correct mgal-range, for use in the inversions. Several seismic profiles over a reflection, proposed to correlate to the Bushveld Complex/Transvaal Supergroup contact, were extracted and interpolated to provide a proposed contact surface. A forward gravity model was calculated using this contact 6 surface, followed by an interactive, iterative inversion process which tests the robustness of the interpreted seismic surface. The Bushveld Complex and Transvaal Supergroup densities are inverted, as well as the geometry of the proposed contact surface, and the densities of IRUPs in the region. An important by-product of this procedure is an inverted 3D geological model, integrating all available geophysical data (i.e. ground and airborne gravity and magnetic data, 3D seismic data, terrain data). The Transvaal Supergroup contact varies from a steeply dipping faulted contact in the SW of the study area, to flattening out, with the Bushveld Complex sub-parallel to the contact in the NE region. 1.4 Dissertation Structure Chapter 2 through to Chapter 5 are dedicated to a literature review of the methods required to collect and enhance data. Chapter 2 covers the ground gravity method, briefly explaining gravity theory, the working of the Scintrex CG-3 gravity meter and Differential GPS theory. Chapter 3 reviews the AGG method, with an explanation of the Air-FTG? system and data processing which provided the FTG gravity dataset. The ground magnetic method is described in Chapter 4, briefly explaining magnetic anomalies and the working of the proton-precession magnetometer. Chapter 5 reviews the data enhancement processes required to convert the AGG data to a free- air Gz dataset which could be used in the 3D gravity inversions. The chapter also explains methods of determining source parameters from magnetic data (including sun-shading and Euler deconvolution). Chapter 6 collates and describes all available geological and geophysical datasets. The data from the ground gravity survey, AGG survey and ground magnetic survey are imaged. An image of the final free-air Gz dataset, converted from the AGG data, is provided, as well as seismic, borehole and rock property data. 7 Chapter 7 and Chapter 8 describe the modelling process carried out in the project. Chapter 7 explains the theory of 3D gravity inversion, discussing the method from basic principles though to the inversion methodology of VPmg. Chapter 8 details the applied modelling processes in the project. These include 2D gravity and magnetic modelling, forward modelling of the 3D geological starting model and interactive, iterative gravity inversions to produce a final model. Chapter 9 uses the modelling results to demonstrate the ability of gravity inversions to model high density contrast contacts at depths. The chapter also presents an integrated geological model, which is consistent with the tectonic history of the region. Chapter 10 concludes the dissertation with a summary of the research and possible future work. 8 CHAPTER 2: METHOD ? GROUND GRAVITY 2.1 Introduction This chapter briefly summarises basic gravity principles, before reviewing the relevant points of operation of the Scintrex CG-3 gravity meter and Trimble Differential GPS (DGPS) that were used in the collection of data. Finally, a discussion of the survey design and methods of gravity data corrections is provided. 2.1.1 Gravity Theory The gravitational force of the Earth is defined by Newton?s Law of Gravitation, which states that two masses (m1 and m2) attract each other with a force of magnitude F r , which is directly proportional to the product of the masses and inversely proportional to the distance between the centres of mass: r r mm GF ?.2 21=r (2.1) where F r is attractive, r? t is the unit vector in a straight line from m1 to m2 and G is the gravitational constan 2g , Telford et al., (1990)). This equation forms the basis for all gravitational theory (especially for two masses with small dimensions relative to the distance between them). ( 211 /10672.6 kmNG ??= ? The gravitational field, f r , is defined as a vector quantity showing the attraction of a unit mass, m1, placed at some point in space: r r Gm f ?.2 2=r (2.2) where m2 is the mass attracting the unit mass, m1 (Telford et al., 1990). Newton?s Second Law of Motion states that the force applied to an object is equal to the object?s mass, m1, multiplied by its acceleration: 9 amF rr 1= (2.3) If one considers the gravitational acceleration of the Earth, g r , and combines Equations 2.1 and 2.3, it is clear that: r r Gm g e e ?. 2 =r (2.4) where me is the mass of the Earth, re is the radius of the Earth and, in this case, r? is pointing towards the centre of the Earth. 2.1.2 Gravitational Potential Gravitational potential, U, is defined as the amount of work done in moving a unit mass from a very distant point of zero potential (i.e. ? ) along any path to a point, p, at a distance, R, from the centre of gravity of a mass, m. Integrating Equation 2.2: R Gm r dr GmdrpfpU R R =?=?= ? ? ? ? 2 )()( r (2.5) The potential is positive by convention and describes the energy of the system in a gravitational field, G. It is independent of the path that the unit mass follows between start- and end-points. It is also independent of time and is therefore considered a conservative field. Hence, for any path taken, when a mass is returned to its original position, the net energy expenditure is zero (Telford et al., 1990). From this definition, it is possible to state properties of an equipotential surface (i.e. a surface formed of points of equal potential): ? The gravitational field is everywhere perpendicular to the equipotential surface (Figure 2.1), ? The intensity of the gravitational field is generally not equal at all points on a equipotential surface, ? No work is done in moving a mass along an equipotential surface, and ? A body of water at rest or a levelled spirit-level define an equipotential surface. 10 a) b) Figure 2.1 ? Equipotential surfaces caused by different substratum. (a) Uniform substratum causes equal gravitational field vectors across the equipotential surface. (b) An ore body of arbitrary shape and with density greater than the substratum (i.e. ?2 > ?1) results in a mass excess causing unequal field lines across an up-warped equipotential surface (after Corner, 1993). Essentially, the gravity method requires a lateral change in density to model geological structure, and measures the density contrast directly. The target body or layer needs to show a significant mass in order to for it to be modelled. The seismic method also measures density changes, but can image horizontally layered targets (i.e. not reliant on lateral changes), and measures this density change through seismic wave travel times. Hence, it can resolve a reflection from a thin layer with a high density contrast (e.g. the high density chromitites of the Merensky Reef and UG2). 2.1.3 Reductions to Ground Gravity Data Latitude, elevation, topography of local terrain, Earth tides and lateral density variations all affect local gravity values. Gravity surveys are primarily concerned with lateral density variations. Hence, collected gravity readings have to be corrected for the first four factors (as well as drift and Earth tide) and reduced to a plane on, or parallel to, the geoid. This ensures that anomalies are only due to lateral spatial and density variations (Fourie, 1998; Telford et al., 1990). 11 2.1.4 Factors Affecting Rock Density Numerous factors can affect the density of rocks, including rock type (igneous, sedimentary, metamorphic), constituent minerals, porosity and depth. Lateral variations in density cause changes in the gravity field; hence an attempt to measure the density of rocks must be made to assist the modelling process. Densities can be measured in a number of ways, including measurements of outcrop or core samples, via borehole logging, or via seismic velocity estimates. Each of these methods only provides an estimate of the bulk density, and not an accurate value. Hence, the density in many field situations cannot be precisely known. Rock densities are highly variable, with differences occurring between rock types, as well as within a single rock; Table 2.1 shows the range of densities in rocks found in the study region. The average crustal density is 2.670 g.cm-3. Table 2.1 ? Densities of rocks in the study area (Mar? et al., 2002; Telford et al., 1990). Stratigraphic Group Rock Type Range (g.cm-3) Average (g.cm-3) Overburden 1.920 Soil 1.200 ? 2.400 1.920 Transvaal Supergroup Pretoria Group Shale 2.530 ? 2.850 Sandstone 2.650 ? 2.810 Bushveld Complex Main Zone Anorthosite 2.790 Norite 2.950 Norite-Anorthosite 2.860 ? 3.230 Gabbronorite 2.780 ? 2.980 Critical Zone Norite-Anorthosite 2.98 Norite 2.910 ? 2.940 Chromite 3.680 ? 4.220 Chromitite Layer (Merensky Reef) 2.770 ? 3.940 Dunite 3.750 Lower Zone Pyroxenite 3.300 Harzburgite 3.010 Norite* 2.930 Intrusives Syenite Dykes 2.630 ? 2.990 *Specific to area south of the Pilanesberg 12 Sedimentary rocks are generally the least dense rock type, relative to igneous and metamorphic rocks. This is because sedimentary rocks, in general, have a higher porosity. Obviously, the degree of consolidation of sedimentary rocks leads to considerable differences in density. Depth and duration of burial both allow for greater consolidation, a reduction in porosity and higher density. Igneous rocks are affected mainly by silica content and depth and time of crystallisation. Density commonly decreases with a decrease in silica content (acidic igneous rocks) and an increase in ferromagnesian content (basic igneous rocks). Rocks that crystallised at depth over a long period of time tend to be denser than extrusive rocks, even if the chemical composition is the same. Finally, if metallic minerals are found in significant proportions, the density can be greatly affected. For example, massive sulphide ore can have a density of twice the host rock. Metamorphic rock densities are affected primarily by the grade of metamorphism. As rocks are heated and compressed, the degree of porosity tends to decrease and mineral recrystallisation occurs. Hence, sedimentary rocks such as limestone and shale are less dense than their metamorphosed equivalents, marble and slate. 2.1.5 Gravity Measurement Gravity meters may be divided into two categories ? those measuring absolute gravity and those measuring relative gravity. Absolute gravity instruments measure the local value of gravity (g ? 980 gal) and include instruments such as Kater?s reversible pendulum and free-fall gravity meters. Relative gravity instruments measure the difference in gravity between stations (field stations are compared to a base station) and are of more use in exploration as they are more portable and take faster readings (Nabighian et al., 2005). 13 Many gravity meters have been invented, each with their own specific design and operation. Relative gravity measurements were first applied to exploration geophysics using portable pendulums, developed by Steneck (1887), where accuracies of 1 mgal were possible. Later, the E?tv?s torsion balance was invented by the Hungarian physicist, E?tv?s, in 1896 (Nabighian et al., 2005). It was also used in exploration from 1918 and continued to be developed into the 1950s. This apparatus measures the horizontal curvature of the gravitational field, not the vertical component, and was the first instrument to measure the gradient of the gravity field directly (discussed further in Chapter 3: Method ? Airborne Full Tensor Gradient Gravity). Spring gravity meters measure the change in equilibrium position of a proof mass as a result of a change in the gravity field between two stations (Nabighian et al., 2005). The two most common spring gravity meters used in the world are the LaCoste-Romberg and the Scintrex CG-3 and CG-5 models. The Scintrex CG-3 was used for this project?s ground gravity survey. 2.2 Scintrex CG-3 Autograv Gravity Meter The Scintrex CG-3 gravity meter (Figure 2.2) is an unstable spring gravity meter which uses a fused quartz zero-length spring and the change in position of a mass to detect variations in gravity. Variations in the mass? position cause a change in impedance of capacitance between the mass and an adjacent plate. An inbuilt micro- processor allows for automated measurements, with a measuring range of 7 000 mgal (without requiring resetting) and an accuracy of 0.005 mgal (Scintrex, 1995). This makes the CG-3 suitable for both detailed local surveys as well as regional surveys. 14 Mechanical X & Y tilt meters Data Output serial port Control Console Figure 2.2 ? Overhead view of the Scintrex Autograv CG-3, showing Control Console, Mechanical Tilt Meters (X and Y axis) and Data Output serial port (Scintrex, 1995). The automated measuring system continuously averages a series of one second samples, providing real-time signal enhancement and statistical analysis which automatically reduces noise. The readings are stored in the CG-3?s memory and automatically corrected for earth-tide variations. Other advantages include: robustness; electronic tilt meters (allowing for automatic adjustments); low machine drift (<0.02 mgal/day); unaffected by ambient temperature (-40 ?C to +45 ?C) and magnetic field (?0.5 mT) changes; 48 kb onboard memory; output to multiple computer devices; and application software (Scintrex, 1995). 15 2.2.1 In-Field Setup & Operation The CG-3 has a set of parameters that need to be setup from the control console before any field readings may be taken. The Autograv Setup is particularly important for an accurate survey. The parameters for the field work were input as shown in Table 2.2. Table 2.2 ? Scintrex CG-3 Autograv Setup parameters chosen for field work, October 2006 and May 2007. Option Selection Mode Field mode Read time 120 seconds Continuous tilt correction Enabled Auto-rejection Enabled Latitude/Longitude 25.5? S, 27? E Time Local time (i.e. GMT + 2 hours) Following the setup, various in-field procedures must be carried out to ensure an accurate survey. These include proper levelling procedures, using the mechanical and electronic tilt meters (i.e. level to within ? 5 arc seconds of the horizontal); safe transporting of the gravity meter (avoiding bumps and knocks); and managing noise- creating factors. Such noise-creating factors include high frequency natural vibrations (e.g. from ocean waves, for surveys near a coastline) and man-made vibrations (e.g. from passing vehicles, excessive movement of the user). The auto-rejection option eliminates noise spikes from a reading by rejecting sample readings four times greater than the standard deviation (Scintrex, 1995). The first five readings are not rejected and are used to calculate the initial standard deviation. Note that in-field processing is automated and completely separate from later office processing (see Chapter 2.6). The LCD screen on the control console shows the calculated standard deviation during measuring intervals. ?Quiet? measurements have a standard deviation of as low as 0.010 ? 0.030 mgal. In the ground survey, readings were taken for a minimum 16 of 60 seconds, after which the measurements were stopped when the standard deviation fell below 0.050 mgal. 2.3 Differential Global Positioning Systems (DGPS) The Global Positioning System (GPS) was initiated by the American Department of Defence in 1978. It is a radio-navigation system that consists of twenty-four satellites orbiting the Earth (Figure 2.3) and seventeen ground stations. These ground-stations monitor the orbit and operation of the satellites, assisting in reducing errors in the system (Merry, 1999). The satellites are arranged in space such that at least four satellites will be ?in view? at any time, all over the world. Figure 2.3 ? Schematic diagram of the orbital paths of the 24-satellite constellation of the GPS system, orbiting the Earth at an altitude of approximately 20 200 km. Note: diagram not to scale (Garmin, 2007). A gravity survey requires highly accurate station positioning, especially height. In order to obtain an in-field accuracy of 0.01 mgal with the CG-3, the horizontal position needs to be accurate to 15 m and the vertical position needs to be accurate to 5 cm. Standard handheld GPSs (e.g. Garmin? GPS 12 XL) are accurate to ~12.5 m in horizon position and up to 20 m in the vertical position. Hence, they are not suitable for any sort of gravity survey. 17 In order to obtain the high accuracy required, a Differential GPS (DGPS) was used for the survey. The Trimble? 5700TM and ReconTM were used to survey the gravity stations, allowing for sub-centimetre accuracies using Real Time Kinematic (RTK) processing. This allows for a horizontal error of ? 10 mm + 1 ppm and a vertical error of ? 20 mm + 1 ppm (Trimble, 2002; Trimble, 2006). An overview of GPS theory is now provided, including factors affecting the accuracy of standard GPSs and how the DGPS overcomes these problems. Finally, the DGPS instrument used in the gravity survey is reviewed. 2.3.1 Theory A standard GPS locates its position by calculating the distance between the receiver and a GPS satellite. The distance is calculated using the difference in time from when a particular code, superimposed on a microwave carrier signal, was emitted from the satellite to when the code was received by the user. The distances between all satellites in the ?field of view? and the receiver are used to determine the position of the receiver more accurately, using a process of trilateration (Merry, 1999). At least four satellites are required to accurately calculate the receiver?s 3D position. However, even with coverage of four or more satellites, there are a number of factors which can reduce the accuracy of standard GPS readings: ? Atmospheric delays: The signal transmitted by the satellite is slowed down and refracted by charged particles in the ionosphere and water vapour in the troposphere. The effect of the ionosphere is commonly a loss of 5 m-accuracy but can be as great as 28 m when satellites appear low on the horizon. The effect of the troposphere can be as great as 0.5 m. The GPS is programmed to calculate and correct for these delays, depending on the positions of satellites. ? Clock errors: The satellites are fitted with atomic clocks (accurate to 3 nanoseconds). However, very small errors can occur in the atomic clocks, which can lead to errors up to 1.5 m. Furthermore, during processing of the 18 received signal, the receiver introduces minor calculation errors and rounding errors, which lead to errors up to 1 m. ? Satellite orbital paths: The paths of satellites vary slightly from their calculated orbit, known as ephemeris errors. This error can be as great as 2.5 m. ? Multi-path errors: If a survey is in an area with large buildings or rugged topography, the satellite signal can be reflected before it reaches the receiver. The increased travel time can lead to errors up to 1 m. Even the more sophisticated GPS receivers are affected by multi-path errors (Merry, 1999). These errors can largely be reduced using a Differential GPS. A DGPS uses two receivers: one as a base station and the other as a ?rover?, which takes readings in the field. The base station is set up at a known point and remains stationary for the duration of the survey. The purpose of the base station is to record deviations from the known point (i.e. comparing the location recorded by the GPS with the known position). This provides a correction factor which is applied to the rover GPS data points. In this way, the DGPS is capable of minimizing the errors introduced by atmospheric effects, clock errors and orbital path errors, although it is still affected by multi-path errors. The correction factor can either be applied to the rover data in real-time or after the survey, depending on the type of DGPS. Real-Time Kinematic (RTK) DGPSs transmit the information from the base station to the rover in real-time, via a radio- signal. This allows processing to be done in real time, during the survey. Post- Processing Kinematic (PPK) DGPSs require data from the base station and rover to be downloaded onto a computer after the survey. Baselines (direct lines from the base station to the rover stations) are then calculated and corrected by the DGPS software. However, PPK surveys may delay effective surveys if there are problems with the data. These problems may only be realised after a day?s work, resulting in a portion or the entire day?s survey being redone. 19 2.3.2 Instrument The DGPS used in the survey consisted of a Trimble? 5700TM receiver (rover), a Trimble? Zephyr receiver (base) and a ReconTM unit. The base station was set up on a levelled tripod at a known point (trig. beacon ?449 PE-E-W?: N -2814562.953 m, E 10295.736 m, height: 1140.370 m) and recorded and transmitted changes in position for the duration of the survey. The rover and was connected to the ReconTM unit, which was set to RTK mode. After the base station has been set up, a radio link must be established between the base station and rover to allow for the Real Time Kinematic processing. Once connected, the rover is ready to survey station positions, ensuring that satellite coverage is always sufficiently maintained. The pre-determined xy-coordinates of each gravity station were programmed into the ReconTM unit. The DGPS user navigated to within ~1 m of the pre-determined position of each station co-ordinate. A DGPS reading was then taken, accurately recording the station position and height, after which a gravity reading was then taken at the same position. 2.3.3 Quality Control The RoverTM unit showed a warning when the Position Dilution of Precision (PDOP) was poor. PDOP is a measure of the quality of the satellite coverage in the region, i.e. the position of each satellite in the constellation as well as the geometry relative to the GPS receiver (Trimble, 1997). Hence, if the majority of satellites are low on the horizon, or if the majority of the satellites are closely grouped together, a poor PDOP value would occur. Poor PDOP values commonly last approximately 15 minutes, during which time the survey had to be stopped to avoid erroneous DGPS readings. 20 2.3.4 Downloading Data The survey data, stored on the ReconTM, are downloaded to computer after each day in the field. The data are transferred via a serial cable and can be saved in Trimble Geomatics OfficeTM as an ASCII file. 2.3.5 Data The final DGPS dataset is an ASCII file with the base station and each gravity station?s name and location in WGS84, LO27, as well as latitude, longitude and elevation. 2.4 Survey Design The gravity survey was designed with three main sections, each with a specific consideration (Figure 2.4): ? A high resolution Grid A (1200 m x 750 m, with a station- and line-spacing of 30 m), primarily to determine overburden thickness as well as any high frequency geological features which would provide a gravity anomaly, ? A larger, wider-spaced Grid B (2600 m x 1800 m, with a station- and line- spacing of 100 m), primarily to determine overburden thickness, and ? Two long, approximately parallel traverses (Line 1 and Line 2, 7000 m and 6500 m respectively, spaced ~800 m apart), to model the relationship of the contact of the Bushveld Complex with the underlying Transvaal Supergroup. The combined area of the grids covers various know geological features, including a large fault, known from seismic data, and numerous IRUPs, known from aeromagnetic data. The traverses cross the sub-cropping contact between the Bushveld Complex and Transvaal Supergroup, and Line 2 also crosses a major IRUP grouping. 21 Grid A Grid B DGPS Base Station Line 2 Line 1 Gravity Base Station Figure 2.4 ? Ground gravity traverses (crosses) and grids (dots) relative to outline of AGG survey area. The parameters of the ground gravity survey (i.e. line direction, line spacing, station spacing) were controlled by Reid?s criteria (Reid, 1980). The line direction of the grids was chosen as the flight line direction of the airborne FTG survey (i.e. bearing 040?). Reid (1980) suggests that, for a ground gravity survey in which the gradient of the field will be computed, station spacing should be equal to the mean height of the sensor from the source (h). Such station spacing would prevent aliasing of the source?s signal. For this survey, h may be considered the thickness of the overburden. Forward modelling from a survey conducted in the area in 2005, showed the weathered overburden to have a thickness of 20 ? 30 m. Hence, a station spacing of 30 m was considered appropriate for Grid A, realising that anomalies with wavelengths less than 30 m would be severely aliased. A line spacing of 30 m was 22 also chosen for Grid A, in order to prevent any directional bias in calculating the FTG components from the vertical gravity component data. A station- and line-spacing of 100 m was chosen for Grid B to model the regional trend of the overburden and to provide a dataset that could be compared with the airborne FTG data (see Chapter 5.4: Merging Airborne and Ground Datasets). A similar process was applied to determining the parameters of the traverses, where the primary target was the Transvaal Supergroup contact. The direction of the traverses followed the same bearing as the grids? lines, where possible (i.e. 040?). The Transvaal contact trends ~125? and it was desirable to cross this contact as close to perpendicularly as possible, allowing for optimum 2D modelling. However, thick vegetation, steep topography and land permission limited the length and direction of the traverse to roads in the region. Hence, a variation in bearing towards the south- west of the traverses was required. A station spacing of 150 m was chosen, giving a line:station spacing ratio of 1:5.33 (line separation of 800 m). The base station, for the entire gravity survey, was set up on solid bedrock in an empty river bed. The location for the base station was chosen as it was easily accessible (close to a road to the north of the grid), it was a quiet area (the road was not well-used), the surface was solid (the gravity meter would remain level during the reading and it was possible to mark the exact location of the station), and the station was easy to survey with the DGPS. 2.5 Quality Control The Scintrex CG-3 was programmed to sample for a maximum of 120 seconds whilst taking a reading. The on-screen standard deviation was observed and if, after 60 seconds, it dropped below 0.050 mgal then the reading was stopped. 23 External and internal repeat readings were carried out to compare readings outside of the current loop and within the current loop, respectively. External repeat readings of stations measured at the end of the previous day were taken at the start of each day. The amount of such repeat readings was determined as 5% of the number of readings taken the previous day (commonly 4-7 stations for Grid A and 2-4 stations for Grid B). Internal repeat readings of two repeat stations in the previous line were also measured. All repeat readings were used for micro-levelling the data. 2.6 Office Processing Once data has been collected, they can be downloaded and processed, applying various gravity data corrections. This produces a dataset with all gravity readings reduced relative to a specific station (e.g. the base station, local government gravity readings). The program Geosoft? Oasis MontaqTM was used to apply all the gravity reductions and the formulae used in the program are discussed below (Whitehead and Musselman, 2006). These gravity reductions are separate to the automatic in-field processing, carried out by the CG-3, and are applied after the survey. 2.6.1 Downloading Data The raw data can be downloaded from the CG-3 to a computer by connecting the CG- 3 via a standard printer cable. The data is transferred using programs such as ?idump.exe? (generic CG-3 software) or Geosoft? Oasis MontaqTM and saved as an ASCII file. The file contains the gravity reading, station number, standard deviation, tilt x, tilt y, sensor temperature, Earth tide correction, reading duration, number of samples rejected and time of the start of the reading. 2.6.2 Drift Correction All gravity meters experience some degree of mechanical drift. Although the CG-3 attempts to automatically correct for this drift, it is common practise to calculate the 24 drift for each loop as a secondary check. Hence, a gravity reading is taken at the base station every 3 hours (on the first few days of the survey) or at the start and end of the day (once satisfied that the drift is linear). Differences in the gravity readings due to mechanical drift are removed from the field readings. Assuming the drift (d) is linear, the base station reading is used at the beginning and end of each loop (gb1 at tb1 and gb2 at tb2, respectively): 12 12 bb bb tt gg d ? ?= (2.6) The value to be removed from the field reading is then calculated using the loop?s drift value and the time that the field station reading was taken after the first base station reading of the loop. 2.6.3 Latitude Correction The gravitational field of the Earth, along with the centrifugal acceleration, causes flattening at the poles and the Earth?s bulge at the equator. The increased distance from the centre of the Earth at the equator results in a decrease the gravity field, a portion of which is countered by the excess mass due to the bulge. The theoretical value of gravity, gt, at any latitude on the mathematically-defined reference spheroid is given by the 1980 International Gravity Formula: gt = 978 032.7(1 + 0.005 302 4 sin2? + 0.000 005 8 sin4? ) (2.7) where ? is the latitude (Whitehead and Musselman, 2006). For regional surveys, exceeding 2 km in latitude from the base station, the full International Gravity Formula has to be applied to the gravity readings in order to provide a full correction. For local surveys, within 2 km latitude of the base station, it is common practise to use an approximation of the 1980 International Gravity Formula to apply the latitude correction, requiring a single latitude value. Differentiating Equation 2.7: 25 ?d dg = 978 037.7 ? 0.005 302 ?2sin = 5 185.56 ?2sin gal.rad-1 (2.8) Now, taking the radius of the Earth at the Equator, re = 6 378 km, and ds as the horizontal distance north or south of the latitude, ? : ?? 2sin812.0 1 == d dg rds dg e lat mgal.km-1 (2.9) The correction, dglat, is added as a survey traverses towards the equator and subtracted as a survey traverses towards the poles. This corrects for the shape of the Earth, with the gravity value decreasing as one traverses towards the equator (further from the centre of the Earth) and increasing as one traverses towards the poles (closer to the centre of the Earth), (Telford et al., 1990). For many present-day gravity meters, a machine accuracy of 0.010 mgal is achievable in the field. In order to obtain such accuracy, the accuracy of the location north or south of the station must be within 15 m. 2.6.4 free-air Correction As seen in Equation 2.1, gravity is inversely proportional to the square of the distance between the source mass and the point of observation (Telford et al., 1990), hence gravity measurements are highly dependent on topography. free-air corrections remove the effect of each station?s difference in height relative to the base station. These corrections do not take into account the difference in material beneath the station, due to the increase or decrease in height. By differentiating the scalar form of Equation 2.4 (Telford et al., 1990): ee e e eFA r g r Gm r Gm dr d dr dg 22 32 ==??? ? ??? ?= (2.10) Hence, 3086.0= dh dg FA mgal.m-1 (2.11) 26 substituting the equatorial radius for re and whe dh is the difference in height above/below the base station. This correction is added levation, compared to latitude. For accuracies of 0.010 mgal, the surveyed station he free-air correction accounted for the change in elevation between a station and a correction accounts for the amount of excess/deficit te, horizontal flat slab of niform thickness and density (Telford et al., 1990). The mass of the theoretical slab re dgFA is the free-air correction and for stations above the base station and subtracted for stations below the base station. As can be seen, the inverse-square relationship leads to a far greater dependency on e elevation should be accurate to within 5 cm. 2.6.5 Bouguer Correction T base station. The Bouguer material, due to a point being above/below the base station. The correction assumes the station is positioned on an infini u is used to approximate the attraction of excess/deficit material. The attraction of the slab is: hGg B ?= ??2 (2.12) which is differentiated by h to give the Bouguer correction: ?04192.0= dh mgal.m (2.13) dg B where dgB is the Bouguer correction, dh is th base station and ? is the density of the rock (g.cm ). The Bouguer correction is in Correction he Bouguer correction assumed a horizontal slab with no terrain irregularities. Obviously, this does not occur in real-life situations and severe topography (mountain -1 e difference in height above/below the -3 subtracted if the station is above the base station and added if the station in below the base station. 2.6.6 Terra T 27 ranges, hills, gorges) should be taken into account. Obviously, hills and valleys in s the digital rrain model into a series of cylinders, cubes and triangles, in order to closely s relatively flat. A ?rule of thumb? is: (2.14) ated, the correction is small enough to discard (Telford et al., 1990). The Styldrift region is a flat landscape, with o quations, providing the free-air gravity value, gFA, and the Bouguer gravity value, close proximity to a gravity meter will result in a mass excess/deficit, respectively, and the effects of these features must be included in the gravity reading. Computer programs (e.g. Geosoft? Oasis MontaqTM) require detailed digital topographical data to run terrain corrections. The program segment te approximate the terrain. The attraction for each segment is then added to give the topographic correction for each station. The correction is always added because the mass of hills, above the gravity station, exert an attraction (upward pull) on the gravity meter. In a similar manner, the mass deficit from valleys, below the gravity station, does not exert an attraction; hence the effect must also be added (Telford et al., 1990). The terrain correction is a computationally intense method, which can be ignored if the terrain i zr 20? where z is the height of the topographic feature and r is the distance from the centre of the segment to the station. If r is greater than st topography of Pilanesberg close to the northern border of the study region. The ratio r/z was ~22 for the survey area; hence no terrain correction was applied to the data. 2.6.7 Gravity Anomalies The corrections used to correct for this project?s data were combined into tw e gB, for the survey: ?? ??? ? ??= FA stationbaseabove stationbasebelowlat south northobsFA dgdggg (2.15) 28 ??)B (2.16) ??? ? ???= ( FA stationbaseabove stationbasebelowlat south northobsB dgdgdggg where gobs is the drift- and tide-corrected value (for the southern hem the latitude correction, is the free-air correction and is the Bouguer var vity Dataset he final product of the office processing for the grid and traverses, was a Geosoft? d processed data (i.e. free-air and Bouguer isphere), latdg is FAdg Bdg correction (Telford et al., 1990). A free-air and Bouguer anomaly map can be produced from this dataset. These maps will show lithological iations as well as regional anomalies. Regional variations may have to be removed in order to show smaller, local anomalies. 2.7 Final Ground Gra T database (?.gdb? file) containing raw an anomalies), locations of stations, elevations and errors in the data. 29 CHAPTER 3: METHOD ? AIRBORNE FULL TENSOR GRADIENT GRAVITY 3.1 Introduction The following chapter provides an overview of the history of the airborne Full Tensor Gradient (FTG) gravity method, as well as the basic theory behind FTG gravity. The two significant Airborne Gravity Gradiometry (AGG) systems in use today (Lane, 2004a), the FALCON? and Air-FTG?, are compared. Of these, the Air-FTG? system was chosen to fly a survey over the study region. Survey parameters and processing procedures are also summarised. 3.1.1 History of Airborne Gravity Instruments Conventional gravity meters and gradiometers have been used in exploration since the late 19th century (Nabighian et al., 2005). While conventional gravity meters (e.g. LaCoste-Romberg, Scintrex CG-3 and CG-5) only measure the vertical component of the gravity field (Gz), gradiometers (e.g. E?tv?s? torsion balance, Air-FTG? and FALCON?) measure one or more of the gradients of the three principle components (Gx, Gy, Gz). The history of AGG instruments can be traced to the E?tv?s torsion balance, invented in 1896, and used commercially in mineral and petroleum exploration from 1918 until as late as 1940 (Nabighian et al., 2005). The torsion balance uses two proof masses suspended from an aluminium balance bar. The bar rotates when placed in a horizontally differential gravitational field, allowing horizontal changes in the gravity field to be mechanically measured. Readings took three to six hours to record, allowing for only a few measurements a day (Nabighian et al., 2005). Advances in conventional gravity meter technology, especially the Lacoste-Romberg gravity meter, saw them gaining favour over the torsion balance after 1940. This was 30 predominantly due to the invention of the zero-length spring (LaCoste, 1934), which made readings faster and gravity meters more robust. Airborne gravity surveys, measuring the vertical component, were introduced in 1983 (Hammer, 1983). An accurate airborne gravity method was highly sought-after, since ground surveys are commonly slow or extremely difficult to carry out (e.g. inaccessible terrain or environments). However, airborne data were very difficult to collect, due to the low signal-to-noise ratio and the problem of measuring accurate accelerations in an accelerating frame of reference (Bell et al., 1991). The introduction of DGPS has made this method increasingly accurate, with latest estimates of airborne gravity errors on a fixed-wing system varying between 0.2 ? 1 mgal (Fairhead and Odegard, 2002) to 2 mgal (Hwang et al., 2007). AGG systems took longer to develop to a suitable standard for accurate surveys. The first airborne gravity gradiometer (AGG), the FALCON?, was developed by BHP Billiton for use in 1999 (Dransfield et al., 2001; Lee, 2001). Bell Geospace?s Air- FTG? system has its history in the 1970?s and 1980?s, when Bell Aerospace (now called Lockheed Martin) developed a gradiometer for the US Navy?s submarine ?stealth? technology (Murphy, 2004). The gradiometer was originally designed to aid submarine navigation in imaging the sea floor. In 1998, this technology was declassified and commercialised by Bell Geospace for use on marine petroleum exploration as the ship-mounted Marine-FTGTM system (Hammond and Murphy, 2003). This system was later developed for airborne use, as the Air-FTG?, which became operational in 2002. These two projects have resulted in the two most successful AGG systems (Lane, 2004a), capable of capturing the high frequency signal associated with near-surface density variations more clearly than conventional vertical component gravity instruments (Murphy, 2004). 31 3.1.2 FTG Theory Gravitational potential, U, is a scalar quantity which describes the potential energy associated with a unit mass in a gravitational field, G. This gravitational field is a vector field, describing the spatial variations in the potential. It may be considered the gradient of the potential and can be represented by the three principle gravity vectors (denoted Gx, Gy, Gz), which are mutually perpendicular. It is common to consider the x-axis trending east-west, the y-axis trending north-south and the z-axis vertical (Murphy, 2004). FTG theory exploits information provided by measuring the rate of change of the gravitational field with respect to the three principle gravity vectors. By taking the derivatives of the three principle components with respect to x, y and z, nine components are formulated, making up the FTG matrix, T r (Figure 3.1). The components of this matrix are Txx, Txy, Txz, Tyy, Tyx, Tyz, Tzz, Tzx and Tzy, where: ? Txx, Tyy and Tzz describe the rate of change of a component in the same direction (e.g. the rate of change of Gz in the z-direction), ? Txy and Tyx describe the rate of change of the x-component in the y-direction and vice-versa, ? Txz and Tzx describe the rate of change of the x-component in the z-direction and vice-versa, and ? Tyz and Tzy describe the rate of change of the y-component in the z-direction and vice-versa. 32 a) b) c) Figure 3.1 ? a) Schematic co-ordinate system showing the three principle components of the gravity field (Gx ? blue, Gy ? green, Gz ? red), with derivatives in the x, y and z direction of each component creating the nine components of the gravity tensor. b) The components of the FTG gravity field, which describe the rate of change of the gravity vector as one moves in the direction of the vector. c) The FTG matrix, Tij, showing the components which are equal (i.e. Txy = Tyx, Txz = Tzx, and Tyz = Tzy), hence one of each of these components can be considered redundant (Bell Geospace, 2006). Heath (2007) shows that the following conditions apply to potential fields: 0 0 =?? =?? T T r r (3.1) which expands to: zyyzzxxzyxxy TTTTTT === (3.2) 0=++ zzyyxx TTT (Laplace?s Equation) (3.3) Hence, due to redundancies and the fact that Tzz can be calculated (since Laplace?s equation must be satisfied, i.e. Tzz = -(Txx + Tyy)), only five of the nine gradients are independent, i.e. Txx, Txy, Txz, Tyy, Tyz, (Hammond and Murphy, 2003). Despite Tzz being a dependent component, it is still commonly mapped together with the five independent components as it reveals geology most tellingly. 33 Figure 3.2 illustrates the Gz and FTG gravity responses for a buried synthetic cube. Bell Geospace (2006) explains that each component represents a different attribute of the geological source. The Txx gradient measures the east-west changes in an east- west direction. As a result, it is suitable for determining the north-south trending edges of a body. Similarly, Tyy gradient data measures the north-south changes in a north-south direction, and determines the east-west trending edges of a body. Txy measures the east-west changes in a north-south direction. It shows apparent ?corners? of mass anomalies and can be used to plot the centre of the causative body in plan view. Txz and Tyz measure the east-west and north-south change in a vertical direction and help determine the central axes of the body, in their respective directions. They also delineate north-south and east-west trending edges (e.g. faults, dykes). Tzz gradient data measures vertical changes in a vertical direction. Since Tzz is a second derivative, it presents targets easily, highlighting all the edges (since Tzz is a combination of Txx and Tyy) and giving the dominant shape of the mass anomaly (Bell Geospace, 2006). 34 Txx Txy Txz Tyy Tyz Eo Gz 4 km mgal Tzz 4 km Figure 3.2 ? Total gravity field (Gz) and FTG gravity components for a buried synthetic cube (outline of cube presented. Parameters: dimensions 1 000 m x 1 000 m x 1 000 m; depth 25 m; density contrast 1.0 g.cm-3). As well as using Laplace?s equation to relate Tzz to Txx and Tyy, it is also possible to calculate various components from other components (e.g. Tzz from Txz and Tyz). Indeed, each of the components can be calculated from a single measured component (Nelson, 1988), using a Fourier transform (further explained in Chapter 5.1: Introduction). Consider a Fourier transform pair ),(),( vuFyxf ? , with the spatial frequency domain variables, u and v (where k2 = u2 + v2). Txz and Tyz can be related to Tzz as follows: ),(),()),(( vuT k iu vuTyxtF zzxzxz == (3.4) 35 ),(),()),(( vuT k iv vuTyxtF zzyzyz == (3.5) Relating Txx to Txz and Tyy to Tyz: ),(),()),(( vuT k iu vuTyxtF xzxxxx == (3.6) ),(),()),(( vuT k iv vuTyxtF yzyyyy == (3.7) Finally, relating Txy to Tyz: ),(),()),(( vuT k iv vuTyxtF yzxyxy == (3.8) where 1?=i . However, if one were to calculate each of the components from a single measured component, the noise from the measured component would be transferred to each of the computed components (Heath, 2007). If one used the calculated data to run FTG gravity inversions, the repeated noise could present itself as a geological anomaly which the algorithm would try to fit. By measuring all five components independently, the noise in each measurement remains independent. Latest developments in FTG gravity inversion algorithms use multiple components in fitting a model. Hence, if noise is related to separate components, it will be discarded by these inversion algorithms and only the geological anomalies will be modelled (Heath, 2007). 3.2 FALCON? The FALCON? system is not a full tensor gradient system, since it only measures Txz and Tyz and uses these to calculate Tzz (Heath, 2007), as in Chapter 3.1.2: FTG Theory. The FALCON? gradiometer takes a reading every 0.1 seconds, hence, if the system is moving at 200 m.s-1 (typical aircraft speed), it takes a reading approximately every 20 m. Filters are applied to the data to suppress high frequency 36 data (i.e. noise) and produce a cleaner signal. Dransfield (1994) developed a filter that suppresses spikes in the data with wavelengths less than 125 m, while flying at 80 m above ground level. Any geological anomalies with wavelengths less than 125 m are deemed to exist near the surface and would not be imaged for exploration (Heath, 2007). The FALCON? has been used in numerous successful surveys around the globe, including kimberlite detection (Hinks et al., 2004), iron ore, iron oxide copper gold (IOCG) and various classes of base metal deposits (Dyke et al., 2002). Forward modelling of kimberlite pipes (Hinks et al., 2004) shows that the smallest pipes which could be targeted have an area of 25 Ha, at a depth of 100 m. Dransfield and Lee (2004) describe the FALCON? as using a gravity gradient instrument (GGI) sensor mounted on an inertial platform. The GGI is composed of four low-noise accelerometers mounted orthogonally on a rotating block. The design differs from Bell Aerospace?s in several ways: increased spacing of the accelerometers (increased by a factor of two), an additional, independent accelerometer on the GGI rotor, and the rotor axis is aligned close to the vertical (Dransfield and Lee, 2004). The gravity gradient signal is provided by summing the output of the accelerometers. The AGG instrument also measures any accelerations and rotations of the system (i.e. aircraft and instrument) and uses this information to model the noise experienced. This noise can be removed via post-processing compensation (Dransfield and Lee, 2004). The FALCON? system is commonly flown in a Cessna Grand Caravan, modified to house the AGG system (Figure 3.3). The latest FALCON? system has been housed in a Eurocopter AS-350 B3 helicopter, producing improved results (Lee et al., 2006). The system also contains a laser altimeter, DGPS and a stinger-mounted caesium- vapour magnetometer. The laser altimeter accurately measures the distance of the instrument from the ground. Combining the altimeter data and accurate positions 37 provided by the DGPS, a Digital Terrain Model (DTM) of the survey region can be created. The DTM is used for terrain corrections which can be directly applied to the gravity gradient data (Dransfield and Lee, 2004), further discussed in Chapter 3.5.2: Terrain Correction. Figure 3.3 ? Cessna Grand Caravan, modified to house the FALCON? system for airborne surveys (BHP Billiton, 2007). 3.3 Air-FTG? Bell Geospace?s Air-FTG? system was the first system capable of measuring FTGs from an airborne platform. It comprises three GGIs, each containing two pairs of orthogonally mounted accelerometers, each pair 10 cm apart, on rotating disks (Figure 3.4). The GGIs axes are mutually perpendicular, with each axis having the same angle with the vertical (35?). Hence, in plan view, the axes are separated by an angle of 120?. The GGIs spin at a set frequency of 0.5 Hz and the entire system of GGIs, mounted on a three gimballed stabilized platform, also rotates about the vertical axis at 300?/hour. These rotations avoid bias in the direction of the primary 38 components (i.e. Gx, Gy, Gz) and reduce noise (Murphy, 2004). The gradient of the gravity field is measured as the difference in readings between the opposing pairs of accelerometers on each disk. Each GGI produces two gradient measurements, in the plane of the rotating disk, dependant on the distance between the accelerometers and the frequency of spin on each disk. In order to obtain the tensor components in the external co-ordinate system, an appropriate linear combination of the six GGI outputs is required (Murphy, 2004). The difference of the gravity field, as sensed by each pair of accelerometers, is used to compensate for aircraft motion and preserve high frequency geological signals (Hammond and Murphy, 2003). A technical paper describing the Air-FTG? gravity gradiometer instrument is provided by Brett (2000). 39 a) b) c) CGIs Figure 3.4 ? a) Schematic diagram of the Air-FTG? instrument, showing the arrangement of the accelerometers within each spinning disk, b) image of the Air-FTG? system, showing positioning of CGIs, and c) geometric arrangement of the GGIs (Murphy, 2004). Using Figure 3.4 as a reference, Zuidweg and Mumaw (2007) show, following Rummel (1986), that by subtracting the opposing accelerometers? readings (e.g. a1 and a2), the disk?s acceleration is cancelled out, leaving the following first-order approximation: jiji dxTda = (3.9) 40 where dai is the acceleration difference vector, dxj is the co-ordinate difference between the two accelerometers and Tij is the gravity gradient tensor. For a horizontal disk, it is possible to define changes in the local co-ordinates: 0 sin cos = ?= ?= dz Ddy Ddx (3.10) with the z-axis vertical, D the distance between the accelerometers and the angle between the x-axis and the line connecting the two accelerometers. Hence, the difference between accelerometers, da, on a horizontal disk (i.e. axis vertical) is: ? tTtTTDda xyxxyy ?? 2cos2sin)(21 +?= (3.11) where the disk is rotating with angular speed ? , at time t. Similarly, for vertical disks in the yz- and xz-plane, respectively: tTtTTDda yzyyzz ?? 2cos2sin)(21 +?= (3.12a) tTtTTDda xzxxzz ?? 2cos2sin)(21 +?= (3.12b) Hence, it is possible to resolve the five independent tensors and Tzz, combining the above equations (Zuidweg and Mumaw, 2007). Similar to the FALCON?, the Air-FTG? system is flown in a Cessna Grand Caravan, together with equipment for terrain measurement, DGPS and magnetometer. In all, the system weighs over 450 kg; hence the choice of aeroplane is critical. Since the Cessna runs on a single engine, the vibrations contributing to noise are more manageable than a twin engine (Murphy, 2004). The GGI is positioned close to the aeroplane?s centre of gravity (i.e. centre of the pitch, roll and yaw), thereby minimizing rotational accelerations. External accelerometers are also installed, which measure any other accelerations (e.g. due to engine activity, varying propeller speeds) which can be removed in post-processing (Hammond and Murphy, 2003). 41 3.4 Air-FTG? Survey Design As with all geophysical surveys, it is important to optimise the survey design in order to extract the most information from the data collected. All priori knowledge of the area to be surveyed must be considered before deciding the station spacing, line spacing, line direction, flight height and area to be covered. The target size, orientation, density contrast and depth of burial all affect the survey design (Murphy, 2004). Forward modelling is commonly used to optimise the design, including factors such as terrain, aircraft climb performance and tie-line spacing (Hammond and Murphy, 2003). Reid?s criteria (Reid, 1980) for aeromagnetic surveys determine optimum survey parameters relative to the distance of the sensor (i.e. magnetometer) from the magnetic source. The FTG gravity signal fall-off is very high because it is proportional to the inverse cube of the distance to the source (i.e. ), (Hammond and Murphy, 2003). However, the same relationship applies to a magnetic signal; hence Reid?s criteria also apply to FTG gravity data. Reid (1980) based a survey?s optimum parameters on the distance, h, of the sensor from the upper surface of the source. Line parameters and flight height will now be discussed with respect to h. 3?? RTij 3.4.1 Flight Height Surveys may either be flown at a constant altitude or in a gentle drape, maintaining a constant height above the ground (Figure 3.5). Due to the high signal fall-off of the FTG signal, it is important to fly as close to the ground as safely possible. A gentle drape with heights of 80 m or above are commonly chosen, although heights of 60 m have been flown (Heath, 2007). This survey was flown at a height of 80 m above the surface of the ground. The survey region is very flat, hence, a draped survey was chosen. The survey?s targets 42 include IRUPs, dykes, faults and geological contacts (i.e. Transvaal Supergroup contact). Each of these features are close to surface or buried by a <15 m thick weathered zone. As an over-estimate, the distance from source to sensor, h, is 100 m. Figure 3.5 ? Comparing surveys flown at a constant height or draped, as well as the actual height flown for a draped survey in regions of steep topography (Syberg, 1972). More closely draped surveys are becoming possible, with improved DGPS systems and laser altimeters, allowing for increasingly accurate detail in airborne surveys. However, in regions of steep topography, planes cannot maintain the constant height above ground and, for safety, have to fly at a different altitude. The advantage of flying a survey at a constant altitude is that the aircraft is always above the highest point in the topography. However, this causes the data quality to vary with altitude, since the measured signal falls off where the altitude is high (i.e. over low-lying areas). Even though it is possible to variably downward continue the data, this process increases the noise in the data. 3.4.2 Production Lines Reid (1980) determined that station- and line-spacing should not exceed 2h. The station spacing of the survey is determined by the frequency which the Air-FTG? system records at (i.e. 1 Hz) and the speed of the aircraft, commonly ~60 m.s-1. This results, effectively, in a ?ground? spacing of ~60 m (i.e. less than 2h). It is important to note that the shortest spatial wavelength, (i.e. the Nyquist wavelength, ?N) that can be determined for a station spacing, ?x, is given by ?N = 2?x (Blackman and Tukey, 43 1959). For the given station spacing, the Nyquist wavelength is 120 m, hence, any wavelengths shorter than 120 m will be aliased. Air-FTG? surveys can be flown with a line spacing of 50 m to 2 000 m, depending on the target and scope of the survey. Detailed surveys (e.g. kimberlite exploration) commonly use line spacings of 50 m to avoid aliasing (Figure 3.6), while regional surveys use line spacings of up to 2000 m (Murphy, 2004). Despite wider line spacings, the system benefits from increased confidence when interpolating from line-to-line, due to the measurement of the horizontal components Txz and Tyz (Schmidt and Clark, 2000). A line spacing of 100 m was chosen for the survey, less than 2h, due to the survey requiring high resolution. a) b) c) Figure 3.6 ? The effect of line-spacing on Tzz data, where lines have been removed from an original survey. a) Original survey using 50 m line-spacing, b) simulated survey using 150 m line-spacing, and c) simulated survey using 250 m line-spacing. The target anomaly (length: 350 m) is still resolved in b) but already shows aliasing in c). The negative anomaly (circled) all but disappears in c), (Murphy, 2004). 3.4.3 Tie Lines Tie lines act as additional flight lines, perpendicular to the flight lines. They serve as checks at cross-points with the flight lines and aid the gridding of data. In aeromagnetic surveys, they are of particular use when correcting for the diurnal variation. However, in AGG surveys they are simply used for micro-levelling purposes. Teskey (1991) recommends an optimum tie line spacing three times the line spacing. However, he also notes that spacings ten to twenty times the line 44 spacing are more commonly used. The FTG gravity survey uses a tie line spacing of 800 m, eight times the line spacing. 3.4.4 Flight Direction Surveys are commonly designed to traverse the target (i.e. fly perpendicular to strike or the length of the body). For multiple targets, targets are prioritised and a survey direction suitable for resolution of the most important bodies is selected. 3.4.5 Styldrift Flight Parameters The flight parameters chosen for the study region were as follows: ? Aircraft: FAS Cessna Grand Caravan ZS/SSa ? DGPS: Fugro Fasdas Digital Acquisition System ? Flight Speed: 60 m.s-1 (~215 km.h-1) ? Data collection frequency: 1 Hz o Resulting in a station spacing of ~60 m ? Line-spacing and trend: 100 m, 040? o Resulting in a station-spacing:line-spacing ratio of 1:1.67 ? Tie-line spacing and trend: 800 m, 130? ? Flight Height: 80 m (draped survey since region is very flat) 3.5 Air-FTG? Processing The collected FTG data is stored on disks, in-flight, and downloaded to a processing computer after each flight. Processing is divided into five key stages: high rate processing, terrain correction, spike and shift removal, levelling, and Full Tensor Processing (FTP). The AGG data was collected and processed by Bell Geospace. A brief overview of the processing stages is presented. Since the processing is all in- house, few published explanations exist to explain each of the stages. 45 3.5.1 High Rate Processing This stage runs algorithms which compensate for the effects of external forces, i.e. self-gradient of the aircraft and the instrument, which include accelerations on the instrument (e.g. turbulence, the rotation of the disks) and mass shifts (e.g. drainage of fuel) (Murphy, 2004). The data is then deconvolved and individual tensor components are calculated (Hammond and Murphy, 2003). 3.5.2 Terrain Correction Since terrain is the closest density contrast to the instrument, it accounts for the largest gradients in the data (Murphy, 2004). Hence, highly accurate terrain measuring equipment is required to construct an accurate DTM. The terrain correction uses a forward model of the terrain?s gravitational effect (using the DTM) which is removed from the data. 3.5.3 Spike and Shift Removal A check for the harmonic fit of the data (i.e. ensuring Laplace?s Equation is solved for by the components, Equation 3.3) and correlation between multiple survey lines is carried out (Murphy, 2004). In this way, the effects of drift, spikes and shifts can be filtered out. Other signals which do not pass these tests are considered noise and are discarded (Hammond and Murphy, 2003). 3.5.4 Levelling A line correction is applied to the data, where DGPS and inertial navigation data correct for heading errors. Standard airborne processing methods are then applied to the data (e.g. filtering, line levelling and micro-levelling). Levelling is applied manually and as a network adjustment (Zuidweg and Mumaw, 2007). These processes remove subtle drift and shift effects. At this stage, the standard deviation of noise is approximately 5 Eo (i.e. ~0.5 mgal.km-1), determined by subtracting 46 measured airborne data from upward continued ground data, converted to FTG components (Murphy, 2004). 3.5.5 Full Tensor Processing This processing algorithm (Murphy et al., 2006) optimises the relationship between the full tensor components and removes any uncorrelated noise. The five independent FTG components are used to produce a single harmonic potential map. The FTG components are then recalculated from this potential data map and compared to the measured data, determining noise which can be removed (Murphy et al., 2006). Improvements in the data due to this processing show a standard deviation of 2 ? 3 Eo over 300 ? 400 m spatial wavelengths. During the processing stages, constant checks are applied to the data, looking for excessive accelerations, calibration errors and repeat differences. Any data that does not meet the standards imposed for the survey are re-flown and merged. Once fully processed, maps for all measured components (Txx, Txy, Txz, Tyy and Tyz) and Tzz can be plotted and interpreted. 3.6 Airborne FTG Gravity Dataset The final product of the office processing was a Geosoft? database (?.gdb?) file, provided Bell Geospace. The file contains xy-locations of stations, raw and processed FTG data, DGPS altimeter elevations, Air-FTG? information and aircraft information (e.g. north- and east-velocity, pitch, roll). 47 CHAPTER 4: METHOD ? GROUND MAGNETICS 4.1 Introduction This chapter covers the theory of geological magnetic anomalies (i.e. spatial anomalies) and variations in the geomagnetic field (i.e. periodic anomalies), before providing parameters used in an aeromagnetic survey of the region. The instruments used in the ground magnetic survey are then described along with parameters of the ground survey. Following this, processing of the magnetic data is explained. Aeromagnetic data and ground magnetic data were used to delineate magnetic geological bodies in the region (i.e. dykes, IRUPs), to aid the interpretation of the gravity data. The aeromagnetic data covers an area larger than the AGG data and the ground magnetic data follow a similar path to the ground gravity traverses. The ground magnetic traverses provide a higher resolution magnetic profile than the aeromagnetic data and also serve to ground-truth the aeromagnetic data. Similar to FTG gravity theory, the magnetic field has an inverse-cubed (B ? R-3) relationship with the distance to source (i.e. for a dipole source). Hence, anomalies measured during aeromagnetic surveys have lower amplitude, longer wavelength responses relative to the same anomalies measured during ground surveys, which have larger amplitude, shorter wavelength responses. Furthermore, subtle variations in the magnetic field are attenuated or aliased with increased distance from source. 4.1.1 Local Magnetic Anomalies The Earth?s geomagnetic field varies approximately with latitude, being strongest in the poles and weakest at the equator. The geomagnetic field strength in South Africa is approximately 28 000 nT. Local magnetic anomalies due to crustal rocks (i.e. less than 40 km depth) can cause significant variations from this value. 48 Magnetic anomalies are most easily understood by considering a magnetic body with an induced dipolar field. The magnetic field of the body and the geomagnetic field interact, creating an anomaly over the body (Figure 4.1), common to the specific latitude. Other factors such the body?s strike direction, dip, geometry and remanence also affect the shape of the anomaly. These properties must be taken into consideration when modelling a magnetic body. N S Figure 4.1 ? South-North schematic profile of a magnetic anomaly over a magnetic dyke with a dipolar field. The dyke is striking E-W, and situated in the southern hemisphere (Roux, 1980). 4.1.2 Magnetic Properties of Rocks & Minerals The magnetic polarisation, M r , of a body is dependent on the body?s magnetic susceptibility, k, and applied magnetic field, B v : BkM vr = (4.1) 49 Magnetic susceptibility is the fundamental rock property in rock magnetism (Telford et al., 1990), similar to rock density in gravity studies. A rock?s magnetic susceptibility is primarily dependent on the amount of ferrimagnetic minerals present in the rock, where minerals divide into sub-domains which are either aligned or opposed to the external magnetic field. Rocks can obtain net positive magnetic alignment with the external field in two situations: when sub-domains have a strong positive alignment (e.g. magnetite), or when large numbers of sub-domains are positively aligned (e.g. pyrrhotite), (Telford et al., 1990). Even though there is a wide variation in rock magnetic susceptibilities, sedimentary rocks generally show the weakest susceptibility while basic igneous rocks show the highest (Table 4.1). Table 4.1 ? Magnetic susceptibilities of minerals and rocks in the study region (Telford et al., 1990). Note minor variations between rocks of the Bushveld Complex and Transvaal Supergroup, but these are small in comparison to the variation with intrusives, which have a prominent magnetic signature. Stratigraphic Group Rock Type Range (x 103 SI) Average (x 103 SI) Minerals Quartz -0.01 Pyrrhotite 1 ? 6000 1500 Magnetite 1200 ? 19200 6000 Transvaal Supergroup Pretoria Group Shale 0.15 ? 0.36 Sandstone 0.02 ? 0.7 Bushveld Complex Main Zone 131.5 (ave.) Norite-Anorthosite 41.0 ? 213.0 Gabbronorite 100.0 ? 450.0 Gabbro 21.5 ? 195.5 Critical Zone Gabbronorite Norite-Anorthosite 35.6 Norite 44.1 Chromite 116.2 Lower Zone Harzburgite-Pyroxenite 55.0 Norite 26.4 Intrusives Syenite Dykes 1 178.1 ? 1 348.1 4.1.3 Induced & Remanent Magnetisation Induced magnetism is shown in rocks containing magnetic minerals, where the Earth?s magnetic field induces a magnetic field. Remanent magnetism occurs when 50 rocks acquire a ?permanent? magnetisation, either intermittently or continuously, during the rock?s formation. This magnetisation is in the direction of the Earth?s magnetic field at the time of the rocks formation. Subsequent movement of the rock (ranging from continental drift to localised folding and faulting) causes the direction of remanent magnetisation to vary widely from the present-day induced magnetic field (Figure 4.2). S N Figure 4.2 ? The effect of remanence on S-N profile: the induced magnetization vector and remanent magnetization vector combine to produce a resultant vector, commonly different to the present-day field (Roux, 1980). The most common causes of such Natural Remanent Magnetisation (NRM) are: cooling of a rock below the Curie point in a magnetic field; magnetic sedimentary grains preferentially aligning themselves in a magnetic field; chemical formation, crystallisation or alteration of magnetic minerals in a magnetic field; reorientation of magnetic grains due to high pressures; and reorientation of the local magnetic field due to lightning strikes (causing very localised remanence), (Telford et al., 1990). 51 4.1.4 Diurnal Variations & Magnetic Storms A constant stream of plasma is issued from the Sun (the solar wind), which impacts the Earth. This plasma causes a daily variation in the magnetic field, known as the diurnal variation. The diurnal has a regular period of 24 hours and is caused by the following factors: interaction of the charged plasma with the ionosphere (the layer of plasma that exists around the Earth between 90 ? 1000 km altitude, Sutcliffe (2004)), heating of the atmosphere from the Sun (causing expansion), and, to a lesser extent, the gravitational attraction of the atmosphere to the moon. The periodic motion of charged particles in the upper atmosphere sets up convective cells, which induce magnetic fields. Essentially, a periodic self-supporting dynamo is generated, creating two cells: one in the sun-lit Northern hemisphere (moving anti- clockwise) and one in the sun-lit Southern hemisphere (moving clockwise, Figure 4.3). Figure 4.3 ? Dynamos in the upper atmosphere generating magnetic fields which causes the diurnal (USGS, 2005). This produces a quiet, regular background variation in the measured magnetic field as the Earth rotates. The variation affects the intensity of the magnetic field at the 52 Earth?s surface in the order of 10?s nT, with the strongest variation occurring at the magnetic poles (up to 30 nT) and decreasing towards the magnetic equator (Figure 4.4). Diurnal variations towards the magnetic equator are largely invariant over hundreds of kilometres, in terms of amplitude and phase (Corner, 1993). During a magnetic survey, the diurnal variation must be recorded and removed from the field data. This is done by setting up a base station, i.e. an immobile magnetometer which measures the magnetic field every minute (or other user-defined period), hence recording the diurnal. 20 nT - 20 nT 0 nT 24 12 (midday) 24 Magnetic equator Intermediate latitudes High latitudes Figure 4.4 ? Schematic diurnal variation, showing the vertical component of the magnetic field through a 24 hour period at varying magnetic latitudes (Roux, 1980). Magnetic storms are caused by solar flares, where increased, irregular amounts of plasma are ejected towards the Earth. This violent interaction with the ionosphere causes irregular variations in the geomagnetic field, with amplitudes up to 1 000 nT (Telford et al., 1990). Variations occur rapidly (minutes to hours) and a storm may last several days before returning to the diurnal (Figure 4.5). Magnetic surveys are 53 discontinued during magnetic storms as the irregularity of the variation cannot be accurately measured and removed from field data. - 300 nT 0 nT 100 nT 24 12 (midday) 24 Figure 4.5 ? Schematic magnetic storm, showing the vertical component of the magnetic field. Note the erratic nature of the field and the change in scale from Figure 4.4. 4.1.5 Aeromagnetic Survey A high-resolution horizontal gradient aeromagnetic survey was flown by Fugro Airborne Surveys for Anglo Platinum in August/September 2002, in order to delineate dykes and faults in the area. The flight parameters chosen for the survey were as follows: ? Aircraft: JetRanger Bell 206 IIIB ? ZS-HWV helicopter ? Magnetometer: 2 x Scintrex cesium vapour, mounted on a horizontal boom with a horizontal gradient separation of 13 m (Figure 4.6). ? DGPS: NovAtel 3151R RealTime DGPS system ? Radar Altimeter: MRA Mk5 dual antenna altimeter ? Flight Speed: 40 m.s-1 (~145 km.h-1) 54 ? Data collection frequency: 10 Hz o Resulting a station spacing of ~4 m ? Line-spacing and trend: 50 m, 055? ? Tie-line spacing and trend: 500 m, 145? ? Flight Height: 20 m ? Draped survey (region is very flat) Figure 4.6 ? JetRanger Bell 206 IIIB, with horizontally mounted magnetometers (Fugro, 2005). The Scintrex cesium vapour magnetometers have an accuracy of 0.05 nT, the radar altimeter is accurate to within 12.5 cm and the DGPS has an accuracy of 1 m. A full review of these instruments, as well as aeromagnetic surveys, is beyond the scope of this dissertation, but is provided by Letts (2004). 55 4.2 Instruments The ground magnetic survey is faster than the ground gravity survey. Magnetic readings, using a proton-precession magnetometer, take less than five seconds and positions accurate to 10 m are sufficient for the ground survey. Hence, a handheld Garmin? GPS was used in place of the Trimble? DGPS. 4.2.1 GeometricsTM G856AX Proton Precession Magnetometer Telford et al. (1990) explains that proton precession magnetometers consist of a container of proton-rich fluid (e.g. kerosene, alcohol or decane) surrounded by a tightly-wound coil (Figure 4.7). Figure 4.7 ? Circuit diagram of a proton precession magnetometer (Telford et al., 1990). When the magnetometer is correctly orientated, the coil is approximately perpendicular to the Earth?s field. When the user takes a reading, current flows through the coil, creating a magnetic field within the coil. The magnetic moments of the protons in the fluid align themselves with this field. The current is then stopped and the proton moments try to realign themselves with the geomagnetic field. They cannot align themselves immediately and start to precess about the Earth?s magnetic field. The frequency of this precession (i.e. the Larmor frequency, fL) is directly proportional to the strength of the geomagnetic field: 56 B G B l m f RL rr ?? 22 == (4.2) where m is the proton magnetic moment, l is the proton angular frequency, B r is the ambient field (in this case, the geomagnetic field) and GR is the gyromagnetic ratio for protons (a constant, G R RR = 0.267513 nT-1s-1, Corner (1993)). Hence: LLR ffGB )002.0487.23(2 ?== ? (4.3) The Larmor frequency is easily measured since the precessing protons induce an alternating current in the surrounding coil, now acting as a pick-up coil. The alternating current continues for 2-3 seconds before protons become disorientated by ordinary molecular motion. The current is measured by comparing the frequency to a high precision oscillator. The G856AX proton precession magnetometer, used in this survey, is accurate to 0.1 nT and only measures the total field intensity. It is important to note that it becomes inaccurate in areas of high magnetic field gradient, since the precession signal decays rapidly (Geometrics, 1995). The magnetometer also cannot be used close to any AC equipment (e.g. power lines). Data is recorded and stored in the instrument?s memory during the survey and downloaded at the end of the survey. 4.2.2 Garmin? GPS 12 XL Theory regarding GPSs has been provided in Chapter 2.3: Differential Global Positioning Systems (DGPS). The handheld Garmin? GPS 12 XL was programmed to direct the user to predetermined waypoints 50 m apart, with a magnetic reading taken every 10 m (estimated by pacing). The position of each magnetic station was determined during office processing by interpolating between the GPS waypoints. Station location using this method has an estimated accuracy of ~2 m, based on handheld GPS location errors, human errors and natural obstacles (e.g. vegetation). 57 4.3 Survey Design Two magnetic traverses were surveyed, following the same route as the ground gravity traverses (i.e. line length ~5500 m, line spacing 800 m). Whilst taking readings, the user of the magnetometer had to be ?magnetically-clean?, i.e. no magnetic or electronic equipment on their person or near them (including watches, cell-phones). Magnetic Base Station Line 1 Line 2 Figure 4.8 ? Ground magnetic traverses, relative to the outline of the AGG data. A base station was set up, in an area of low magnetic field gradient, to automatically record the diurnal variation. The Line 2 was unfortunately shortened as the road 58 which the gravity survey traverses runs parallel to an electric power line. The electromagnetic field generated by the power line produced spurious readings in the magnetic data, hence the line was abandoned. 4.4 Quality Control Base station data were checked to ensure no magnetic storms had occurred during the survey and that the diurnal was smooth and regular (i.e. <2 nT/min, Figure 4.9). Figure 4.9 ? Plot of base station data for the two days spent on the ground magnetic survey. Note: surveys were approximately one year apart, accounting for the ~100 nT shift from Day 1 to Day 2. This shift was corrected for in the survey data. In the field, repeat readings of sections of the traverses were taken each day. 200 m of the previous day?s readings were recorded before the start of each day?s survey. A 200 m-repeat line was also done at the end of each day, repeating the first part of the day?s survey. Due to the surveys occurring approximately one year apart, these 59 repeat lines corrected for instrument drift (difference ~100 nT) as well as levelling of the day?s survey (due to the diurnal variation). If the traverse happened to pass any magnetic objects (e.g. fences, vehicles, power- lines), a note of the station number was made in the field book and the reading could be manually removed if it caused any spikes or problems in the data. 4.5 Processing Once the dataset has been collected, it is downloaded onto computer and the diurnal is removed to produce the final data set. 4.5.1 Downloading Data Magnetic data are downloaded from the field magnetometer and base station magnetometer to computer via a serial connection. The DOS-based program Magloc (?magloc.exe?) was used to import the readings. The base station data and field data are given different file extensions, in order to keep the two datasets separate. The output ASCII file contains columns showing line, station number, time and magnetic reading. The GPS data used to navigate to the 50 m-stations were used for positioning. 4.5.2 Corrections & Interpolation of Field Readings During the downloading process, Magloc uses the base station data file to automatically remove the diurnal variation from the field readings. The global magnetic field, originating from the Earth?s roughly dipolar field, can also be removed from the data. The model for the International Geomagnetic Reference Field (IGRF) is used to calculate the main component of the magnetic field at a user- 60 defined position on the Earth, which is then removed from the magnetic readings (Maus and MacMillan, 2005). The corrected magnetic data and GPS data were imported into Excel and each GPS station was correlated with the corresponding magnetic station. The co-ordinates of the magnetic stations in-between GPS readings were interpolated. Once completed, the stations? positions could be plotted. 61 CHAPTER 5: DATA ENHANCEMENT FOR INTERPRETATION 5.1 Introduction A full geophysical suite of data was used to build a 3D geological model of the region, covering the AGG survey area. Once created, it would serve as a starting model upon which gravity inversions could be run to improve the model?s density and layer geometry. In order to run inversions using vertical component (Gz) data across the largest possible area, the airborne FTG gravity dataset had to be converted to a Gz dataset. This chapter explains how the ground Gz data and airborne FTG data were processed and merged to form a single Gz dataset. The method of vertical continuation of data (ground and airborne) is explained, followed by the method for converting the FTG data to Gz. Then, the method of merging airborne and ground data to a single Gz dataset is described. The aeromagnetic data is also used to interpret the position of geological and man- made features, to determine if any of these features are related to gravity anomalies. An explanation of sun-shading filters is provided and applied to enhance features (e.g. lineaments, geological boundaries). Finally, Euler deconvolution is described as a method of automatically locating magnetic sources using the field data. 5.1.1 Geophysical Filters Geophysical filters are routinely applied to data to manipulate the signal in order to extract different kinds of information. The two common means of filtering data are convolution and Fourier transforms. Filters that attenuate high frequency data and accentuate low frequency data are termed low-pass filters, while filters that attenuate low frequency data and accentuate high frequency data are called high-pass filters (Heath, 2007). 62 In the convolution method, data are filtered in the space domain. A moving window, with an associated filter response, w(x,y), is convolved with the data, D(x,y). The output of the convolution is a new dataset D?(x,y). In order to more fully explain this, a smoothing filter is now described as an example. Following Figure 5.1, the data point in the centre of the 3 x 3 window (i.e. ?4?) undergoes a mathematical process, relating to the window and surrounding data points. In this case, each of the nine points are multiplied by the corresponding value in the window (i.e. 1), summed and then divided by the number of points in the window (i.e. 9). This produces an average of the surrounding points. The window moves on to the next point (i.e. ?7?), repeats the mathematical averaging, and continues through the dataset, thus producing a newly filtered dataset. Spatial filters are capable of altering the data in many ways, e.g. smoothing data, increasing detail, and highlighting lineaments and circular features. D(x,y) w(x,y) D?(x,y) Fourier transforms require the dataset to be converted into the frequency domain, where a filter operation is applied before converting the dataset back to the spatial domain. This method is capable of enhancing or suppressing different frequencies. The Fourier transform of a function, f(x,y), is defined as: ?? ?? ?= dxexfuF uxi ?2)()( (5.1) 9 4 8 6 8 3 4.9 6.2 6.9 5.6 1 4 7 9 2 1 4.0 1 3 7 8 7 6 3 5 5 9 0 2 Figure 5.1 ? Schematic example of a spatial 3 x 3 smoothing filter kernel. Note, the 6 x 4 input grid is reduced to a 4 x 2 filtered grid, after Cooper (2004c). 6.3 6.0 4.9 1 1 1 1 1 1 1 1 1 63 while the inverse Fourier transform is: ?? ?? = dueuFxf uxi ?2)()( (5.2) In two dimensions, the Fourier transform and inverse transform become: ? ?? ?? ? ?? +?= dxdyeyxfvuF vyuxi )(2),(),( ? (5.3) ? ?? ?? ? ?? += dxdyevuFyxf vyuxi )(22 ),(4 1 ),( ?? (5.4) In this manner, a Fourier transform pair is created, represented as (Cooper and Olavsdottir, 2004). Once in the frequency domain, the function can be multiplied by a number of operators to, e.g. smooth data, increase detail, apply derivatives and vertically continue data. ),(),( vuFyxf ? 5.1.2 Processing Software Packages Two programs were used to enhance the data for interpretation: ? Geosoft Oasis Montaq? allows viewing of multiple overlying datasets, enabling powerful integration of various geophysical datasets. It is also capable of database management, co-ordinate transformations, gridding, filtering (in the spatial and frequency domain) and mapping the data. ? Intrepid has produced new processing techniques and visualisation of FTG data, including gridding, filtering and interpretation using multiple FTG components in the process. 5.2 Vertical Continuation Continuation is a filter which calculates the potential field at heights varying from the altitude it was measured. This proves very useful when comparing ground potential field data to airborne data, smoothing data, and merging airborne surveys that have 64 been flown at different heights. The filter operator for profile data in the frequency domain is: uzeuFuF ?= )()(' (5.5) where u is the frequency and z is the distance above or below the point of measurement (z being positive or negative, respectively), (Cooper and Olavsdottir, 2004; Dean, 1958). Obviously, upward continuation, i.e. moving the grid further from the potential source, smoothes the data. Conversely, downward continuation is a high-pass filter, which increases both detail and noise. The noise may either be a result of increased noise of the original dataset or introduced by artefacts from Fast Fourier Transform (FFT) operations (Cooper, 2004d). Due to the increased noise levels, a low-pass filter is commonly applied after downward continuation. 5.2.1 Lowering Noise Levels during Downward Continuation Downward continuation is highly unstable and the result quickly diverges below a given height for a dataset. However, a number of methods may be applied to the dataset to reduce noise levels and downward continue the potential field data in a stable manner (Cooper, 2004d). The first method involves minimising the noise introduced by the FFT by performing parts of the downward continuation process in the space domain (since continuation in the space domain is less susceptible to noise). This is most effectively done by calculating the horizontal derivative of the data signal, f, in the space domain (Cooper and Olavsdottir, 2004): x ff dx df xxxx ? ?= ???+ 2 (5.6) The downward continuation of this derivative and simultaneous integration of the data may then be carried out in the frequency domain to obtain a smoother downward continued dataset. 65 The horizontal integration in the frequency domain is: niuuFuF ?= ))(()(' (5.7) where n is the order of horizontal integration (Blakely, 1995). Hence, in order to simultaneously downward continue and horizontally integrate the data, Equations 5.5 and 5.7 are combined and the following operator is applied: n uz iu e uFuF )( )()(' ? = (5.8) By downward continuing the horizontal derivative, frequencies near the Nyquist frequency (fN = 1/2?x) are not as enhanced and edge effects are reduced. In applying this operation to map data, artefacts are introduced at 90? to the derivative direction. To avoid this, the second derivative of the data should be calculated in the space domain (using Laplace?s equation) and downward continued (Cooper, 2004d). The second approach compares the data as downward continued in the space domain and in the frequency domain. By subtracting the space domain downward continued dataset from the frequency domain downward continued dataset, the amount of FFT- induced noise may be measured. The difference between the two datasets can be upward continued and added to the original dataset, thereby compensating it. This compensated dataset is then used in all subsequent continuations (Cooper, 2004d). Finally, Cooper (2004d) suggests using inversion to control FFT-induced noise and any other noise forming part of the original data. Considering three functions fa, fb and fc of a dataset y: )())(( yfyff cba = (5.9) where fa and fc are known functions, while fb has to be determined. So, if fa represents the upward continued data, by a distance h, and fc represents the unchanged data, then fb must be the downward continued data by the same distance h. 66 The downward continued data may be found by solving (5.10) here A is the gradient matrix, k is the damping factor, I is identity matrix and e is e misfit between the original data y and the upward continued data. field, as e noise levels are much lower. In other words, noise levels of directly measured Tzz In order to calculate the horizontal derivative and v frequency domain, the following operators are app for the following using least- squares: ekyfb TT AIAA 1)()( ?+= w th 5.3 Integrating FTG data to Gz data Derivatives of datasets are useful in defining detail and edges in the dataset, since they act as high-pass filters. However, they also increase the level of noise. This is one of the benefits of directly measuring the FTG components of a potential th data would be less than a computed vertical derivative of measured Gz data. ertical derivative of a signal in the lied, respectively: niuuFuF ))(()(' = (5.11) n uuFuF )()(' = (5.12) From the derivatives, it is clear that the operator fo where u is the frequency and n is the order of derivative (Blakely, 1995). r vertical integration of a dataset in the frequency domain is: n uuFuF ?= )()(' (5.13) here n is negative. The integration acts a low pass filter, smoothing data and w emphasising deeper sources (Cooper and Olavsdottir, 2004). Standard integration of Tzz to Gz is possible, providing a suitable result. However, it is possible to exploit some or all of the tensor components during processing. 67 Fitzgerald and Holstein (2006) follow Vassiliou (1986)?s work in obtaining a transform function in the frequency domain to obtain the free-air vertical component, Gz. Vassiliou (1986) shows that any one of the following combinations of components may be used in the integration: itzgerald and Holstein (2006) implemented these procedures into their work. This is roprietary software, hence the only details published are those provided. avity data are direct easurements of the gravity. Hence they are intrinsically more accurate than the od does handle high radients, but may not appear as smooth as the cubic method. Areas lying outside of the correction width are not altered in any way (Geosoft, 2006). )(aTzz )( )( c b TTT TT zzyzxz yzxy ?? ? (5.14) F p 5.4 Merging Airborne and Ground Datasets Although the preceding procedures result in a Gz dataset, the highly accurate ground gravity data are not actually included. The ground gr m calculated Gz data and should be included in the final dataset. The upward continued ground gravity dataset and the airborne Gz dataset were merged using a ?grid knitting? function, called suturing. The suture method was applied where a line is defined, along which the two grids will be joined. The mismatch between the grids is corrected for by adjusting the grids on either side of this line. The suture method uses a circle, radius equal to a user-defined correction width (e.g. 25 m), which moves along the path, interpolating between the two grids and producing a smooth transition (Figure 5.2). The interpolation method may be linear, cubic or Akima spline. The linear method uses a basic linear variation between the two grids and the cubic method fits a minimum curvature surface, but does not handle areas of high gradient. The Akima spline meth g 68 R Figure 5.2 ? The suture method, with circle (radius: R) traversing the boundary of the ground gravity survey, interpolating between the ground gravity and calculated airborne Gz data. Note: this figure does not show the final sutured image. The amount of correction applied relatively to each grid is determined by weighting between each grid. A value of 0 to 1 can input, with 0.5 being the default value and allowing equal corrections between the grids. A value of 0 applies all corrections to the first chosen grid, while a value of 1 applies all corrections to the second chosen grid. This method allows one grid to remain unchanged in the merging process (Geosoft, 2006). 5.4.1 Final Gz Dataset The Tzz data was converted to Gz data, with a statistical approach applied to obtain an equal mgal range between the ground gravity data and airborne converted Gz data. The dataset was further improved in the region of the ground gravity survey by 69 merging the upward continued ground gravity dataset with the airborne Gz data. This is the final Gz dataset used in the inversion process. 5.5 Sun-shading To aid interpretation of the gravity data, aeromagnetic data was analysed to see if features in both datasets correlated. The aeromagnetic data was provided already processed, so only two data enhancement techniques were applied to it: sun-shading and Euler deconvolution. Sun-shading is a filter which enhances or suppresses lineaments. It can be applied to any gridded data, see topographic data (Figure 6.3, Figure 6.23) and gravity data (Figure 6.14). Pelton (1987) considers any 2D data as topography (e.g. positive amplitudes represent hills and negative amplitudes represent depressions) with the sun placed at a particular height and azimuth. The reflectance is calculated for each data/gridded point: 2 0 2 0 22 00 11 ..1 QPQP QQPP R +++++ ++= (5.15) with the sun at inclination ? (angle from the vertical) and azimuth ? (angle anticlockwise from east). P and Q are the horizontal gradients of the data and ?? tancos0 ?=P and ?? tan0 sin?=Q . Linear features approximately perpendicular to the azimuth are highlighted, while those approximately parallel to the azimuth are suppressed. When the sun?s elevation is taken as 90?, all the edges around objects are enhanced; hence sun-shading may be used as an edge detector. 70 5.6 Euler Deconvolution Euler deconvolution is a semi-automatic interpretation method for potential field data. It was first described by Thompson (1982) along magnetic profiles and later extended to gridded magnetic data by Reid et al.(1990). The method quickly estimates the location and depth of potential field bodies using Euler?s homogeneity relation. 5.6.1 Theory Thompson (1982) considered a function in the Cartesian co-ordinate system, f(x,y,z), with the plane of observation z = 0 and z being positive downwards. The function is homogenous of degree n if the following two equations are satisfied: ),,.(),,( zyxfttztytxf n= (5.16) nf z f z y f y x f x =? ?+? ?+? ? (5.17) where Equation 5.17 is known as Euler?s equation. Next, consider a function f(x,y,z) having the general form: Nr Q zyxf =),,( (5.18) where Q is independent of (x,y,z), and N = 1,2,3,? (also known as the structural index). To make the above equation homogeneous, n = -N. Many simple magnetic point sources have this form of equation (Thompson, 1982). 2/1222 )( zyxr ++= For a point source (e.g. point mass, magnetic dipole) located at a point (x0,y0,z0), the total magnetic intensity, ?T, is given as: ( ))(),(),(),( 000 zzyyxxfyxT ???=? (5.19) Hence, Euler?s equation can be presented as: TN z T zz y T yy x T xx ??=? ???+? ???+? ??? )()()( 000 (5.20) where the derivatives may be directly measured or calculated. The magnetic field, T, may be separated into the source field (?T) and regional field (B) as follows: 71 TBT ?+= (5.21) and so, Equation 5.20 can be written as: )()()()( 000 TBNz T zz y T yy x T xx ?=? ??+? ??+? ?? (5.22) This equation changes in the case of a magnetic contact, where the structural index is zero and an offset of A is introduced (Reid et al., 1990): A z T zz y T yy x T xx =? ??+? ??+? ?? )()()( 000 (5.23) The basic Euler convolution has been developed significantly since 1982. Programs exploiting the ?tightness? of the vertical derivative are now common (Hsu, 2002), allowing for improved source resolution. Cooper (2004a) provides an improved method for calculating the horizontal and vertical derivatives, using the second horizontal derivatives in the space domain and Laplace?s equation, to improve the source results. In the same paper, he recommends a computationally efficient method to determining invalid Euler solutions. Euler deconvolution of gravity FTG data has also been developed (Zhang et al., 2000), which exploits all the components in constraining the Euler solutions. Finally, the Euler method has been extended to obtain source parameters (e.g. dip, susceptibility) of simple magnetic objects (Cooper, 2006; Mushayandebvu et al., 2001). 5.6.2 Structural Index The structural index (SI) is defined as the measure of the rate of change with distance of the field, for a given geometry (Reid et al., 1990), e.g. the field surrounding a vertical pipe has fall-off proportional to the inverse-square, hence its structural index is two. The ?rule-of-thumb? for choosing structural indices is presented in Table 5.1. 72 Table 5.1 ? Structural indices for simple potential field bodies (after Thompson (1982) and Geosoft (2006)). Structural Index (N) Basic magnetic model Basic gravity model 0 Contact, step Dyke, sill, step 1 Dyke, sill Pipe 2 Pipe Sphere 3 Sphere N/A The rejection of invalid Euler solutions is commonly associated with maximum user- defined tolerances, i.e. solutions with error estimates smaller than, say, 15% are allowed. A smaller tolerance with the correct structural index will result in more reliable results but the fewer solutions. Hence, structure may be poorly delineated. However, high tolerance and poor structural index choice will result in poorly defined solutions, masking the reliable solutions. This trade-off results in a high number of reliable solutions to occur with some spurious results. 5.6.3 Application of Euler Deconvolution Euler deconvolution takes the following steps to obtain solutions for a gridded dataset (Reid et al., 1990): 1. Perform clean vertical and horizontal calculations of the gradients ( xT ?? , yT ?? , zT ?? ), or use the measured gradients. 2. Pass a square window over the gradient dataset. This window should be at least 3 x 3 grid points, although Reid et al. (1990) recommends a window of 10 x 10 to produce efficient, reliable results. Obviously, the smaller the window, the higher the resolution of the results. 3. (a) For non-zero structural indices: Use all the points in the window to solve for Equation 5.22. The resulting estimate gives a source position (x0,y0,z0), over the central point of the window, and background value B, using Moore- Penrose inversion. E.g. for a 10 x 10 window, 100 equations must be solved which gives the standard deviation of the four unknown variables. If the standard deviation is less than a user-defined tolerance, the solution is allowed. 73 3. (b) For a structural index of zero: Use Equation 5.23 in solving for the source position and offset, A. Continue as for step 3 (a). 4. Move the window and repeat steps (2) and (3) for all possible window positions. 5. Plot a map of the solutions, with each solution?s (x,y) plotted with a symbol/colour proportional the depth, z (Reid et al., 1990). 5.6.4 Euler Deconvolution Applied to Profile Data Although theory for Euler convolution over gridded data has been presented, it is useful to demonstrate Euler solutions over a synthetic magnetic profile across two dykes (Table 5.2, Figure 5.3). Table 5.2 ? Parameters of synthetic dykes used in Euler deconvolution. DYKE A DYKE B Depth 20 m 40 m Width 15 m 30 m Mag. Susc. Contrast 0.05 SI 0.05 SI Dip 75? 85? Strike 1000 m 1000 m Depth Extent 1000 m 1000 m 74 nT nT/m S N A B m Figure 5.3 ? Euler solutions over a synthetic dyke S-N profile using computer program Euler (Cooper, 2004b). Upper panel: Original data (solid black line) and Reduction-to-Pole (RTP, dashed red line) calculated from original data. Middle panel: Horizontal gradient (solid black line) and vertical gradient (dashed red line) data calculated from original data. Lower panel: Euler solutions for a structural index of 1 (i.e. green dots), with overlying position of synthetic dykes. Note the spread of the solutions, both horizontally and vertically. Cooper (2004b)?s program presents the original data over the two dykes together with the Reduction-to-Pole data. RTP transforms data at mid- to low-latitudes to what the profile would look like if measured at a magnetic pole. It is a common method of removing the dipole effect of magnetic bodies, allowing the user to see the approximate position of the causative body. The horizontal and vertical gradients are also calculated (as per Step 1, above). The Euler solutions are presented, showing a large vertical and horizontal spread, even though the data is perfect and noise-free. In spite of the spread, the solutions provide an estimate of position and depth of the dykes. 75 CHAPTER 6: GEOPHYSICAL DATASETS A number of geological and geophysical datasets, ranging from regional (>100 km) to local (<10 km) in extent, were obtained in order to create a 3D geological starting model for the area covered by the AGG survey. A second 3D gravity model of the overburden/bedrock contact was also created for the ground gravity region. Datasets from the ground gravity survey (grid and traverses), AGG survey (grid) and magnetic survey (traverses) are presented, as well as regional gravity and magnetic data, and local geological, topographic, QuickBird imagery, borehole, seismic and rock property data. 6.1 Geological Map The 1:250 000 ?2526 Rustenburg? geological map was used to show lithological changes, structural features (especially faults) and the traced outcrop of the Merensky Reef and UG2 chromitite layer (Figure 6.1). 76 f Aeromagnetic survey outline Rustenburg Fault Rustenburg Fault: conjugate set Airborne gravity survey outline Suboutcrop: Merensky Reef Suboutcrop: f f f f LEGEND Pilanesberg Complex Bushveld Complex Transvaal Supergroup UG2 f Figure 6.1 ? Geological map (after Walraven (1981)) showing ground gravity survey (grid and traverses), AGG survey and aeromagnetic survey. The position of the ground gravity, AGG and aeromagnetic surveys are also visible on the map. The aeromagnetic data overlaps most of the AGG data, except for two <6 km2 areas, along the south-western edge and the eastern corner of the AGG region. The geological map shows the Rustenburg Fault trending NNW-SSE past the AGG region and a conjugate set of the Rustenburg Fault, trending NW-SE and extending ~7.5 km across the AGG region. This conjugate fault presents a contact between the Bushveld Complex and Transvaal Supergroup which is different to the more common Marginal Zone contact. 77 6.2 Regional Gravity Data Regional Bouguer ground gravity data (Figure 6.2) are obtained from the South African Council for Geoscience. The data were collected with random spacing (i.e. no grid alignment) and with stations ranging 1-5 km apart (average ~2 km). The dataset shows the position of the eastern, western and northern limbs of the Bushveld Complex and the Pilanesberg Complex. Northern Lobe Eastern Lobe 100 km Pilanesberg Complex mgal Figure 6.2 ? Regional Bouguer gravity data over the Bushveld Complex. (right) Northern, western and eastern limbs of the Bushveld Complex are clearly visible as gravity high anomalies. (left) Zooming in on the western limb of the Bushveld Complex, showing the outline of the AGG survey area (solid box), and the approximate position of the Transvaal Supergroup contact (dashed line) and Pilanesberg Complex (dotted line). Data courtesy of the South African Council for Geoscience. The high density rocks of the mafic phase of the Bushveld Complex are clearly represented by the gravity highs. The complex is surrounded by the less dense sediments of the Transvaal Supergroup and Karoo Supergroup cover sequences, represented by gravity lows. Another observed gravity low is the circular outcrop of Pilanesberg Complex alkaline rocks. 78 6.3 Ground Gravity Dataset The ground gravity grid and traverses produced a high resolution Gz dataset. High frequency anomalies in the Gz field (e.g. variations in overburden thickness) were not measureable in the airborne survey, since the data are collected at a flight height of 80 m above the ground ( ). The free-air and Bouguer gravity data were gridded to produce two maps of the data ( 3?? RTij Figure 6.3a and Figure 6.3b, respectively). The DGPS data for each each gravity station in the grid was used to produce a Digital Terrain Model (DTM) of the ground gravity survey region (Figure 6.3c). The elevation descends from ~1060 m in the south to ~1020 m along the Elands River in the north of the survey region. The steepest gradients are seen next to the dry river beds, acting as tributaries to the Elands River after heavy rainfall. Besides these river beds, the land is very flat with a low gradient of ~75:1 traversing N-S. 79 a) b) c) Elands River Dry river beds Figure 6.3 ? a) free-air gravity map of the ground survey grid, with 100 m x 100 m grid and 30 m x 30 m grid in the eastern corner. b) Bouguer gravity map of the ground survey grid, with similar grid positions. c) DTM of the ground gravity survey region. The free-air gravity data shows a strong E-W regional trend, from highs of ~17 mgal to lows of ~10 mgal. The effect of elevation is significantly reduced by the free-air correction (i.e. gravity signal is not correlated to elevation), except in areas of steeply changing elevation around the Elands River and dry river beds. The Bouguer gravity data has a strong ENE-WSW regional trend, from highs of ~-100 mgal in the ENE to lows of ~-107 mgal in the WSW. The effect of elevation changes relating to the dry 80 river beds is largely removed or significantly reduced. This effect cannot be completely removed as it is based on a theoretical infinite slab, which does not occur in nature. It is common practise in potential field modelling to remove the long wavelength features of the regional field and model the short wavelength features of the residual field. However, the modelling program used in this project uses the original free-air gravity data in the inversions (discussed in Chapter 7.4: VPmg Inversion Methodology). The accuracy of the gravity survey was determined from the Root Mean Square (RMS) error of the repeat station readings (Table 6.1). Table 6.1 ? Statistics of ground gravity survey. Total number of stations 1302 Number of repeat stations 182 Percentage of stations repeated 14.0% Total number of readings 1503 Number of repeat readings 389 Percentage of repeat readings 25.9% Maximum repeat error 0.210 mgal Mean repeat error 0.012 mgal RMS error 0.015 mgal A high proportion of repeat readings were recorded (25.9% of the readings were repeats) obtaining a RMS error of 0.015 mgal. This is sufficiently close to the advertised error of 0.005 mgal (Scintrex, 1995) to be acceptable. The maximum repeat error occurred on an extremely windy day over ploughed, unconsolidated land, with the gravity meter experiencing vibrations during recording. 6.3.1 Overburden Thickness Borehole data in the study region were used to provide the overburden thickness in the hole. The bedrock contact surface was approximated by joining these points in 81 the borehole data. By subtracting this contact surface from the DTM, a map of the approximate overburden thickness is created (Figure 6.4). Figure 6.4 ? Thickness of the overburden approximated from the DTM and borehole data. Borehole positions represented by crosses. The region is generally very well constrained, with thicknesses ranging from 0 m (i.e. outcrop) to ~12 m. Along the SW edge, the data is poorly constrained by boreholes and relies on interpolation with boreholes outside of the ground gravity study region. 6.3.2 Gravity Traverses Bouguer values of the two gravity traverses (Line 1 and Line 2) were plotted with corresponding elevations to produce two profiles (Figure 6.5). The profiles changed orientation, as one traversed from SW-NE, to avoid steep topography and dense vegetation. The predominant direction of each of the traverses lies at a bearing of 040?. In order to model the data in 2D, the data had to be projected onto a straight line. The data were projected onto an extended line bearing 040?, with perpendicular 82 vector lines projecting each data point. Although this distorts the data, the traverses were mainly used to model long wavelength features which would not be significantly affected by this distortion. Gravity high - IRUP SW NE Figure 6.5 ? Bouguer gravity profiles (station spacing: 150 m) for Line 1 (solid black) and Line 2 (solid red), and corresponding elevations for Line 1 (dashed black) and Line 2 (dashed red), running from NE to SW. The gravity anomaly in Line 2 is related to a traversed IRUP. The trend seen in the gridded ground data is also present in the traverses (i.e. a change from high to low from ENE to WSW. Line 1 shows a relatively constant gradient, while Line 2 shows a gravity high anomaly as it traverses an IRUP. Although the gravity high also traverses an elevation high, the two do not correspond, as the gravity high clearly begins ~300 m SW of the start of the elevation high and ends ~750 m to the NE of the elevation high. 6.3.3 Gravity Traverses ? Vertical Continuation Airborne profiles of the free-air gravity data (see Chapter 6.4.4: Final Gz Dataset), matching Line 1 and Line 2, were extracted with 150 m station spacing. The airborne survey was flown at 80 m, hence, in order to compare ground and airborne data, the 83 airborne profiles were downward continued by 80 m and smoothed in Cooper (2000)?s program, SignProc. This program uses an FFT to transform the data in the frequency domain (Equation 5.5). Downward continuation accentuates high frequency data as well as noise; hence the program automatically applies a low-pass filter to prevent data getting too noisy. Artefacts from the FFT are seen as edge effects in the downward continued profiles. The ground, airborne and downward continued free-air Gz gravity data for each line are compared in Figure 6.6 and Figure 6.7. SW NE Figure 6.6 ? Line 1: Profiles of free-air ground gravity data (black), airborne gravity data flown at 80 m (red), and airborne gravity data downward continued 80 m (blue). 84 FFT artefact SW NE Figure 6.7 ? Line 2: Profiles of free-air ground gravity data (black), airborne gravity data flown at 80 m (red), and airborne gravity data downward continued 80 m (blue). The anomaly in the downward continued profile in the NE is a result of edge effects from FFTs in the downward continuation. The gradient variation in the gravity data is relatively low for the majority of each profile. Hence, downward continuation only provides minor improvements to the data. As one traverses towards the SW, the correlation between the ground and airborne data degrades. This may be attributed to: ? Aliasing problems relating to the extraction of the profiles from the airborne data (airborne stations and ground stations not coinciding, variations in ground-/flight-line direction, and differences from projection to a straight line), and ? Small variations related to the statistical DC shifting of the converted Gz data (Chapter 5.3: Integrating FTG to Gz Data). 85 6.3.4 Upward Continued Ground Gz Data The ground Gz data was upward continued 20 m, 40 m, and 80 m (the flight height of the AGG data) using Equation 5.5 (Figure 6.8). The smoothing effect of the upward continuation is clearly evident. a) Ground b) z = 20 m c) z = 40 m d) z = 80 m Figure 6.8 ? a) Original ground free-air Gz data, upward continued to b) 20 m, c) 40 m and d) 80 m (units: mgal). Upward continuation is commonly used in airborne gravity and AGG surveys, comparing a ground gravity survey with the corresponding portion of an airborne 86 gravity grid (Dransfield and Lee, 2004; Hatch et al., 2006; Hinks et al., 2004; Murphy, 2004). Ground gravity surveys are more accurate and closer to the geological source. Hence, comparing two equivalent datasets allows an estimate of the noise and error in the airborne data. The filter can also be used to approximate the regional trend in a dataset, an alternative to a smoothing convolution filter. 6.4 Airborne FTG Gravity Dataset The 10 km x 10 km airborne FTG gravity survey (Figure 6.9) was flown with previously provided parameters. The five independent components and the Tzz data are presented in Figure 6.10. The equivalent ground station spacing is approximately 60 m, hence, any wavelengths shorter than 120 m are aliased into longer wavelength features (Reid, 1980). Transvaal Supergroup contact IRUPs Eo Figure 6.9 ? Final levelled airborne FTG gravity grid, showing Tzz data and xy-recording positions (showing lines and tie lines). Data courtesy of Anglo Platinum. 87 Txx Txy Txz 10 km Tyy Tyz E?tv?s Tzz Figure 6.10 ? The individual tensor components of the AGG data (units: E?tv?s). The Tzz data clearly delineates the Transvaal Supergroup contact with a prominent low gradient anomaly. The Tzz data also defines several high gradient anomalies in the Bushveld Complex which correlate to magnetic highs and are interpreted as high density IRUP bodies. Finally, the gradient highs close to the SW edge of the Tzz data are seen as topographic highs in the Transvaal Supergroup. 88 6.4.1 Downward Continued AGG Data The airborne Tzz data was downward continued 10 m, 20 m, and 40 m using Equation 5.5 (Figure 6.11). a) z = 80 m b) z = 70 m Eo c) z = 60 m d) z = 40 m Figure 6.11 ? Original AGG Tzz data, flown at 80 m above ground level downward continued to b) 70 m, c) 60 m, and d) 40 m (units: E?tv?s). It is clear how quickly the data becomes unstable when using a conventional FFT downward continuation algorithm. Due to instability in downward continued data, even when applying smoothing filters, it was decided to upward continue the ground 89 Gz data to 80 m above ground level and merge it with the AGG data at this height. Hence, the inversions would use a merged Gz dataset at this altitude. 6.4.2 Integrated AGG Data Fitzgerald and Holstein (2006)?s process, using Equation 5.14 (c), was used to integrate the AGG dataset (Figure 6.12a). However, the integration does not result in a final product ? the z-values in the dataset must be divided by 10 000 to obtain mgal units. Following this, a statistical method is applied to enforce a DC shift on the data. The DC shift was calculated using the upward continued ground gravity data (to 80 m) and fitting the equivalent portion of the integrated grid using least-squares (Figure 6.12b). a) b) 10 km mgal Figure 6.12 ? a) The integrated airborne data (units: arbitrary), and b) the integrated airborne data, divided by 10 000 and with a DC shift applied to fit the upward continued ground data (units: mgal). Note the position of the ground gravity grid is outlined in b). 6.4.3 Error in Integrated Airborne Data The error in the newly created airborne Gz dataset can be calculated by finding the difference between the upward continued Gz ground data and the equivalent section of the airborne Gz dataset (Figure 6.13). 90 a) b) -2.0 0.0 1.5 mgal Figure 6.13 ? a) Difference between the upward continued ground data and the shifted airborne Gz data, and b) histogram of the grid provided in a). Note the trend in the difference not resolved by the DC shift using least-squares. Qualitatively, the map of the difference between two datasets shows a strong regional gradient, from low in the west to high in the east. This indicates that the statistical least-squares fit of the two grids was not completely accurate. Quantitatively, the distribution of the misfit can be seen in the histogram of the data (Figure 6.13b). This shows a population approximately symmetrically centred about zero and with the majority of the positive and negative values extending to ? 1.500 mgal. The statistics of the histogram are provided in Table 6.2. Table 6.2 ? Statistics of histogram for the difference between upward continued ground Gz data and the equivalent section of the shifted airborne Gz data. Statistic Value (mgal) Minimum -2.062 Mean 0.010 Maximum 1.583 Standard deviation 0.870 The mean shows that the distribution is centred close to zero. Although the minimum and maximum are both greater than ? 1.500 mgal, the standard deviation is 0.870 mgal. The standard deviation of the histogram provides the accuracy of the airborne 91 Gz dataset, i.e. 0.870 mgal, compared to the accuracy of the ground data, 0.015 mgal. Although the airborne gravity accuracy is significantly greater than the ground gravity accuracy, it is still within presented estimates of airborne gravity errors on a fixed-wing system. These estimated accuracies vary between 0.2 ? 1 mgal (Fairhead and Odegard, 2002) to 2 mgal (Hwang et al., 2007). 6.4.4 Final Gz Dataset The vertically continued ground gravity data set (to 80 m height) was combined with the airborne Gz dataset to create a final Gz dataset, which would be used in the inversion process (Figure 6.14). The two Gz gravity grids were knitted together using the suture method. The suture parameters included Akima spline interpolation, a correction width of 250 m and default weighting of 0.5, which produced a smooth merge of the two grids. The upward continued ground data is smoother than the converted airborne Gz data. This is due to the high accuracy of the ground data relative to the converted airborne data, and artefacts from the conversion of FTG data to Gz data. 92 a) Gravity high - topography Figure 6.14 ? a) Final merged Gz dataset, 80 m height (units: mgal), and b) zooming in to the ground gravity grid section, with sun-shading from the NE applied. Note how much smoother the upward continued data is relative to the integrated AGG data. The contact between the Bushveld Complex and Transvaal Supergroup is clearly delineated from the rapid change in gravity. This is seen in the change from high b) 10 km mgal Gravity high - IRUPs BC/Tvl Supergroup contact 93 gravity values (>10 mgal) over the high density mafic intrusives of the Bushveld Complex to lower gravity values (<8.5 mgal) over the low density sediments of the Transvaal Supergroup towards the west. A number of clear gravity high anomalies exist. Some of these can be related to IRUPs from the aeromagnetic data, while one is specifically related to topography, with no magnetic signature. 6.5 Aeromagnetic Dataset The regional aeromagnetic dataset was obtained from the South African Council for Geoscience, flown with a line spacing of 1 km and flight height of 150 m (Figure 6.15). Major features in the data include the magnetite layers in the Upper Zone of the Bushveld Complex, the Pilanesberg Complex (positive anomaly), the NW trending Pilanesberg dykes (negative anomalies, believed to be associated with the Pilanesberg Complex intrusion) and the late-stage E-W trending dykes (positive anomalies) which cross-cut other features. Although regional features are seen, the wide line spacing and high flight height have caused the data to be aliased. The local aeromagnetic data has more defined resolution, showing positive anomalies relating to Iron Rich Ultramafic Pegmatoids (IRUPs), NW trending Pilanesberg dykes and the E-W trending dykes (both more clearly resolved). 94 NW Pilanesberg dykes Late-stage E-W dykes Pilanesberg Complex Magnetite layers in Upper Zone NW Pilanesberg dykes IRUPs Faulting Late-stage E-W structure Figure 6.15 ? (top) Regional aeromagnetic map of the Pilanesberg region, data courtesy of the South African Council for Geoscience. (bottom) Aeromagnetic map of the study region, sun-shaded from NE. Data courtesy of Anglo Platinum. 95 The local aeromagnetic map shows two polarities of dykes in the region as well as cross-cutting relationships. The NW trending Pilanesberg dykes show a negative polarity while the E-W trending late-stage dykes show a positive anomaly. This indicates a magnetic reversal between the intrusion of the Pilanesberg Complex and the dyke swarm intrusion, but is beyond the scope of this project. The cross-cutting relationships of the magnetic anomalies indicate age relationships between the magnetic bodies. It is clear that the IRUPs are cross-cut by the Pilanesberg dykes and E-W trending dykes, indicating they are the oldest intrusive bodies. The Pilanesberg dykes are cross-cut by the E-W trending dykes; hence they are older than the late- stage E-W trending dykes. Again, accurate age constraints are beyond the scope of this project. The E-W trending anomaly is actually a complicated structure, appearing as a series of parallel dykes, at times offset by near-parallel faults. 6.5.1 Cultural Noise Combining aeromagnetic data and QuickBird orthophotographs, it was also possible to interpret a number of man-made magnetic anomalies: ? The local village produced a fair amount of cultural noise, causing the dappled magnetic signature on the aeromagnetic map (Figure 6.16a), and ? Other infrastructure such as power-lines and dumps caused straight-lined noise in the aeromagnetic data (Figure 6.16b). 96 a) Village Power Line b) Power Line Village Dump Figure 6.16 ? Sun-shaded colour aeromagnetic maps zooming in on areas overlaying QuickBird orthophotographs. a) Village causing the dappled features in the magnetic data, b) Power lines and a dump causing other man-made noise in the aeromagnetic data. This data was compared to Euler solutions in deciding if solutions relate to geological anomalies or cultural noise. 6.5.2 Sun-shaded Aeromagnetic Data Sun-shading grey-scale aeromagnetic data from different orientations (Figure 6.17) is a fast method of realising varying trends of magnetic anomalies. 97 a) b) c) Figure 6.17 ? a) Grey-scale aeromagnetic data, b) the same data, sun-shaded from inclination 45?, azimuth 45? (origin N, clockwise rotation), highlighting features perpendicular to 45?, and c) sun- shaded from inclination 90?, highlighting the edges of all magnetic sources. Sun-shading from inclination 45?, azimuth 45? highlights dykes trending NW-SE as well as emphasising the prominent E-W trending structure. Sun-shading from 90? delineates the edges of dykes and structures and also clearly highlights the positions of IRUPs in the region (seen as circular ?pock-marks? in the data). This is a useful dataset to combine with Euler solutions to relate the solutions to magnetic sources. 6.5.3 Euler Deconvolution Applied to Aeromagnetic Dataset Euler deconvolution solutions were calculated for the aeromagnetic dataset, for a structural index of 1, in order to help delineate and provide depth resolution on dykes in the area. A tolerance of 10% was applied for a 10 x 10 window and no depth limit was applied (i.e. solutions could occur at any depth). 98 Over four million solutions were presented, with a minimum depth of 5 m and a maximum of 135 m. Depth and location uncertainties were calculated for each solution. The solutions were divided into four categories, according to realistic depths of dykes and with a limit on the spatial uncertainties (Table 6.3). Solutions below 45 m were discarded, based on the high vertical spread of Euler solutions and that the goal was to image dykes close to surface. The magnetic data, sun-shaded from 90? inclination and zoomed into an appropriate area, is presented with the four solution depth ranges (Figure 6.18). Table 6.3 ? Parameters for the four initial solution categories, based on depth and with user-defined spatial uncertainty restrictions. Depth Range (m) Location Uncertainty, dxy (%) Depth Uncertainty, dz (%) 5 ? 15 0 ? 5 0 ? 5 15 ? 25 0 ? 5 0 ? 5 25 ? 35 0 ? 5 0 ? 5 35 ? 45 0 ? 5 0 ? 5 99 Power lines Local village Dump Figure 6.18 ? Euler solutions (structural index: 1) overlaying aeromagnetic data (sun-shaded from 90? inclination) with solutions for varying depth ranges: green (5-15 m), blue (15-25 m), red (25-35 m) and black (35-45 m). Inset: Zoom of the solutions over two prominent dykes, with depths dominated by 5- 15 m solutions. The high spread in Euler solutions is clearly evident. However, over-laying the solutions on vertically sun-shaded aeromagnetic data allows dykes, especially the NW trending Pilanesberg dykes, to be delineated and depth-matched. Interestingly, the prominent E-W structure shows poor Euler solutions. Solutions for the circular IRUP anomalies are also present. This is because Euler deconvolution is an automatic process and models all anomalies in the data. However, a structural index of 1 was used in the calculations, i.e. for magnetic dykes or sills. More accurate estimates of the tops of pipes could be obtained by running Euler deconvolution with a structural index of 2. However, this was not done as the 100 small pipe-like IRUP anomalies are commonly too small to be imaged in the gravity data. Larger IRUP anomalies are seen in the gravity data, but have poor Euler solutions. This is because their shape does not approximate to a dyke (structural index: 1) or a pipe (structural index: 2). Cultural noise also provides Euler solutions, which need to be interpreted and discarded. This includes the local village, power lines and the dump. 6.6 Ground Magnetic Dataset Aeromagnetic data, ground-truthed by two ground magnetic traverses (Figure 6.19), were used to delineate dykes and IRUPs in the region. Line 1 Line 2 Figure 6.19 ? Ground magnetic profiles locations (black), with associated ground profile data (red), and gravity stations (black crosses), overlaying aeromagnetic data in the study region. 101 The magnetic traverses were planned to follow the same path as the gravity traverses. However, towards the south-western end of the lines, the gravity traverses followed roads that had electric power lines running parallel. The electromagnetic signal from the power lines caused false magnetic measurements, so the magnetic lines were curtailed. The lines cross two NW trending Pilanesberg dykes, both having a low- amplitude negative signature. They also both cross the dominating E-W trending structure, which has a dominating high-frequency, high-amplitude positive signal. 6.6.1 Magnetic Traverses Two magnetic profiles were created from the ground magnetic traverses (Figure 6.20). E-W trending structure IRUPs Pilanesberg dykes SW NE Figure 6.20 ? Magnetic profiles of Line 1 (black) and Line 2 (red), running SW to NE. A 100 nT shift was applied to Line 2, allowing for easier viewing. The significant anomalies in the SW relate to the E-W trending structure. Data for Line 2 was collected along a straight line. However, Line 1 changed orientation, following the corresponding gravity traverse. In the profile above, the 102 stations were not projected onto a straight line. However, the sections of the traverses which were modelled were extracted and projected onto a line bearing 040? (Chapter 8.1.2: Magnetic Modelling), similar to the gravity traverses. The high-frequency, high-amplitude signal of the E-W trending structure is apparent in Line 1 and 2. Line 1 shows an anomaly >7 000 nT, while Line 2 has a >4 000 nT anomaly. The high-frequency signal seen within the broader anomaly indicates complicated structure (i.e. multiple parallel dykes, faulting). North of these large anomalies, minor positive magnetic signals (~100 ? 200 nT) appear. They relate to circular anomalies on the aeromagnetic map, indicating the presence of small IRUP bodies. Finally, the negative Pilanesberg dyke anomalies (100 ? 1 000 nT) are clearly apparent. The data is magnetically quiet, with no cultural noise, seen by the flat gradient traversing towards the north-east. 6.6.2 Magnetic Traverses ? Vertical Continuation The aeromagnetic profiles, corresponding to Line 1 and Line 2, were extracted with 10 m station spacing. The airborne survey was flown at 20 m, hence the data are high resolution and downward continuation does not show major improvements. The data were downward continued 18 m (i.e. to 2 m above ground level, the equivalent of the height of the ground magnetic sensor) and smoothed, as for the downward continued gravity profiles, in SignProc (Cooper, 2000). The ground, airborne and downward continued data for each line are compared in Figure 6.21 and Figure 6.22. Since all three datasets correspond closely. A positive shift of 150 nT was applied to the downward continued airborne data and a positive shift of 300 nT was applied to the unchanged airborne data, allowing for easier viewing of the data. 103 Improved ground resolution SW NE Figure 6.21 ? Line 1: Profiles of ground magnetic data (black), airborne magnetic data downward continued 18 m, offset 150 nT (blue), airborne magnetic data flown at 20 m (red), offset 300 nT. Improved resolution Improved ground resolution SW NE Figure 6.22 ? Line 2: Profiles of ground magnetic data (black), airborne magnetic data downward continued 18 m, offset 150 nT (blue), and airborne magnetic data flown at 20 m, offset 300 nT (red). Note: change of y-axis scale from Figure 6.21. 104 Despite the already high resolution and the predominantly flat gradient, the downward continued profiles do show increased detail which corresponds to the ground data. The data has no cultural noise and is predominantly magnetically quiet. Three small positive anomalies (~150 nT amplitude change) are visible in Line 1?s ground data, relating to circular magnetic anomalies in the aeromagnetic map (possibly IRUPs). A major positive anomaly (~1 500 nT amplitude change) relating to the E-W trending dyke is also aliased in the aeromagnetic data in Line 1. Line 2 shows two similar small anomalies (<100 nT) and improved resolution towards the E- W trending dyke. 6.7 Topographic Data Accurate elevations are critical to gravity surveys, both in processing and modelling. The Digital Terrain Model (DTM) was provided by the aeromagnetic survey, using DGPS and radar altimeter readings (Figure 6.23). Elands River Tributaries Transvaal Supergroup quartzite hills Figure 6.23 ? DTM of the study region, determined from aeromagnetic survey (data courtesy of Anglo Platinum). Outline of AGG survey, ground gravity stations and geological boundaries included. Horizontal locations and elevations reported at cm level accuracy. Station spacing is ~4 m and line spacing ~50 m. 105 The DTM clearly shows the flat plains resulting from the easily weathered Main Zone of the Bushveld Complex. Other obvious features are the Elands River (and associated tributaries) in the north of the region, and the hills of the resistant Transvaal Supergroup quartzites in the south-west. 6.8 Seismic Data Davidson and Chunnett (1999) showed the importance of using high-resolution seismic data to map the Merensky Reef and UG2 in the Bushveld Complex. Their seismic profiles were recorded close to this study region (locations not published), and show a flat Transvaal Supergroup contact at approximately 2 000 m. Two 3D seismic surveys were carried out over the Styldrift region in 2001 and 2005, using Vibroseis trucks which provided seismic energy from 30-250 Hz over a 16 second period. Interpretations of the seismic data have shown faults accurate to the data resolution (i.e. 7 m) and have imaged layers as thin as 17 m (Gibson et al., 2006). Seismic profiles were used to pick the Transvaal Supergroup contact horizon, constrained by borehole data close to the Transvaal contact sub-outcrop (Figure 6.24). 106 Inline profiles Cross profiles Figure 6.24 ? Final Gz gravity map showing position of seismic NE-SW inline profiles and perpendicular cross profiles. The inline profile in bold is shown (left) with the trace of the Transvaal Supergroup contact (yellow). Data courtesy Anglo Platinum. Eight inline profiles were picked and joined by seven cross profiles, after which a 3D surface was interpolated and used to define the Transvaal Supergroup contact. This surface was used, together with the borehole markers, to build the 3D geological starting model in the vicinity of the seismic data. An important observation in the picked profiles is the significant change in slope of the contact. This ?bottoming-out? is confirmed in Davidson and Chunnett (1999)?s seismic interpretation and is also seen in Odgers et al. (1993)?s seismic interpretation in the southern Bushveld. 107 6.9 Borehole Data A large number of boreholes (277) were used to constrain the 3D geological starting model and inversion processes. Due to boreholes containing significant confidential information, the only data provided were the borehole collar locations, borehole deviations for the length of each borehole, and overburden and Transvaal Supergroup contact geology information. The target of the boreholes is the economic layers of the Critical Zone. Boreholes commonly extend <200 m deeper than the target, hence only 27 boreholes intersect the Transvaal Supergroup contact. These shallow boreholes exist close to the contact, in the south-western portion of the grid (Figure 6.25). a) 108 b) Boreholes intersecting basement contact Figure 6.25 ? a) Converted Gz free-air gravity map with borehole positions and outline of the ground gravity survey (borehole data courtesy of Anglo Platinum and PTM), b) 3D view of Transvaal Supergroup contact (interpolated from seismic profiles) and boreholes in the study region, with only 27 shallow boreholes intersecting the contact in the south-western region. Gravity inversions were carried out over two models: the region covered by the AGG data (a two-layer model consisting of the Bushveld Complex and underlying Transvaal Complex) and the region covered by the ground gravity data (also a two- layer model consisting of the overburden and underlying Bushveld Complex bedrock). Hence, the boreholes? geological data were edited, with geological markers for bedrock contact and Transvaal Supergroup contact. From the overburden markers, a simplified surface was modelled in the ground gravity grid (as shown). The borehole markers of the overburden and Transvaal Supergroup contact, as well as the bottom of the boreholes, were used as constraints in the inversion process. 6.10 Rock Property Data Densities of the two stratigraphic units, as well as simple realistic density ranges, were estimated from the South African Geophysical Atlas: Physical Properties of South African Rocks (Mar? et al., 2002) and density logs in Ashwal et al. (2005), (Table 6.4). 109 Table 6.4 ? Starting densities of the regional model (after Mar? et al. (2002) and Ashwal et al. (2005)). Stratigraphic Unit Starting Density (g.cm-3) Minimum Density (g.cm-3) Maximum Density (g.cm-3) Bushveld Complex (Main Zone) 2.950 2.750 3.150 Transvaal Supergroup 2.600 2.400 2.800 IRUPs 3.200 3.000 3.400 The starting density contrast between the Bushveld Complex and Transvaal Supergroup is 0.350 g.cm-3. The Bushveld Complex and Transvaal Supergroup show a density contrast of 0.250 g.cm-3 and 0.600 g.cm-3 with IRUPs, respectively. 110 CHAPTER 7: GEOPHYSICAL MODELLING 7.1 Theory In this chapter, 3D methods of modelling potential field data are reviewed. Forward and inverse modelling methods are explained, with reference to 3D modelling using rectangular prisms. This leads to an understanding of the inversion methodology of used in this project, which is fully described. Once all data has been collected, corrected and processed, it can be used to interpret the underlying geology. Potential field modelling is an inherently non-unique process (i.e. several different models can provide the same response). For this reason, interpretation of the source bodies must incorporate all priori geological and geophysical knowledge of the area (e.g. borehole information, outcrop maps). This is required in order to constrain the model to a realistic construction. There are two methods of modelling geophysical data: forward modelling and inversion. Each method attempts to determine the model parameters (e.g. depth, thickness, density, etc.) via a different approach (Figure 7.1). 111 Figure 7.1 ? Comparing the forward and inverse method. A represents the measured anomaly, and A0 the calculated anomaly. Parameters p1, p2,? are the source attributes (e.g. depth, thickness, density), Blakely (1995). The forward modelling approach is heavily based on the user?s knowledge of the geology of the area and the likely model parameters. As seen in Figure 7.1, the method is based on a three-stage process. An initial model is constructed and the calculated anomaly is compared to the measured anomaly. From this comparison, the parameters can be adjusted by the user until the misfit is sufficiently small. Forward modelling may be done in 2-, 2.5-, 2.75- or 3-dimensions. In 2D modelling, bodies are considered to have an infinite strike length, with the profile being perpendicular to the body?s strike. 2.5D modelling algorithms allow a limited strike length of the body, with the profile also perpendicular to the body?s strike length. 2.75D modelling is the same as 2.5D modelling, with the added feature of allowing the profile to be at any angle relative to the body. 3D modelling uses algorithms which take the 3D geometry into account and produce realistic responses. 112 Inversion modelling uses mathematical methods to automatically determine the body parameters from the observed anomaly, with suitable geological constraints (Blakely, 1995). The inversion method can also be adapted for 2-, 2.5-, 2.75- or 3-dimensions. The mathematical approaches used in inversion theory are based on forward modelling. In iterative inverse methods, a computational method automatically applies a similar three-step process until the difference between the measured and calculated anomalies fall below a particular cut-off value. This project concentrates on 3D modelling of the study area and, for this reason, various methods underlying 3D gravity forward modelling are briefly explained. Inversion software such as Grav3D and VPmg (Vertical Prism Magnetic and Gravity inversions) are becoming the industry standard in the oil and mining industries. As such, the modelling methodologies presented here lead to the methodology underlying VPmg. 7.2 3-Dimensional Gravity Forward Modelling Before considering individual forward modelling techniques, it is necessary to explain how gravity anomalies can be calculated (Blakely, 1995). Consider the gravitational potential (U) and gravitational attraction ( g r ) at point P due to a volume of mass V with density ? : ?= V dv r Q GPU )( )( ? (7.1) ??=?= V dv r r QGUPg 3 ? )()( ?v (7.2) where r is the distance from P to an element of the body Q, and G is the gravitational constant (Figure 7.2). 113 ^ Figure 7.2 ? Arbitrary 3D body with density ?(x?,y?,z?) and unit vector r pointing from an element of the body to the calculation point P(x,y,z), (Blakely, 1995). ? If one uses a Cartesian co-ordinate system (as in Figure 7.2) with z-axis positive downwards and x- and y-axes arranged according to a right-hand system, the vertical gravitational attraction (g) at point P(x,y,z) is: ''' )'( )',','(),,( ''' 3 dzdydx r zz zyxG z U zyxg xyz ??? ??=??= ? (7.3) where 222 )'()'()'( zzyyxxr ?+?+?= . Equation 7.3 has the general form: ??? ???= ''' ''')',','()',','(),,( xyz dzdydxzzyyxxzyxzyxg ?? (7.4) where 2/3222 )( ),,( zyx z Gzyx ++?=? . ),,( zyx? is a Green?s function and, in this case, it is the gravitational attraction at (x,y,z) of a point mass located at (x?,y?,z?). 114 The method requires the geology be sub-divided into N simple geometric shapes. The gravity response of each sub-divided shape may then be calculated individually, using Equation 7.3, and all parts summed to give the overall modelled response. Common forms of sub-divided geometrical shapes include rectangular prisms and stacked lamina (Blakely, 1995), or using a series of Fast Fourier Transforms (FFTs) to obtain the gravitational response of a layer bound by two surfaces (Parker, 1972). Once this has been done, Equation 7.3 equates to: ? = = N n mnnmg 1 ?? (7.5) where gm is the vertical attraction at the mth observational point, ?n is the density of part n, and ?mn is the gravitational attraction at point m due to part n with unit density. 7.2.1 Rectangular Prisms Rectangular prisms can be used to sub-divide a geological model (Figure 7.3). Prisms provide a simple way to approximate the volume of the mass, provided they are small enough to assume constant density within each prism. The gravitational anomaly of the body at any observational point may be calculated by summing the effects of all the prisms. Obviously, the larger the prisms, the less accurate the approximation of the body will be. Conversely, the smaller the prisms, the more computational time required. The computational time is proportional to the number of inputs and the number of outputs, N. 115 Figure 7.3 ? Schematic diagram of a 3D approximation of a body, using rectangular prisms (Blakely, 1995). The gravitational attraction of any single prism may be calculated by integrating Equation 7.3 over the limits of a single rectangular prism. Consider a prism of constant density ? with dimensional limits x1 < x < x2, y1 < y < y2 and z1 < z < z2. The vertical gravitational attraction at the origin is: ? ? ? ++= 2 1 2 1 2 1 ''' )'''( ' 2/3222 z z y y x x dzdydx zyx z Gg ? (7.6) The gravitational attraction of each prism can be calculated and summed. This approach shows numerous impracticalities. Firstly, it is difficult to model geological features using rectangular prisms (e.g. dipping faults, tight folds, thin stratigraphic units) and, secondly, the computation does not simplify the model by seeing two or more adjacent bodies of equal density as a single body. Although computer processing speed is increasing, this method is still computationally inefficient. It will be shown how VPmg?s inversion methodology uses an improved form of this approach. 116 7.3 3-Dimensional Gravity Inversion As discussed, the forward method is a three step process, where the parameters (e.g. dimensions, depth of burial, physical properties) of the causative body are estimated through trial-and-error calculations. The inverse method differs in that the source parameters are calculated directly from the gravity measurements. The basic principle of inversion can be explained using the generalised equation for the vertical attraction of gravity (gz, from Equation 7.4), which can be written as: ? ? ??= V i i iz dv rr zz QGPg 3)()( ? (7.7) where V is the volume of the causative source, Pi is the ith observation point located at (x,y,z) and always outside R, Q is the point of integration (xi,yi,zi) within R, and ?(Q) is the density function (Blakely, 1995). This can be presented in the generalised form: ?= V dvQPQsPf ),()()( ? (7.8) where f(P) is the potential field at observation point P, s(Q) describes the physical quantity (in this case, density, but also holds for magnetization) at source point Q, and ?(P,Q) is a function dependant on the geometric placement of point P at Q (Blakely, 1995). In the forward method, f(P) is calculated from user-input functions s(Q) and ?(P,Q) and volume V, where s(Q), ?(P,Q) and V are repeatedly adjusted until f(P) fits the measured data. Hence, f(P) is determined from a solution of s(Q), ?(P,Q) and V (Blakely, 1995). 117 In the inverse method, the measured data are inserted as f(P) and the algorithm solves for s(Q) or V (Blakely, 1995). Solving for s(Q) is known as a linear inverse problem and solving for V is a nonlinear inverse problem. 7.3.1 Linear Inversion The dependence of a gravity field on the sub-surface density is a basic linear inverse problem (e.g. doubling the density of the causative body doubles the gravity anomaly). In the simplest form, assume a known basic shape of homogeneous density. The gravity anomaly may be calculated, at N discrete locations, as: ,...,...,2,1 Nig ii == ?? (7.9) where ? can be determined by linear regression (Blakely, 1995). If the body is heterogeneous in nature, it can be divided into smaller, simple bodies (e.g. rectangular prisms, stacked lamina), with each small body being assigned an individual density. The least-squares method can effectively solve for the density of each body (Blakely, 1995). 7.3.2 Non-Linear Inversion Non-linear inverse problems relate to geometric parameters of the body, contained in the ?(P,Q) function. This includes depth of body, thickness and shape. In order to create inversion algorithms for a geological body, it is necessary to simplify the body. The validity of the inversion depends on whether the assumptions of the simplified body are geologically sensible. Following the development of the methodology, as outlined in Blakely (1995), two examples of simplifying a body for non-linear inversion are now provided (Bott, 1960; Cordell and Henderson, 1968). This leads to an explanation of the methodology used by Grav3D and VPmg. A simple method of estimating the cross-sectional shape of a sedimentary basin was provided by Bott (1960). The basin is assumed to extend infinitely and have a 118 uniform density contrast, ??, to the background density. The basin is divided into N rectangular blocks of equal width, extending infinitely perpendicular to the profile. Each block continues to a depth tj, where j = 1, 2, ? , N. There are N field points, gj, where j = 1, 2, ? , N, along a profile, with each point centred above a block (Figure 7.4). gi tj Figure 7.4 ? Cross-section model of a sedimentary basin. Basin is approximated by rectangular blocks, extending infinitely perpendicular to the profile and with one block per gravity station (Bott, 1960). The initial thickness is calculated for each block, by assuming each block is an infinite horizontal slab: ?? ?= G g t jj 2 1 j = 1, 2, ? , N (7.10) The level of iteration is provided by the superscript. Following the calculation to estimate the initial thickness, a three-step process is followed to iteratively modify the block thickness (with the superscript k indicating the level of iteration): 1. The gravitational field gkj is calculated at each field observation point, taking into account the effect of all the blocks (using the thickness provided from the previous iteration). 2. The residual (gj - gkj, i.e. observed minus calculated) is calculated at each point. 3. The infinite slab approximation (Equation 7.10) is used to estimate a new set of thicknesses. A correction is applied to each block, assuming the block is an 119 infinite slab, with a thickness required to accommodate the original. Hence, the new thickness is: k j k jjk j tG gg t +? ?=+ ??2 1 (7.11) These three steps are repeated until the user is satisfied that convergence is obtained (Bott, 1960). Cordell and Henderson (1968) translated the previous method to 3-dimensions. Gravity data points are measured on, or gridded to, a regular rectangular grid and the source is divided into a grid of rectangular blocks, with one block per data point (Figure 7.5). Observation surface Reference surface Figure 7.5 ? 3D model, defined by Cordell and Henderson (1968), for iterative inversions. Block thicknesses are relative to a reference surface and the calculated gravity data occur on an observation surface. The thickness of the blocks, tj (where j = 1, 2, ? , N), is defined relative to a reference surface (representing the top or bottom of the blocks). As for Bott (1960), 120 the initial block thickness is estimated by calculating the gravity response for an infinite slab (Equation 7.10). The three-step iterative process (calculation, comparison, adjustment) is followed as for Bott (1960), however, the new block thickness, , is adjusted using the following ratio (Cordell and Henderson, 1968): 1+kjt k j j k j k j g g t t = +1 (7.12) 7.3.3 3-Dimensional Inversion Algorithms The development of inversion algorithms from 2D to 3D has allowed for improved 3D inversion algorithms, particularly by incorporating priori knowledge as constraints in the inversion. This produces models that fit the acquired data as well as the constraints from other sources (e.g. geological and geophysical borehole data, surface mapping). Two industry-standard inversion programs were considered for this project: ? VPmg (Fullagar et al., 2000; Fullagar and Pears, 2007; Fullagar et al., 2004): Uses a 3D geological model, built in the 3D modelling software package Gocad. The geological model is made up of geological surfaces (e.g. topography, lithological contacts, faults). These surfaces are exported to VPmg, which divides the Earth into vertical rectangular prisms which extend from the surface of the Earth, hence honouring the topography, infinitely down into half-space. The prisms are subdivided to fit the geological surfaces. In this way, geological units can be defined in terms of lithology, physical property (e.g. density, magnetic susceptibility) and geometry. The inversions alter the geometry and physical properties to fit the data. ? Grav3D (Li and Oldenburg, 1998): Divides the Earth into rectangular cells of equal density, forming an orthogonal mesh. The inversion alters the density of the cells to fit the data. From these inverted densities, the geology (e.g. lithological units, faults) can be inferred to create a 3D geological model. The 121 user can set specific sections of the model to within defined densities, if there is increased confidence of the region. VPmg was chosen for this project as it provided a more applied approach to creating a 3D starting model from various geological and geophysical datasets, constrained the model more precisely and was less computationally intensive. 7.4 VPmg Inversion Methodology The inversion methodology of VPmg is now explained, first providing a brief overview of the functions it can perform. The model parameterisation is then discussed, showing how the vertical prisms are divided such that the lithological model is more closely represented (compared to conventional cell division) and also leading to faster computational times. The inversion algorithm is then summarised showing how the method of steepest descent obtains the best-fit calculated model. Finally, a description of geometrical and physical property constraints is provided. 7.4.1 Overview VPmg uses a pre-defined 3D starting model, containing all priori knowledge in order to better constrain the inversion processes. This 3D model honours: borehole data (imposing constraints on geological units? topography); density ranges of geological units and the basement; the topography; and the original free-air gravity data (Fullagar et al., 2000; Fullagar and Pears, 2007; Fullagar et al., 2004). The model discretization process divides the earth into closely-packed vertical rectangular prisms, which have internal contacts. A starting model is created in a 3D modelling package (i.e. Gocad), which comprises litho-stratigraphic contacts (i.e. surfaces) and structures (e.g. faults), dividing the model into rock-type domains. The VPmg prisms honour the topography implicitly and the internal contacts, 122 approximated from the imported geological surfaces, create the geological unit topography. If the internal contacts are ?pierced? by a drill-hole they can be held fixed, thus constraining the model (Fullagar et al., 2000). During the inversion process, the shape and density of each geological unit can change, by inverting the geometry and physical property, respectively. However, the geological or topological integrity of the model is maintained (Fullagar et al., 2000). Topography plays a significant role in gravity data reductions. The effect of direct modelling of the topography, by the vertical prisms, has a significant advantage in the modelling process. The two corrections that use elevation in gravity reductions are the free-air and Bouguer corrections. The free-air correction is a relatively straightforward application that uses elevation, a fundamentally unambiguous measured quantity (Fullagar et al., 2000). However, the Bouguer correction operates on the assumption of an infinite horizontal plate, with a single user-input density. These assumptions do not hold in areas of rapidly varying topography or quickly changing near-surface densities. By directly modelling the terrain effects, VPmg avoids the Bouguer correction and uses the free-air gravity data in the inversion process. The geometry of the model can be constrained by borehole data. The contacts within the vertical prisms may be fixed, if pierced by a borehole, bounded or free to move during the inversion process. If the contact exists below a borehole, it will be bounded above, and if the contact exists above a borehole, it is bounded below. All contacts are bounded from above by the topography (Fullagar et al., 2000). Physical properties of geological units can also be constrained by user-defined upper and lower bounds. 7.4.2 Model Parameterisation Surfaces are imported into VPmg, dividing the vertical prisms into cells which create geological units. Individual cells is then assigned the lithology and physical property 123 of the geological unit (e.g. density, magnetic susceptibility), (Fullagar and Pears, 2007). By assigning the cells a lithology as well as a physical property, the lithology does not have to be inferred from the rock properties (as for pure property inversions, e.g. Grav3D). It also allows greater control in the inversion process (e.g. allowing inversions to be carried out on a particular lithological unit while keeping other units constant). Once individual cells have been classified according to lithology, their physical properties become characteristic of the entire geological unit. Hence, if a geological unit is homogeneous, all the cells share the same value. If a geological unit is heterogeneous, all the cells share equal upper and lower bounds as well as an appropriate probability density function (Fullagar and Pears, 2007). In geometry inversion algorithms, it is possible to vary geological boundaries by two methods. Firstly, cells can be reclassified from one lithology to another, within a fixed mesh (Figure 7.6a). Here, the basic model mesh is invariant and boundaries have to change in discrete shifts (e.g. Grav3D). This requires that cells sharing a geological boundary must use the physical property of the lithology with the greatest weight in the cell. In this case, entire units can be lost in the model (e.g. the thin red layer). The second option uses an adaptive mesh which consists of closely-packed vertical prisms (e.g. VPmg). These prisms have tops honouring the topography and internal horizontal contacts that divide the prisms into cells. In this method, the horizontal cell boundaries can move vertically while lithology of the unit remains constant (Figure 7.6b). Here, the mesh deforms, with arbitrarily small vertical boundary adjustments. The vertical prisms can be divided again into sub-cells which allow for geological units to be homogeneous or heterogeneous (Fullagar and Pears, 2007). 124 a) b) Figure 7.6 ? Comparing geometry variation methods for two sets of similar geological units (lithological boundaries represented by red and blue lines). a) Conventional fixed mesh (black lines), with invariable horizontal cell divisions and cells assigned according to the lithology occupying the greatest percentage of the cell. b) Adaptive mesh with variable horizontal cell divisions (solid horizontal lines) and sub-cells (dotted horizontal lines), after Fullagar and Pears (2007). Note how the adaptive mesh represents the geology more closely. Various advantages are evident from the section view in Figure 7.6: ? The adaptive mesh allows for greater detail in the geological model, especially for thin units, ? Surfaces are represented more accurately, including the topography, ? Fewer cells are required compared to a conventional mesh, especially for homogeneous units, which reduces inversion run times, and ? For heterogeneous units, the size of the sub-cells can be varied from unit to unit, allowing for increased intra-unit detail (Fullagar and Pears, 2007). 125 The total number of VPmg model parameters depends on the number of prisms and the number of lithological boundaries per vertical prism. This decreases the number of parameters relative to a conventional mesh, usually by an order of magnitude, with a similar effect on the inversion run time (Table 7.1). Table 7.1 ? Comparing the forward modelling run time of VPmg and Grav3D for a model with the same number of data points, using a 1 GHz Pentium III with 512 MB RAM, after Fullagar (2007b). Program Contacts/Cells Data Points Run Time VPmg ~30 000 contacts 216 45 s Grav3D ~106 cells 216 239 s 7.4.3 Inversion Algorithm Once a 3D adaptive mesh has been created, using the 3D geological starting model, four inversion styles are possible: ? Basement physical properties (density, magnetic susceptibility) variable (homogeneously or heterogeneously), with geological units? physical properties and geometries held fixed, ? Geological units? physical properties (density, magnetic susceptibility) variable and homogeneous, geometries and basement physical properties fixed, ? Geological units? geometries variable, all physical properties fixed, and ? Geological units? physical properties variable and heterogeneous, geometries and basement physical properties fixed. These inversion styles are applied sequentially, in any user-defined order and number of times. The sequential style of inversion shows numerous advantages over simultaneous inversion styles, where, for example, geometry and property inversions are performed together. These include enhanced flexibility and control for the user, increased speed in inversions (fewer inversion parameters allow faster inversions) and reduced demands on computer memory (Mira Geoscience, 2006). 126 The inversion process alters the model according to the user-defined style, with associated user-defined parameters (number of inversion iterations, limit on the amount of change per iteration). For geometry inversions, the vertical motion of the cells? horizontal divisions is limited by a percentage of the depth of the division (called ?relative elevation perturbation per iteration?). For property inversions, a ?maximum density change per iteration? is allowed, with upper and lower density bounds for a geological unit. In VPmg, these adjustments are carried out stochastically with random adjustments made, tested and then accepted or rejected (Fullagar and Pears, 2007). In order to calculate the gravitational response of the model, the vertical gravitational responses of all rectangular prisms in the model are summated. For a review of the gravitational response of a vertical prism, see Appendix B. This calculated response is then compared to the observed response and accepted or rejected depending on the improvement to the model. This improvement is defined using the method of steepest descent, now discussed. Since no matrix inversion is required, the steepest descent method is very fast. Fullagar (2007a) shows the steepest descent method is calculated by considering B r , the vector of data residuals: n nn n co B ? ?= (7.13) where and are the observed and calculated data vectors, respectively, and no r nc r n?r is the standard deviation, associated with the nth measurement of no r . The chi-squared, , misfit is defined as: 2? ( )? = = N n nBN 1 22 1? (7.14) where N is the number of data points and (Fullagar and Pears, 2007). 12 ?? 127 The steepest descent perturbation, p? , is anti-parallel to the gradient of the chi- squared misfit. Hence, the steepest descent perturbation vector is defined as: 2?? ?= rr fp (7.15) where f is the scaling factor. The scaling factor, f, is defined by various means, e.g. if the user wishes to halve the chi-squared misfit for each iteration, then: 2 2 2D f ??= (7.16) where f must be negative in order for the chi-squared misfit to decrease (i.e. for the data to improve), and D2 is the squared magnitude of the gradient of the chi-squared misfit, (? ??= k kpD 222 ? ) , (Fullagar, 2007a). Note, kp?? 2? are the elements of the gradient vector: ? ???=?? n knnnk p cB Np ? ? 22 (7.17) where pk denotes the kth inversion parameter. The inversion process ceases when a user-defined fit is achieved or when successive iterations do not produce a significant improvement in the misfit. This is, when the Root Mean Squared (RMS) misfit between the calculated and observed response falls below a user defined value, or when the value of the steepest descent perturbation falls below a certain value (i.e. showing no significant improvement to the model). 7.4.4 Geometrical Constraints Finally, a description of constraints placed on the geometry and physical properties is provided. Besides geometry constraints based on borehole data, a number of general, seemingly obvious, constraints exist: geological boundaries cannot lie above the 128 topographical surface, they cannot pass through each other, and the inversion algorithm dampens large geometric inversion changes at shallow depths (preventing near-surface changes dominating improvement in misfit), (Fullagar and Pears, 2007). Hence, even unconstrained models maintain some form of geological realism. A number of constraints also exist in the VPmg framework. These include so-called hard constraints (borehole pierce points) and soft constraints (weighted changes), (Fullagar and Pears, 2007). In honouring borehole data (i.e. hard constrains), the position of the pierce point in a vertical prism boundary must be considered (Figure 7.7). For a pierce point lying close to the centre of a vertical prism, with gently dipping boundaries, it is clear which cell boundary must be fixed (Figure 7.7a). For a pierce point close to the edge of a prism, it is necessary to introduce an ?activation distance? to fix interfaces which lie within a certain range of the pierce point. This is necessary as the pierce point may be far above or below the nearest cell boundary. In Figure 7.7b, boundary A is within activation distance of pierce point P and is fixed, while boundary B can move. a) b) Figure 7.7 ? a) Gently dipping geological boundaries with pierce points near the prism centre, boundary B is fixed and boundary A can vary, b) Steeply dipping geological boundary with pierce point near the prism edge, contact A is within range of pierce point P (after Fullagar and Pears (2007)). 129 Borehole traces also act as constraints on the model, with cell boundaries limited above and below by the position of boreholes, termed ?bounds? (Fullagar and Pears, 2007). In Figure 7.7a, a continuous geological unit (in green) has been logged between the two pierce points P. During inversion, the interpreted contact A cannot move higher than B. In Figure 7.7b, the interpreted contact at B is bounded above by pierce point P (Fullagar and Pears, 2007). Although pierce points fix the position of cell boundaries, the allowed geometrical change in surrounding cell boundaries must also be constrained to maintain geological sense of the model. These soft constraints are performed by a distance weighting algorithm. The radius of influence of a pierce point is defined as its depth below the topography or the distance to the nearest pierce point, depending on which is smaller (e.g. R1 and R2 in Figure 7.8), (Fullagar and Pears, 2007). Weights are applied to the sensitivity of each data point with respect to changes in elevation of a free interface. This sensitivity is the z-derivative of Equation B.4 (Appendix B), which simplifies to: ??? ? ??? ? +++ ?=? ? ? ))(( )( tan2 21 2 121 mmmm mmz xuxuz uuz G z g ? (7.18) where nnmmn yzyxu ?++= 222 , one for each y-parallel prism face (Fullagar et al., 2004). 130 m m2 n Figure 7.8 ? Radius of influences, with RR1 defined by the depth of the pierce point below the topography and R2R defined as: defined by the distance to nearest pierce point. Geometry changes within the radius of influence are damped during inversion (after Fullagar and Pears (2007)). Using Figure 7.8 as reference, consider a jth free boundary at a distance rjk from the kth pierce point (e.g. boundary m a distance rm2 from pierce point P2). The dampening weight of vertical motion, wj, is ?= k k jk j R r w (7.19) where RRk is the radius of influence of the k pierce point, and provided rth jk < RkR . Using this weight, unconstrained cell boundaries within a radius of influence undergo dampened vertical changes during inversions. For example, in Figure 7.8, boundary m is further from pierce point P2 than boundary n. Hence, boundary n?s vertical changes will be dampened more than boundary m. 7.4.5 Physical Property Constraints As for geometry constraints, hard constraints and soft constraints can be applied to the physical properties of cells. Users can also control which cells? physical properties remain fixed, based on their geological knowledge of the area. 131 Within a homogeneous unit, properties remain constant during inversion, allowing for inter-unit variations. Each homogeneous unit can be classified as active or inactive during an inversion, with upper and lower density bounds applied to the active units, retaining geological sense in the model. The homogeneous inversion problem has a small number of parameters and, even in a large model, the method is fast (Fullagar and Pears, 2007). In the case of heterogeneous units, physical properties vary between cells. This allows any down-hole density logs or core measurements to be honoured for pierced cells (similar to geometry pierce points). The user can decide to hold the value of a pierced cell constant or allow it to vary according to a weight proportional to the standard deviation. The weight is applied as described in the previous section, using a radius of influence (Fullagar and Pears, 2007). 7.5 Conclusion VPmg was chosen for the project, based on availability, ease-of-use, compatibility with Gocad and its ability to invert large amounts of data rapidly in a sequential fashion, with various styles of inversion possible. It is one of the few programs that can invert for property and geometry, while honouring the known constrains from borehole data. Geology is not inferred from a density model, rather it is discretized into geological units from the starting model, and the geological units are not blurred during the property inversions. 132 CHAPTER 8: APPLIED GEOPHYSICAL MODELLING A 3-dimensional gravity model of the contact region of the Bushveld Complex and underlying Transvaal Supergroup was developed from ground and airborne gravity and magnetic data. The chapter shows how the 3D geological starting model was constructed using a host of geophysical data (predominantly, the interpolated surface from the seismic data). Aeromagnetic data was used to qualitatively determine the position of magnetic IRUPs in the study region which were modelled in the geological starting model as vertical pipe-like bodies. The iterative inversion process is explained, using the methodology described in Chapter 7: Geophysical Modelling. Initially, a forward model was created to check the validity of the starting model. This was followed by iterative inversion steps, first modelling long wavelength features in the gravity data and proceeding to model short wavelength features. The Bushveld Complex and Transvaal Supergroup?s bulk densities are modelled, followed by heterogeneous density variations in the Transvaal Supergroup, the geometry of the Transvaal Supergroup contact, heterogeneous density variations in the Bushveld Complex, and homogeneous density variations of IRUPs. The second 3D geological model over the ground gravity survey area is then discussed, attempting to improve the model of overburden thickness. The regional field was removed using the large-scale 3D regional model, after which inversions of the overburden/bedrock contact geometry were carried out. 8.1 2.5-Dimensional Modelling The ground gravity traverses were used to create two 2.5D geological models, which were used in the 3D geological starting model. 2.5D magnetic models were also 133 created to constrain magnetic parameters of the NW trending Pilanesberg dykes. The advantages of 2.5D modelling is that the modelled bodies have a limited strike length, as opposed to 2D modelling which models bodies with infinite strike length. This helps to produce a more geologically accurate model. 8.1.1 Gravity Modelling The Bouguer data of the gravity traverses were modelled in Grav2DC (Cooper, 2003). The Bushveld Complex was modelled perpendicularly to the gravity profiles in 2.5D, with the Transvaal Supergroup as the background model (Table 8.1). The results of modelling Line 1 and Line 2 are presented in Figure 8.1 and Figure 8.2. Table 8.1 ? Modelling parameters for Line 1 and Line 2 gravity profiles. Background density i.e. Transvaal Supergroup 2.600 g.cm-3 Bushveld Complex/Transvaal Supergroup starting density contrast 0.350 g.cm-3 Bushveld Complex strike length 2 000 m 134 13 8 mgal 4 0 0 1000 m 2000 SW NE Figure 8.1 ? 2.5D gravity model of Line 1, with measured data (green line) and calculated data (black line). Inversions produced a density contrast of 0.333 g.cm-3. Line 1 produced a model extending to a depth of ~2 000 m and with a Bushveld Complex density of 2.933 g.cm-3. 135 13 8 mgal 4 0 0 1000 m 2000 SW NE Figure 8.2 ? 2.5D gravity model of Line 2, with measured data (green line) and calculated data (black line). Inversions produced a density contrast of 0.334 g.cm-3. An IRUP was modelled with an inverted density contrast of 0.338 g.cm-3. Modelling of Line 2 produced a similarly shaped model to Line 1. Again, the model extended to a depth of ~2 000 m, with a Bushveld Complex density of 2.934 g.cm-3. The location of the IRUP was postulated from coincident magnetic and gravity anomalies. The dimensions of the IRUP were estimated from the magnetic grid, with a profile length of 900 m and a strike-length of 600 m. The depth to the top of the IRUP, providing a best-fit to the data, was 30 m. Finally, the inverted density of the IRUP was 2.938 g.cm-3. The gravity profiles show a general shape of the Transvaal Supergroup contact (i.e. steeply dipping then flattening out at ~2 000 m depth), consistent with the seismic interpretations. However, the seismic data provided a fully 3D, consistent contact surface which was used for the 3D geological starting model. The above results, i.e. 136 density of the Bushveld Complex and IRUP and geometry of the IRUP, provided additional constraints for the 3D starting model. 8.1.2 Magnetic Modelling Line 1 and Line 2 traversed two dykes which were modelled in 2.5D, using constraints from the Euler solutions, in ModelVision Pro. The results of modelling Dyke 1 and Dyke 2 are presented in Figure 8.3, with a range of solutions from forward modelling shown in Table 8.2. The dykes were modelled with a strike length of 2 000 m and assumed to be sub-vertical. Remanent inclinations and declinations were taken from Gough (1957), while the remanent intensity was modelled in the program. Furthermore, the E-W trending feature is a complicated geological structure consisting of a major dolerite dyke and multiple faults, fractures and IRUPs, hence it is not well suited to modelling in 2.5D. 137 c) b) a) Dyke 1 Dyke 2 Figure 8.3 ? Profiles of a) Line 1 across Dyke 1, b) Line 2 across Dyke 1, and c) Line 1 across Dyke 2. The 2.5D modelling provided a suitable fit to the magnetic data. The dykes are thin and shallow, with intense remanent magnetism (Table 8.2). 138 Table 8.2 ? Range of parameters for 2.5D profiles presented in Figure 8.3. First row shows range, second row (italics) shows value used in presented model. Parameter a) b) c) Magnetic Susceptibility contrast (SI) 0.31 ? 0.40 (0.35) 0.33 ? 0.39 (0.35) 0.27 ? 0.34 (0.32) Depth (m) 5.5 ? 11.0 (8.3) 10.5 ? 13.0 (11.7) 8.0 ? 10.5 (9.4) Thickness (m) 5.0 ? 6.5 (6.0) 4.0 ? 5.5 (4.7) 7.0 ? 8.5 (7.7) Dip (degrees) 87.5 ? 90 (88.5 SW) 86.0 ? 90.0 (88.1 SW) 85.0 ? 89.0 (87.0 NE) Remanent Intensity (nT) 2345.0 ? 2420.0 (2390.8) 310.0 ? 345.0 (330.9) 1505.0 ? 1580.0 (1543.5) Remanent Declination (degrees) 63 ? 72 (69) 64 ? 78 (69) 62 ? 75 (69) Remanent Inclination (degrees) 14 ? 39 (24) 12 ? 42 (24) 10 ? 40 (24) The tops of dykes occur at a depth of 8-12 m, which corresponds to the depth of overburden seen in the region (Figure 6.4) and the majority of solutions provided by Euler deconvolution (Figure 6.18). The dykes are thin, ranging from 4-8 m, and sub- vertical. The remanent intensity shows high values, which is unlikely, and need to be further studied to constrain the remanence (beyond the scope of this project). A palaeomagnetic study of these dykes has not been carried out and would provide this information. 8.2 3-Dimensional Starting Model A simplified 3D geological starting model was constructed in Gocad, which was used as an input for the gravity inversions. The 3D model was simplified to two layers: the Bushveld Complex and the Transvaal Supergroup and includes vertical IRUP bodies cutting through both units. An overburden was not modelled for the airborne survey area as the data were acquired at 80 m height, too far from the source to measure accurate changes in overburden thickness. However, the final inverted 139 model was used as a background model for overburden thickness inversions, using the ground gravity data. The starting model covers the extent of the AGG survey area in plan view (i.e. 10 km x 10 km) and extends to a depth of 2500 m. Data used in constructing the 3D starting geological model included: 1:250 000 geological map, DTM from aeromagnetic survey, seismic profile interpretations, published rock property data, and qualitatively interpreted aeromagnetic data. The geological map, within the AGG survey area, was simplified to two stratigraphic units (i.e. Bushveld Complex and Transvaal Supergroup). The surface contact of the Transvaal Supergroup was simplified, removing the complicated structure in the eastern corner and simplifying the structure in the southern corner (Figure 8.4). The geology was simplified as these areas were too structurally complicated, which the time frame of the project did not allow. Also, the ultimate goal of the project was to investigate long wavelength features, which would not have been affected by the increased detail in these areas. 140 Simplified Transvaal Supergroup Simplified Transvaal Supergroup contact Simplified Bushveld Complex LEGEND Pilanesberg Complex Bushveld Complex Transvaal Supergroup Figure 8.4 ? Simplified geological map (after Walraven (1981)) of the study area, within the AGG survey area (i.e. the white box). The Transvaal Supergroup contact was simplified to a single contact (dashed white line), separating the Bushveld Complex (green) and Transvaal Supergroup (vertical white lines). All other structure was removed. The simplified Transvaal Supergroup sub-outcrop contact and interpreted seismic surface was used to create the contact surface covering the extent of the study area. The starting bulk densities assigned to the Bushveld Complex and Transvaal Supergroup were approximated as 2.950 g.cm-3 and 2.600 g.cm-3, respectively. These values were taken from the averages provided by Mar? et al. (2002) and Ashwal et al. (2005). The starting density of each of the IRUPs was 3.200 g.cm-3. This was considered a reasonable estimate as IRUPs are high density bodies with densities that vary laterally and vertically. The positions of the IRUPs were estimated through a qualitative analysis of the aeromagnetic data (Figure 8.5). 141 11 12 10 7 13 8 6 5 14 9 3 1 4 2 Figure 8.5 ? Plan view of the study area, showing the Transvaal Supergroup contact with IRUPs labelled 1-14. IRUPs were chosen, as a first-pass, based on the following parameters: approximate diameter greater than 150 m, magnetic anomaly showed correlation with Gz gravity high anomaly, and no correlation to cultural or man-made magnetic noise. The outlines of the fourteen interpreted IRUPs were projected down to 2 500 m, piercing through the Transvaal Supergroup contact. The tops of the IRUPs were modelled at 30 m below the surface, based on the 2.5D gravity model presented previously. The geophysical units (i.e. air, Bushveld Complex, Transvaal Supergroup, IRUPs) were separated by the DTM, the Transvaal Supergroup contact and the surface of the IRUPs (Figure 8.6). 142 y = 10 000 m x = 10 000 m DTM IRUPs z = 2 500 m Transvaal Supergroup contact 3D model bounding box N Figure 8.6 ? 3D geological starting model of the study region, showing topographic surface (i.e. surface separating air from the Bushveld Complex and Transvaal Supergroup), IRUPs, and Transvaal Supergroup contact. Dimensions of the bounding box are shown and apply to all subsequent images of the regional model. Vertical exaggeration: 3. The starting model?s surfaces were imported into VPmg which discretized the model into 200 m x 200 m vertical prisms, bounded by the surfaces, i.e. 51 N-S prisms, 51 E-W prisms, hence 2601 prisms in the active area (Figure 8.7). The dimensions of the prisms were chosen in order to provide a short inversion run-time, while ensuring that long wavelength features and IRUPs would still be modelled. 200 m x 200 m allowed a suitable trade-off between time and resolution. 143 N b) a) N Figure 8.7 ? a) 3D voxet of the starting model, with the transparent Transvaal Supergroup contact and cross-sections showing the prisms defining the geophysical units: air (blue), Bushveld Complex (orange), Transvaal Supergroup (lime) and IRUPs (red). b) 3D voxet model showing the discretized 200 x 200 m prisms of the Transvaal Supergroup (lime) and IRUPs (red). Vertical Exaggeration: 3. IRUP 13 and 14?s small diameter caused them to be lost in the discretization process. Hence, only IRUPs 1-12 were modelled. This exclusion process is suitable as the small diameter IRUPs would not significantly affect the gravity data. 8.2.1 Forward Modelling A gravity forward model of the geological starting model was run and compared to the measured data to determine the initial root mean-squared (RMS) misfit. The starting model produced a RMS misfit of 2.431 mgal (Figure 8.8). 144 a) b) Edge effects c) Transvaal Supergroup IRUPs Simplified structure Figure 8.8 ? a) Observed data, b) calculated data from the starting model, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. Qualitatively, various areas of significant difference exist between the measured gravity data and forward model?s gravity response. The positive difference in the Transvaal Supergroup shows the forward gravity response is too low, indicating the estimated bulk density is too low. The southern corner also shows a positive difference, attributed to the simplified structure in that region. Positive and negative difference anomalies are seen along the edges of the study area. This may be due to edge effects related to the calculation of the free-air Gz data from FTG data, or from problems modelling data at the edge of the model. Finally, the negative difference 145 relating to the IRUP-dominated region indicates that the initial bulk density is too high. These regions of positive and negative differences are shown quantitatively in the histogram of the residual data (Figure 8.9). Positive differences (Transvaal Supergroup, simplified structure, edge effects) Negative differences (edge effects, IRUPs) Figure 8.9 ? Histogram of residual data for the starting model, bin size: 0.25 mgal, with normal Gaussian overlay. The histogram shows how the residual data is distributed. Although one population of data is centred about zero, there is also a negative population (relating to differences from edge effects and IRUPs) and a positive population (relating to differences from the Transvaal Supergroup, simplified structure and edge effects). Despite the significant differences in data distribution, the starting model is a reasonable fit, seen in the histogram statistics (Table 8.3). Table 8.3 ? Statistics for the histogram of the residual data, starting model (units: mgal). Samples 10 000 Minimum -7.394 25th Percentile -1.717 Median -0.054 75th Percentile 1.731 Maximum 5.257 Mean 0.000 Standard Deviation 2.431 146 The minimum and maximum values do not necessarily define the histogram distribution. Although, one would expect the minimum and maximum to decrease through the inversions, two important statistical descriptions are the 25th and 75th percentile. These show the residual value below and above which 25% and 75%, respectively, of the observations fall. Assuming an approximately Gaussian distribution centred about zero, the absolute value of these two quantities should be approximately equal. As inversions improve the misfit, these values should also get closer to zero. The median is the value separating the lower 50% from the upper 50% of the observations (also known as the 50th percentile). A residual misfit dataset should show a median close to zero. The mean is the average of all the numbers in the dataset. Again, residual misfit data should have a mean close to zero. Finally, the standard deviation is a measure of the spread about the median. The more closely a model fits the calculated data (i.e. the smaller the residual), the lower the standard deviation. This value can be related to the error in the gravity data. Hence, for the final Gz dataset with an error of 0.870 mgal, a standard deviation close to, or less than, this value indicates a satisfactory model. From the above table, the minimum and maximum residual values are both large. However, the 25th and 75th percentiles (both within ?1.750 mgal and centred about the median, -0.054 mgal), demonstrate that the starting model has a reasonable fit. The standard deviation of 2.431 mgal gives the starting error which must be improved. 8.3 VPmg Inversion: Regional Model An iterative interactive inversion methodology was emplaced using a systematic approach to create the final 3D geological model. VPmg?s four styles of inversion (Chapter 7.4.3: Inversion Algorithm) allows the user to model from long wavelength to short wavelength features, i.e. deep to shallow features (Figure 8.10). 147 STARTING 3D GEOLOGICAL MODEL Inversion 5: Homogeneous density inversion IRUPs Inversion 4: Heterogeneous density inversion Bushveld Complex Inversion 3: Geometry inversion Transvaal Supergroup contact Inversion 2: Heterogeneous density inversion Transvaal Supergroup Inversion 1: Homogeneous density inversion Bushveld Complex & Transvaal Supergroup FINAL 3D GEOLOGICAL MODEL In cr ea si ng w av el en gt h in ve rs io ns D ec re as in g R M S m is fit Figure 8.10 ? Schematic diagram showing the inversion steps carried out in the project to create a final model. As the features undergoing inversion become shallower, the RMS misfit increases. Five inversion steps were carried out using the final free-air Gz data, with the model of the previous inversion used as input to the following inversion. This process allowed the user to reduce the RMS misfit to a suitable value, whilst maintaining geological sense of the model (Fullagar et al., 2000; Fullagar and Pears, 2007; Fullagar et al., 2004). Each inversion was defined using certain input parameters, including the number of inversion iterations, maximum density change per iteration (for density inversions) or relative elevation perturbation per iteration (for geometry inversions), and RMS cut- off (iterations stop after the RMS reaches this value). These criteria were commonly chosen to account for bulk changes in the data. Hence, the number of iterations were 148 relatively low (5-10) and large changes per iteration were not allowed. The RMS cut- off was always input at 0.870 mgal, i.e. the accuracy of the free-air Gz data. The inversion was allowed run until all iterations were completed or until the RMS improvement ceased (Chapter 7.4.3: Inversion Algorithm). 8.3.1 Inversion 1: Homogeneous Density Inversion The initial inversion optimised the bulk densities of the two geological units in the model (i.e. Bushveld Complex and Transvaal Supergroup). Since the inverted property was the bulk density, in a large volume of rock, only ten iterations were allowed with a maximum density change of 0.005 g.cm-3 per iteration. The density values were also bound above and below the starting density value, based on estimated variations in Ashwal et al. (2005)?s density log of the Main Zone. The inversion stopped after six iterations, due to the RMS misfit not improving in the seventh iteration (this is termed a stalled inversion). The average density contrast between the Bushveld Complex and the Transvaal Supergroup was reduced from 0.350 g.cm-3 to 0.303 g.cm-3 (Table 8.4). Table 8.4 ? Starting densities of the regional model (after Mar? et al. (2002) and Ashwal et al. (2005)), with inverted densities. Unit Starting Density (g.cm-3) Minimum Density (g.cm-3) Maximum Density (g.cm-3) Inverted Density (g.cm-3) Bushveld Complex (Main Zone) 2.950 2.750 3.150 2.926 Transvaal Supergroup 2.600 2.400 2.800 2.623 The inversion shows no significant bulk density change in the two stratigraphic units. Although the inversion was restricted to a small number of iterations, with small property changes per iteration, the iterations stopped short or the allowed number. This indicates that the initial density values were a reasonable approximation. The small percentage change in the Bushveld Complex (0.814%) and Transvaal 149 Supergroup (0.885%) further demonstrate that the starting density values were close to the real bulk value. The inversion improved the RMS misfit from 2.431 mgal to 2.091 mgal (Figure 8.11). Hence, a misfit improvement of 0.340 mgal was achieved. a) b) Edge effects c) Transvaal Supergroup IRUPs Simplified structure Figure 8.11 ? a) Observed data, b) calculated data after inversion 1, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. As seen in the qualitative assessment of the forward model, there are several areas where significant differences exist between the measured gravity data and inverted 150 model?s gravity response. The differences in the Transvaal Supergroup, area of simplified structure, regions relating to edge effects and IRUPs are still present. Quantitatively, these differences are shown in the histogram of the residual data (Figure 8.12). Figure 8.12 ? Histogram of residual data after inversion 1, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.9. The histogram does not show significant improvements from the histogram of the forward model, still not matching a Gaussian distribution. The negative population, relating to differences from edge effects and IRUPs, and the positive population, relating to differences from the Transvaal Supergroup, simplified structure and edge effects, still exist. Inversion 1 shows a standard deviation of 2.091 mgal, seen in Table 8.5?s statistics of the histogram. Table 8.5 ? Statistics for the histogram of the residual data after inversion 1 (units: mgal). Samples 10 000 Minimum -7.344 25th Percentile -1.529 Median 0.185 75th Percentile 1.387 Maximum 5.472 Mean 0.000 Standard Deviation 2.091 Positive differences (Transvaal Supergroup, simplified structure, edge effects) Negative differences (edge effects, IRUPs) 151 The minimum and maximum values of the residual data did not change significantly from the forward model, even though the model improved. The improvement is reflected in the 25th and 75th percentiles, which improved by 0.188 and -0.344 mgal, respectively. The median increased to 0.184 mgal, indicating that the bulk of the residual data points are positive. However, the average of all the values is still 0.000 mgal, with a standard deviation of 2.091 mgal. The decrease in standard deviation also indicates the improvement in the model. The residual map shows a number of areas of poor misfit, which is quantitatively confirmed in the histogram and corresponding statistics. This indicates that bulk homogeneous changes are not sufficient to model the geology. Therefore, lateral density changes must be accounting for differences in the measured and calculated data. 8.3.2 Inversion 2: Heterogeneous Transvaal Supergroup Density Inversion In order to remove the long wavelength misfit in the data, the density of the Transvaal Supergroup was inverted, allowing the unit to be heterogeneous (with the Bushveld Complex density and Transvaal Supergroup contact kept constant). This follows the approach of inverting for long wavelength to short wavelength features. Inversion 1?s model was used as an input, with inversion parameters of five allowed iterations, a maximum density change of 0.005 g.cm-3 per iteration and inversion 1?s upper and lower density bounds. Again, the small number of iterations and limit on density change per iteration were used to account for bulk changes. The inversion ran for all five iterations, producing an inverted density range from 2.595 g.cm-3 to 2.630 g.cm-3 (Figure 8.13). 152 Modelled Transvaal Supergroup contact N g.cm-3 Figure 8.13 ? Plan view of the density of the inverted heterogeneous Transvaal Supergroup (units: g.cm-3). Density properties are constant in each vertical prism, extending from the Transvaal Supergroup contact to half-space. The inverted density range show a small percentage change (-1.067% to 0.267%) from the input density of 2.623 g.cm-3. Again, this small change suggests that the input model shows an accurate approximation of the Transvaal Supergroup?s bulk density. This variation is seen in the density map, relating to the conjugate set of the Rustenburg Fault. The offset of this fault within the Transvaal Supergroup is expressed as a change from high to low density. Inversion 2 improved the RMS misfit from 2.091 mgal to 1.858 mgal (Figure 8.14). Hence, a misfit improvement of 0.233 mgal was achieved. 153 a) b) c) Transvaal Supergroup (improvement) Simplified structure Improvements Figure 8.14 ? a) Observed data, b) calculated data after inversion 2, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. The density inversion of the Transvaal Supergroup led to a general improvement in the residual map. This is seen particularly in the anomaly over the Transvaal Supergroup, as well as smaller improvements labelled above. The histogram (Figure 8.15) shows that these qualitative improvements are reflected quantitatively. 154 Small population ? positive differences Small population ? negative differences Figure 8.15 ? Histogram of residual data after inversion 2, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.12. For the first time, the histogram shows an approximately Gaussian distribution. Both the positive and negative populations have been significantly reduced. These improvements indicate that the model is approaching a fit of the measured data, and that the inversion methodology, inverting from long wavelength to short wavelength features, is correct. The histogram shows a standard deviation of 1.858 mgal, seen in the statistics (Table 8.6). Table 8.6 ? Statistics for the histogram of the residual data after inversion 2 (units: mgal). Samples 10 000 Minimum -6.980 25th Percentile -1.174 Median 0.006 75th Percentile 1.140 Maximum 5.276 Mean 0.000 Standard Deviation 1.858 Both the minimum and maximum residual values improved, as did the 25th percentile, by 0.355 mgal, and the 75th percentile, by -0.247 mgal. The histogram clearly shows an approximate Gaussian distribution, reflected in the approximately equal absolute values of the 25th and 75th percentiles, which are centred about 0.006 mgal (i.e. even 155 amount of negative and positive residual values). The decrease in the standard deviation also indicates the improvement in the model. The lateral density changes seen in the second inversion produced a satisfactory development in the model, improving part of the long wavelength misfit. 8.3.3 Inversion 3: Transvaal Supergroup Contact Inversion The next step in modelling long wavelength features would result from an improvement to the geometry of the Transvaal Supergroup contact. This geometry inversion was carried out, keeping the physical properties of the Bushveld Complex, Transvaal Supergroup and IRUPs constant. The model from inversion 2 was used as input, with five iterations and a relative elevation perturbation of 0.005% (of depth) per iteration allowed. These parameters were employed as the IRUPs caused problematic geological models if large numbers of iterations and changes per iteration were allowed. That is, the contact surface showed large changes close to the IRUPs, where the inversion attempted to adjust the geometry to fit the data, but produced unrealistic results. Away from the IRUPs, the surface remained largely unaltered. The Transvaal Supergroup contact was poorly constrained, with only 27 boreholes intersecting the contact, all of which were in close proximity to the outcrop of the contact (i.e. poor spatial distribution). The bottoms of non-intersecting boreholes were used as upper bounds for the surface. The inversion ran for all five iterations, producing an inverted surface which varied from -10.00 m to 10.05 m (i.e. above and below the original surface, respectively, Figure 8.16). 156 N Figure 8.16 ? 3D model showing the minor differences between the original Transvaal Supergroup contact (yellow) and inverted Transvaal Supergroup contact (green). Vertical exaggeration: 3. These changes were small (<2%), relative to the depth of the Transvaal Supergroup contact. Unfortunately, the inversion did not improve the RMS misfit; hence the original surface was maintained for the rest of the modelling process. However, these minor changes indicate a high confidence in the interpreted seismic boundary. As a test, this inversion step proves the robustness of the seismic interpretation. 8.3.4 Inversion 4: Heterogeneous Bushveld Density Inversion Continuing with the process of inverting for long wavelength to short wavelength features, the fourth inversion allowed for a heterogeneous Bushveld Complex. The densities of prisms, comprising the Bushveld Complex unit, were allowed to vary to create lateral density changes. However, the density of each prism was kept vertically constant; hence there were no vertical density changes. The Transvaal Supergroup contact geometry, Transvaal Supergroup density and IRUP density did not change in the inversion. The model from inversion 2 was used as the input. The inversion parameters were input as ten iterations, a maximum density change of 0.005 g.cm-3 per iteration and 157 using the density bounds of inversion 1 (i.e. 2.750 g.cm-3 to 3.150 g.cm-3). Similar to inversion 1 and 2, these constraints were used to account for the bulk density changes in the geological unit. The inversion stalled after seven iterations, showing no improvement in the eighth iteration. The inverted densities varied from 2.870 g.cm-3 to 2.980 g.cm-3 (Figure 8.17). E-W trending feature Simplified structure Transvaal Supergroup (?<2.870) N g.cm-3 Figure 8.17 ? Plan view of the inverted densities of the heterogeneous Bushveld Complex (units: g.cm-3). The heterogeneous Transvaal Supergroup appears blue as it falls below the density range (?<2.870 g.cm-3), and the IRUPs appear red as they fall above the density range (? = 3.200 g.cm-3). The inversion leads to density range with a small percentage change (-1.914% to 1.846%) from the input density of 2.926 g.cm-3. Again, this small variation reflects an accurate approximation of the Bushveld Complex?s input bulk density. The density map shows a distinct change from positive to negative values, relating to the E-W trending magnetic feature. This is a complicated structure, consisting of parallel 158 dykes and faulting. The offset of the Bushveld Complex is expressed as a density variation. Inversion 4 improved the RMS misfit from 1.858 mgal to 1.364 mgal (Figure 8.18). Hence, a misfit improvement of 0.494 mgal was achieved. b) a) c) Edge effects (minor improvement) Improvements Simplified structure (improvement) Figure 8.18 ? a) Observed data, b) calculated data after inversion 4, and c) residual data, taking the difference between the observed and calculated data. Note, the change in mgal range. The heterogeneous density inversion of the Bushveld Complex led to a significant improvement in the RMS misfit. This is seen as a general improvement in the 159 residual map (some specific areas are highlighted above). Although the anomaly over the region of simplified geology is still a prominent feature, the gravity response has been improved. The edge effects also show minor improvements, but the residual is not fully resolved. Quantitatively, these improvements are seen in the residual map?s histogram (Figure 8.19). Improvement in negative differences Single population centred about 0 mgal Improvement in positive differences Figure 8.19 ? Histogram of residual data after inversion 4, bin size: 0.25 mgal, with normal Gaussian overlay. Note, change in x-axis range from Figure 8.15. The histogram has improved beyond a Gaussian distribution, with a single defined population centred about zero. All positive and negative populations have been removed. This inversion significantly reduced the misfit of the calculated data relative to the measured data, maintaining a geologically realistic model. Again, this indicates that the inversion methodology is following the appropriate steps. The standard deviation of the histogram is 1.364 mgal (Table 8.7). Table 8.7 ? Statistics for the histogram of the residual data after inversion 4 (units: mgal). Samples 10 000 Minimum -6.543 25th Percentile -0.622 Median -0.054 75th Percentile 0.629 Maximum 3.777 160 Mean 0.000 Standard Deviation 1.364 The statistics show further improvements in the minimum and maximum residual values. The 25th and 75th percentiles indicate the tightening of residual values about the median: large improvements in the 25th percentile (by 0.552 mgal) and 75th percentile (by -0.511 mgal) make the absolute values approximately equal. The median shows a slight offset to zero, but is close enough to indicate equal amounts of negative and positive values. The standard deviation decreased to 1.364 mgal, another indication of the enhancement to the misfit. As mentioned, the E-W trending structure is a complicated feature which was not modelled in the geological model due to lack of information and constraints. Modelling this feature would clearly improve the model and may be considered for future work. The model could be improved further by vertically sub-dividing each prism (e.g. forming 200 x 200 x 200 m cells) and running another heterogeneous density inversion on the model. This process was not carried out due to the lack of physical property constraints in the model. 8.3.5 Inversions 5-7: Homogeneous IRUP Density Inversion The only features which remain to be inverted are the IRUP bodies, showing the shortest wavelength anomalies. This final step required three separate inversions of the homogeneous IRUPs and the heterogeneous Transvaal Supergroup, with the Transvaal Supergroup contact geometry and Bushveld Complex density remaining unchanged in each inversion. Inversion 5: Homogeneous density inversion of IRUPs, allowing three iterations with a maximum density change of 0.200 g.cm-3 per iteration, and the IRUPs? density bounded at 2.700 g.cm-3 to 4.000 g.cm-3. A large value for maximum density change per iteration was chosen, since the IRUPs occupy a small volume of the model. Data 161 far from the IRUPs would therefore be insensitive to the IRUP density changes, resulting in a rapidly stalled inversions (Pears, 2008). Inversion 6: Heterogeneous density inversion of the Transvaal Supergroup, with two iterations, a maximum density change of 0.005 g.cm-3 and the density range bound between 2.400 g.cm-3 to 2.800 g.cm-3. This was run on the Transvaal Supergroup density, with negligible effect. The reason for running this inversion is related to the reason provided above, i.e. to change the model by a tiny amount, allowing another inversion on the IRUP densities. Inversion 7: A final homogeneous density inversion was run on the IRUPs. Again, three iterations were allowed with a maximum density change of 0.200 g.cm-3 per iteration and the IRUPs? density bounded from 2.700 g.cm-3 to 4.000 g.cm-3. This inversion stalled after two iterations, indicating negligible improvements in further inversions. The inverted bulk densities varied between 2.770 g.cm-3 to 3.185 g.cm-3 (Table 8.8, Figure 8.20) Table 8.8 ? IRUPs inverted densities, including change in density (using the starting density of 3.200 g.cm-3), average inverted density and average change in density. IRUP number Inverted Density (g.cm-3) Change in Density (g.cm-3) 1 3.180 -0.020 2 3.060 -0.141 3 3.002 -0.198 4 2.982 -0.218 5 2.774 -0.427 6 3.144 -0.057 7 3.179 -0.022 8 3.051 -0.149 9 3.172 -0.028 10 3.185 -0.015 11 3.180 -0.020 12 3.184 -0.016 Ave. 3.091 -0.109 162 Transvaal Supergroup (?<2.870) 12 9 10 8 4 11 1 2 5 3 6 7 5 000 m Bushveld Complex (2.870 2.05 Ga Intrusion of RLS controlled by far-field stresses. Possible feeder fissures occurring in western BC (Kruger, 2005) in zones of weakness created by pre-BC far-field stresses (Du Plessis and Walraven, 1990). > 2.20 Ga Right-lateral wrenching along TML, causing development (or possible reactivation) of extensional, NW-orientated features in the western BC (focus of stress-field), (e.g. Rustenburg Fault) and short-wavelength folds, axes trending NE, in the TVL, (Du Plessis and Walraven, 1990). ~2.30 Ga NE-SW compressive force Resulting in TVL folding, axes trending NW (Bumby et al., 1998) 2.72-2.57 Ga NNW-SSE compressive force Limpopo orogeny incl. emplacement of Great Dyke (Zimbabwe) and Ventersdorp Supergroup (SA), i.e. collisional rifts (Silver et al., 2004). 2.95-2.85 Ga NW-SE compressive force Collision of terranes along northern and western margins of Kaapvaal Craton, resulting in NE trending mantle anisotropy (Silver et al., 2004). ~2.96 Ga (Good and de Wit, 1997) TML present on Kaapvaal Craton, early suture feature (de Wit et al., 1992; Eglington and Armstrong, 2004). 184 The above table provides a summary of the published tectonic forces which had an influence on the study area and surrounding regions. There are clearly conflicting far-field stress directions, especially around the intrusion of the Bushveld Complex (~2.05 Ga). Du Plessis and Walraven (1990) show NE-SW compressive stress, while Bumby et al. (1998) and Silver et al. (2004) agree on a general NW-SE compressive stress, relating to orogenic processes on the northern and western edges of the Kaapvaal Craton. The NW-trending open folds in the region necessitate post-Bushveld intrusion NE- SW compressive stress, in agreement with Du Plessis and Walraven (1990). A significant aspect of the TML?s history is the changes in orientation of the primary far-field stress. The change from NW-SE in pre-Bushveld Complex to NE-SW in post-Bushveld Complex allows for compressive features (e.g. open folds, axes trending NW) and extensional features (e.g. Rustenburg Fault) to run parallel to each other in the region. Du Plessis and Walraven (1990) note that the strike of the TML remains constant in both strain ellipsoid orientations, and that reactivation of pre- Bushveld NW-trending extensional features resulted from the intrusion of the Bushveld Complex and later Pilanesberg Complex. 9.1.3 Intrusion of the Pilanesberg Complex The nearby Pilanesberg Complex, found at the northern edge of the study area, is the remnants of a ~1.3 Ga volcano. Major pre-Bushveld NE-SW extensional faults in the area (e.g. Rustenburg Fault) controlled the emplacement of this alkaline intrusion, acting as pathways for the intrusive magma (Du Plessis and Walraven, 1990). Subsequent uplift and erosion has exposed the roots of this ancient volcano. Vermaak (1976) suggests that the Pilanesberg intrusion caused reactivation of pre- existing faults and structures in the surrounding region, while two further sets of fractures are interpreted to be associated with the intrusion. The first set occurs within the complex itself and consists of concentric faults which have been filled with 185 so-called ?ring dykes? (Cawthorn, 1988). The second is a regional set of NW to NNW faults and joints running through the complex from Botswana in the north to the Witwatersrand gold-fields in the south, running parallel to the Rustenburg Fault. Many of these faults and joints were filled with magma during a period of stress- relaxation (Walraven and Darracott, 1976), forming characteristic negative polarity dykes. These dykes are assumed to be Pilanesberg age (Gough, 1957; Rompel et al., 2002) due to their spatial association; however, the Pilanesberg shows a positive magnetic polarity. This indicates that a magnetic reversal must have occurred between the emplacement of the Pilanesberg and the dykes. 9.2 Interpretation of Regional Modelled Features The geological features interpreted from the inverted 3D gravity model and corresponding suite of geophysical data include: the Bushveld Complex/Transvaal Supergroup contact (and related folding and faulting), faulting relating to the prominent E-W trending magnetic anomaly (and related dykes), the Pilanesberg dykes, IRUPs and the thickness of the overburden (associated with recent weathering of faults and dykes). 9.2.1 Contact of Bushveld Complex/Transvaal Supergroup The depth of the Bushveld Complex?s economic layers in the Critical Zone is a significant factor for mining, including feasibility studies and mine planning. The conjugate set of the Rustenburg Fault poses a geological problem: does the Critical Zone onlap the Transvaal Supergroup or does the Bushveld Complex occur ?conformably? above the Transvaal Supergroup, sub-outcropping close to surface? The conjugate set of the Rustenburg Fault was simplified in the 3D model, as the underlying geology of the study region is very complicated. The mapped conjugate fault actually terminates within the study region, indicating the possibility of the both 186 the above scenarios occurring. 2D gravity-gradient forward modelling of each scenario shows the likely Tzz response, compared to the mapped measured Tzz response. For a layered model (Figure 9.3), the Critical and Lower Zones were combined to form a ~200 m thick layer with a mean density of 3.100 g.cm-3 (range: 2.800 ? 3.300 g.cm-3). The Critical and Lower Zones occur between a 0 ? 3500 m thick Main Zone and 0 ? 2000 m thick Transvaal Supergroup, with all layers dipping at 15? (Antoine, 2004). Referring to the simplified geological map (Figure 6.1), the southern corner of the study region shows sub-outcrop of the Critical Zone (hosting the Merensky Reef and UG2). The associated simplified geological cross-section shows a similar geological situation, with a modelled Tzz response similar to that mapped across the sub-outcropping Critical Zone (approximate position indicated by the sub-outcrop of the UG2, Figure 9.3). Eo UG2 NE fault Gravity highs associated with the BC/TVL contact SW Figure 9.3 ? (left) Synthetic 2D Tzz response across a Bushveld Complex layered with the Transvaal Supergroup (after Antoine (2004)) using Grav2DC (Cooper, 2003) and SignProc (Cooper, 2000), (right) Gridded Tzz map over the study area, showing the response across a similar geological scenario (black line), from AGG data. The 2D model shows a clear positive Tzz anomaly over the sub-outcrop of the layered Critical/Lower Zone, with an increase of ~40 Eo from background values (range: 20- 60 Eo, with varying density of the Critical/Lower Zone). The southern corner of the 187 study region shows a similar positive response (shape and amplitude) across the sub- outcropping Critical Zone. The alternative model shows the fault-bounded Critical/Lower Zone onlapping the Transvaal Supergroup (Figure 9.4). The model has similar parameters to that used above, except with the Transvaal Supergroup dipping at 30?, simulating the dip of the fault (range: 25-45?), (Antoine, 2004). Again referring to the simplified geological map (Figure 6.1), the conjugate set of the Rustenburg Fault can be modelled, with no apparent sub-outcrop of the Critical/Lower Zone. Figure 9.4 ? (left) Synthetic 2D Tzz response from a Bushveld Complex onlapping the Transvaal Supergroup (after Antoine (2004)) using Grav2DC (Cooper, 2003) and SignProc (Cooper, 2000), (right) Gridded Tzz map over the study area, showing the response across a similar geological scenario (black lines), from AGG data. The 2D model shows an obvious negative Tzz response to the onlapping contact between the Bushveld Complex and Transvaal Supergroup, showing a decrease of ~100 Eo (range: 85-115 Eo, with varying density of the Critical/Lower Zone and dip of the Transvaal Supergroup) from background values. The negative response is apparent in the Tzz grid of the AGG data, with a smaller amplitude change (~65 Eo), and ignoring the high Tzz values of the IRUPS. The gravity gradient response of the Gravity lows associated with the BC/TVL contact SW fault NE Eo Response of onlapping Critical/Lower Zones masked by regional response IRUP UG2 188 Critical/Lower Zone onlapping at depth is masked by the regional response of the Main Zone/Transvaal Supergroup contact. The 2D models provide dramatically different results, indicating that both the layered contact and fault-bounded onlapping contact are likely to occur in the structurally complicated south-east region of the study area. The onlapping contact changes to a layered contact as the conjugate fault terminates and sub-outcrop of the Critical Zone occurs. A detailed model of the layered region in the southern corner would improve the present model and is one aspect to be explored in future work. 9.2.2 Folds The contact of the Bushveld Complex with the Transvaal Supergroup, derived from seismic data, shows two sets of open fold structures. Firstly, the contact shows a clear change in gradient, with a steep gradient (20-25?) towards the surface changing to a shallower gradient (5-10?) at the base of the fold, as traversed SW-NE. This change in gradient can be related to the open fold structures, axes trending NW (wavelength: ~20 km), which dominate the area. Shallowing of this contact has considerable economic implications. Secondly, folds with axes trending NE (wavelength: ~5 km) are clearly imaged in the Transvaal Supergroup contact from the seismic data. The folds, axes trending NW, as interpreted by Walraven and Darracott (1976), have a half-wavelength (i.e. from peak of anticline to trough of syncline) of 7.5 ? 12 km. In the study region, the axis of an interpreted syncline (S3) has been traced close to the north-eastern edge of the AGG data (Figure 9.5). It is reasonable to assume that a corresponding anticline may be present to the south-west of syncline 3, i.e. anticline 3 (A3), forming another syncline-anticline pair. In this case, the top of the anticline would be present close to the south-western edge of the AGG data. 189 LEGEND Figure 9.5 ? Geological map (after Walraven (1981)) of the study area, showing the axis of syncline 3 (S3), NE of the study area, and the inferred axis of anticline 3 (A3), SW of the study area. The conjugate set of the Rustenburg Fault exists close to the edge of the study area, causing the Bushveld Complex to onlap the Transvaal Supergroup (Figure 9.5). Hence, only one limb of A3 (i.e. limb dipping NE) may have occurred in the Bushveld Complex, with the other limb (i.e. dipping SW) occurring in the Transvaal Supergroup. Anticline 3 could also, possibly, have been bounded by the Rustenburg Fault or may not have extended into the Transvaal Supergroup, just as the axes of all folds die out before the Transvaal Supergroup (Walraven and Darracott, 1976). If A3 does exist, the study region falls above the ?steepest? section of anticline 3, which can be explained by a series of folding, tilting and erosion (Figure 9.6). Cawthorn (1998) and Bumby (1998) both suggest ~10? of tilting of the Bushveld Complex, also imposed on the underlying Pretoria Group, towards its centre due to thermal subsidence. If the folded region was subjecting to such tilting and erosion, this would Pilanesberg Complex Bushveld Complex Transvaal Supergroup Aeromagnetic survey outline S3 Airborne gravity survey outline Rustenburg Fault Rustenburg Fault: conjugate set S3 A3 Suboutcrop: Merensky Reef Suboutcrop: UG2 A3 190 lead to steeper dips (20-25?) at surface, which decrease (to 5-10?) towards syncline 2 in the north-east region of the study area, as observed. In this way, the Transvaal Supergroup contact shallows towards the north-east. Study Area a) SW NE 10? b) Far-field stress 20-25? c) 5-10? Figure 9.6 ? Schematic cross-section (trending from SW to NE) representing possible deformational history of the Bushveld Complex in the study region. a) Intrusion of the Bushveld Complex (layers originally horizontal, Letts (2007)) into the Transvaal Supergroup (dipping ~10?, Cawthorn, 1998) adjacent to the Rustenburg Fault, b) NE-SW far-field stress causes gentle open folds, forming a syncline-anticline pair (i.e. S3-A3), followed by ~10? tilting and erosion (Cawthorn, 1998 and Bumby, 1998), and c) Present tilted fold, overlaying interpreted seismic profile, with Transvaal Supergroup contact drawn in yellow. 191 This sequence of events leads to flattening of the Transvaal Supergroup contact, as one traverses SW to NE. Unfortunately, the seismic profile and gravity data do not extend far enough NE to determine the fold?s continuity. Hence, the model is seen as an extension of the Walraven and Darracott (1976) fold model. The model also incorporates a folded Bushveld Complex which onlaps the Transvaal Complex in the region of the fault, as suggested by 2D gravity gradient modelling (Antoine, 2004). The interpretation of the seismic data also suggests the presence of a second set of folds, with axes trending NE, in the Transvaal Supergroup contact (Figure 9.7). Again, the seismic interpretation was not altered by the geometric gravity inversions, proving the robustness of the seismic data. The half-wavelength of these folds are ~2 ? 2.5 km, which is much tighter than the long wavelength folds, axes trending NW. These folds are possibly related to a pre-2.2 Ga compressive force and, as it is unclear if the folds are restricted to the Transvaal Supergroup. 192 LEGEND Figure 9.7 ? Geological map (after Walraven (1981)) of the study area, showing short wavelength, NE-trending folds (wavelength: ~5 km), in the Transvaal Supergroup, as seen in the Bushveld Complex/Transvaal Supergroup contact. Again, if this change in the Transvaal Supergroup contact topography is inferred in the Critical Zone of the Bushveld Complex, there are economic implications for mining. 9.2.3 Faults A number of faults have been interpreted in the study region, either from geological maps, seismic data or interpreted magnetic data. Two known faults were delineated in the gravity modelling process: the conjugate fault forming the contact between the Bushveld Complex and Transvaal Supergroup, and the fault, or series of faults, related to the major E-W trending dyke (Figure 9.8). Aeromagnetic survey outline Pilanesberg Complex Airborne gravity survey outline Bushveld Complex Transvaal Supergroup Rustenburg Fault Rustenburg Fault: conjugate set Suboutcrop: Merensky Reef Suboutcrop: UG2 193 a) b) N N g.cm-3 g.cm-3 Figure 9.8 ? a) Plan view of the inverted Transvaal Supergroup heterogeneous density, delineating the Rustenburg Fault conjugate set (white dashed line), and b) Plan view of the inverted Bushveld Complex heterogeneous density, delineating the E-W-trending fault (white dashed line). The Transvaal Supergroup (? ? 2.620 g.cm-3) is imaged in blue along the SW edge as it falls below the colour range for the density scale, while IRUPs (? ? 3.091 g.cm-3) occur in red as they fall above the colour range. A fault model of the conjugate set of the Rustenburg Fault was imaged, due to the high density contrast between the Bushveld Complex and Transvaal Supergroup rocks. The model showed a fault through the heterogeneous Transvaal Supergroup, delineated as a small density change (?? ? 0.035 g.cm-3), along the distinctive gravity gradient between the Bushveld Complex and Transvaal Supergroup (?g ? 5 mgal). Gravity modelling is not necessarily well-suited to faults within the Bushveld Complex, due to the low density contrast between Bushveld units and thickness of the units. However, the fault relating to the E-W trending dyke was delineated as a small change in density in the heterogeneous Bushveld Complex (?? ? 0.095 g.cm-3). This result was unexpected as there is no clear gravity gradient in the data to suggest that the fault would be imaged. 194 The method of delineating faults, using inverted densities, provides a first pass at indicating major structure in the study region. The presence of faults in the Bushveld Complex has economic implications which must be considered in mine planning. 9.2.4 Magnetic Interpretation Three prominent magnetic features are seen in the aeromagnetic data: the high intensity (>28 800 nT) E-W trending dyke and the negative (<28 250 nT) NW- to NNW-trending Pilanesberg dykes. Both these features traverse large portions of South Africa, as seen in the regional aeromagnetic map. Finally, numerous highly magnetic IRUP bodies can be interpreted in the data. The E-W trending dyke is related to a late-stage, possibly Karoo, event (Rompel et al., 2002), with dykes cross-cutting most geological and geophysical features in the study area, including the Pilanesberg dykes and IRUPs. The dyke is structurally complex (i.e. faulting and jointed), which is characteristic of this swarm of E-W trending dykes (Rompel et al., 2002). These features affect mining in the Bushveld Complex and are prominent as far south as the gold mines of the Witwatersrand Supergroup. The Pilanesberg dykes show a characteristic negative signature, indicating intense remanence. They cross-cut IRUPs but are cross-cut by the E-W trending dykes, indicating they are younger than the IRUPs but older than the E-W dykes. The dykes appear to be spatially associated with the Pilanesberg and trend NW to NNW. They can be traced from Botswana through to the gold fields of the Witwatersrand, where they also complicate mining. Although the dykes are assumed to be of Pilanesberg age (Gough, 1957; Rompel et al., 2002), it should be noted that the Pilanesberg Complex has a positive magnetic signature, indicating a magnetic reversal occurred from the time of emplacement of the complex to the intrusion of the dykes. Magnetic 2.5D modelling suggests the dykes are sub-vertical, consistent with Campbell?s 195 (2006) observation that proximal Pilanesberg dykes to the west of the complex dip sub-vertically eastwards. Geological magnetic anomalies were interpreted using a combination of aeromagnetic data, Euler solutions, geological maps and QuickBird orthophotographs. Aeromagnetic data, sun-shaded from 90? inclination to enhance the edge of magnetic bodies, was overlaid with Euler solutions to delineate dykes, with associated depth information (Figure 9.9). Figure 9.9 ? Vertically sun-shaded aeromagnetic data with lines of Euler solutions delineating dykes. Colours represent the depth of dykes (green = 5-15 m, blue = 15-25 m, red = 25-35 m). Pilanesberg dykes trend NW-SE, with a predominant depth of 5-15m, while the E-W trending dyke shows a depth of 25-35m. Note: the depth of dykes varies within the specified range. A positive magnetic anomaly sub-parallel to the UG2 sub-outcrop is represented by the yellow line. 196 Euler solutions for the E-W trending dyke were unexpectedly poor, possibly related to its structurally complex nature. Hence, the feature was mapped in Figure 9.9 as a single dyke with the dyke tops at a depth of 25-35 m. This is very deep (boreholes and cuttings indicates depth to top of dykes <20 m, Rompel et al. (2002)) and indicates a possible error in the Euler solutions. However, this probably relates to the highly variable magnetic signal resulting from the complicated structure associated with the dykes. The NE-SW striking Pilanesberg dykes (mostly green and blue in Figure 9.9) show a more realistic depth of 5-15 m depth, with some possibly extending to 15-25 m, concurring with the depth of overburden (Rompel et al., 2002). Finally, a positive WNW trending anomaly, with poor Euler solutions and running parallel to sub-outcrop of the UG2 chromitite, was highlighted in the vertically sun- shaded aeromagnetic data (shown in yellow in Figure 9.9). This anomaly may relate to a dyke running parallel to the UG2 sub-outcrop, or to a harzburgite layer below the UG2, where weathering has caused olivine alteration to magnetite. If the anomaly is related to a harzburgite layer, it is a useful marker horizon of the UG2 (Campbell, 2006). If the anomaly is a dyke, it is not related to the Pilanesberg dyke swarm, due to its positive magnetic anomaly and the change in orientation (i.e. WNW trending, not NW-NNW trending). No magnetic analysis of the IRUP bodies was performed. However, Viljoen et al. (1986) note that IRUPs are often spatially related to potholes, another loss-of-ground feature common to the Merensky Reef and UG2. Potholes show an excursion of the economic chromitite units into the footwall stratigraphy. Excursions may be 1 m to 100 m down and the diameter of potholes can be up to 500 m. Magnetic surveys can model shallow potholes, but gravity surveys are not well suited to imaging these features, as the affected high density chromitite layers are commonly only 40 ? 60 cm thick (Kiefer and Viljoen, 2006). Hence, the gravitational effect is negligible or falls below the accuracy of the instrument. Despite only imaging IRUPs qualitatively by vertical sun-shading, the area may also be susceptible to potholes due to their established spatial relationship with IRUPs. 197 9.3 Interpretation of Local Modelled Features 9.3.1 Ground Gravity Survey Results The overburden thickness of the ground gravity survey was determined from 3D inversion, constrained by 56 boreholes, in the 100 m x 100 m station- and line- spacing survey grid. The thicknesses determined correlate to structure in the bedrock (i.e. faults and dykes). These structures were independently mapped from geological field data, and the locations were imported to overlay the overburden thickness (Figure 9.10). The structure includes dextral faults, joints and dykes delineated from aeromagnetic data (Rompel et al., 2002). Dextral Faults Elands River Outcrop in dry riverbed Fault ?Saddles? related to minor faults ?Valley? Village Pilanesberg dyke Figure 9.10 ? Overburden thickness (in colour, transparent) based on boreholes, with warm colours for thick overburden and cool colours for thick overburden. Overlays QuickBird orthophotograph and with plotted dextral faults (pink lines), unspecified faults (black lines) and Pilanesberg dykes (blue lines) based on geological mapping from boreholes and field exposure. 198 A dry riverbed, tributary to the Elands River and containing outcrop, follows a dextral fault. The fault had a likely control on the tributary, since faults allow the flow of water relative to the surrounding bedrock, and also weather faster as the rock immediately surrounding the fault is not as consolidated. Other subtle features are also distinguishable in the overburden, relating to lineaments from geological mapping and remote sensing (Rompel et al., 2002). These features include two ?saddles?, where overburden thickness increases (~1 m) in relation to minor faults, and a ?valley? in the western section, where the overburden thickness is slightly thicker relative to surroundings, which is spatially correlated to another minor fault. 9.3.2 High Resolution Ground Gravity Survey Results The 30 m x 30 m station- and line-spacing ground gravity grid was also interpreted for high frequency anomalies. For the given station spacing, the Nyquist wavelength is calculated at 60 m, i.e. the minimum determinable wavelength response is 60 m. A regional gradient was removed from the Bouguer gravity data, producing a residual map. This residual data showed two distinctive anomalies: a negative anomaly (?g ? -0.400 mgal) running parallel to a dextral fault, and a positive anomaly (?g ? 0.350 mgal) over a small IRUP (Figure 9.11), (Rompel et al., 2002). 199 Figure 9.11 ? Residual Bouguer gravity grid (colour) from 30 m x 30 m line- and station-spacing survey, overlaying aeromagnetic data (greyscale). The surrounding outline of the 100 m x 100 m line- and station-spacing ground gravity survey (black line) also plotted. Features mapped in geological field survey and from aeromagnetic data: dextral faults (pink lines), IRUPs (red polygons) and delineated dykes (blue). (inset) Profile across the linear gravity low relating to dextral fault (density values from Mar? et al. (2002) and density inversion results). The NW-SE striking linear gravity low in the north-east of the grid corresponds to a dextral fault mapped during a geological field survey (Rompel et al., 2002). Faults act as good transporters of water; hence the surrounding rock is more susceptible to weathering, creating an increase in the overburden thickness, from 2 m to 10 m thick. The IRUP shows an expected gravity high anomaly due to the high density of the body. The advantages of closely spaced stations are clearly seen in this region: the above change in overburden thickness could not have been so well mapped if the line- and Dextral faults Modelled line IRUP -0.1 -0.2 -0.3 -0.4 2.0 4.0 6.0 8.0 10.0 m mgal 100 200 300 Distance (m) 200 station-spacing had not been so closely spaced. High frequency gravity anomalies are apparent and relate to geological features in the region. The two prime examples are the weathered channel, showing an 8 m change in overburden thickness, and the IRUP anomaly which would be aliased with larger station spacing. This is station- and line-spacing is commonly too detailed for most areas: the trade-off for such high frequency data is that a field survey with similar grid parameters is very time consuming, leading to increased cost. However, accurate mapping of regions of rapidly changing overburden can be critical in seismic static corrections, in order to properly process and interpret the seismic data. Such a survey would be highly valuable to seismic surveys over important areas (e.g. shaft positions). 9.4 Gravity Modelling as a Seismic Tool Results of the gravity inversion provide a number of insights into the region and of the application of the inversion method. The initial seismic interpretation of the Bushveld Complex/Transvaal Supergroup contact shows a distinct change in slope, from steeply dipping close to the Rustenburg Fault to flattening out towards the north-east. Geometry inversions of the Transvaal contact are consistent with this interpolated seismic surface; hence the flattening is also supported by the gravity inversions. This consistency between the gravity inversions and seismic interpretations suggests that the densities and velocities used in the seismic interpretation are reasonable. This is an important result as in the deeper regions, where the dip flattens out, no borehole information was available to constrain either the seismic or gravity interpretations. In tests where the dip of the Transvaal contact was systematically different from the seismic interpretation, the gravity data could not be fitted. The flattening of the Bushveld Complex indicates the possibility that the complex is closer to surface, in areas proximal to the study region, than previously expected 201 (Figure 9.12). Similar flattening is seen in a southern region of the Bushveld Complex from 2D seismic lines (Odgers et al., 1993), which shows that this is not an isolated occurrence. Hence, the Bushveld Complex may be shallower towards its centre, hidden beneath the Lebowa Granite Suite and cover sequences (e.g. Karoo Supergroup). This is especially true in areas which have experienced folding, up- thrusting, doming or other tectonic uplift events. Distance (km) Depth (km) 0 2 4 6 8 10 a) 0 BC (UZ) 1 BC (MZ) 2 BC (CZ/LZ) 3 TVL 4 Basement 5 Figure 9.12 ? Comparison of a) dipping Bushveld Complex, and b) proposed flattening Bushveld Complex. The dipping Bushveld Complex shows the economic Critical Zone dipping into the Earth, beyond the depths of current shaft technology (i.e. >4 km). However, the flattening Bushveld Complex shows the dip of the Critical Zone decreasing down-dip, making it shallower than previously expected and within minable range. Note abbreviations: Bushveld Complex (BC), Upper Zone (UZ), Main Zone (MZ), Critical Zone (CZ), Lower Zone (LZ), Transvaal Supergroup (TVL). b) 0 BC (UZ) BC (MZ) 1 BC (CZ/LZ) 2 TVL 3 202 If the Bushveld Complex is shallower than previously assumed in these areas, it may be possible that the depth of the Critical Zone is within range of modern-day shaft technology. Flattening also indicates that the Bushveld Complex is not part of a steeply dipping sheet, which in turn supports the theory of a continuous Bushveld Complex (Webb et al., 2004). Airborne gravity surveys and the gravity inversion process show various advantages, both in terms of cost and modelling capability. The airborne method is cheaper, taking considerably less time to fly a 10 km x 10 km area than a seismic survey over the same area. The inversion method is capable of modelling vertical bodies (e.g. pipe-like structures, IRUPs) in 3D which are difficult to image, even with 3D seismic data. These vertical bodies? parameters, including dimensions, density and depth-to- top, can be optimised to fit the data, as well as provide a range of possible solutions. AGG surveys present the most potential, as long wavelength features can be modelled using converted Gz data, and detailed short wavelength features can be modelled using the FTG data. With continuous improvements to 3D FTG inversion packages, software will soon be readily available to model near-surface features. However, the method also has intrinsic disadvantages, related to potential field modelling. The user must be fully aware that potential field data can be fitted by an infinite number of solutions (i.e. geological models). In order to limit solutions, all available geological and geophysical data must be interpreted to create the geological starting model. The model also shows significantly less resolution, especially compared to seismic interpretations, but may still be useful in adjacent areas without detailed seismic coverage. In honouring the seismic interpretation, it is clear that the presented gravity modelling method can be carried out in blocks adjacent to the seismic survey. This model would provide a first-pass at the structural geometry of the high density contrast contact, aiding in the planning of future seismic surveys. The model would also 203 provide an approximate depth of the Bushveld Complex/Transvaal Supergroup contact, which would aid in picking the correct reflector during the seismic interpretation. 9.5 Conclusions The final inverted 3D gravity model shows features which agree with the historical far-field stresses along the TML. The model of open folds in the region (axes trending NW) allows for an anticline along the south-western edge of the study area. Tilting of the region, towards the centre of the Bushveld Complex during thermal collapse, and erosion explains the flattening of the Transvaal Supergroup contact towards the north-east of the region. Together with 2D gravity gradient modelling, a change in contact between the Bushveld Complex and Transvaal Supergroup can be inferred: from a fault-bounded, onlapping contact in the west of the study area, to being layered with the Transvaal Supergroup in the southern portion. Loss-of-ground features such as faults, dykes and IRUPs (with spatially inferred potholes) are modelled from the geophysical data. The overburden thickness agrees with the depth of dykes in the region (ranging from 0-14 m), showing changes across various faults and dykes. The gravity inversion method, using Gz data, has provided a model that supports the seismic interpretation. As a first-pass model in adjacent blocks, it could be capable of approximating the depth of the Transvaal Supergroup contact. Used together with available geological and geophysical data, it also aids in planning future seismic surveys. 204 CHAPTER 10: CONCLUSION AND SUGGESTIONS FOR FURTHER WORK Vertical component gravity (Gz) data, converted from AGG data and incorporating ground Gz data, were used to run property and geometry inversions over a portion of the western Bushveld Complex. The final inverted model was used, together with a comprehensive suite of geophysical and geological data, to create an integrated 3D model of the region. The main features of the 3D model include the Transvaal Supergroup contact, with associated folds, and various loss-of-ground features, particularly IRUPs, dykes and faults. A 3D geological starting model of the Bushveld Complex and Transvaal Supergroup contact was created from interpreted seismic data, borehole data, geological and rock property data. While the upper portion of the model is well constrained by boreholes and seismic data, the deeper sections have no borehole control. Hence the densities and velocities determined from seismic data are not well constrained. Thus the gravity data also provided constraints in the geometrical and property inversions of the deeper sections of the data. An initial forward model, comparing the gravity response of the starting geological model with the measured data, showed a RMS misfit of 2.431 mgal. Seven inversions steps were carried out, producing a final model with an unchanging RMS misfit of 1.076 mgal. Key features of the final model are: a steeply dipping Transvaal Supergroup contact in the south-west, which flattens to the north-east; a heterogeneous Transvaal Supergroup (density ranging from 2.595 g.cm-3 to 2.630 g.cm-3, consistent with measured data); and a heterogeneous Bushveld Complex (density ranging from 2.870 g.cm-3 to 2.980 g.cm- 3, also consistent with measured data). Twelve IRUPs were identified and modelled, with inversions producing homogeneous bulk densities which ranged from 2.773 g.cm-3 to 3.185 g.cm-3, and an average of 3.091 g.cm-3. 205 The gravity model was also used to approximate and remove the regional response from the ground gravity survey. This was used to determine the thickness of the overburden (ranging from 0 ? 14 m thick, constrained by boreholes). Several regions of rapidly varying overburden thicknesses were identified in the 100 m x 100 m station- and line-spacing survey, relating to mapped faults and joints. The detailed 30 m x 30 m station- and line-spacing survey also revealed a small IRUP and a channel of increased overburden thickness, related to a fault. This channel would not have been identified without such close station spacing. This is an important result for seismic static corrections, which need accurate overburden thicknesses, especially in regions requiring highly detailed interpretations (e.g. over future shaft locations). The gravity model of the study area has been integrated with other geophysical data to provide a 3D geological model of the Bushveld Complex/Transvaal Supergroup contact (which changes from a fault-bounded, onlapping Bushveld Complex, to a Bushveld Complex layered with the Transvaal), density variations in the Bushveld Complex and Transvaal Supergroup and loss-of-ground features in the region (e.g. dykes, faults, IRUPs). Aeromagnetic data were enhanced using a variety of sun- shading angles (with vertical inclination highlighting the edges of dykes and IRUPs) to delineate magnetic features. Euler deconvolution was run on the dataset to automatically determine the approximate depth of the tops of dykes. There are two major sets of dykes: E-W trending, structurally complex dyke from late-stage activity, and NW-trending dykes, assumed to be related to the Pilanesberg intrusion. Basic 2.5D profile modelling of the negative Pilanesberg dykes showed the Euler depth solutions to be compatible (commonly 5-15 m), with the dykes being sub-vertical and varying from 5-9 m thick. The dykes were also, clearly, strongly remanent. A positive magnetic lineament running parallel to the UG2 could be a dyke or could relate to an underlying harzburgite, where olivine has altered to magnetite. Major geological features in and around the study region relate to tectonic activity and associated motion along the TML. The primary features include extensional 206 NW-trending faults (e.g. Rustenburg Fault) and a series of open folds, axes trending NW with wavelengths of ~20 km. The Rustenburg Fault and associated extensional faults were a result of NW-SE far-field stress, causing dextral motion along the TML, pre-Bushveld emplacement. The series of open folds in the Bushveld Complex and Transvaal Supergroup is likely a response to NE-SW far-field stress, causing sinistral motion along the TML, syn- to post-Bushveld emplacement. By including an anticline to the south-west of the study area, and accounting for tilting and erosion, it is possible to explain the shallowing of the contact between the Bushveld Complex and Transvaal Supergroup. If similar folding applies to the economic Critical Zone, it is possible that these layers occur at shallower depths than previously thought, especially around regions of tectonic uplift events (e.g. up-thrusting, doming). Secondary faults, axes trending NE and wavelength ~5 km, are also apparent in the Transvaal Supergroup, relating to a NW-SE far-field stress, possibly pre-Bushveld intrusion. Future work can investigate testing these gravity inversions outside of the region constrained by seismic data. The region should, however, contain seismic data in order to compare the gravity inversion. The change in Bushveld Complex/Transvaal Supergroup contact, from fault-bounded, onlapping to layered, could also be investigated in order to improve the model. 207 APPENDIX A Units Note on Gravity Units All vertical gravitational component values, Gz, provided in the text are in gals (1 gal = 1 cm.s-2), named after Galileo (the first person to investigate gravity). The most commonly used unit in exploration is the milligal (1 mgal = 0.001 gal = 10-5 m.s-2). Hence, the Earth?s gravitational acceleration is ~ 980 gal or ~ 980 000 mgal. The FTG gravity data are presented in units of E?tv?s (1 Eo = 10-4 mgal.m-1 = 0.1 mgal.km-1 = 10-9 s-2). The unit is named after Baron Roland von E?tv?s, who invented the torsion balance, which accurately measured the curvature of the gravitational field (Heath, 2007). Note on magnetic units The international unit of magnetic field intensity is the Tesla (T). The Tesla is, however, too large for exploration purposes and, instead, the nanotesla (1 nT = 10-9 T) is used. Geophysical anomalies are commonly between 10?s and 10 000?s nT (Roux, 1980). Note on Datum and Projection Unless otherwise stated, the datum and projection used in the presented maps throughout this dissertation are WGS 1984, LO 27. However, in certain cases, units have been removed or altered to preserve the confidentiality of the location of the data. 208 APPENDIX B Gravity Response of a Vertical Prism Instead of computing the gravitational response of a rectangle, the inversion algorithm calculates the response of a vertical prism. Consider a point outside an arbitrary body of uniform density, ?. The vertical component of gravity, gz, can be written as a surface integral (Fullagar et al., 2004), where: ds r nz Gg z ? ??= r ??? (B.1) where is the unit vector in the vertical direction, is the outward normal to the body surface, and z? n? r r is the distance from the observation point to a point on the surface of the body. For a vertical rectangular prism, with sides parallel to the x- and y-axes, it is only the top and bottom faces of the prism which contribute to gz. Coggon (1976) showed that the vertical gravity component, in this case F(x,z,y1,y2), of an individual prism is reduced to a summation of four definite line integrals: ? += 2 1 )ln(),,,( 21 y y dyrxGyyzxF ? (B.2) along the y-parallel edges of the top and bottom of the prism, where ? is the density of the prism. Fullagar (1975) derived the analytic expression for F: 2 1 1 21 tan.2)ln(.)ln(.),,,( y yz yrx zryxyrxyGyyzxF ?? ??? ? ?? ??? ? ?+?++?+= ?? (B.3) which simplifies to: [ ] ?? ? ??? ? ??? ? + +++= 11 22 21 ln.)ln(.),,,( 2 1 ry ry xrxyGyyzxF yy? ?? ?????? ? ??? ? +++ ?? ? ))(( )( tan.2 21 2 121 xuxuz uuz z (B.4) where u = r ? y (Fullagar et al., 2004). 209 REFERENCES Antoine, G., 2004. Gravity gradiometer modelling over Styldrift: Memo for Anglo American, 13 pp. Ashwal, L.D., Webb, S.J. and Knoper, M.W., 2005. Magmatic stratigraphy in the Bushveld Northern Lobe: continuous geophysical and mineralogical data from the 2950 m Bellevue drillcore. South African Journal of Geology, 108(2): 199-232. Bell Geospace, 2006. 3D FTG Geophysical Surveys - Airborne and Marine Gravity Gradiometry. Bell Geospace, available at http://www.bellgeo.com/tech/overview.html (26 July 2006). Bell, R.E., Coakley, B.J. and Stemp, R.W., 1991. Airborne gravimetry from a small twin engine aircraft over Long Island Sound. Geophysics, 56(9): 1486-1493. BHP Billiton, 2007. Orion Operations: FALCONTM Instrumentation. BHP Billiton, available at http://falcon.bhpbilliton.com/falcon/ (8 November 2007). Blackman, R.B. and Tukey, J.W., 1959. The measurement of power spectra. Dover publications, New York. Blakely, R.J., 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge University Press, Cambridge. Bott, M.H.P., 1960. The use of rapid digital computing methods for direct gravity interpretation of sedimentary basins. Geophys J. R. astr. Soc., 3: 63-67. Brett, J., 2000. The gravity gradiometer instrument: Theory of FTG measurements. Bell Geospace, available at http://www.bellgeo.com/tech/measurements.html (26 July 2006). Buick, I.S., Maas, R. and Gibson, R., 2001. Precise U-Pb titanite age constraints on the emplacement of the Bushveld Complex, South Africa. Jour. Geol. Soc., London, 158: 3-6. Bumby, A.J., Eriksson, P.G. and van der Merwe, R., 1998. Compressive deformation in the floor rocks to the Bushveld Complex (South Africa): evidence from the Rustenburg Fault Zone. Journal of African Earth Sciences, 27(3): 307-330. Cairncross, B. and Dixon, R., 1995. Minerals of South Africa. Geological Society of South Africa, Johannesburg, 290 pp. Campbell, G., 1990. The seimic revolution in gold and platinum prospecting. SAGA Yearbook, 1(1): 37-45. Campbell, G., 2006. High resolution aeromagnetic mapping of "loss-of-ground" features at platinum and coal mines in South Africa. South African Journal of Geology, 109: 439-459. Cawthorn, G. (Editor), 1988. The Geology of the Pilanesberg. Information Series, 1. National Parks Board, Mafikeng, South Africa, 24 pp. Cawthorn, R.G., 1998. Geometrical relations between the Transvaal Supergroups, the Rooiberg Group and the mafic rocks of the Bushveld Complex. S. Afr. J. Geol., 101(4): 275-279. 210 Cawthorn, R.G., Cooper, G.R.J. and Webb, S.J., 1998. Connectivity between the western and eastern limbs of the Bushveld Complex. S. Afr. J. Geol., 101(4): 291-298. Cawthorn, R.G. and Walraven, F., 1998. Emplacement and Crystallization time for the Bushveld Complex. Jour. Petrology, 39(9): 1669-1687. Coggon, J.H., 1976. Magnetic and gravity anomalies of polyhedra. Geoexploration, 14: 93-105. Cooper, G.R.J., 2000. SignProc for Windows, Version 1.56 (Cooper, G.R.J., University of the Witwatersrand, Johannesburg, South Africa). Cooper, G.R.J., 2003. Grav2DC, Version 2.01 (Cooper, G.R.J., University of the Witwatersrand, Johannesburg, South Africa). Cooper, G.R.J., 2004a. Euler deconvolution applied to potential field gradients. Exploration Geophysics, 35: 165-170. Cooper, G.R.J., 2004b. Euler, Version 1.15 (Cooper, G.R.J., University of the Witwatersrand, Johannesburg, South Africa). Cooper, G.R.J., 2004c. Numerical Methods in Geophysics: Space Domain Filtering, Geophysics Honours Lecture Notes, 97 pp. Cooper, G.R.J., 2004d. The stable downward continuation of potential field data. Exploration Geophysics, 35: 260-265. Cooper, G.R.J., 2006. Obtaining dip and susceptibility information from Euler deconvolution using the Hough transform. Computer and Geosciences, 32(10): 1592-1599. Cooper, G.R.J. and Olavsdottir, B., 2004. Advanced potential field data processing (course notes, March 2004), Uppsala, Sweden. Cordell, L. and Henderson, R.G., 1968. Iterative three-dimensional solution of gravity anomaly data using a digital computer. Geophysics, 33(4): 596-601. Corner, B., 1993. Lecture notes for geophysics honours students on gravity and magnetic potential fields, University of the Witwatersrand, Johannesburg. Davidson, G.E. and Chunnett, G.K., 1999. Seismic exploration for Merensky Reef: The way ahead. South African Journal of Geology, 102: 261-267. de Wit, M.J. et al., 1992. Formation of an Archaean Continent. Nature, 357: 553 - 562. Dean, W.C., 1958. Frequency analysis for gravity and magnetic interpretation. Geophysics, 23: 97-127. Dransfield, M.H., 1994. Airborne Gravity Gradiometry, PhD thesis, University of Western Australia, 254 pp. Dransfield, M.H., Christensen, A., Rose, M., Stone, P. and Diorio, P., 2001. FALCON test results from the Bathurst Mining camp. Exploration Geophysics, 32: 243-246. Dransfield, M.H. and Lee, J.B., 2004. The FALCON(R) airborne gradient gravity survey systems. In: R.J.L. Lane (Editor), Airborne Gravity 2004 - Abstracts from the ASEG-PESA Airborne Gravity 2004 Workshop. Geoscience Australia Record, Sydney, pp. 15-20. 211 du Plessis, A. and Levitt, J.G., 1987. On the structure of the Rustenburg Layered Suite - insight from seismic data, Abstracts Indaba on the Tectonic Setting of Layered Intrusives, Pretoria, South Africa. Du Plessis, C.P. and Walraven, F., 1990. The tectonic setting of the Bushveld Complex in Southern Africa, Part 1. Structural deformation and distribution. Tectonophysics, 179: 305-319. Dyke, A.L., Harman, P. and Mahanta, A.M., 2002. Falcon spreads its wings. Preview, 99: 25-28. Eales, H.V. and Cawthorn, R.G., 1996. The Bushveld Complex. In: R.G. Cawthorn (Editor), Layered Intrusions. Elsevier, Amsterdam, pp. 181-229. Eglington, B.E. and Armstrong, R.A., 2004. The Kaapvaal Craton and adjacent orogens, southern Africa: a geochronological databaseand overview of the geological development of the craton. S.Afr.J.Geol, 107: 13-32. Fairhead, D. and Odegard, M.E., 2002. Advances in gravity survey resolution. The Leading Edge, January: 26-27. Ferguson, J., 1973. The Pilanesberg Alkaline Province, Southern Africa. Trans. Geol. Soc. S. Afr., 87: 249-270. Fitzgerald, D.J. and Holstein, H., 2006. Innovative data processing methods for gradient airborne geophysical data sets. The Leading Edge, 25(1): 87-94. Fourie, C.J.S., 1998. Spreadsheet reductions of Gravity Data. Computers and Geoscience, 24(5): 495-499. Fugro, 2005. Services > Survey Technology > Aeromagnetics > MIDASTM II Helicopter-borne Horizontal Gradiometer. Fugro Airborne Surveys, available at http://www.fugroairborne.com/service/midas.php (8 March 2007). Fullagar, P.K., 1975. Interpretation for underground gravity profiles from the North Broken Hill Mine, B.Sc. (Hons), Monash University (unpubl.). Fullagar, P.K., 2007a. Email: Outline of the inverse formulation in VPmg (personal communication, 2 June 2007). Fullagar, P.K., 2007b. VPmg 3D potential fields inversion (Powerpoint presentation). Fullagar Geophysics Pty. Ltd., pp. 45. Fullagar, P.K., Hughes, N.A. and Paine, J., 2000. Drilling-constrained 3D gravity interpretation. Exploration Geophysics, 31(1): 17-23. Fullagar, P.K. and Pears, G., 2007. Towards geologically realistic inversion, Exploration 2007, Toronto. Fullagar, P.K., Pears, G., Hutton, D. and Thompson, A., 2004. 3D gravity and aeromagnetic inversion for MVT Lead-Zinc exploration at Pillara, Western Australia. Exploration Geophysics, 35(2): 142-146. Garmin, 2007. Garmin: What is GPS? GarminTM, available at http://www.garmin.com/aboutGPS/ (13 March 2007). Geometrics, 1995. G-856 Memory-MagTM Proton Precession Magnetometer Operator's Manual, California, USA, 61 pp. Geosoft, 2006. Oasis Montaq, Standard Edition v6.4 (9M) Help Notes. Gibson, M., Pretorius, C. and Lionnet, M., 2006. The 2005 Styldrift 3D seismic survey, Anglo Technical Division, Johannesburg, South Africa. 212 Good, N. and de Wit , M.J., 1997. The Thabazimbi-Murchison Lineament of the Kaapvaal Craton, South Africa: 2700 Ma of episodic deformation. J. geol. Soc. Lond., 154: 93-97. Good, N. and de Wit, M.J., 1997. The Thabazimbi-Murchison lineament of the Kaapvaal craton, South Africa: 2700 Ma of episodic deformation. Jour. Geol. Soc. London, 154: 93-97. Gough, D.I., 1957. Study of the paleomagnetism of the Pilanesberg dykes. Monographs At. Royal Australian Society Geophysics.(Supplement 7): 193- 213. Hammer, S., 1983. Airborne gravity is here! Geophysics, 48(2): 213-223. Hammond, S. and Murphy, C., 2003. Air-FTGTM: Bell Geospace's airborne gravity gradiometer - a description and case study. ASEG Preview, 105: 24-26. Hatch, D., Murphy, C., Mumaw, G.R. and Brewster, J., 2006. Performance of the Air-FTG? system aboard an airship platform, AESC 2006, Melbourne, Australia. Heath, P., 2007. Analysis of potential field gradient tensor data: forward modelling, inversion and near-surface exploration, University of Adelaide, Adelaide, 206 pp. Hinks, D., McIntosh, S. and Lane, R.J.L., 2004. Comparison of the Falcon(R) and the Air-FTGTM airborne gravity gradiometer systems at the Kokong Test Block, Botswana. In: R.J.L. Lane (Editor), Airborne Gravity 2004 - Abstracts from the ASEG-PESA Airborne Gravity 2004 Workshop. Geoscience Australia Record, Sydney, pp. 125-134. Hsu, S.K., 2002. Imaging magnetic sources using Euler's equation. Geophysical Prospecting, 50: 15-25. Hwang, C. et al., 2007. Geodetic and geophysical results from a Taiwan airborne gravity survey: Data reduction and accuracy assessment. J. Geophys. Res., 112. Kiefer, R. and Viljoen, M.J., 2006. PGE exploration targets to the west of the Pilanesberg, South Africa. South African Journal of Geology, 109: 459-474. Kruger, F.J., 2005. Filling the Bushveld Complex magma chamber: lateral expansion, roof and floor interaction, magmatic unconformities, and the formation of giant chromitite, PGE and Ti-V-magnetitite deposits. Mineralium Deposita, 40: 451-472. LaCoste, L.J.B., 1934. A new type long period vertical seismograph. Physics, 5(July 1934): 178-180. Lane, R.J.L. (Editor), 2004a. Airborne Gravity 2004 - Abstracts from the ASEG- PESA Airborne Gravity Workshop 2004. Geoscience Australia Record. Lane, R.J.L., 2004b. Integrating ground and airborne data into regional gravity compilations. In: R.J.L. Lane (Editor), ASEG-PESA Airborne Gravity 2004 Workshop. Geoscience Australia Record, pp. 81-98. Lee, J.B., 2001. FALCON gravity gradiometry technology, ASEG 15th Geophysical Conference and Exhibition, August 2001, Brisbane, Australia. Lee, J.B. et al., 2006. First test survey results from the FALCONTM helicopter-borne airborne gravity gradiometer system, AESC 2006, Melbourne, Australia. 213 Letts, S.A., 2004. The development of aeromagnetic mapping and dyke interpretation methodologies with application to a high resolution aeromagnetic survey in the eastern Bushveld Complex, University of the Witwatersrand, Johannesburg, 217 pp. Letts, S.A., 2007. Palaeomagnetic significance of the Bushveld Complex and related 2 Ga magmatic rocks in ancient continental entities, University of the Witwatersrand, Johannesburg, 228 pp. Li, Y. and Oldenburg, D.W., 1998. 3-D inversion of gravity data. Geophysics, 63(1): 109-119. Mar?, L.P., Oosthuizen, B.C. and Tabane, L.R., 2002. South African Geophysical Atlas Volume IV: Physical Properties of South African Rocks, Council for Geoscience, Pretoria. Maus, S. and MacMillan, S., 2005. 10th generation International Geomagnetic Reference Field. EOS Trans. AGU, 86(159). Merry, C.L., 1999. Fundementals of the global positioning system, South African Geophysical Association, Rondebosch. Mira Geoscience, 2006. Constrained inversion using Gocad and VPmg. Mira Geoscience Limited, Vancouver, Canada, 62 pp. Murphy, C., 2004. The Air-FTG? airborne gravity gradiometer system. In: R.J.L. Lane (Editor), Airborne Gravity 2004 - Abstracts from the ASEG-PESA Airborne Gravity 2004 Workshop. Geoscience Australia Record, Sydney, pp. 7-14. Murphy, C., Brewster, J. and Robinson, J., 2006. Evaluating Air-FTG(R) survey data: bringing value to the full picture, Australian Earth Sciences Convention 2006, Melbourne, Australia. Mushayandebvu, M., van Driel, P., Reid, A.B. and Fairhead, J.D., 2001. Magnetic source parameters of 2D sturctures using extended Euler deconvolution. Geophysics, 66(3): 814-823. Nabighian, M.N. et al., 2005. Historical development of the gravity method in exploration. Geophysics, 70(6): 63-89. Nelson, J.B., 1988. Calculation of the magnetic gradient tensor from total field gradient measurements and its application to geophysical interpretation. Geophysics, 53(7): 957-966. Odgers, A.T.R., Hinds, R.C. and Von Gruenwaldt, G., 1993. Interpretation of a seismic reflection survey across the southern Bushveld Complex. S. Afr. J. Geol., 96: 205-212. Parker, R.L., 1972. The rapid calculation of potential anomalies. Geophys. Jour. Roy. Astro. Soc., 31: 447-455. Pears, G., 2008. Email: Final VPmg request! (personal communication, 4 January 2008). Pelton, C., 1987. A computer program for hill shading digital topographic data sets. Computer and Geosciences, 13(5): 545-548. Reczko, B.F.F., Antoine, L.A.G. and Eriksson, P.G., 1997. Three-dimensional computer-assisted basin modelling to generate exploration target areas: an 214 example from the late Archaean-early Proterozoic Transvaal Supergroup, South Africa. Mineralium Deposita, 32: 392-400. Reid, A.B., 1980. Aeromagnetic survey design. Geophysics, 45(5): 973-976. Reid, A.B., Allsop, J.M., Granser, H., Millett, A.H. and Somerton, I.W., 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55(1): 80-91. Reid, D.L. and Basson, I.J., 2002. Iron-rich ultramafic pegmatite replacement bodies within the Upper Critical Zone, Rustenburg Layered Suite, Northam Platinum Mine, South Africa. Miner. Mag., 66(6): 895-914. Rompel, A.K.K., du Plessis, S.J., Luck, E. and Gibson, M., 2002. Integration and Interpretation of Geophysical Datasets for the Styldrift JV, Western Bushveld. Geosciences Resource Group, Anglo Technical Division, pp. 65. Roux, A.T., 1980. The Magnetic Method. Geophysical Field Manual for Technicians, No. 1. South African Geophysical Association, Johannesburg, South Africa, 71 pp. Rummel, R., 1986. Satellite Gradiometry, Mathematical and Numerical Techniques in Physical Geodesy. Springer Verlag, pp. 317-363. Schmidt, P.W. and Clark, D.A., 2000. Advantages of measuring the magnetic gradient tensor. ASEG Preview, April, 2000. Scintrex, 1995. CG?3/3M Autograv Automated Gravity Meter - Operator Manuel, Ontario, Canada, 198 pp. Scoon, R.N. and Mitchell, A.A., 2004. Petrogenesis of Discordant Magnesian Dunite Pipes from the Central Sector of the Eastern Bushveld Complex with Emphasis on the Winnaarshoek Pipe and Disruption of the Merensky Reef. Econ. Geol., 99: 517-541. Silver, P.G., Fouch, M.J., Gao, S.S., Schmitz, M. and Group, K.S., 2004. Seismic anisotrpy, mantle fabric and the magmatic evolution of Precambrian Southern Africa. South African Journal of Geology, 107: 45-58. Sutcliffe, P.R., 2004. Lecture notes for the Hermanus magnetic winter school, Hermanus Magnetic Observatory, Hermanus. Syberg, F.J.R., 1972. a Fourier Method For The Regional-Residual Problem Of Potential Fields. Jeophysical Prospectors, 20: 47-75. Telford, W.M., Geldart, L.P. and Sheriff, R.E., 1990. Applied Geophysics. Cambridge University Press, New York, 770 pp. Teskey, D., 1991. Guide to aeromagnetic specifications and contracts, Geological Survey of Canada, Open File 2349, pp. 73. Thompson, D.T., 1982. EULDPH: a new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47(1): 31-37. Trimble, 1997. 4600LS Surveyor: Operation Manual. Trimble Navigation Limited, Sunnyvale, California, 108 pp. Trimble, 2002. 4600 LS: Economical, fully integrated single-frequency GPS survey unit, Trimble, Dayton, Ohio. Trimble, 2006. Trimble 5700 GPS system: Datasheet, Trimble, Dayton, USA. 215 216 USGS, 2005. USGS National Geomagnetism Program: A brief introduction to geomagnetism. USGS, available at http://geomag.usgs.gov/intro.php (5 July 2007). Vassiliou, A.A., 1986. Comparison of methods for the processing of gravity gradiometer data, 14th Annual Gravity Gradiometry Conference, Colorado, USA. Vermaak, C.F., 1976. The Merensky Reef - thoughts on its environment and genesis. Econ. Geol., 71: 1270-1298. Viljoen, M.J. et al., 1986. The Amandelbult Section of Rustenburg Platinum Mines Limited, with reference to the Merensky Reef. In: C.R. Anhaeusser and S. Maske (Editors), Mineral Deposits of South Africa. Vol. II. Geological Society of South Africa, pp. 1107-1134. Viljoen, M.J. and Reimold, W.U., 2002. An Introduction to South Africa?s Geological Heritage. Mintek in association with the Geological Society of South Africa, Randburg, 193 pp. Walraven, F., 1974. Tectonism during emplacement of the Bushveld Complex, South Africa and the resulting fold structures. Trans. Geol Soc. S. Afr., 77(323-328): 323-328. Walraven, F., 1976. Notes of the late-stage history of the western Bushveld Complex. Trans. Geol Soc. S. Afr., 79(1): 13-21. Walraven, F., 1981. 2526 Rustenburg, 1:250 000 Geological Series. Dept. of Mineral & Energy Affairs, Geological Survey, Pretoria, South Africa. Walraven, F. and Darracott, B.W., 1976. Quantitative interpretation of a gravity profile across the western Bushveld Complex. Trans. geol. Soc. S. Afr., 79: 22-26. Webb, S.J., Nguuri, T.K., Cawthorn, R.G. and James, D.E., 2004. Gravity modelling of Bushveld Complex connectivity supported by southern African seismic experiment results. South African Journal of Geology, 107(1/2): 207-218. Whitehead, N. and Musselman, C., 2006. montaj Gravity and Terrain Correction (Tutorial and User Guide). Geosoft, Toronto, Canada, 50 pp. Zhang, C., Mushayandebvu, M., Reid, A., Fairhead, D. and Odegard, M.E., 2000. Euler deconvolution of gravity tensor gradient data. Geophysics, 65(2): 512- 520. Zuidweg, K. and Mumaw, G.R., 2007. Airborne gravity gradiometry for exploration geophysics: The first five years, International Gravity Field Service General Assembly 2007, Istanbul, Turkey.