Gigayear stability of cratonic edges controls global distribution of sediment-hosted metals

Karol Czarnota*, Mark J. Hoggard*, Fred D. Richards, David L. Huston & A. Lynton Jaques 1. Geoscience Australia, GPO Box 378, Canberra ACT 2601, Australia. 2. Research School of Earth Sciences, Australian National University, Canberra, ACT 0200, Australia 3. Department of Earth and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA 02138, USA. 4. Lamont-Doherty Earth Observatory, Columbia University, 61 Rte 9W, Palisades, NY 10964, USA. *karol.czarnota@ga.gov.au; mark hoggard@fas.harvard.edu ***Add proper acknowledgment of ACS and define tau prime eta and tau s in the anelasticity section.*** 1 Sustainable development and transition to a clean-energy economy is placing ever-increasing de2 mand on global supplies of base metals (copper, lead, zinc and nickel). Alarmingly, this demand 3 is outstripping the present rate of discovery of new deposits, with significant shortfalls forecast in 4 the coming decades. Thus, to maintain growth in global living standards, dramatic improvements 5 in exploration success rate are an essential goal of the geoscience community. Significant quantities 6 of base metals have been deposited by low-temperature hydrothermal circulation within sedimen7 tary basins over the last 2 billion years. Despite over a century of research, relationships between 8 these deposits and geological structures remain enigmatic. Here, for the first time, we show that 9 85% of sediment-hosted base metals, including all giant deposits (> 10 megatonnes of metal), occur 10 within 200 km of the edges of thick lithosphere, mapped using surface wave tomography and a 11 parameterisation for anelasticity at seismic frequencies. This remarkable observation implies long12 term lithospheric edge stability and a genetic link between deep Earth processes and near-surface 13 hydrothermal mineral systems. This result provides an unprecedented global framework for iden14 tifying fertile regions for targeted mineral exploration, reducing the search-space for new deposits 15 by two-thirds on this lithospheric thickness criterion alone. 16 Consumption of base metals over the next ∼25 years is set to exceed the total produced in human history to 17 date. Moreover, trace metals (e.g. cobalt, indium and germanium) are often produced as by-products of base 18 metal mining and are essential in many high-tech applications. A growing concern is that the rate of exploitation 19 of existing reserves is outstripping discovery of new deposits, despite exploration expenditure tripling during the 2

In summary, this work illustrates a new and remarkably clear link between giant sediment-hosted base metal 140 mineral systems and the edges of thick lithosphere. Approximately 55% of the world's lead, 45% of its zinc and 20% 141 of known copper is found within ∼200 km of this edge. We have demonstrated the value of regional seismic arrays to 142 better resolve this edge and enhance the mineral exploration efforts required to sustain ongoing global development. 143 Significantly, deposit ages indicate that, following rifting, edges of thick-lithosphere are generally stable over billion- 144 year timescales. The far-reaching geodynamic and societal implications of our observations highlight the urgent 145 need for further research. To improve resolution of mapped lithospheric structure, higher fidelity seismic imaging 146 must be coupled with enhanced mantle xenolith coverage and tighter constraints on seismic anelasticity from 147 mineral physics experiments. More generally, these maps need to be integrated with models of basin dynamics, 148 surface processes and reactive transport modelling, and bench-marked against additional geological information, 149 such as sedimentary facies variations, tectonic structures and alteration zones. These multiple research strands 150 will yield fundamental new insights into sediment-hosted mineral systems and lead to substantial improvements in 151 exploration success rates. 152 Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure 1: Distribution of sediment-hosted and iron-oxide-copper-gold base metal deposits as a function of lithospheric thickness in Australia. (a) Carpentaria Zinc Belt; red/blue = variably reduced to pole aeromagnetic intensity data 13 ; grey polygons = generalised outcrop of Cretaceous marine sediments in Eromanga and Karumba Basins 28 ; black dashed contour = 170 km LAB thickness; symbols = deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 2 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clasticdominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed); stars = iron-oxide-copper-gold (IOCG). (b) LAB mapped by converting FR12 tomography 16 to temperature using an anelasticity parameterisation 15 calibrated on local paleogeotherms (Supplementary Material) and illuminated by free-air gravity anomalies 13 ; black/green crosses = geotherms used as constraints/tests in anelasticity calibration; box = location of panel (a).
Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 where ρ is the density and J 1 and J 2 represent real and imaginary components of the complex compliance, J * , 193 which is a quantity describing the sinusoidal strain resulting from the application of a unit sinusoidal stress. J 1 194 represents the strain amplitude in phase with the driving stress, whilst the J 2 component is π 2 out of phase, resulting 195 in dissipation. These terms are expressed as where A B = 0.664 and α B = 0.38 represent the amplitude and slope of background stress relaxation and J U is the 200 unrelaxed compliance. Parameters A P and σ P represent the amplitude and width of a high frequency relaxation 201 peak superimposed on this background trend such that

209
where µ 0 U is the unrelaxed shear modulus at surface pressure-temperature conditions, the differential terms are 210 assumed to be constant and the pressure, P , in GPa is linearly related to the depth, z, in km by z 30 . The Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 and τ M = η µ U is the normalised Maxwell relaxation timescale. τ P represents the normalised shear-wave period 213 associated with the centre of the high frequency relaxation peak, assumed to be 6 × 10 −5 . The shear viscosity, η, is where d is the grain size, m the grain size exponent (assumed to be 3), R the gas constant, E a the activation 216 energy and V a the activation volume. Subscripts [X] r refer to reference values, assumed to be d r = d = 1 mm, 217 P r = 1.5 GPa and T r = 1200 • C for the upper mantle. A η represents the extra reduction of viscosity due to an 218 increase in E a near the solidus, expressed as where T η is the homologous temperature above which activation energy becomes E a + ∆E a and γ = 5 is the factor 221 of additional reduction. λφ describes the direct effect of melt on viscosity, assumed to be negligible here. The where K 0 = 130 GPa is the bulk modulus at zero pressure and K T = 4.8 is the pressure-derivative of the isothermal 234 bulk modulus. The associated isothermal density change with pressure, ρ(P ), is given by

236
where ρ 0 = 3.33 Mg m −3 is the density of mantle at surface pressure and temperature. The effect of pressure on 237 Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 thermal expansivity is included according to where δ T = 6 is the Anderson-Grüneisen parameter. Thus, the final density, ρ(P, T ), can be calculated using

241
where T 0 = 273 K is temperature at the surface. In a similar manner to Equation (1), the shear-wave attenuation, 242 Q −1 S , can be defined as Xenolith and xenocryst thermobarometry. Temperature estimates across a range of depths are required 245 to generate a series of V S -T-P tie points in order to calibrate the regional seismic tomography models. We therefore

251
In these samples, we use a thermometer 39 that exploits exchange of calcium and magnesium between orthopy-252 roxene and clinopyroxene and a barometer 40 based upon aluminium exchange between orthopyroxene and garnet, 253 given by equation (5) of Nickel & Green (1985). This approach therefore requires compositions of garnet, diopside 254 (clinopyroxene) and enstatite (orthopyroxene) for each xenolith. This barometer and thermometer pair both also 255 depend upon the temperature and pressure, respectively. These two equations are therefore solved simultaneously 256 by iteration to obtain equilibration P-T conditions. Samples are discarded if they fail more than one of the eight 257 oxide, cation and equilibration checks. 41 Analyses from locations on thicker lithosphere are predominantly obtained 258 from heavy mineral concentrates generated during diamond exploration (plus rare diamond inclusions and occa-259 sional whole peridotite xenoliths), where the association of one mineral grain with any other has been lost. Thus, 260 the approach outlined above using multiple phases is unavailable, and we instead turn to single grain combined 261 thermobarometers for deriving equilibration P-T conditions. For these samples, we use the chrome-in-diopside 262 barometer 42 that exploits the exchange of chromium between clinopyroxene and garnet (Equation (9) , 2000). It uses only diopside compositions, but requires that garnet was also present in the source region.

264
The associated thermometer 42 exploits enstatite-in-diopside, again using only diopside compositions but requiring 265 that orthopyroxene was present within the source. The temperature is given by Equation (17) 2000). Again, these two equations must be solved by iteration to obtain P-T conditions for each diopside grain.

267
Calibration on laboratory experiments has shown that this thermobarometer may become innacurate at low pres-268 sures and at temperatures <700 • C. 41 We therefore only use P-T estimates derived from this thermobarometer that 269 yield depths >60 km and pass both of the clinopyroxene cation and oxide checks.

270
Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Fitting a geotherm to P-T estimates. For each locality, P-T estimates derived from thermobarometry and T S (z) must be independently determined by inverting real-Earth observational constraints on temperature, where V o ij are observed shear-wave velocities with associated standard deviation σ ij , V c ij is the prediction from 289 Equation (1), M is the number of age bins at a given depth and N is the number of depth slices. A second 290 suite of tie-points is created by assuming that temperatures are isentropic at depths well below the upper thermal 291 boundary layer. We calculate average V S as a function of depth over oceanic regions in the global model, and 292 over the whole spatial domain in regional models. Over the depth range 250-400 km, beyond which the resolving 293 power of surface waves drops significantly, these values are combined with an isentrope calculated for pyrolite with 294 a potential temperature of 1334 • C using Perple X. Misfit for the isentrope, H 2 , is It has been observed that over the depth range 150-400 km, both V S and Q −1 S are relatively consistent for oceanic 297 ages ≥ 100 Ma. Over this range, we stack the QRFSI12 attenuation model 50 , generating a suite of Q −1 S to V S 298 tie-points as a function of depth. Equations (1) and (15) are coupled such that average temperature is obtained Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 observed and predicted attenuation is We also adopt the bulk viscosity measurement 51 of η ref = 3 × 10 20 Pa s for upper mantle (∼ 100-670 km) and 303 compare it to the mean predicted value for 225-400 km depths obtained from Equation (7). Misfit, H 4 , is calculated where η c i is predicted viscosity. Finally, for calibration of regional tomography models, we take the better con- where M is the number of paleogeotherms, N is the number of tie-points associated with each geotherm and σ ij 319 reflects uncertainty in the V S measurement, assumed to be a constant 0.1 km s −1 which captures typical variations 320 between different tomography models at a given location. Combined misfit, H, is given by where w represents weighting applied to each misfit constraint. H is minimised in two steps. Initially, a pa-323 rameter sweep is performed to identify the approximate location of the global minimum. µ 0 U is varied between 333 For the regional model FR12 we constrain the calibration using the paleogeotherms. All weights are set to zero 334 except for w 2 = 1 and w 5 = 10, yielding minimum misfit H = 0.578 when µ 0 Mapping the lithosphere-asthenosphere boundary. A recent study 18 on the thermal structure of oceanic 337 lithosphere found that the 1175 ± 50 • C isotherm provides a good match to seismological observations of the 338 lithosphere-asthenosphere boundary (LAB), such as peak variation in the orientation of azimuthal anisotropy. In 339 this study, we therefore adopt this isotherm as a proxy for lithospheric thickness beneath the continents. T (z) is   The SL2013sv model contains only 6 seismometers in Australia, so has limited resolution within this continent.

575
Therefore, we also investigate three regional seismic tomography models to generate high resolution maps for the to the computationally intensive methodology, ∼ 3, 000 waveforms are used in this inversion.

590
The third and final regional model considered in this study is the radially anisotropic Y14. 58 It combines 591 Rayleigh waves (8000 fundamental, ∼ 2500 higher mode) and Love waves (approximately two-thirds as many) with 592 periods ∼ 25-200s, corrected for local crustal structure using a fixed crustal model. It adopts the same three-step 593 inversion procedure as YK04. 54 All three models are plotted alongside the global SL2013sv model in Figures S1, S2 Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S1: 100 km depth slice through Australian seismic tomography models. Black/green crosses = paleogeotherms used as constraints/tests in anelasticity calibration. (a) FR12 = regional isotropic V S 16 . (b) AuSREM = regional V SV 53 . (c) Y14 = regional Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S2: 175 km depth slice through Australian seismic tomography models. Black/green crosses = paleogeotherms used as constraints/tests in anelasticity calibration. (a) FR12 = regional isotropic V S 16 . (b) AuSREM = regional V SV 53 . (c) Y14 = regional Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S3: 250 km depth slice through Australian seismic tomography models. Black/green crosses = paleogeotherms used as constraints/tests in anelasticity calibration. (a) FR12 = regional isotropic V S 16 . (b) AuSREM = regional V SV 53 . (c) Y14 = regional Temperature estimates across a range of depths are required to generate a series of V S -T-P tie points in order 597 to calibrate the regional seismic tomography models. We therefore assemble a suite of Australian paleogeotherms 598 derived from thermobarometric analysis of mantle xenoliths and xenocrysts from fifteen locations in thick and thin 599 lithosphere ( Figure S4). The resulting P-T estimates are entered into FITPLOT to generate the palaeogeotherms 600 shown in Figure S5 (Methods).

601
The results of regional calibration using the paleogeotherms are shown in Figures S6 and S7. Note that the 602 global model SL2013sv yields good fits to paleogeotherms away from south Australia (Monk Hill, Orroroo and 603 Cleve), despite being lower resolution than the local models and being calibrated completely independently of this 604 information (red lines in Figure S7). Conversely, regional models often provide a poorer fit to the full range of 605 the paleogeotherms and can exhibit substantial crustal bleeding artefacts at depths shallower than ∼ 125 km.

606
Generally amongst the regional models, FR12 performs the best, followed by AuSREM and then Y14.
607 Figure S4: Location of Australian xenolith and xenocryst suites. Labels give site name and age (in million years); black crosses = locations used to constrain anelasticity calibration, green crosses = locations used to visually test validity of results; red/blue colours = lithospheric thickness (from Figure 1b), derived from FR12 seismic tomography model. 16 Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S5: Australian paleogeotherms derived from xenolith and xenocryst thermobarometry. Labels give site name and age (in million years) from Figure S4; red circles = P-T estimates derived from multiphase thermobarometry 40,39 ; blue circles = P-T estimates derived from single chrome diopside thermobarometry 42 ; dashed line = crustal thickness from AusMoho 47 ; solid line = FITPLOT optimal paleogeotherm.
Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S6: V S as a function of depth at sites of fifteen Australian paleogeotherms. Labels give site name (locations in Figure S4); red = global SL2013sv model 17 ; purple = regional FR12 model 16 ; blue = regional AuSREM model 53 ; orange = regional Y14 model 58 .
Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0 Figure S7: Calibration of anelasticity parameterisation on Australian paleogeotherms. Labels give site name and inferred age of paleogeotherms in million years (locations in Figure S4); sites Argyle to Wandagee are used to constrain calibration; sites Bow Hill to Sapphire Hill are used to visually check output; dashed line = crustal thickness from AusMoho 47 ; solid line = optimal FITPLOT geotherm from Figure S5; purple = regional FR12 model 16 ; blue = regional AuSREM model 53 ; orange = regional Y14 model 58 ; red = global SL2013sv model 17 , for comparison, calibrated independently of palaeogeotherm constraints.
Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0

608
For each of the individually calibrated seismic tomography models in this study, we have mapped out the LAB in 609 a consistent manner. The resulting maps for Australia are shown in Figure S8, whilst in Figure S9 we compare our 610 preferred FR12 regional model to previously published maps of LAB depth beneath Australia. Manuscript submitted to Nature Geoscience on 1st April 2019. c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license: https://creativecommons.org/licenses/by-nc-nd/4.0

622
In order to test the statistical significance of real deposit locations, test suites of random points on a sphere are 623 generated. Example test suites of 100, 1000 and 10,000 points are shown in Figure S12. For each Kolmogorov-Smirnov test, a number of random points are generated that is equivalent to the number 625 of real deposits of that type (109 for PbZn-CD, 147 for PbZn-MVT and 139 for sedimentary copper). Given the 626 low sample size for some of the deposit classes, the distribution of this random set can vary somewhat from the true 627 average distribution of continental locations. We therefore draw a test set in this manner 100 times ( Figure S13).

628
These random CDFs are relatively consistent but have some outliers. The D-value and Kolmogorov-Smirnov 629 statistics between each random CDF and the real one is calculated and reported within a histogram ( Figure S14).