Global distribution of sediment-hosted metals controlled by craton edge stability

Sustainable development and the transition to a clean-energy economy drives ever-increasing demand for base metals, substantially outstripping the discovery rate of new deposits and necessitating dramatic improvements in exploration success. Rifting of the continents has formed widespread sedimentary basins, some of which contain large quantities of copper, lead and zinc. Despite over a century of research, the geological structure responsible for the spatial distribution of such fertile regions remains enigmatic. Here, we use statistical tests to compare deposit locations with new maps of lithospheric thickness, which outline the base of tectonic plates. We find that 85% of sediment-hosted base metals, including all giant deposits (>10 megatonnes of metal), occur within 200 kilometres of the transition between thick and thin lithosphere. Rifting in this setting produces greater subsidence and lower basal heat flow, enlarging the depth extent of hydrothermal circulation available for forming giant deposits. Given that mineralization ages span the past two billion years, this observation implies long-term lithospheric edge stability and a genetic link between deep Earth processes and near-surface hydrothermal mineral systems. This discovery provides an unprecedented global framework for identifying fertile regions for targeted mineral exploration, reducing the search space for new deposits by two-thirds on this lithospheric thickness criterion alone. Major sediment-hosted base metal deposits are located within 200 km of the border between thick and thin lithosphere, according to statistical comparisons between global lithospheric thickness and known deposit locations.

Sustainable development and the transition to a clean-energy economy drives ever-increasing demand for base metals, substantially outstripping the discovery rate of new deposits and necessitating dramatic improvements in exploration success. Rifting of the continents has formed widespread sedimentary basins, some of which contain large quantities of copper, lead and zinc. Despite over a century of research, the geological structure responsible for the spatial distribution of such fertile regions remains enigmatic. Here, we use statistical tests to compare deposit locations with new maps of lithospheric thickness, which outline the base of tectonic plates. We find that 85% of sediment-hosted base metals, including all giant deposits (>10 megatonnes of metal), occur within 200 km of the transition between thick and thin lithosphere. Rifting in this setting produces greater subsidence and lower basal heat flow, enlarging the depth extent of hydrothermal circulation available for forming giant deposits. Given that mineralisation ages span the last 2 billion years, this observation implies long-term lithospheric edge stability and a genetic link between deep Earth processes and near-surface hydrothermal mineral systems. This discovery provides an unprecedented global framework for identifying fertile regions for targeted mineral exploration, reducing the search-space for new deposits by two-thirds on this lithospheric thickness criterion alone.
Consumption of base metals (copper, lead, zinc and nickel) over the next ∼25 years is set to exceed the total produced in human history to date (Ali et al., 2017;Schodde, 2017). Moreover, critical minerals (e.g. cobalt, indium and germanium) are often produced as by-products of base metal mining and are essential in many high-tech applications (Nassar et al., 2015;Mudd et al., 2018;IRENA, 2019;Dominish et al., 2019;Sovacool et al., 2020). A growing concern is that the rate of exploitation of existing reserves is outstripping discovery of new deposits, despite exploration expenditure tripling during the 2005-2012 minerals boom (Ali et al., 2017;Schodde, 2017). To reverse this trend and supply the resources necessary to comply with policies such as the Paris Climate Agreement and United Nations' Sustainable Development Goals, improved techniques for locating new deposits are required, particularly those buried under shallow sedimentary cover or ice (The Uncover Group, 2012).

Narrowing the search-space for new deposits
In mineral exploration, initial area selection at continental scales is arguably the most important step, as successful identification of fertile regions can compensate for many subsequent errors (McCuaig et al., 2010). Over the last two decades, the search for analogues of known deposits has progressed towards a more holistic determination of factors controlling deposit generation and preservation (Wyborn et al., 1994;Bierlein et al., 2006;McCuaig & Hronsky, 2014;Dentith et al., 2018;Skirrow et al., 2019). Mineral systems analysis has resulted in a growing acceptance that the spatial distribution of deposits associated with magmatic processes is controlled by lithospheric-scale structure (McCuaig et al., 2010;Begg et al., 2010;Griffin et al., 2013). For example, porphyry copper deposits are generated by wet melting in the mantle wedge above a subducting slab, followed by emplacement of these melts into the shallow overlying crust and subsequent concentration by high-temperature hydrothermal circulation (Griffin et al., 2013). Thus, by combining the plate tectonic setting with geological constraints on the location of key mineral system ingredients, the search-space for new magmatic deposits can be substantially reduced (Rosenbaum et al., 2005;Butterworth et al., 2016;O'Reilly et al., 2017).
In the case of sediment-hosted deposits, most assessments to date have focused on their genesis within the context of Earth's secular evolution, as well as past tectonic and geographic settings (McCuaig et al., 2018). The majority are found in failed rift and passive margin settings, and it is generally agreed that basin-scale hydrothermal circulation is required to scavenge sufficient metals to form giant deposits ( Figure 1a; Leach et al., 2010;Hitzman et al., 2010;Manning & Emsbo, 2018). Metals are mobilised and transported by oxidised brines with moderate temperatures (80-250 • C) and moderate-to-high salinity (10-30 wt.% NaCl), limiting their maximum age to the Great Oxidation Event at 2.4 Ga (Leach et al., 2010;Hitzman et al., 2010). These fluids are sourced from evaporites at low latitudes and remain buffered as they pass through voluminous oxidised terrestrial sediments, allowing them to scavenge lead from arkosic sandstones and felsic volcanics, as well as copper and zinc from mafic rocks (Leach et al., 2010;Hitzman et al., 2010). Transport along faults, either during rifting or basin inversion, focuses these fluids into oxidation-reduction interfaces, such as distal-facies black shales, where metals precipitate (Figures 1b and 1c; Huston et al., 2016). Sediment-hosted base metal deposits are desirable due to their greater quantity of contained metal (in contrast to volcanogenic massive sulfides) and higher grades of ore (in comparison to porphyry copper), resulting in lower environmental degradation during extraction (Dominish et al., 2019;Azadi et al., 2020). However, narrowing the search-space for new sediment-hosted deposits has been less successful than for magmatic mineral systems. Sedimentary basins cover ∼75% of the continental surface, and the key ingredients of evaporites associated with brine formation, felsic and mafic volcanic rocks for sourcing metals, and organic rich shale precipitation sites, are widespread and do not substantially reduce this search-space. The first-order geological control that localises their spatial distribution throughout the continents remains unknown, severely limiting predictive power for identifying new targets. A classic example comes from the Carpentaria Zinc Belt in northern Australia, which contains several world class clastic-dominated lead-zinc deposits formed between 1.8-1.4 Ga (Figure 2a). These deposits lie along an arcuate trend that runs oblique to mapped geology and crustal geological boundaries, as demonstrated by gravity and magnetic datasets (Geoscience Australia, 2018).
Crucially, despite the absence of a clear crustal relationship, the linear distribution of sediment-hosted deposits in the Carpentaria Zinc Belt hints at an underlying regionalscale control. An important advance in understanding the genesis of magmatic mineral systems has come from probing their relationship with major crustal and lithospheric structures (Griffin et al., 2013). Given that sedimentary basins are themselves the result of lithospheric scale processes, we therefore investigate both regional and global-scale links between sediment-hosted base-metal deposits and the most fundamental shallow mantle structure -the lithosphere-asthenosphere boundary (LAB).

Relationship with lithospheric structure
We begin by collating global inventories of six major basemetal mineral systems from published sources (Methods). Three are magmatic and three are sediment-hosted, which include sedimentary copper (Cu-sed), clastic-dominated leadzinc (PbZn-CD, commonly also referred to as sedimentary exhalative), and Mississippi Valley-type lead-zinc (Pb-Zn-MVT). We next refine a method developed by Priestley & McKenzie (2013) for mapping the thermal LAB from seismic tomography, taking into consideration recent laboratory experiments concerning the effect of anelasticity on shear-wave velocities (Yamauchi & Takei, 2016;Methods). This benchmarking procedure is necessary in order to increase consistency between LAB maps obtained for different tomography models, which can image surprisingly variable seismic velocities. A high resolution regional LAB map over Australia is obtained from the FR12 model (Fishwick & Rawlinson, 2012) and is calibrated using nine local paleogeotherms derived from thermobarometry of mantle peridotite xenoliths and xenocrysts. To expand our analysis to other continents, a global LAB is also produced using the SL2013sv model (Schaeffer & Lebedev, 2013) and calibrated using multiple constraints, including the latest thermal structure of cooling oceanic lithosphere (Richards et al., 2018). This global LAB exhibits a bi-modal thickness distribution, with peaks at 90 km and 190 km, separated by a minimum at 150 km (Supplementary Information).
Inspection of the Australian model reveals a striking correlation between major sediment-hosted mineral deposits and the edge of thick lithosphere, defined here by the 170 km thickness contour (Figure 2b). Major PbZn-CD and sedimentary copper deposits in the Carpentaria Zinc Belt overlie this contour, which runs obliquely to geological boundaries, such that intersections between these two features con-2 Preprint accepted at Nature Geoscience on 15 th May 2020 and made available under the CC-BY-NC-ND 4.0 license. c 2020 Figure 2: Distribution of sediment-hosted and iron-oxide-copper-gold base metal deposits as a function of Australian lithospheric thickness. (a) Carpentaria Zinc Belt; red/blue = variably reduced to pole aeromagnetic intensity data (Geoscience Australia, 2018); grey polygons = generalised outcrop of Cretaceous marine sediments in Eromanga and Carpentaria Basins (Raymond, 2018); black dashed contour = 170 km LAB thickness; symbols = deposit locations; symbol area proportional to estimate of total contained mass of metal (either lead plus zinc, or copper) in mega-tonnes (Mt); unknown deposit size given 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley-type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed); stars = iron-oxide-copper-gold (IOCG). (b) Australian LAB mapped by converting FR12 tomography (Fishwick & Rawlinson, 2012) to temperature using an anelasticity parameterisation (Yamauchi & Takei, 2016) calibrated on local paleogeotherms (Methods) and illuminated by free-air gravity anomalies (Geoscience Australia, 2018); black/green crosses = geotherms used as constraints/tests in anelasticity calibration (Supplementary Information); box = location of panel (a); large star in South Australia shows location of the Olympic Dam IOCG deposit.
sistently coincide with deposit locations. This behaviour is particularly useful for highlighting new prospective regions for exploration. Other observables that correlate with this lithospheric thickness change include variations in lead isotopes from Proterozoic galena and pyrite minerals (Huston et al., 2019), long-wavelength gravity anomaly gradients (Hobbs et al., 2000), changes in crustal character in deep reflection seismic profiles (Kennett et al., 2016), a topographic ridge, and the western extent of Cretaceous marine sediments ( Figure 2a). These latter two associations demonstrate the post-Proterozoic stability of this edge and its influence on local geology and topography. There is also a strong relationship with iron-oxide-copper-gold deposits, including the Olympic Dam mine in South Australia (84 Mt of copper, world's largest known uranium resource; Heinson et al., 2018;Skirrow et al., 2018;Curtis & Thiel, 2019). However, a lack of consensus on global classification schemes means that we have limited analysis of this deposit type to only Australia. Extending our analysis globally further confirms the strength of this relationship ( Figure 3). The link between the 170 km lithospheric thickness contour and location of all large sediment-hosted deposits holds regardless of deposit age, which spans the last 2 billion years. Given the 180-220 km cluster of LAB thicknesses is likely to represent standard cratonic lithosphere, the 170 km contour therefore demarcates the outer boundary of cratons. Within the PbZn-CD deposit class, those more strongly associated with abundant mafic rocks systematically occur on the thinner lithosphere side of the contour compared to their carbonaterich counterparts (e.g. Carpenteria Zinc Belt and northwest North America). These observations are consistent with an extensional origin of the host basins.
To quantify these visual relationships, the shortest distance is calculated between each deposit and the 170 km LAB thickness contour and results are plotted in a cumulative distribution function (CDF). Weighting deposits by the mass of contained metal and substituting the Australian LAB from the global model with our high-resolution regional version substantially improves the correlation for PbZn-CD ( Figure 4a). Globally, we observe that ∼95% of sedimentary copper, ∼90% of clastic-dominated lead-zinc and ∼70% of Mississippi Valley-type lead-zinc resources are located within a 200 km-wide corridor either side of the 170 km LAB thickness contour (Figure 4b). This region corresponds to only ∼35% of continental surface area. Given that this swath width is similar to the ∼280 km node spacing in SL2013sv, tighter constraints are only possible with higher resolution tomography models. For example, all giant deposits in Australia lie within 100 km of the 170 km contour for the higher resolution FR12 model. The significance of the relationship is globally examined using the two-sample Kolmogorov-Smirnov test, which estimates that the probability of these sedimenthosted deposits representing random continental locations is less than 1 in 10 12 (Kolmogorov, 1933;Methods). This relationship holds true for three other global surface wave tomography models (Extended Data Figure 1).
All >10 mega-tonne sediment-hosted deposits are located along the boundary of thick lithosphere, but amongst the smaller deposits, there are some notable exceptions. Minor PbZn-CD outliers occur in Europe, the Caribbean, Indonesia and east China. Anomalous PbZn-MVT deposits are found in Ireland, east China and along the Tethys sub-3 Preprint accepted at Nature Geoscience on 15 th May 2020 and made available under the CC-BY-NC-ND 4.0 license. c 2020 Figure 3: Global distribution of sediment-hosted base metal deposits as a function of lithospheric thickness. LAB derived from SL2013sv tomography model (Schaeffer & Lebedev, 2013) using a calibrated anelasticity parameterisation (Yamauchi & Takei, 2016;Methods). Symbols = deposit locations; symbol area proportional to estimate of total contained mass of metal (either lead plus zinc, or copper) in mega-tonnes (Mt); unknown deposit size given 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed).
duction zone across Europe, whilst small sedimentary copper deposits occur in southwestern North America and southern South America. This observation indicates that minor sediment-hosted mineral systems can develop in a variety of basins, whilst giant deposits form only at the edges of cratonic lithosphere. PbZn-MVT deposits are generally more widely distributed, which likely reflects that a subset of these deposits are linked with orogenic processes along basin margins (Leach et al., 2010), resulting in longer hydrothermal fluid flow paths. Nevertheless, not all sediment-hosted outliers were necessarily anomalies at the time of ore formation. The majority now occur in accretionary terranes, whereby plate tectonic processes may have rifted segments off thick lithosphere and transported them into subduction zone settings. Other areas, such as east China, are known to have undergone lithospheric thinning some time after deposit formation, based on thermobarometric constraints (Menzies et al., 2007).
Regardless of age, sediment-hosted base-metal deposits predominantly cluster on the edges of present-day thick lithosphere. Therefore, many of these lithospheric steps appear to be remarkably robust on billion-year timescales, despite the assembly and disaggregation of several supercontinents, the impacts of large igneous provinces and the possible erosional effect of edge-driven convection (Currie & van Wijk, 2016). Deposit ages in northwestern North America span ∼1.5-0.5 Ga, pointing to the stability and importance of this boundary in localising multiple deformational and oreforming processes.

Mineral system implications
Our results indicate that the edges of thick lithosphere place first-order controls on the genesis of extensional basins and their associated mineral systems (Figure 1). Rifting causes localised thinning and produces a lateral transition from oxidising terrestrial environments into marine settings. This transition provides the optimal juxtaposition of the ingredients necessary for deposit formation. The adjacent unstretched cratons provide a bountiful source of oxidised sediments and extensive low-elevation platforms, which enhance evaporite formation. Proximal land masses also aid the development of restricted marine settings that are ideal for accumulating thick evaporite sequences (through periodic cycles of evaporation and flooding), and promote euxinic water conditions that are favourable for deposition of reducing shales that have high organic carbon content. Thinning of the lithosphere in the centre of the basin causes decompression melting, providing mafic and felsic volcanic rocks from which metals are scavenged. Intercalation of proximal and distal facies components is further modulated by transient vertical motions, generally thought to be associated with edge driven convection across lithospheric steps (Davies & Rawlinson, 2014). Nevertheless, these mineral system components are common to both rifts in thick lithosphere and regular passive margins, and the question remains -what is favourable about rifting cratonic lithosphere for formation of the shallow hydrothermal systems necessary to produce giant deposits?
From a geodynamic perspective, these lithospheric edges represent rheological contrasts that focus strain and localise repeated cycles of extensional deformation and basin contraction, thereby controlling both the spatial distribution of required lithologies and the focusing of mineralising fluids (Sloan et al., 2011;Gibson et al., 2016;O'Reilly et al., 2017). Thick cratonic lithosphere is colder than standard continental lithosphere and has a larger seismogenic thickness, resulting in the development of deeper, longer, and more widely spaced normal faults during rifting (Biggs et al., 2010). This architecture increases the horizontal aspect ratio of hydrothermal cells, providing greater sediment volumes for fluid-rock 4 Preprint accepted at Nature Geoscience on 15 th May 2020 and made available under the CC-BY-NC-ND 4.0 license. c 2020 interaction. These faults are active for longer periods of time, and the entire syn-rift phase of basin formation can last 50-100 Myr, in contrast to standard continental rifts that typically last ∼25 Myr, yielding a more extensive time window for mineralisation (Allen & Armitage, 2012). Furthermore, basins forming in continental interiors are more likely to be restricted, which is not only favourable for high siliciclastic supply and the formation of evaporites and reductants, but also results in laterally continuous sediment packages that trap hydrothermal fluid systems and potentially promote recirculation of brines (Hitzman et al., 2010).
A key observation is that metal precipitation in sedimenthosted base metal deposits is generally driven by oxidationreduction reactions, which become ineffective when brine temperatures exceed ∼ 250 • C (Figures 1b and 1c; Huston et al., 2016). As hydrothermal fluid temperatures are buffered by conditions towards the base of the sediment pile (often where the mafic metal source rocks are located), the basal temperature of the sedimentary pile must not substan-tially exceed this threshold value. Total extension in a basin can be estimated using a stretching factor, β, which is the ratio of original to final crustal thickness. Failed rifts on standard continental lithosphere such as the North Sea typically have β ≈ 2, and simple thermal modelling assuming pureshear rifting indicates that this produces ∼ 4 km of syn-rift sediment with basal temperatures cooler than ∼ 250 • C (Figure 5a;Methods). Given that all the necessary ingredients occur within basins, the likelihood of developing a successful mineral system is higher for a larger sediment pile, which can be achieved by increasing the stretching factor. However, more extreme rifting causes the asthenosphere to upwell to substantially shallower depths, producing elevated basal heat flow that heats the sediment pile above this threshold and so inhibits metal precipitation (Figure 5b). Two important differences occur during rifting of cratonic lithosphere. First, the greater initial thickness of the lithosphere results in a lower geothermal gradient, such that the basal heat flow spike during rifting is substantially lower than for standard continental rifts (Supplementary Information). Secondly, the density of cratonic lithosphere is reduced by up to ∼50 kg m −3 by chemical depletion compared to standard lithosphere (Jordan, 1978). This increased buoyancy reduces the dampening effect during syn-rift subsidence that is associated with replacing cold continental lithosphere with lower density asthenosphere, resulting in substantially larger thicknesses of syn-rift sediments for any given stretching factor. Assuming unlimited sediment supply, β ≈ 2 yields a 8-9 km syn-rift sediment pile, the base of which stays cooler than the threshold ∼250 • C ( Figure 5c). Thus, rifting cratonic lithosphere produces more than twice the volume of mineral system ingredients without exceeding the thermal conditions necessary for successful precipitation, over a duration of time that can be up to a factor of four more extensive (Figure 5d and 5e). This mechanism explains why smaller deposits can occur in any extensional setting (e.g. Irish PbZn-MVT deposits), but giant deposits requiring the largest volumes of fluid-rock interaction are restricted to rift basins at the edges of the thickest lithosphere.
A final consideration is that a setting on the edge of thick lithosphere enhances the preservation potential of deposits through subsequent orogenic events and supercontinent cycles. For example, the Boleo Copper District in Baja California formed only three-million years ago and sits in shallow crust on thin lithosphere, resulting in poor long-term preservation potential. In contrast, the 1.7 Ga Broken Hill deposit in Australia (world's largest known lead deposit) has been metamorphosed to amphibolite-granulite facies, yet survives on the edge of the Curnamona part of the South Australian Craton.
Surprisingly, given the results of previous studies (e.g. Griffin et al., 2013), deposits associated with magmatic systems generally exhibit a weaker association with the edge of cratonic lithosphere than sediment-hosted systems (Supplementary Information). Porphyry copper deposits are predominantly Cenozoic in age and are generally positioned on thin lithosphere (≤100 km). Their formation in subduction zone settings at shallow crustal depths leads to poor preservation potential within the geological record, making this association unsurprising. Volcanogenic massive sulphides have an episodic age distribution from 3.5 Ga to present. Their generation is thought to require moderate-degree partial melting of hydrated mantle in back-arc settings (Huston et al., 2010). We observe that they spatially occur randomly on thick and thin lithosphere, but exhibit systematic temporal ordering, with the oldest positioned over thick lithosphere 5 Preprint accepted at Nature Geoscience on 15 th May 2020 and made available under the CC-BY-NC-ND 4.0 license. c 2020 that are rimmed by progressively younger deposits, consistent with growth of cratons by accretion. Finally, magmatic nickel deposits are mostly Archean and Proterozoic in age and commonly occur on thick lithosphere (≥150 km). Unlike other base metal deposits, their distribution is associated with edges of even thicker lithosphere (∼200 km), broadly consistent with previous studies showing major lithospheric structural controls on these deposit locations (Begg et al., 2010;Regis et al., 2017;Alghamdi et al., 2018). Their generation requires large fraction partial melting of peridotite, indicative of high mantle temperatures (more prevalent in a early, hotter Earth) and decompression melting at shallow depths (Arndt et al., 2005). Therefore, their present distribution suggests lithospheric thickness must have locally increased since formation, simultaneously enhancing preservation potential.
In summary, this work illustrates a new and robust link between giant sediment-hosted base metal mineral systems and the edges of thick lithosphere. Approximately 55% of the world's lead, 45% of its zinc and 20% of known copper is found within ∼ 200 km of this boundary. We have demonstrated the value of regional seismic arrays to better resolve this edge and enhance the mineral exploration efforts required to sustain ongoing global development. Importantly, deposit ages indicate that, following rifting, the edges of thicklithosphere are generally stable over billion-year timescales. The far-reaching geodynamic and societal implications of these observations highlight the need for extensive further research. To improve resolution of mapped lithospheric structure, higher fidelity seismic imaging needs to be coupled with enhanced mantle xenolith coverage and tighter constraints on seismic anelasticity from mineral physics experiments. More generally, these maps should be integrated with models of basin dynamics, surface processes and reactive transport modelling, and benchmarked against additional geological information, such as sedimentary facies variations, tectonic structures and alteration zones. These multiple research strands will yield fundamental new insights into sedimenthosted mineral systems and lead to substantial improvements in exploration success.
For each deposit, we include the type (based on established classification schemes), location, age (direct measurement or inferred based on geological relationships) and total resource size by combining historical production with estimated resources. Our Cu-sed deposit dataset follows the classification scheme and compilation of Hitzman et al. (2005), cross-checked against Cox et al. (2007). Where these two compilations disagree on deposit size, the larger value has been used. Our PbZn-CD and PbZn-MVT deposit compilations extensively revise and build on the work of Taylor et al. (2009). References for each deposit type were manually checked and additional references have been included. We exploit the compilation of Sillitoe (2010) for Cu-por deposits. Our magmatic Ni-Cu-PGE compilation follows Hoatson et al. (2006), with deposit location populated from disparate sources. Our catalogue of VMS deposits is an extensive revision of the compilation by Franklin et al. (2005). Australian information for all the above deposit types, with the addition of 25 iron-oxide-copper-gold deposits, was updated using the authors' own knowledge building from the Geoscience Australia OZMin database (Sexton, 2011).
We have endeavoured to assemble the most complete deposit dataset possible by revising and extending pre-existing compilations. Our database can be found in the online Supplementary Datasets. Importantly, patchy or absent reporting of mineral deposit information from some countries inevitably means our global database is incomplete, but we do not believe that this will impact the veracity of our main conclusions.
Metal scavenging window. Following Cooke et al. (2000), we have selected a precipitation gradient from ∼100 ppm to 1 ppm for the metal scavenging window in Figure 1. Estimates of the metal concentration of mineralising brines are obtained from fluid inclusion studies contained within ore and other gangue minerals. Care must be taken that host minerals formed during the period of ore deposition and that analyses are not contaminated by metals contained within the mineral lattice itself, which is a particular problem in the early stages of LA-ICP-MS analysis when the laser initially vaporises a mixture of mineral and fluid inclusion to produce a mixed signal. Generally, inclusions from European sedimenthosted deposits and Irish and US MVTs hosted by quartz, calcite, and dolomite have values at the lower end of our range, whilst those in sphalerite can sometimes reach several hundred ppm of Pb (Appold et al., 2004;Kostova et al., 2004;Wilkinson et al., 2005;Kotzeva et al., 2011;Fusswinkel et al., 2014). Values in the several thousands of ppm have been inferred in rare cases for single fluid inclusions (e.g. Wilkinson et al., 2009;Davey, 2019) However, these high values are often found for inclusions with high homogenisation temperatures, whilst those fluids thought to be cooler than ∼ 250 • C generally return concentrations an order of magnitude or more below 1000 ppm (Davey, 2019). While there remains debate on complexity and exact composition of mineralising fluids (Fusswinkel et al., 2013;Schlegel et al., 2018), for ores deposited via the reduction mechanism, our inference of an enlarged operating window in cratonic rift basins holds regardless of the exact fluid metal concentrations.
Choice of seismic tomography model. Our LAB maps are based on recent shear-wave velocity (V S ) models that contain a lot of surface wave data and have nominal vertical resolution on the order of 25-50 km (Priestley & McKenzie, 2006). For the global map, we use SL2013sv, which is an upper mantle-only model built from a combination of body and surface waves, including fundamental and higher modes (Schaeffer & Lebedev, 2013). Periods considered are 11-450 s, ∼750,000 seismograms are included, and misfits are calculated between synthetics and the full waveform up to the 18 th overtone. Crucially, simultaneous inversion for the crustal model results in reduced smearing of slow crustal velocities down into the upper mantle in comparison to other models, thereby allowing us to use more depth slices in our V S to temperature calibration. Checkerboard resolution tests indicate that features ∼600 km in diameter at lithospheric depths are generally well resolved. Finer features should be resolvable in regions with dense ray path coverage, such as North America, Europe and southeast Asia. The SL2013sv model utilises data from only 6 seismometers in Australia, so has limited resolution within this continent. Therefore, we also investigate the FR12 regional seismic tomography model to generate a high resolution map for the Australian continent (Fishwick & Rawlinson, 2012). FR12 is a radially isotropic V S model derived from Rayleigh wave travel times (Fishwick et al., 2008). Periods considered are 50-120 s and the fundamental and first four higher modes have been used where possible, leading to good sensitivity down to ∼250 km depths. It contains a greater number of source-receiver paths (>13,000) compared to other Australian models. However, it uses an a priori crustal model that remains fixed throughout the inversion, resulting in noticeable smearing of crustal velocities into the upper mantle. Checkerboard tests indicate that features ∼300 km in diameter at lithospheric depths are well resolved, and where higher mode information is included.
For comparison, we include LAB maps derived from three additional upper mantle seismic tomography datasets with global coverage in Extended Data Figure 1. These include the 3D2015-07Sv model of Debayle et al. (2016), the CAM2016 model of Ho et al. (2016) and Priestley et al. (2018), and a version of SL2013sv into which we have blended the regional updates SL2013NA in North America (Schaeffer & Lebedev, 2014), AF2019 in Africa (Celli et al., 2020a), and SA2019 in South America and the South Atlantic Ocean (Celli et al., 2020b) to produce a combined model SLNAAFSA. For the continental analyses in Australia, we also consider regional model AuSREM of Kennett et al. (2013) and Y14 of Yoshizawa (2014). Deposit locations are compared with all seven of our new LAB maps in the Supplementary Information, in addition to eleven previously published LAB maps derived from a mixture of heat flow data, seismic tomography datasets, and potential field data. Many giant sediment hosted mineral deposits lie along LAB edges defined by these other studies, testifying to the robustness of the observed relationship.
Parameterising anelasticity. Seismic tomography models provide images of the upper mantle and have been extensively used to constrain its thermomechanical structure, composition, and the depth of the lithosphere-asthenosphere boundary (Priestley & McKenzie, 2006;An & Shi, 2006;Goes et al., 2012;Afonso et al., 2016;Cammarano & Guerri, 2017;Klöcking et al., 2018;Afonso et al., 2019). For accurate mapping from V S into temperature, it is essential to include the effect of anelasticity on this conversion (Karato, 1993;Cammarano et al., 2003). When a viscoelastic material such as the mantle is cold, deformation associated with passage of acoustic energy is predominantly elastic, yielding a linear dependence of V S on temperature referred to as the anharmonic velocity. As temperature increases, a special case of viscoelastic deformation known as anelasticity becomes increasingly important and gives rise to a strongly non-linear relationship between V S and temperature. This behaviour has been extensively studied in laboratory experiments on silicates and organic analogues of mantle rocks, revealing that the strength of the anelastic regime varies with both the frequency of seismic waves and as a function of material properties, such as melting temperature and grain size (Jackson et al., 2002;Sundberg & Cooper, 2010;McCarthy et al., 2011;Takei et al., 2014;Faul & Jackson, 2015). Several studies have attempted to parameterise these complex dependencies, and have been regularly updated as forced oscillation and creep experiments in the laboratory have been pushed towards increasingly realistic frequencies, pressures, temperatures, grain sizes and strain rates (Faul & Jackson, 2005;Jackson & Faul, 2010;Takei, 2017). In this study, we adopt the parameterisation of Yamauchi & Takei (2016), which includes effects of anelasticity in pre-melt conditions (temperatures above ∼ 90% of melting temperature) and is outlined in full in the Supplementary Information.
Xenolith and xenocryst thermobarometry. Temperature estimates across a range of depths are required to generate a series of V S -temperature-pressure tie points in order to calibrate the regional seismic tomography models. We therefore assemble a suite of fifteen Australian paleogeotherms derived from thermobarometric analysis of mantle xenoliths and xenocrysts (Supplementary Information). These come from a range of settings within both thick and thin lithosphere. Localities with thin lithosphere tend to have data obtained from whole xenolith samples, typically hosted in basaltic volcanic products. For these cases, the compositions of multiple phases (garnet, clinopyroxene, orthopyroxene and olivine) can be obtained that all equilibrated under the same pressure-temperature (P-T) conditions. In these samples, we use a thermometer that exploits exchange of calcium and magnesium between orthopyroxene and clinopyroxene (Taylor, 1998), and a barometer based upon aluminium exchange between orthopyroxene and garnet, given by Equation (5) of Nickel & Green (1985). This approach therefore requires compositions of garnet, diopside (clinopyroxene) and enstatite (orthopyroxene) for each xenolith, and we only use samples with all three of these minerals present. This barometer and thermometer pair both also depend upon the temperature and pressure, respectively. These two equations are therefore solved simultaneously by iteration to obtain equilibration P-T conditions. Samples are discarded if they fail more than one of the eight oxide, cation and equilibration checks suggested by Nimis & Grütter (2010).
Despite all samples containing garnet, a small number return depths as shallow as ∼25 km (see Bullenmerri, Monaro, Mt St Martin, and Sapphire Hill). The presence of garnet in xenoliths from shallow depths is well documented. The garnet-spinel transition can occur at pressures as low as 1 GPa (∼30 km depth) in pyroxenite and 1.5 GPa (∼45 km depth) in lherzolite, with the exact pressure of the transition depending on relative abundance of Cr and Al in each assemblage (Gasparik, 1984;Klemme, 2004;Nimis & Grütter, 2010). Our shallow samples are dominantly pyroxenites and mostly give pressures larger than the 1 GPa lower limit. Of these four sites with shallower samples, we select only Bullenmerri and Monaro for the anelasticity calibration, as these geotherms also contain samples at greater depths. In both cases, the deeper samples are consistent with the shallow results.
Analyses from locations on thicker lithosphere are predominantly obtained from heavy mineral concentrates generated during diamond exploration (plus rare diamond inclusions and occasional whole peridotite xenoliths), where the association of one mineral grain with any other has been lost. Thus, the approach outlined in the preceding paragraphs using multiple phases is unavailable, and we instead turn to single grain combined thermobarometers for deriving equilibration P-T conditions. For these samples, we use the chrome-in-diopside barometer that exploits the exchange of chromium between clinopyroxene and garnet (Equation (9) of Nimis & Taylor, 2000). It uses only diopside compositions, but requires that garnet was also present in the source region. The associated thermometer exploits enstatite-in-diopside, again using only diopside compositions but requiring that orthopyroxene was present within the source. The temperature is given by Equation (17) of Nimis & Taylor (2000). Again, these two equations must be solved by iteration to obtain P-T conditions for each diopside grain. Calibration on laboratory experiments has shown that this thermobarometer may become inaccurate at low pressures and at temperatures <700 • C (Nimis & Grütter, 2010). We therefore only use P-T estimates derived from this thermobarometer that yield depths >60 km and pass both of the clinopyroxene cation and oxide checks.
There are two sources of error to consider for each suite of P-T estimates. The first is uncertainty in the microprobe analyses of elemental oxide concentrations in each of the mineral samples. For the three-mineral thermobarometer, this introduces uncertainty of ±30 • C and ±10 km at low temperatures (∼700 • C), reducing to ±10 • C and ±3 km by ∼1200 • C (Mather et al., 2011). For the diopside-only thermobarometer, uncertainties are larger at ±70 • C and ±12 km for low temperatures (∼600 • C) and ±15 • C and ±3 km for higher temperatures (∼1200 • C) (Mather et al., 2011). However, these uncertainties in pressure and temperature are positively correlated, such that samples broadly move up and down the geothermal gradient, with limited effect on the best fitting geotherm. The second and more significant source of uncertainty arises from error in the thermobarometers themselves, which are calibrated on laboratory samples over a range of pressure-temperature conditions and do not necessarily trade-off in the same manner. Quoted uncertainties are ±50 • C and ±15 km for the three-mineral, and ±100 • C and ±15 km for the diopsideonly thermobarometer (Nimis & Taylor, 2000;Nimis & Grütter, 2010;Mather et al., 2011).
Fitting a geotherm to P-T estimates. For each locality, P-T estimates derived from thermobarometry are entered into FITPLOT to constrain the best-fitting paleogeotherm (McKenzie et al., 2005;Mather et al., 2011;Supplementary Information). Within the crust, we adopt a constant conductivity of 2.5 W m −1 • C −1 , whilst a pressure-and temperaturedependent parameterisation is used within the mantle (Osako et al., 2004). Bulk crustal radiogenic heat production is assumed to be 0.7 µW m −3 , with a standard deviation of 0.2 µW m −3 (Jaupart et al., 2007). Crustal thickness at each location is ob-tained from the AusMoho model with standard deviation assigned as 10% of the total thickness (Kennett et al., 2011). We assume a potential temperature of 1330 ± 50 • C, which is consistent with both seismological observations and the thickness and geochemistry of mid-ocean ridge basalts, assuming a dry lherzolite source using a corner-flow melting parameterisation (Dalton et al., 2014;Katz et al., 2003;Shorttle et al., 2014). Kinematic viscosity of the mantle is set to 2 × 10 16 m 2 s −1 , with a standard deviation of 0.7 orders of magnitude, which is consistent with constraints from glacial isostatic adjustment (Lau et al., 2016). Thermal parameters that are consistent with the melting parameterisation are used to calculate the adiabatic gradient, including a reference density of ρ 0 = 3.3 Mg m −3 , thermal expansivity of α = 3 × 10 −5 • C −1 and specific heat capacity of C P = 1187 J kg −1 • C −1 . Uncertainty in the crustal thickness, radiogenic heat production, mantle potential temperature, and kinematic viscosity are propagated through FIT-PLOT using a Monte Carlo approach. 1000 combinations of these four parameters are randomly drawn assuming Gaussian distributions of the uncertainties. Geotherms are strongly consistent in the vicinity of P-T constraints, but can vary by ±50 • C when greater than ∼30 km from a P-T estimate (Supplementary Information).
Calibrating V S -to-temperature conversion. Some of the anelasticity parameters have been directly constrained by forced oscillation experiments on borneol (Yamauchi & Takei, 2016). Others, however, are material properties that must be independently determined. A widely adopted approach is to fix these parameters for a given mineral assemblage, often calculated using mineral physics tables and a thermodynamic Gibbs energy minimisation algorithm (Stixrude & Lithgow-Bertelloni, 2005;Connolly, 2009;Stixrude & Lithgow-Bertelloni, 2011;Holland & Powell, 2011;Cottaar et al., 2014). In this manner, an anelastic conversion can be used in a forward sense to map between V S and temperature (An & Shi, 2006;Cammarano et al., 2009;Goes et al., 2012;Dannberg et al., 2017;Cammarano & Guerri, 2017). However, inferred temperature structures are variable as a result of uncertainty in the mantle's chemical composition and grain size, and differences in absolute V S between tomography models arising from different reference models and regularisation schemes.
An alternative approach to constraining these material properties is to invert real-Earth observations of the relationship between temperature, shear-wave velocity, attenuation and viscosity in the upper mantle (Priestley & McKenzie, 2006;Afonso et al., 2013a,b;Yamauchi & Takei, 2016). The general philosophy is that there are certain properties that are 'known' about Earth, including the typical thermal structure of oceanic lithosphere (Richards et al., 2018), the average adiabatic gradient within the convecting mantle (Connolly, 2009), the attenuation structure of the upper mantle beneath old oceanic lithosphere (Dalton et al., 2009), and the bulk diffusion creep viscosity of the upper mantle from studies of glacial isostatic adjustment (Lau et al., 2016). Thus any thermal model inferred from shear-wave velocities should be compatible with these observations. This general approach was pioneered by Priestley & McKenzie (2006) and Priestley & McKenzie (2013). The approach has been further refined by (Richards et al., 2020, under review) and we adopt their approach for calibration of global tomography models (see Supplementary Information).
For the Australian regional tomography models, we cannot use oceanic observations for calibration due to insufficient offshore coverage, and instead use the better constrained paleogeotherms derived from thermobarometry on mantle xenoliths (in combination with the adiabatic gradient). Away from three close together sites in South Australia in the vicinity of the Gawler Craton, it is notable that the global SL2013sv model provides a surprisingly good fit to the Australian paleogeotherms, despite being calibrated independently (Supplementary Information). This observation is unexpected for two reasons. First, the nominal resolution of the global model is lower than the local models. There are only six seismometers in Australia (located in the far west, north and east of the continent, with none in South Australia), and the density of crossing ray paths is much lower than in Europe, Asia, North, and South America (Schaeffer & Lebedev, 2013). Secondly, the Australian geotherms occur in continental lithosphere that is thought to be chemically depleted by melt extraction, reducing the quantity of garnet and clinopyroxene with respect to more fertile oceanic mantle. Nevertheless, the global model calibrated on fertile mantle constraints provides a good match to independent V S -T-P observations in depleted continental lithosphere. This result implies that temperature plays the dominant role in controlling variations in seismic wave speed in the shallow mantle, whilst the effects of compositional variation are substantially smaller (Goes et al., 2000;Priestley & McKenzie, 2006;Schutt & Lesher, 2006).
Mapping the lithosphere-asthenosphere boundary (LAB). A recent study on the thermal structure of oceanic lithosphere found that the 1175 ± 50 • C isotherm provides a good match to seismological observations of the lithosphere-asthenosphere boundary (Richards et al., 2018). In this study, we therefore adopt this isotherm as a proxy for lithospheric thickness beneath the continents. The temperature (T ) as a function of depth (z) is extracted from the V S model and ∂T ∂z is calculated over 25 km increments. Starting from the surface and progressing downwards, when temperature passes the 1175 • C threshold, LAB depth is calculated using linear interpolation, with one important exception. In locations of thick crust, low V S values at shallow depths arising from crustal bleeding can lead to lithospheric mantle being erroneously interpreted as hot. In the regional seismic tomography models, this crustal bleeding can be observed down to ∼125 km in some locations (see calibration Figure S7 in Supplementary Information). To circumvent this issue, when an inverted temperature gradient is found at shallow depths, we move on to deeper levels until temperature starts to increase with depth. This crustal bleeding is only considered down to 200 km. Maximum LAB depth is limited to 350 km or the deepest slice in the seismic tomography model. Our 1175 • C isotherm LAB proxy is shallower than used in some other studies that define the LAB using the intersection of conductive and adiabatic temperature gradients in the thermal boundary layer (typically occurring at temperatures 1350-1450 • C; McKenzie et al., 2005;Mather et al., 2011;Priestley & McKenzie, 2013). However, in addition to matching oceanic observations, the 1175 • C isotherm corresponds to lower homologous temperatures, where uncertainty in anelastic parameters has a smaller impact on the recovered LAB.
As in previous studies using seismic tomography (e.g. Priestley & McKenzie, 2013;Steinberger & Becker, 2018;Priestley et al., 2018;Afonso et al., 2019), our LAB map exhibits regions of thick lithosphere in some subduction zones (e.g. west coast of South America, south Alaska and Japan). Many of these features are likely to represent subducting slabs rather than cratonic lithosphere. None of the giant (>10 Mt of contained metal) sedimenthosted deposits is found in these settings, although some minor sedimentary copper deposits do occur, particularly in the Andes. These deposits may well represent distal components of porphyry coppers, but we have left them in our sedimentary copper dataset in line with pre-existing classification schemes. It is possible to manually exclude potential slab-related features from the analysis utilising the Slab2 model (Hayes et al., 2018;Supplementary Information). Doing so actually improves the results of statistical tests, with the chances of the relationship between sediment-hosted deposits and the edge of cratonic lithosphere being random reducing by a factor of three. This occurs because the continental area within 200 km of the 170 km LAB contour decreases from 34.3% to 31.0%, while only marginally increasing the proportion of small outlier deposits. Nevertheless, we have deliberately retained these regions in the main manuscript in order to avoid introducing subjectivity and bias into our LAB maps, as opinions are likely to differ on which features to exclude. Furthermore, some studies argue that over long periods of time, thick lithosphere may actually be generated at subduction zones by thrust stacking . Thus, exclusion of these features is potentially unwarranted.
Test suites of random continental locations. In order to test the statistical significance of real deposit locations, a test suite of random points on a sphere have been generated by randomly selecting two variables, a and b, in the range 0-1 and converting into longitude, θ, and latitude, φ, using area-normalised relationships These are subsequently filtered to select only those points that lie onshore (Supplementary Information). For each location, the closest approach to the 170 km lithospheric thickness contour is calculated and the resulting distances are plotted in a cumulative distribution function (CDF).
Kolmogorov-Smirnov statistical tests. We use the twosample Kolmogorov-Smirnov test to examine whether the difference between two cumulative distribution functions is significant, given their respective population sizes (Kolmogorov, 1933). The D-value is the maximum magnitude of the difference between two CDFs at any point. The test calculates the probability that a D-value of this magnitude might accidentally occur, had the two CDFs been randomly selected from the same underlying population. The probability, P , is approximated using where p and q are the number of samples in each CDF and D is the D-value expressed as a fraction between 0 and 1. For each Kolmogorov-Smirnov test, a number of random points are generated that is equivalent to the number of real deposits of that type (109 for PbZn-CD, 147 for PbZn-MVT and 139 for sedimentary copper). Given the low sample size for some of the deposit classes, the distribution of this random set can vary somewhat from the true average distribution of random continental locations. We therefore draw a test set in this manner 100 times and report the Kolmogorov-Smirnov statistics associated with each separate test within a histogram. For PbZn-CD deposits, the D-value between the real non-weighted, regionally enhanced CDF and each random CDF is individually calculated, yielding a mean and standard deviation of D = 0.36 ± 0.05, with extremes of 0.26-0.46. For the combined sediment-hosted deposits in Figure 3, the equivalent values are D = 0.27 ± 0.02 with extremes of 0.22-0.32. A D-value of 0.27 for the 395 combined sedimentary-hosted deposits suggests that the probability this CDF is drawn from randomly distributed continental points is less than 1 in 10 12 (Supplementary Information).
Thermal modelling of lithospheric rifting. Rifting of continental lithosphere causes subsidence of the surface to form a basin that progressively infills with sediments. An initial syn-rift subsidence phase occurs during lateral extension and vertical thinning of the crust and lithospheric mantle, which is contemporaneous with normal faulting. Following cessation of extension, faulting stops and post-rift thermal subsidence occurs as hot, upwelled asthenospheric mantle conductively cools back to an equilibrium lithospheric thickness (McKenzie, 1978). To predict the subsidence and basal heat flow of the basin, we model the 1-dimensional thermal evolution of the lithosphere during rifting using a finite difference scheme. Following McKenzie (1978), we assume thinning occurs by pure shear and that vertical heat transfer dominates. For each rift scenario, we select an initial lithospheric template. For regular continental lithosphere, the crustal thickness is set to 30 km and the total lithospheric thickness to 140 km, which matches results from plate cooling models of oceanic lithosphere (Richards et al., 2018) and places the 1175 • C isotherm at ∼ 120 km. Radiogenic heat production in the mantle is set to zero, whilst the crustal value is tuned to 1.0 µW m −3 such that the steady state geotherm yields a surface heat flow of ∼63 mW m −2 , which is the average for Phanerozoic continental lithosphere (Lucazeau, 2019) For cratonic lithosphere, we assume an initial crustal thickness of 50 km, lithospheric thickness of 280 km (1175 • C isotherm at ∼ 240 km), and crustal radiogenic heat production of 0.57 µW m −3 , which yields an initial surface heat flux consistent with the average of ∼48 mW m −2 for Archean and cratonic areas (Lucazeau, 2019) Based on the typically low paleowater depth of sediments found in proximal portions of these basins and the high supply of clastic material from adjacent cratons, we assume the basin is constantly filled by sediments. We subsequently predict the temperature of the sediment pile using the basal heat flux and a constant sediment conductivity of 2.5 W m −1 K −1 , assuming a steady state conductive geotherm and negligible internal heat generation. Details of the finite difference scheme, thermal parameterisations, and individual model runs are given in the Supplementary Information.  Supplementary Information for "Global distribution of sediment-hosted metals controlled by craton edge stability " 1 Summary of selected seismic tomography models Our LAB maps are based on recent, high-resolution shear wave tomography models that are constructed from large quantities of surface wave data. For the global map, we use vertically polarized shear wave velocity from SL2013sv upper mantle-only model, which represents the isotropic component of their azimuthally anisotropic velocities (Schaeffer & Lebedev, 2013 The SL2013sv model contains only 6 seismometers in Australia, so has limited resolution within this continent. Therefore, we also investigate three regional seismic tomography models to generate high resolution maps for the Australian continent. The main model used throughout this paper is the radially isotropic V S model FR12 of Fishwick 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 VS (Fishwick & Rawlinson, 2012). (b) AuSREM = regional VSV (Kennett et al., 2013). (c) Y14 = regional VSV (Yoshizawa, 2014). (d) SL2013sv = global VSV (Schaeffer & Lebedev, 2013). 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 VS (Fishwick & Rawlinson, 2012). (b) AuSREM = regional VSV (Kennett et al., 2013). (c) Y14 = regional VSV (Yoshizawa, 2014). (d) SL2013sv = global VSV (Schaeffer & Lebedev, 2013).
& Rawlinson (2012), which is derived from Rayleigh wave travel times (Fishwick et al., 2008). Periods considered are 50-120 s and the fundamental and first four higher modes have been used where possible, leading to good sensitivity down to ∼ 250 km depths. It contains a greater number of source-receiver paths (> 13, 000) compared to other Australian models. However, it uses an a priori crustal model that remains fixed throughout the inversion, resulting in noticeable smearing of crustal velocities into the upper mantle. Checkerboard tests indicate that features ∼ 300 km in diameter at lithospheric depths are well resolved. The second regional model is AuSREM and is a hybrid model constructed by linear combination of several previous studies (Kennett et al., 2013). It combines FR12 with YK04 and AMSAN.19 (Yoshizawa & Kennett, 2004;Fichtner et al., 2010). YK04 is a radially anisotropic Rayleigh wave model using > 8000 ray paths for the fundamental mode and ∼ 2000 for the first three higher modes, yielding a maximum period range of 40-150 s. It includes off-great circle and finite frequency effects, but also uses a fixed crustal model. AMSAN.19 is a radially anisotropic, 3D waveform, spectral element model that uses an inversion scheme based on the adjoint approach (Fichtner et al., 2009a,b). Periods considered are 30-200 s and a fixed crustal model is used. Due to the computationally intensive methodology, ∼ 3, 000 waveforms are used in this inversion.
combines Rayleigh waves (8000 fundamental, ∼ 2500 higher mode) and Love waves (approximately two-thirds as many) with periods ∼ 25-200s, corrected for local crustal structure using a fixed crustal model. It adopts the same three-step inversion procedure as YK04 (Yoshizawa & Kennett, 2004). All three models are plotted alongside the global SL2013sv model in Figures S1, S2 and S3. At any given location within the continent, V S varies between models by ∼ 0.1 km s −1 . 4 2 Parameterisation for anelasticity at seismic frequencies As introduced in the Methods section, we have adopted the parameterisation of Yamauchi & Takei (2016) to convert shear wave velocity into temperature, which includes effects of anelasticity in pre-melt conditions (temperatures above ∼ 90% of melting temperature). V S is defined as where ρ is the density and J 1 and J 2 represent real and imaginary components of the complex compliance, J * , which is a quantity describing the sinusoidal strain resulting from the application of a unit sinusoidal stress. J 1 represents the strain amplitude in phase with the driving stress, whilst the J 2 component is π 2 out of phase, resulting in dissipation. These terms contain a high temperature background absorption band and and an additional low temperature absorption peak, expressed as where J U is the unrelaxed compliance and the third term on the right of Equation (3) represents viscous relaxation. A B = 0.664 and α B = 0.38 represent the amplitude and slope of the high temperature background stress relaxation, whilst A P and σ P represent the amplitude and width of the relaxation peak superimposed on the background trend and are given by and where T ′ = T Ts is homologous temperature, with T the temperature and T s the solidus temperature, both in Kelvin. φ m is the melt fraction and β(φ m ) describes the direct poroelastic effect of melt, which is assumed to be not important within the upper mantle where only very low volumes of melt are expected to be retained (∼ 0.1%; McKenzie, 2000;Takei, 2017). For this case, J U is the inverse of the unrelaxed shear modulus, µ U (P, T ), such that where µ 0 U is the unrelaxed shear modulus at surface pressure-temperature conditions, the differential terms are assumed to be constant and the pressure, P , in GPa is linearly related to the depth, z, in km by P = z 30 . The normalised shear wave period, τ ′ S , in Equations (2) and (3) is equal to τ S 2πτ M , where τ M = η µ U is the normalised Maxwell relaxation timescale and τ S = z 1.4 is the Rayleigh wave period most sensitive to ambient velocity structure at that depth (Forsyth, 1992). τ ′ P represents the normalised shear-wave period associated with the centre of the low temperature relaxation peak, assumed to be 6 × 10 −5 . The steady-state diffusion creep viscosity, η, is where d is the grain size, m the grain size exponent (assumed to be 3 for this diffusion creep deformation mechanism), R the gas constant, E a the activation energy and V a the activation volume. Subscripts [X] r refer to reference values, assumed to be d r = 1 mm, P r = 1.5 GPa and T r = 1200 • C for the upper mantle. In this study, we make the simplifying assumption that d = d r , which indicates an endmember scenario whereby lateral changes in V S within the upper mantle arise purely from variations in temperature rather than grain size. It is also possible that grain size may vary significantly within the shallow mantle, but remains poorly constrained (Behn et al., 2009;Dannberg et al., 2017). A η represents the extra reduction of viscosity due to an increase in E a near the solidus, expressed as where T ′ η = 0.94 is the homologous temperature above which the effective activation energy increases beyond its original value and γ = 5 is the factor of additional viscosity reduction. λφ describes the direct effect of melt on viscosity, also assumed to be negligible at low melt volumes. The solidus temperature, T s , is fixed to a value of 1326 • C at 50 km equivalent to the dry peridotite solidus of Hirschmann (2000), and linearly increases below this depth according to where ∂Ts ∂z is the solidus gradient. We use a temperature-dependent, compressible density, ρ(P, T ), following the approach of Grose & Afonso (2013). First, we define a linear temperature-dependence on thermal expansivity, α(T ), such that where α 0 = 2.832 × 10 −5 K −1 and α 1 = 0.758 × 10 −8 K −2 are constants calibrated from mineral physics experiments (Bouhifd et al., 1996). To include pressure-dependence, the isothermal volume change, (V 0 /V ) T is calculated at each pressure using a Brent minimisation algorithm and the third-order Birch-Murnaghan equation of state where K 0 = 130 GPa is the bulk modulus at zero pressure and K ′ T = 4.8 is the pressure-derivative of the isothermal bulk modulus. The associated isothermal density change with pressure, ρ(P ), is given by where ρ 0 = 3.33 Mg m −3 is the density of mantle at surface pressure and temperature. The effect of pressure on 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 where T 0 = 273 K is temperature at the surface. In a similar manner to Equation (1), the shear-wave attenuation, Q −1 S , can be defined as

Thermobarometry of mantle xenoliths and xenocrysts
In order to calibrate the anelasticity parameterisation for Australian regional seismic tomography models, a series of local V S -T -P tie points are required. For temperature estimates, we have assembled a suite of Australian paleogeotherms derived from thermobarometric analysis of mantle xenoliths and xenocrysts from fifteen locations in thick and thin lithosphere ( Figure S4; Nickel & Green, 1985;Nimis & Taylor, 2000;Taylor, 1998). As outlined in the methods, the resulting P-T estimates are entered into FITPLOT to generate palaeogeotherms ( Figure S5; Mather et al., 2011). 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 2b), derived from FR12 seismic tomography model (Fishwick & Rawlinson, 2012). 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 (Nickel & Green, 1985;Taylor, 1998); blue circles = P-T estimates derived from single chrome diopside thermobarometry (Nimis & Taylor, 2000); dashed line = crustal thickness from AusMoho (Kennett et al., 2011); solid line = FITPLOT optimal paleogeotherm (Mather et al., 2011); coloured band = spread of 1000 geotherms from Monte Carlo FITPLOT analysis (Methods).

Inverse calibration of anelasticity parameterisation
Anelasticity parameters A B , α B , τ ′ P , β(φ m ), γ, T ′ η and λφ have been directly constrained by forced oscillation experiments on borneol (Yamauchi & Takei, 2016). However, µ 0 U , ∂µ U ∂T , ∂µ U ∂P , η r , E a , V a and T S (z) are material properties that must be independently determined. As outlined in the Methods, we adopt the methodology of Richards et al. (2020) to invert a suite of observational constraints. First, the SL2013sv global V S model is stacked in oceanic regions to calculate average V S as a function of depth and lithospheric age (Schaeffer & Lebedev, 2013). The age grid and optimal thermal model for a cooling oceanic plate are taken from Richards et al. (2018).
Due to the limited vertical resolution of tomographic models (25-50 km), we extract a suite of V S versus temperature tie-points by taking the average V S between 75 km and 100 km and cross-plotting it with temperature obtained from the thermal model at 87.5 km, which marks the midpoint of this depth range. A second curve is obtained using the where V o ij are observed shear-wave velocities with associated standard deviation σ ij , V c ij is the prediction from Equation (1), M = 76 is the number of age bins at a given depth and N is the number of depth slices (two in this case).
A second suite of tie-points is created by assuming that temperatures are isentropic at depths well below the upper thermal boundary layer. We calculate average V S as a function of depth over oceanic regions in the global model, and over the whole spatial domain in regional models. Over the depth range 225-400 km, beyond which the resolving power of surface waves drops significantly, these values are combined with an isentrope calculated for a potential temperature of 1333 • C and an adiabatic gradient determined using a reference density of ρ 0 = 3.3 Mg m −3 , thermal expansivity of α = 3 × 10 −5 • C −1 and specific heat capacity of C P = 1187 J kg −1 • C −1 . 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 ages ≥ 100 Ma (Adenis et al., 2017). Over this age range, we stack the QRFSI12 attenuation model of Dalton et al. (2009), generating a suite of Q −1 S to V S tie-points as a function of depth. Equations (1) and (15) are coupled such that average temperature is obtained from the average V S , rather than assuming isentropic temperatures extend up to 150 km. Misfit, H 3 , between observed and predicted attenuation is We also adopt a bulk diffusion creep viscosity of η ref = 3×10 20 Pa s for the upper mantle (∼ 100-670 km) obtained from glacial isostatic adjustment studies (Lau et al., 2016), and compare it to the average predicted value for 225-400 km depths obtained from Equation (7). Misfit, H 4 , is calculated using where η c i is predicted viscosity and the viscosity uncertainty, σ i , is assumed to be one order of magnitude. Finally, for calibration of regional tomography models where these global oceanic observations are unavailable, we take the better constrained paleogeotherms derived from thermobarometry on mantle xenoliths. Argyle, Boowinda Creek, Bullenmerri, Ellendale, Merlin, Monaro, Monk Hill, Orroroo and Wandagee are used to directly constrain each anelasticity model.
None of these paleogeotherms show evidence of having been perturbed by heating events immediately prior to xenolith entrainment, and we therefore take the calculated P-T conditions to represent ambient mantle conditions immediately where M is the number of paleogeotherms, N is the number of tie-points associated with each geotherm and σ ij reflects uncertainty in the V S measurement, assumed to be a constant 0.1 km s −1 which captures typical variations 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 parameter sweep is performed to identify the approximate location of the global minimum. µ 0 U is varied between 69-82 GPa (in increments of 1 GPa), ∂µ ∂T between −20 and −8 MPa • C −1 (2 MPa • C −1 increments), ∂µ ∂P between 1.5-2.9 (0.2 increments), η r between 10 17 -10 23 Pa s (10 0.5 Pa s increments), E a between 100-1000 kJ mol −1 (100 kJ mol −1 increments), V a between 0-30 cm 3 mol −1 (2 cm 3 mol −1 increments) and ∂Ts ∂z between 0-4.5 • C km −1 (0.25 • C km −1 increments). These extrema are chosen to exceed the range of previous parameter estimates obtained from laboratory experiments and other studies (Priestley & McKenzie, 2013;Yamauchi & Takei, 2016;Jain et al., 2018). Secondly, Powell's conjugate gradient algorithm is used to further minimise H using best-fitting parameters from the initial sweep as the starting point (Press et al., 1992). There are two broad sets of parameter trade-offs. The first occurs between parameters controlling anharmonic velocities (µ 0 U , ∂µ ∂T , and ∂µ U ∂P ). The second is between parameters controlling anelastic effects (η r , E a , V a , and ∂Ts ∂z ). The change in misfit value, H, along each of these trade-offs is minor, such that multiple parameter combinations yield similar conversions between V S and temperature (Richards et al., 2020). This effect means that whilst the uncertainty on each anelastic parameter may be high for a given tomography model, the resulting mantle temperature structures remain largely consistent.
For calibration of the four global models, we set w 1 = 10, w 2 = 1, w 3 = 2, w 4 = 2 and w 5 = 0. In the case of SL2013sv, the minimum misfit is H = 0.568 and occurs when µ 0 67, η r = 4.45 × 10 22 Pa s, E a = 400 kJ mol −1 , V a = 0.092 cm 3 mol −1 and ∂Ts ∂z = 0.919 • C km −1 . These parameters are then used to convert the full three-dimensional V S model to temperature. Optimal anelasticity parameters for each of the calibrated seismic tomography models are listed in Table 1. For comparison, the optimal misfit for CAM2016 is H = 0.929 and for 3D2015-07Sv is H = 0.887, which is partly why we prefer the SL2013sv model. Our blended global model SLNAAFSA achieves the lowest optimal misfit of H = 0.533. For Australian regional models, we constrain the calibration using the nine most reliable paleogeotherms. All Figure S6: VS as a function of depth at sites of fifteen Australian paleogeotherms. Labels give site name (locations in Figure S4); red = global SL2013sv modelSchaeffer & Lebedev (2013); purple = regional FR12 modelFishwick & Rawlinson (2012); blue = regional AuSREM modelKennett et al. (2013); orange = regional Y14 modelYoshizawa (2014).
Note that the global model SL2013sv yields good fits to paleogeotherms away from south Australia (Monk Hill, Orroroo and Cleve), despite being lower resolution than the local models and being calibrated completely independently of this information (red lines in Figure S7). Conversely, regional models often provide a poorer fit to the full range of the paleogeotherms and can exhibit substantial crustal bleeding artefacts at depths shallower than ∼ 125 km. 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 AusMohoKennett et al. (2011); solid line = optimal FITPLOT geotherm from Figure S5; purple = regional FR12 modelFishwick & Rawlinson (2012); blue = regional AuSREM modelKennett et al. It is important to note that of the nine geotherms used to calibrate the anelasticity parameterisation for the regional FR12 model, two are only constrained by a single P-T estimate (Orroroo and Boowinda Creek). We have therefore tested the effect of removing these two sites from the calibration scheme. As Figure S8 shows, there is no discernible effect on the inferred temperature structure. Given this result and that these two samples pass all of the thermobarometry cation and oxide tests, we have chosen to keep Orroroo and Boowinda Creek within the set of nine geotherms used in calibration of local tomography models. Figure S8: Calibration of anelasticity parameterisation on Australian paleogeotherms. Labels give site name and inferred age of paleogeotherms in million years (locations in Figure S4); red line = FR12 model calibrated using sites Argyle through Wandagee; green line = same but excluding Boowinda Creek and Orroroo from the calibrations set.

Global lithospheric thickness maps
We have calibrated two other recent surface wave tomography models that each have global coverage (3D2015-07Sv and CAM2016; Debayle et al., 2016;Priestley et al., 2018). In addition, we have generated a higher resolution version of SL2013sv, referred to as SLNAAFSA, by blending in three regional updates including SL2013NA in North America, AF2019 in Africa, and SA2019 in the South Atlantic Ocean (Schaeffer & Lebedev, 2014;Celli et al., 2020a,b). The lithospheric thickness maps for all four global models are visually similar and are shown in the extended data figure of the manuscript and available in ASCII, GeoTIFF and GMT NetCDF formats in the online Supplementary Information.
All four exhibit robust relationships with sediment-hosted deposits.
A histogram of global lithospheric thickness derived from the SL2013sv model is shown in Figure S9 and reveals a bi-modal population, with peaks at ∼ 90 km and ∼ 190 km, separated by a minimum at 150 km. There is also a noticeable drop-off deeper than 200 km, which we attribute to a change in the gradient of V S with depth in the initial starting profile used to construct the tomography model. Figure S9: Area-weighted histogram of global LAB depths. LAB derived from the SL2013sv tomography model (Schaeffer & Lebedev, 2013); black bars = oceanic regions; red bars = continental regions.
For comparison, we also provide seven previously published global lithosphere-asthenosphere boundary (LAB) maps derived from a mixture of heat flow data, seismic tomography datasets, and potential field data. Interestingly, many giant sediment hosted mineral deposits lie along LAB edges defined by these other studies, testifying to the robustness of the observed relationship. Figure S10: Previously published global maps of depth to the lithosphere-asthenosphere boundary. Symbols = sediment-hosted deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed). (a) LAB derived from surface heat flow measurements (Artemieva, 2006); (b) LAB derived from surface wave tomography (Priestley & McKenzie, 2013); (c) LAB derived from vertical shear-wave travel time anomalies in the continents (Bird et al., 2008); (d) LAB derived by Steinberger & Becker (2018)

Australian lithospheric thickness maps
For each of the Australian regional seismic tomography models considered in this study, we have individually calibrated and mapped out the LAB in a consistent manner. The resulting maps for Australia are shown in Figure S11, whilst in Figure S12 we compare our preferred FR12 regional model to previously published maps of LAB depth beneath Australia.

Kolmogorov-Smirnov Statistical Tests
In order to test the statistical significance of real deposit locations, test suites of random points on a sphere are generated (Methods). Example test suites of 100, 1000 and 10,000 points are shown in Figure S13.  For each Kolmogorov-Smirnov test, a number of random points are generated that is equivalent to the number of real deposits of that type (109 for PbZn-CD, 147 for PbZn-MVT and 139 for Cu-sed). Given the low sample size for some of the deposit classes, the distribution of this random set can vary somewhat from the true average distribution of continental locations. We therefore draw a test set in this manner 100 times ( Figure S14). These random CDFs are relatively consistent but have some outliers. The D-value and Kolmogorov-Smirnov statistic between each random CDF and the real one is calculated and reported within a histogram ( Figure S15).

Effect of subducting slabs on deposit statistics
The relationship between seismic wavespeed and temperature results in a tomography model at depths > 150 km imaging cold subducting slabs as fast velocities, in addition to those associated with thick, cold cratonic lithosphere.
Thus it is possible that some of the features imaged in our LAB maps are not related to cratonic lithosphere. Examples include features along the Andean, Caribbean, Aleutian, Japanese, Philippines, Indonesian, and eastern Mediterranean subduction zones ( Figure S16). None of the giant (> 10 Mt contained metal) sediment-hosted deposits are located at these subduction zones.
Nevertheless, some of the smaller deposits can be, such as those in South America, Indonesia and Turkey. We have therefore manually excised 170 km LAB contours that may potentially be related to slabs (red polygons on Figure S16) and repeated the statistical tests.
The CDF for all sediment-hosted base metal deposits is essentially unchanged by this procedure, with ∼ 85% of total metal still found within 200 km of the 170 km contour ( Figure S17c). However, the deposit statistics are actually improved, with the D-value increasing from 0.270 ± 0.020 to 0.276 ± 0.022, changing the probability of the relationship occurring by random chance from 1 in 3.35 trillion to 1 in 10.6 trillion. The reason for this is that the reduction in contour extent results in fewer random continental locations falling near the line, with the percentage of total continental area within 200 km of the 170 km LAB thickness contour dropping from 34.3% to 31.0%.
Nevertheless, we have chosen to use the full, non-excised LAB in the paper. Manual identification of potential slabs is a subjective process, with results dependent upon an individual's opinion. Furthermore, it is possible that rifted cratonic fragments may be transported into subduction zone settings , and thus not all of these subduction zone features are necessarily anomalous. We therefore prefer to keep our testing of the veracity of the observed relationship between LAB thickness and ore deposit locations as objective as possible.  Figure S16); dotted line = simple count of number of deposits with increasing distance from the 170 km contour; solid black line = deposits weighted by mass of contained metal; grey line/bounds = mean and standard deviation of 100 sets of equivalent number of randomly drawn continental locations. (b) Histogram of D-values for ensemble of 100 random CDFs calculated for a random test set of continental points compared with the non-mass-weighted CDF. (c) and (d) same but using 170 km LAB thickness contours with potentially anomalous subduction zone features removed (only green polygons in Figure S16).

Rift modelling of continental lithosphere
To predict the subsidence and basal heat flow of a basin formed by rifting, we model the thermal evolution of the lithosphere during extension. Following McKenzie (1978), we assume thinning occurs by pure shear and that vertical heat transfer dominates. We start with the one-dimensional heat flow equation where t is time, z is depth, T is temperature, P is pressure, X is composition, ρ is density, C P is the isobaric specific heat capacity, k is the thermal conductivity, and H is the internal radiogenic heat production.
We solve Equation (22) numerically using an unconditionally stable time-and space-centered Crank-Nicholson finitedifference scheme with a predictor-corrector step (Press et al., 1992). Equation (22) is recast as where ∆t is the time step, ∆z is the depth spacing between nodes, and n and j are the time and depth indices, respectively. Equation (23) is solved by tridiagonal elimination (Press et al., 1992). For the initial predictor phase of each time step, m = n, whilst in the subsequent corrector phases, m = n + 1. We use a Lagrangian reference frame, whereby ∆z is initially set to 1 km and updates using the strain rate for each timestep. Timesteps are calculated using a Courant-Friedrichs-Lewy condition with the Courant number set equal to five, such that and T n+1 typically convergences to within a tolerance of 0.001 • C after two corrector phases. The strain rate is constant during rifting and is set by rift duration and a stretching factor, β, which is the ratio of initial to final crustal thickness.
For the crustal layer, we adopt constant thermal parameters of C P = 750 J kg −1 K −1 and k = 2.5 W m −1 K −1 .
Crustal density is ρ = 2800 kg m −3 for regular continental lithosphere and is ρ = 2900 kg m −3 for cratonic lithosphere, in accordance with expectations for an increase in bulk density with increasing crustal thickness (Christensen & Mooney, 1995). For the mantle, we use temperature and pressure dependent parameters that have been derived from experimental data on olivine. Density is calcuated according to Equation (14), with ρ • = 3330 kg m −3 in regular lithospheric mantle, and 3280 kg m −3 in cratonic lithosphere, which has been chemically depleted by melt extraction. The effect of pressure on heat capacity is minimal throughout the upper ∼ 300 km of the mantle (Hofmeister, 2007). We therefore use the Korenaga & Korenaga (2016) purely temperature-dependent parameteristaion where c 0 = 1580, c 1 = 12230 and c 2 = 1694 × 10 6 . Mantle conductivity is pressure-and temperature-dependent and includes contributions from both phonons (lattice conductivity, k lat ) and photons (radiative conductivity, k rad , which is poorly constrained by experimental data). Total conductivity is expressed using where ∂ln(k) ∂P = 0.05 W m −1 K −1 GPa −1 (Hofmeister, 2007). κ lat (T ) is the temperature-dependent lattice diffusivity and is given by where κ 0 = 0.565, κ 1 = 0.67, κ 2 = 590, κ 3 = 1.4, κ 4 = 135, and T 0 = 273 K (Pertermann & Hofmeister, 2006). The parameterisation for k rad is taken from Grose & Afonso (2013), based on earlier work by Hofmeister (2005), and is given by and d is the grain size in centimeters, assumed to be d = 0.5 cm.
For the finite difference scheme boundary conditions, we fix the surface node to have T n 0 = T 0 = 273 K, whilst the initial basal node has an adiabatic value of (1606 + 0.44z) K, equivalent to a potential temperature of 1333 • C. In cratonic areas, the lithospheric mantle is thicker than standard continental lithosphere and has been chemically depleted.
During the rift phase, this basal node shallows through time and non-depleted asthenospheric mantle rises adiabatically beneath. If this basal node becomes shallower than the initial thickness of standard continental lithosphere, we update the index at which this lower boundary condition is applied to the node closest to this standard depth. Heat flow through the top of the crust, Q H (t), and subsidence, S(t), are calculated according to where J is the index of the node at the depth of the original lithospheric thickness, ρ J is the adiabatic density of standard mantle at this depth, and ρ inf ill = 2200 kg m −3 is the density of material that infills the basin, which we assume to be 24 sediments. Figures S18-S20 show the results of thermal modelling of rifting continental lithosphere on basin subsidence and temperature of the sedimentary pile.
10 Regional deposit maps  Figure S21: Distribution of sediment-hosted base metal deposits as a function of lithospheric thickness in Africa. Global LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016); 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 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed). Figure S22: Distribution of sediment-hosted base metal deposits as a function of lithospheric thickness in Europe. Global LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016); 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 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed). Figure S23: Distribution of sediment-hosted base metal deposits as a function of lithospheric thickness in Asia. Global LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016); 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 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey; circles = clastic-dominated lead-zinc (PbZn-CD); triangles = Mississippi Valley type lead-zinc (PbZn-MVT); squares = sedimentary copper (Cu-sed). Figure S24: Distribution of sediment-hosted and IOCG base metal deposits as a function of lithospheric thickness in Australia.

31
Figures S30-S32 show deposit locations, age distributions with respect to lithospheric thickness, and Kolmogorov-Smirnov statistical test results for each individual deposit type compiled in this study (Methods). These are followed by tabulated information and associated references for each individual deposit, which is also available as a spreadsheet through the online Supplementary Information. Figure S28: 109 clastic-dominated lead-zinc deposits from Table 2. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey.   Table 3. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 1 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey.   Table 4. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = 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.  : 108 magmatic nickel-copper-platinum group element deposits from Table 5. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 0.5 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey.   Table 6. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = deposit locations; area proportional to estimate of total contained mass of metal (Mt = megatonnes); unknown deposit size given 0.5 Mt symbol; colour = ore body formation age (billion years); unknown age plotted in grey.   Table 7. (a) LAB derived from SL2013sv tomography model using a calibrated anelasticity parameterisation (Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). Circles = 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.   Figure S28, building extensively on Taylor et al. (2009). Ga = billion years; Mt = million tonnes; ND = no data.