The Hudson Bay Lithospheric Experiment (HuBLE): insights into Precambrian plate tectonics and the development of mantle keels

Abstract Hudson Bay Lithospheric Experiment (HuBLE) was designed to understand the processes that formed Laurentia and the Hudson Bay basin within it. Receiver function analysis shows that Archaean terranes display structurally simple, uniform thickness, felsic crust. Beneath the Palaeoproterozoic Trans-Hudson Orogen (THO), thicker, more complex crust is interpreted as evidence for a secular evolution in crustal formation from non-plate-tectonic in the Palaeoarchaean to fully developed plate tectonics by the Palaeoproterozoic. Corroborating this hypothesis, anisotropy studies reveal 1.8 Ga plate-scale THO-age fabrics. Seismic tomography shows that the Proterozoic mantle has lower wavespeeds than surrounding Archaean blocks; the Laurentian keel thus formed partly in post-Archaean times. A mantle transition zone study indicates ‘normal’ temperatures beneath the Laurentian keel, so any cold mantle down-welling associated with the regional free-air gravity anomaly is probably confined to the upper mantle. Focal mechanisms from earthquakes indicate that present-day crustal stresses are influenced by glacial rebound and pre-existing faults. Ambient-noise tomography reveals a low-velocity anomaly, coincident with a previously inferred zone of crustal stretching, eliminating eclogitization of lower crustal rocks as a basin formation mechanism. Hudson Bay is an ephemeral feature, caused principally by incomplete glacial rebound. Plate stretching is the primary mechanism responsible for the formation of the basin itself.

Much of the geological record can be interpreted in the context of processes operating at the presentday plate boundaries. While the plate tectonic paradigm works well to explain processes and products during the Phanerozoic era, during Precambrian times, when the oldest rocks were forming, conditions on the younger, hotter, more ductile Earth were probably very different, making analogies with modern-day tectonics less certain. The precise onset of 'continental drift' is disputed; it has, for example, been estimated to be as early as c. 4.1 Ga (e.g. Hopkins et al. 2008), or as late as c. 1 Ga (e.g. Stern 2005) -a time span covering approximately two-thirds of Earth history (Fig. 1). Gathering evidence preserved deep within the plates in stable Precambrian regions (shields) is thus essential to improve our understanding of early Earth processes.
Cratonic regions are readily identified in global tomographic images where they are characterized by their deep-seated, fast wavespeed lithospheric keels. Roots can extend to depths of ≥250 km into the upper mantle, in contrast to the oceans and Phanerozoic continents where the plates are usually ≤100 km thick ( Fig. 2; e.g. Ritsema et al. 2011). The tectosphere, or lithospheric mantle beneath a craton (e.g. Jordan 1988), is thus thought to have a thermochemical signature that differs from average lithospheric mantle, and keel formation is commonly associated with Archaean processes, such as the extraction of komatiitic magmas (e.g. Griffin et al. 2003), to explain the intrinsic low density of the tectosphere. It is estimated that the chemical lithosphere must be more viscous than normal mantle by a factor of c. 20 (Sleep 2003), enabling it to survive thermal and mechanical erosion during multiple Wilson cycles over billions of years. Precisely how keels formed remains poorly understood, however.
The geological record of northern Canada spans more than 2 billion years of early Earth history (c. 3.9-1.7 Ga; Hoffman 1988), making it an ideal study locale for the analysis of Precambrian Earth processes. In the heart of the Canadian Shield lies Hudson Bay, a vast inland sea, which masks a significant portion of the Trans-Hudson Orogen (THO), a Palaeoproterozoic collision between the Archaean Superior and Western Churchill cratons. The THO signalled the final stages of assembly of Laurentia at c. 1.8 Ga (Fig. 3: Hoffman 1988;St-Onge et al. 2006;. The region is underlain by one of the largest continental keels on Earth (Fig. 2) and is also the site of one of the largest negative geoid anomalies (e.g. Hoffman 1990). The cratonic keel beneath the Hudson Bay region extends beneath both Archaean and Proterozoic terranes. How these various lithospheric subdivisions of the Canadian shield were assembled in Precambrian times, how the deep cratonic keel formed beneath them, and how the Hudson Bay basin subsequently developed within the shield, have each been difficult questions to address until recently. Hypotheses have been based principally on evidence from the geological record and potential field maps (e.g. gravity and magnetics), and seismological studies of the deeper mantle were limited by the availability of data from only two or three permanent seismic stations.
To address outstanding research questions exemplified by the Hudson Bay region, a broadband seismograph network in the Hudson Bay region was deployed by the Hudson Bay Lithospheric Experiment (HuBLE); stations have operated since 2003. The purpose of this contribution is to review the results of HuBLE, and to synthesize their implications for Precambrian processes, the formation of the Laurentian keel, and the reasons for the Phanerozoic development of the Hudson Bay basin. In addition to reviewing already-published HuBLE constraints (e.g. Thompson et al. 2010;Bastow et al. 2011a, b;Pawlak et al. 2011Pawlak et al. , 2012Steffen et al. 2012;Snyder et al. 2013;Darbyshire et al. 2013), this contribution presents the results of a new P-wave relative arrival-time tomographic study of mantle wavespeed structure beneath the region.

Tectonic setting
The Canadian Shield lies in the heart of Precambrian North America. Comprising several Archaean terranes brought together during a series of orogens during the Palaeoproterozoic, the region is one of the largest exposures of Precambrian rocks on Earth (Hoffman 1988). The major phase of mountain building in the region was believed to have occurred at c. 1.8 Ga during the THO (Hoffman 1988;St-Onge et al. 2006). The Superior craton formed the indenting lower plate to the collision; the Churchill plate, presently lying to the north of the Superior craton, was the upper plate. Extensive geological mapping in recent years has added considerable detail to this two-plate picture (Fig. 3).
It is now thought that northern Hudson Bay comprises seven distinct crustal blocks, including the Superior, Rae, Hearne, Chesterfield, Meta Incognita, Sugluk and Hall Peninsula ( Fig. 3; see e.g. Snyder et al. 2013;Berman et al. 2013, for a more detailed review). The Superior craton to the south is a collage of Meso-to Neoarchaean microcontinents and oceanic terranes amalgamated 2.72 -2.66 Ga. The Ungava Peninsula bounds eastern Hudson Bay and exposes 3.22 -2.65 Ga felsic orthogneiss and plutonic rocks that underlie c. 2.0-1.87 Ga volcanic and sedimentary rocks of the Cape Smith fold belt (e.g. St-Onge et al. 2009, and references therein).
The Hearne craton is characterized by c. 2.7 Ga mafic to intermediate volcanic rocks (e.g. central Hearne supracrustal belt; Davis et al. 2004) cut by c. 2.66 Ga granitic plutons associated with regional metamorphism (Davis et al. 2004). The Hearne is distinguished from Rae by its relatively young (Mesoarchaean) volcanic rocks, and by its lack of Archaean magmatism or tectonism post 2.66 Ga. The Chesterfield block (C.B. in Fig. 3 Wodicka et al. 2011b), overlain by Palaeoproterozoic clastic carbonates of the c. 1.9 Ga Lake Harbour Group (Scott 1997;Wodicka et al. 2011b). The Sugluk block (e.g. Corrigan et al. 2009, and references therein) includes Mesoarchaean crust that crops out on SW Baffin Island and nearby islands in Hudson Bay. It is inferred from evolved Nd signatures to underlie the northern tip of Quebec (e.g. Corrigan et al. 2009) where it is heavily intruded by younger c. 1.86 -1.82 Ga plutons (St-Onge et al. 2006). Corrigan et al. (2009) suggested that it may extend west under Hudson Bay to merge with a gravity high referred to as the 'Hudson protocontinent' (Roksandic et al. 1987;Berman et al. 2005). Hall Peninsula, southeastern Baffin Island, exposes Mesoarchaean crust (2.92-2.8 Ga zircon crystallization ages; Scott 1997) which, given its spatial distribution, is tentatively considered a separate crustal block (e.g. Whalen et al. 2010). Amalgamation of Meta Incognita and the Sugluk block is postulated to have occurred in an intra-oceanic setting shortly before c. 1.9 Ga (Wodicka et al. 2011b). The Hall Peninsula block was also considered to be a part of this composite block, which is referred to by  Berman et al. 2010a, for a detailed summary and metamorphic maps). The oldest event, the 2.56-2.50 Ga MacQuoid orogeny, affected the Chesterfield block and adjacent southeastern flank of the Rae craton. The 2.5-2.3 Ga Arrowsmith orogeny (Berman et al. 2005(Berman et al. , 2010aSchultz et al. 2007) affected much of the western side of the Rae craton from northern Saskatchewan through Boothia Peninsula to northern Baffin Island (Berman et al. 2010a. The 2.0-1.9 Ga Thelon orogeny (Hoffman 1988) affected the westernmost Rae craton, including northern Boothia Peninsula (Berman et al. 2010a). The most widespread tectonic/ metamorphic reworking of the Hudson Bay region occurred during the 1.9-1.8 Ga Hudsonian orogeny (e.g. Hoffman 1988;Berman et al. 2010b) when accretion of microcontinents to the SE flank of the Rae craton occurred c. 1.9-1.87 Ga. These tectonic events are known collectively as the Snowbird phase of the Hudsonian orogeny (Berman et al. 2007).
Proposed under-thrusting of the Superior craton beneath the Churchill collage and the Cape Smith belt was initiated by c. 1.82 Ga (St-Onge et al. 2006). Palaeozoic events include deposition of carbonate rocks, presently estimated to be up to 1.5 km thick, on crystalline basement rocks of Hudson Bay, western Southampton Island, and nearby Coats and Mansel islands. At least five separate, potentially diamondiferous, kimberlite fields are known to have erupted during the Mesozoic and Cenozoic, providing mantle xenoliths and garnet xenocrysts that sample the mantle below the study area (Fig. 3).

Previous seismological studies of the Canadian Shield
The first constraints on the seismic structure of the Hudson Bay region were presented in the early 1960s by Brune & Dorman (1963) using surfacewave dispersion computed from two-station paths typically thousands of kilometres long. Hobson (1967) and Ruffman & Keen (1967) subsequently studied crustal structure using controlled source data recorded by both land and sonobuoy recorders. Since then, several studies have mapped wavespeeds beneath the Canadian Shield using surface waves in regional to global-scale models. In general the global models use data from permanent seismograph networks, resulting in reduced resolution beneath the Bay region; the two or three seismograph stations contributing to the inversions yield a limited number of surface wave paths across the shield (e.g. Megnin & Romanowicz 2000;Shapiro & Ritzwoller 2002;Lebedev & van der Hilst 2008). They indicate that Hudson Bay is underlain by a c. 200-300 km deep fast wavespeed mantle keel (e.g. Lebedev & van der Hilst 2008;Li et al. 2008;Nettles & Dziewonski 2008;Ritsema et al. 2011). Continent-scale models such as those of Bedle & van der Lee (2009) and Yuan et al. (2011) also indicate a deep-seated, fast wavespeed mantle keel beneath the Canadian Shield.
Some of the tomographic inversions take seismic anisotropy into account, but the results are variable. The model of Debayle et al. (2005), for example, yields east-west fast anisotropic directions beneath the Bay region, while the model of Marone & Romanowicz (2007) indicates NW-SE fast directions in the same depth interval, rotating to NE-SW at 300 km. Yuan et al. (2011) presented a continent-scale model of seismic heterogeneity within the North American craton, using joint inversions of long-period waveforms and SKS splitting data. They found evidence for a two-layered lithosphere. The 150-200 km-thick upper layer has fast polarization directions correlated with trends in surface geology, and was interpreted as the chemically depleted 'lid' of the shield. The lower lithospheric layer, with different anisotropic characteristics, was interpreted by Yuan et al. (2011) as a significantly younger, less-depleted thermal boundary layer. The Yuan et al. (2011) model attributes anisotropic fast polarization directions at asthenospheric depths to anisotropic fabric development owing to absolute plate motion. This model places the lithosphere-asthenosphere boundary (LAB) at a depth of c. 200 km depth everywhere beneath the Bay region and throughout the entire Canadian Shield, based on the anisotropic stratification of the upper mantle model. S to P receiver function studies (Rychert & Shearer 2009;Abt et al. 2010;Miller & Eaton 2010;Kumar et al. 2012) provide support for an abrupt change in seismic wavespeed at c. 80 -150 km depth across much of North America. If interpreted in the context of the two-layered lithosphere hypothesis of Yuan et al. (2011) or Darbyshire et al. (2013), this feature probably represents a mid-lithospheric discontinuity beneath cratonic North America.
The c. 200 -260 km lithospheric thickness in central and northern Canada inferred in the surface wave studies is corroborated by estimates of the thermal boundary layer, as inferred from joint interpretation of surface heat flow and S-wave travel time delays (e.g. Lévy et al. 2010). Although based on sparse data coverage, heat flow data from northernmost Ontario and central-northern Quebec are the lowest anywhere across the Canadian Shield. The thickness of the thermal lithosphere beneath the Canadian Shield varies regionally by up to c. 100 km (Artemieva 2006;Lévy et al. 2010).
Previous studies of the mantle transition zone beneath the Canadian Shield were limited to imaging the receiver-side structure using a small number of permanent stations located large distances apart (Bostock 1996;Chevrot et al. 1999), geographically restricted temporary networks (Li et al. 1998) or global studies utilizing mid-point reflections from distant source -receiver combinations (e.g. Gossler & Kind 1996;Gu et al. 1998;Flanagan & Shearer 1998;Gu & Dziewonski 2002;Lawrence & Shearer 2006). Several of these studies showed that the transition zone structure was uniform across the region (less than +5 km; Bostock 1996), yet sparse geographical coverage of receiver functions or lower resolution (PP and SS precursors) meant that it was unclear whether this pattern was true beneath large swathes of the shield.
On a more local scale, the Teleseismic Western-Superior Transect (TWIST; see Kendall et al. 2002 for a review) experiment deployed 11 short-period and 14 broadband seismometers along a 600 km line in the Superior region of SW Hudson Bay in May-November 1997. A further three broadband stations were deployed further north on the Hudson Bay coastline (e.g. Kay et al. 1999). TWIST data were used in regional seismic tomographic imaging of upper mantle wavespeed structure  in an SKS study of mantle anisotropy (Kay et al. 1999) and a receiver function study of crustal structure (Angus et al. 2009). The tomographic study did not yield evidence for differences between Archaean and Proterozoic mantle wavespeed. In contrast, the splitting study revealed marked differences between the two domains, with null measurements characterizing the Proterozoic coastal areas, while moderate to large splitting was found to parallel Archaean Superior trends.

The HuBLE seismograph networks
Preliminary discussions began in the early 1990s concerning a Hudson Bay Lithospheric Experiment designed to study crust -mantle dynamics beneath Hudson Bay, to determine the reason for this intracratonic basin, and to understand the Precambrian processes that shaped the central Canadian Shield. Four co-ordinated activities were discussed: a network of three-component broadband stations deployed for several years, deep seismic reflection profiling along 600 km profiles, ocean-bottom seismometers and single-component land stations to record the large air-gun array to be used in the profiling. International co-ordination and funding of the marine profiling has not yet been achieved, but the network of three-component broadband stations has now been completed by the Canadian and UK components of HuBLE.
The UK component of the HuBLE project was deployed in the summer of 2007 by personnel from the University of Bristol, in collaboration with the Geological Survey of Canada. These 10 stations were situated mainly around the islands of northern Hudson Bay, providing excellent azimuthal coverage around the Bay when combined with the existing POLARIS and a few permanent Canadian National Seismograph Network stations (Fig. 3). Figure 4 shows a completed HuBLE-UK Natural Environment Research Council seismograph station in northern Hudson Bay. Each site was equipped with a Güralp CMG-3TD broadband seismometer, recording at 40 samples per second (sps). Güralp data collection modules were used at the stations, which were powered by up to six solar panels (providing 100-140 W power) and three 100 A h deep-cycle batteries. Each remote site was equipped with an iridium satellite modem that provided access to state-of-health data from the stations. Utilizing modems over the iridium network provides pole-to-pole global coverage of both short message and short data burst services. In the case of the Güralp data collection module, this allowed the UK base-station to access weekly reports of remote-system state of health. Where problems were diagnosed, low-latency two-way communication for reconfiguration of the remote systems was also utilized via simple terminal interaction. Further details of the UK component of the HuBLE field campaign can be found in Bastow et al. (2011a).
Nunavut, the homeland of the Inuit, is the most sparsely populated region of Canada, so centres of population and infrastructure are few and far between. Wherever possible, seismograph stations were deployed in secure compounds such as airports and weather stations with mains power supply. Elsewhere, vast tracts of wilderness meant that remote, independently powered recording sites had to be designed. Transportation to these locations was by chartered helicopters or light aircraft with large tundra tyres that permitted landing on relatively flat and well-drained glacial deposits (Fig.  4). Some stations shared logistics with exploration companies by co-location at their camps. All HuBLE-UK stations have now been removed and, with the exception of a few vaults left at airports for potential future re-occupation, no trace of the stations' existence remains at the sites.

HuBLE: the salient results
This section reviews the findings of the major phases of HuBLE. The 'Crustal structure' section presents the crustal studies of Thompson et al. (2010) and Pawlak et al. (2011Pawlak et al. ( , 2012, which adopted receiver function and ambient noise tomographic methods to constrain fundamental parameters such as crustal thickness and seismic wavespeed structure across the Bay region. The 'Seismicity in northern Hudson Bay' section looks at seismicity in northern Hudson Bay and its implications for the state of crustal stress in the region. The section entitled 'Mantle seismic anisotropy: evidence from SKS splitting' reviews the seismic anisotropy studies of Bastow et al. (2011b) and Snyder et al. (2013), who performed shear wave splitting of SKS phases to study mantle anisotropy. The 'Surface wave tomography' section reviews various regional surface wave studies conducted across the Bay region (e.g. Darbyshire 2005; Darbyshire & Eaton 2010; Darbyshire et al. 2013). Finally, the section entitled 'Mantle transition zone structure' summarizes the work of Thompson et al. (2011), who performed a receiver function study of the mantle transition zone with a view to improving understanding of its thermal structure.

Crustal structure
Receiver function analysis captures P-to S-wave conversions that occur at velocity contrasts in the crust and mantle recorded in the P-wave coda from distant (teleseismic) earthquakes (Fig. 5a;Langston 1979;Helffrich 2006). The abrupt change in wavespeed usually encountered at the Moho makes receiver function analysis a particularly powerful means of studying the properties of the crust. The arrival-times of P to S conversions and subsequent reverberent phases from the Moho can be used to constrain bulk crustal properties: crustal thickness (H ) and V p /V s ratio (k) (Zhu & Kanamori 2000). These parameters can be related subsequently to bulk crustal composition via Poisson's ratio (e.g. Christensen 1996). Thompson et al. (2010) carried out such a study of crustal structure in northern Hudson Bay, providing for the first time detailed constraints on crustal architecture across the boundaries between many of the tectonic subdivisions that comprise the Canadian Shield.
Receiver functions from within the Palaeoarchaean Rae domain reveal remarkably simple crust, with high-amplitude, impulsive Moho conversions and reverberations (Fig. 5b). The crust is also seismically transparent, with little evidence for internal architecture. The Rae domain has the thinnest crust (c. 37 km; Fig. 6a) and the lowest V p /V s ratios ( Fig. 6b; ≤1.73) in the region. More northerly stations display slightly thicker crust, up to c. 42 km. In contrast to the Rae, the crust and Moho of the Hearne domain crust is more complex (thus making Hk analysis more challenging), with evidence for internal architecture (Fig. 5c). Crustal thickness within the Hearne is c. 38 km and the mean V p /V s ratio is c. 1.76 (Fig. 6b). Uncertainties in the Thompson et al. (2010) study are of the order +2 km for crustal thickness, and +0.03 for V p /V s ratio.
Within the Palaeoproterozoic Quebec-Baffin segment of the THO, bulk crustal properties are more variable than in the neighouring Archaean Hearne and Rae domains. Higher V p /V s ratios (.1.75) are found around the Hudson Strait than further north within Baffin Island and further west towards NW Hudson Bay (c. 1.73: Fig. 6b). The thickest crust (c. 43 km) underlies stations in central and southern Baffin Island (Fig. 6a). Thompson et al. (2010) also found evidence, particularly in the northern Rae domain, for a Hales discontinuity in the shallow lithospheric mantle at c. 60-90 km depth (Hales 1969).
Ambient-noise tomography uses the crosscorrelation of diffuse wavefields (e.g. ambient noise, scattered coda waves) to estimate the Green's function between pairs of seismograph stations (e.g. Shapiro et al. 2005). This technique is a popular tool for crustal studies and was particularly useful in the HuBLE project to glean information on crustal structure beneath Hudson Bay. Pawlak et al. (2011Pawlak et al. ( , 2012 used 21 months of continuous recording at 37 stations around Hudson Bay to measure dispersion characteristics of fundamentalmode Rayleigh wave group velocity. The signals extracted contain preferred azimuths that are indicative of stationary coastal noise sources near southern Alaska and Labrador. Tomographic methods are then used to reconstruct a velocity model that is best resolved in areas of dense, crossing path coverage.
An important feature of the tomographic models is a prominent low-velocity region beneath Hudson Bay (Fig. 7). At mid-crustal depths (i.e. longer Tr an s-Hu ds on Or og en Fig. 6. Crustal thickness (a) and V p /V s ratio (b) in the Hudson Bay region determined using the method of Zhu & Kanamori (2000). Black regions are 2.6-2.7 Ga greenstone belts. BaS, Baffin Suture; BeS, Bergeron Suture; SRS, Soper River Suture. Modified after Thompson et al. (2010).  (Hanne et al. 2004). The contrast between slow wavespeeds beneath the Trans Hudson Orogen and the neighbouring Archaean Superior craton persists to at least a 40 s period (see Pawlak et al. 2011, figs 11 & 12). Modified after Pawlak et al. (2011Pawlak et al. ( , 2012. periods than shown in Fig. 7), the Pawlak et al. (2011) study reveals that fundamental-mode Rayleigh velocity within the Superior craton (3.18 + 0.03 km s 21 ) is significantly greater than the velocity within the Trans Hudson Orogen beneath Hudson Bay (3.10 + 0.03 km s 21 ). Rayleigh-wave anisotropy inferred from azimuthal analysis of ambient noise (Pawlak et al. 2012) reveals an arcuate pattern of fast directions interpreted to be indicative of the double-indenter geometry of the Superior craton. At most periods, their results suggest a significant change in anisotropic direction across the inferred primary suture beneath Hudson Bay. The region of lowest velocity beneath Hudson Bay (Fig. 7) corresponds remarkably well with the pattern of lithospheric stretching proposed by Hanne et al. (2004): c. 3 km of crustal thinning.

Seismicity in northern Hudson Bay
Northern Hudson Bay is a region of moderate intraplate seismicity, mainly within the Boothia Uplift -Bell Arch structure (Basham et al. 1977), where earthquakes up to Mw ¼ 6.2 have been documented. Using data from the HuBLE network, focal mechanisms for five moderate earthquakes in the time period 2007-2009 have been determined by fitting surface waveforms (Fig. 8: Steffen et al. 2012). These small earthquakes, of depth 3 -17 km, have thrust-fault mechanisms, consistent with previously determined focal mechanisms for this region as well as model predictions for glacial isostatic adjustment. Steffen et al. (2012) performed a stress inversion using available focal mechanisms and found a maximum compressive stress direction, SH max , oriented approximately NNW-SSE. This orientation is surprising, as it is rotated by about 908 from previous glacial isostatic adjustment model estimates, which assume a background stress field in North American with SH max oriented NE -SW. These results indicate that existing crustal fault zones exert a strong influence on the local stress field.

Mantle seismic anisotropy: evidence from SKS splitting
Seismic anisotropy -the directional dependence of seismic wavespeed -can be measured from the waveforms of teleseismic earthquakes via analysis of shear-wave splitting (e.g. Silver & Chan 1991). When a horizontally polarized shear wave, such as SKS, enters an anisotropic medium it will split into two orthogonally polarized waves. Splitting can be quantified by the time delay (dt) between the two shear waves, and the orientation (w) of the fast shear wave. These splitting parameters can subsequently be used to understand the preferential alignment of minerals in the crust and/or mantle, or the preferential alignment of fluid or melt (e.g. Blackman & Kendall 1997). Many processes can lead to such anisotropy, including flow of the asthenosphere parallel to absolute plate motion (e.g. Bokelmann & Silver 2002;Assumpção et al. 2006), mantle flow around deep continental keels and slabs (e.g. Fouch et al. 2000;Di Leo et al. 2012), and a pre-existing fossil anisotropy frozen in the lithosphere (e.g. Bastow et al. 2007). Bastow et al. (2011b) and Snyder et al. (2013) used shear wave splitting of SKS phases at HuBLE and POLARIS stations in the Hudson Bay region to infer patterns of seismic anisotropy in the Laurentian crust and mantle (Fig. 9).
In the northern part of Hudson Bay, Bastow et al. (2011b) found no significant variation in splitting parameters at most stations across the HuBLE network. One exception was permanent station FRB on Baffin Island, for which more than a decade of data were available. Here, significant backazimuthal variation in w and dt was found, raising the possibility that complex patterns of anisotropy exist beneath the region (Fig. 10). In a parallel analysis of SKS splitting, Snyder et al. (2013) proposed that two layers of anisotropy exist beneath the Bay, one paralleling near-surface tectonic trends and an underlying fabric paralleling presentday plate motion.

Surface wave tomography
Studies of fundamental-mode Rayleigh wave dispersion were carried out for two-station paths across the Hudson Bay region by . The dispersion curves clearly indicated a thick, fast lithospheric keel beneath the region, with phase velocities significantly greater than those associated with global reference models. The phase velocities were also systematically higher than the Canadian Shield average (Brune & Dorman 1963), corroborating global tomographic models that place the centre of the high-wavespeed lithospheric keel beneath Hudson Bay (e.g. Nettles & Dziewonski 2008). Each two-station dispersion curve was inverted for a 1D shear wavespeed profile representing the average structure along the inter-station path. The period ranges recovered from the dispersion analysis allowed for reliable models from lower-crustal depths (c. 35-40 km) to mantle depths of c. 300 km. Each path-averaged model showed a prominent high-wavespeed anomaly, interpreted as the lithospheric 'lid' in the upper mantle.
Different proxies for lithospheric thickness exist in the literature because fundamental-mode ig. 9. Shear wave splitting parameters (arrows) and null results (solid lines) from HuBLE-UK and neighbouring POLARIS stations (triangles) in northern Hudson Bay from the study of Bastow et al. (2011b). STZ, Snowbird Tectonic Zone. Solid black lines are sutures. B.S., Baffin Suture. Inset: Back-azimuth and distance distribution of earthquakes used in study. Concentric circles indicate 308 intervals from centre of network at 758W, 638N. APM-absolute plate motion from the HS3-Nuvel-1A model (Gripp & Gordon 2002) in both hotspot reference frame (white APM arrow) and no-net rotation reference frame (black APM arrow). surface wave models are essentially insensitive to first-order discontinuities (see e.g. Eaton et al. 2009). In some studies, negative gradients in the models are used, although these can be difficult to constrain. Other authors choose (somewhat arbitrarily) a contour of positive V s anomaly with respect to the global reference (e.g. PREM, ak135). The 1.5-2.0% anomaly range is commonly used (see  for a more detailed discussion). Using a proxy based on the depth to the 1.7% fast anomaly, the depth to the base of the lithosphere was interpreted to vary from c. 190 km beneath the Hudson Strait area to at least 240 -280 km beneath central Hudson Bay. Within the lithospheric lid, maximum seismic wavespeed anomalies varied from c. 4 to 7%, which is fast relative to the global reference (Fig. 11).
The lithosphere-asthenosphere boundary depth variations inferred from the surface-wave analysis, and from interpretation of previous studies in the Canadian Shield, indicate that the lithosphere is thickest beneath central Hudson Bay (and significantly thicker than in the Superior craton to the south). There appeared to be no spatial correlation between lithospheric thickness, path-averaged seismic wavespeed anomalies and surface ages (Archaean v. Proterozoic). This result is similar to interpretations of the Fennoscandian (Bruneton et al. 2004) and Australian (Simons et al. 1999;Fishwick et al. 2005) lithosphere, but contrasts with systematic changes in lithospheric thickness between Archaean and Proterozoic domains in central Asia (Lebedev et al. 2009) and southern Africa (Li & Burke 2006).
The dispersion curves from a total of 172 interstation paths across the Hudson Bay region were combined in a tomographic inversion by Darbyshire et al. (2013), solving simultaneously for isotropic phase velocity heterogeneity and azimuthal anisotropy. The resulting phase velocity maps were subsequently used as the basis for a full 3D model of shear wave velocity and anisotropy beneath the region (Fig. 12). The model confirmed the earlier findings of a thick, fast lithospheric keel, but highlighted internal velocity variations and stratification of seismic anisotropy within the cratonic lithosphere.
In the uppermost mantle, probably associated with the intracratonic basin, seismic velocities immediately beneath the Bay are relatively low compared with those beneath the surrounding Archaean landmass, and anisotropic fast directions wrap around the Bay. In contrast, in the 70 -160 km depth range, the model shows two highvelocity cores beneath the central Superior and Churchill cratons, separated by a near-vertical curtain of lower-velocity material. In this depth range, no large-scale coherency in anisotropic pattern is evident. In contrast, the model also images a basal lithospheric layer with a significantly more homogeneous velocity structure than the mid-lithosphere. Anisotropy within this layer is coherent, with a pattern strongly reminiscent of the inferred geometry of the THO structure (Darbyshire et al. 2013).

Mantle transition zone structure
The seismic discontinuities observed at depths of c. 410 and 660 km define the mantle transition zone (TZ), and are commonly attributed to phase transitions in the olivine system (olivine to wadsleyite, the '410', and ringwoodite to perovskite + magnesiowüstite, the '660', respectively; Bina & Wood 1987;Helffrich 2000;Ita & Stixrude 1992). Described by their Clapeyron slopes, the '410' and '660' will vary in depth if temperature conditions are varied. The opposite Clapeyron slopes of these discontinuities result in a thickened transition zone in cold regions; hot regions will act to thin it (e.g. Helffrich 2000). By analysing the nature of transition zone seismic discontinuities using receiver function analysis, Thompson et al. (2011) explored the thermal and chemical state of the mantle beneath the Hudson Bay region. Thompson et al. (2011) showed that, beneath a significant portion of cratonic North America, the TZ seismic discontinuities are unperturbed compared with the global mean, with the implication that there is no seismically detectable thermal variation at TZ depths beneath one of the deepest and most laterally extensive continental keels on the planet (Fig. 13).

A new relative arrival-time upper mantle tomographic model for northern Hudson Bay
Overview Using teleseismic data recorded by HuBLE and surrounding POLARIS seismograph stations since 2007 a new study of P-wave mantle velocity structure using least-squares tomographic inversion of relative arrival-times is presented using the method of VanDecar et al. (1995). Body wave tomographic methods such as this benefit from particularly good lateral resolution of velocity anomalies. The study thus presents an excellent opportunity to examine whether variations in seismic wavespeed exist between Phanerozoic and Archaean mantle juxtaposed during the THO. This study, which samples only Precambrian geology, is particularly valuable because most studies of cratons using the method include data from stations deployed within younger (Phanerozoic) terranes. In Tanzania, for example, the high wavespeed Tanzania craton contrasts with the neighbouring East African Rift (Ritsema et al. 1998

Method of relative arrival-time determination
Manual picking of the first arriving P-wave phase identifiable across the seismograph network was performed on waveforms that were filtered with a zero-phase two-pole Butterworth filter with corner frequencies 0.4-2 Hz. Subsequently, phase arrivals and relative arrival-time residuals were more accurately determined using the multi-channel crosscorrelation technique (MCCC) of VanDecar & Crosson (1990); we selected a 3 s window containing the initial phase arrival and typically one or two cycles of P-wave energy to cross-correlate. The chosen bandwidth is similar to that used in other teleseismic tomographic studies in both tectonically active (e.g. Bastow et al. 2008) and shield (e.g. Frederiksen et al. 2007) regions. The filter bandwidth is designed to retain as high a frequency as possible since our inversion procedure adopts ray theory (the infinite frequency approximation). All waveforms with cross-correlation coefficients ,0.85 were eliminated from the analysis. The MCCC method also provides a means of quantifying the standard error associated with each arrival time. In this study relative arrival times determined in this way have mean standard deviation of 0.02 s. In line with other studies using the MCCC method (Tilmann et al. 2001;Bastow et al. 2005), we regard the MCCC-derived estimates of timing uncertainty as optimistic.
Relative arrival-time residuals t RES for each station are given by: where t i is the relative arrival time for each station i; t ei is the expected travel time based on the IASP91 travel time tables (Kennett & Engdahl 1991) for the ith station; and t e is the mean of the IASP91 predicted travel times associated with that particular event. The final travel time dataset comprises 3682 P-wave travel times (Fig. 14).

Analysis of travel-time residuals
International Seismological Catalogue travel-time data for permanent station FRB (Fig. 15) in Iqaluit, Baffin Island (Fig. 3) show that the mean absolute delay time with respect to the IASP91 travel-time tables for P-wave arrivals is amongst the earliest of all permanent stations; mantle seismic wave-speeds beneath FRB are therefore amongst the fastest worldwide (Poupinet 1979). All P-wave velocity anomalies shown in this tomographic study should thus be considered markedly fast compared with normal mantle, with red low-velocity regions slower and blue high-velocity regions faster than the background mean of the shield, not the global average.

Model parameterization and inversion procedure
Upper mantle wavespeed structure is imaged using regularized, least-squares inversion of Canadian relative arrival-time residuals following the method of VanDecar et al. (1995). The parameterization  scheme consists of 23 knots in depth between 0 and 1800 km, 50 knots in latitude between 50 and 788N and 85 knots in longitude between 46 and 1148W: a total of 97 750 knots parameterizing slowness. Knot spacing is 35 km in the innermost resolvable parts (58-718N, 63-1008W, 0-350 km depth). Outside this region, knot spacing increases to 50 km between 350 and 500 km, to 100 km between 500 -1000 km, and then to 200 km at 1800 km depth. Structural interpretations are thus limited to features of spatial wavelength ≥70 km in the upper part of the model. The parameterization scheme continues outside the area of interest so that an unwarranted and spurious structure is not mapped into the region where structural interpretations will be made: an Occam's Razor approach (VanDecar et al. 1995). In the regularized least-squares inversion procedure, slowness perturbations, source terms and station terms are determined simultaneously (VanDecar et al. 1995). The source terms are free parameters used in the inversion procedure to account for small variations in back-azimuth and incidence angle caused by distant heterogeneities and source mis-locations. The station terms account for traveltime anomalies associated with the region directly beneath the station where the lack of crossing rays prevents the resolution of crustal structure. Since the inverse problem is under-determined (more unknowns than observations), even in the absence of errors, a unique solution cannot be found. A model is therefore chosen that contains the least amount of structure necessary to satisfy the arrivaltime data (e.g. Constable et al. 1987).
By investigating the trade-off between the root mean square (RMS) residual reduction (the percentage difference between the initial and final RMS misfit to the travel-time equations) and RMS model roughness, a preferred model is selected that fits the data well but does not account for more relative arrival-time residual reduction than can be justified by the a priori estimation of data noise levels. All models in this study account for 95% (from 0.384 to 0.02 s) of the RMS of the relative arrival-time residuals. Estimates of RMS timing uncertainty (0.02 s) are thus treated as optimistic bounds when fitting the data. Subtracting the station terms from the delay times reduces the RMS of the relative arrival-time residuals from 0.384 to 0.369 s; these corrected residuals reflect more accurately the proportion of the delay-time anomalies that will be mapped into the region of the model where tectonic interpretations are drawn.

Resolution
The resolving power of the inversion technique in this study is assessed by analysing the ability of  Fig. 15. P-Wave travel-time delays for permanent seismograph stations. Note the extremely early arrivals at station FRB in Iqaluit. Modified after Poupinet (1979). the ray geometry to retrieve a checkerboard model using raypaths through a 1-D Earth. In the checkerboard test, positive and negative slowness-anomaly (V p ¼ +5%) spheres described by Gaussian functions across their diameter are placed in two layers at 175 and 400 km depth (Fig. 16a, b). The checkerboard approach permits assessment of model sensitivity by highlighting areas of good ray coverage, and the extent to which smearing of anomalies occurs. These synthetic velocity structures are inverted to use identical model parameterization and inversion regularization as during the inversion of the observed data. A Gaussian residual time error component with a standard deviation of 0.02 s is added to the theoretical P-wave travel times (the standard deviation of noise estimated for the observed data). Lateral resolution of the anomalies in Figure 16b and c is good, with c. 30% amplitude resolution not uncommon in regularized, under-determined tomographic inversions such as this. Vertical resolution is good in some parts of the model (e.g. beneath the central Rae domain; Fig. 16e), but moderate in others. For example, vertical smearing in areas such as southern Baffin Island (Fig. 16f ) is considerable, meaning that the model cannot be used to place tight depth constraints on wavespeed variations.

The HuBLE P-wave velocity model
Cross sections and depth slices through the new P-wave velocity model are shown in Figure 17. Depth sections at lithospheric depths (Fig. 17a, b) indicate relatively little heterogeneity, with peakto-peak amplitudes of ,1%. This compares with peak-to-peak anomalies of c. 3% in active areas such as the Ethiopian Rift (e.g. Bastow et al. 2005). Across the Snowbird Tectonic Zone (Rae-Hearne boundary), for example, there is almost no discernible variation in seismic wavespeed (Fig.  17a, b, e), although it is acknowledged that this could be partly the result of poor resolution at the edge of the model (Fig. 16b).
One area where relatively strong variations in seismic wavespeed can be observed is the Quebec -Baffin segment of the THO (Fig. 17a, b, e). Here, fast wavespeeds associated with the southern edge of the Churchill plate (i.e. Baffin Island) contrast with slower wavespeeds beneath the Baffin Strait, with a subvertical to south-north-dipping boundary between the two (Fig. 17e). This region corresponds to a region of active seismicity and Palaeozoic faulting (e.g. Steffen et al. 2012). The depth extent of the low-wavespeed anomaly beneath the Baffin Strait (Fig. 17e) is c. 300 km, implying that it may extend through most or all of the lithosphere (allowing for a finite amount of vertical smearing in the model; Fig. 16). It must be noted, however, that resolution in this part of the velocity model is moderate: amplitude recovery is low and vertical smearing is clearly evident in the resolution tests (Fig. 16).

Causes of seismic heterogeneity in tomographic models
A number of factors can affect seismic velocity in the mantle, including temperature, partial melt and composition (Karato 1993;Goes et al. 2000;Goes & van der Lee 2002). Variations in seismic anisotropy may be manifest as velocity heterogeneities. However, temperature is often cited as the main source of mantle heterogeneity in the upper mantle (e.g. Goes & van der Lee 2002).
The markedly fast wavespeeds observed in global and continent-scale tomographic images of the Canadian Shield mantle indicate strongly that the lithosphere beneath the Hudson Bay region is cold relative to the global average (e.g. Darbyshire et al. 2007;Nettles & Dziewonski 2008;. Composition is also an important contributing factor, however: shields are generally associated with cold, depleted, mechanically strong mantle material characterizing their tectospheric 'roots' (e.g. Jordan 1978). It is often commonly assumed that the Archaean cores of the continents are characterized by the seismically fastest and thickest lithosphere on Earth, with thinner and slower wavespeed structure characterizing Proterozoic regions (Durrheim & Mooney 1994). More recent seismic studies of mantle structure in Australia (Simons et al. 1999;Fishwick et al. 2005;Fishwick & Reading 2008), Canada  and Fennoscandia (e.g. Bruneton et al. 2004) question this simple correlation between surface geology age and underlying mantle structure and we revisit this issue in the 'Implications for the formation of the Laurentian keel' section.
As indicated in the 'Analysis of travel-time residuals' section, when interpreting tomographic images derived from relative arrival-time data (Fig.  17), it is essential to appreciate the background velocity structure of the region, which is represented in the tomographic images by the dV p ¼ 0% contour. Studies of Precambrian lithosphere commonly include data from Phanerozoic terranes as well (e.g. Tanzania: Ritsema et al. 1998;SE Canada: Villemaire et al. 2012). These models illuminate high-/ low-velocity anomalies that are sometimes genuinely fast/slow compared with the global mean. However, our study area's mean wavespeed structure is all fast compared with the global average (Fig. 15), with no stations lying off-shield. Both  'low'-and 'high'-velocity anomalies presented in Figure 17 are thus fast compared with the global mean.
The c. 1.8 Ga elapsed since the last major tectonic event (the THO) means that lithospheric thermal anomalies associated with magmatism in the region will have long since dissipated. The velocity anomalies in Figure 17 imaged are thus probably sensitive to variable composition across the network.
Implications for the assembly of the Canadian Shield: evidence for modern-style plate tectonics?
Whether or not modern-style plate tectonics operated throughout Precambrian times is debated: its onset has been estimated to be as early as the Hadean (e.g. Hopkins et al. 2008) or as late as c. 1 Ga (e.g. Stern 2005) -a time span of approximately two-thirds of Earth history (Fig. 1). At crustal depths, the receiver function study of Thompson et al. (2010) revealed that crustal character has a strong age dependence, with the implication that crustal formation processes have evolved over time.
The low bulk crustal V p /V s ratios of the Palaeoarchaean Rae domain crust (Fig. 6b) are indicative of felsic-to-intermediate composition to rocks in the mid to lower crust (e.g. Christensen 1996). Crustal formation hypotheses involving island arc accretion are thus not easily applied to the Rae domain because arcs are believed to have a basaltic bulk composition (Rudnick 1995). Nair et al. (2006) explain low bulk crustal V p /V s ratios in the Kaapvaal craton through delamination of the basaltic lower crust during collision in an island-arc setting, thereby favouring the uniformitarian view that modern plate tectonics can be used to describe crustal formation in Archaean times. In the Rae domain, however, there is a paucity of evidence to support subduction-related hypotheses. For example, extensive linear orogenic trends in surface geology and potential field maps, which are usually associated with island-arc accretion, are not in evidence (e.g. . The absence of greenstone terranes associated with collisional tectonics in the Rae domain (e.g. Hartlaub et al. 2004) also makes accretionary crustal formation models difficult to invoke. The felsic, uniform thickness (c. 37 km; Fig. 6a) Rae domain crust, with its remarkably sharp Moho and paucity of evidence for intra-crustal reflectivity (Fig. 5b) thus lacks evidence for modern-style plate tectonics. Models favouring vertical tectonic processes, such as crustal delamination or plume activity, were thus considered better suited to the observations by Thompson et al. (2010).
The Meso-to-Neoarchaean Hearne domain has a slightly more complex crust (Fig. 5c) with higher V p /V s ratios than the Rae. These properties were interpreted by Thompson et al. (2010) as a potentially transitional period between non-plate tectonic and plate tectonic processes. Receiver function results from the Quebec -Baffin segment of the THO, however, map out the first-order shape of the Superior plate under-thrusting Meta Incognita, the southern tip of the Churchill collage ( Fig. 3; Thompson et al. 2010). The elevated V p /V s ratios in this region (Fig. 6b) are considered to be representative of the rifted margin of the Superior craton. Markedly thicker crust (c. 43 km; Fig. 6a) compared with that observed in the Archaean domains is coincident with widespread medium to high-grade metamorphic geology outcropping at the surface (e.g. St-Onge et al. 2006) and can be explained by crustal thickening owing to stacking of accreted terranes during continent-continent collision, analogous to the present-day Tibetan Plateau, followed by erosion (Thompson et al. 2010), consistent with the hypothesis of St-Onge et al. (2006).
From SKS splitting analysis of seismic anisotropy, Bastow et al. (2011b) and Snyder et al. (2013, with support from receiver functions) showed that fast polarization directions in northern Hudson Bay parallel THO structural trends at the surface (Fig. 9). Relatively large splitting delay times (dt up to 1.45 s) require a relatively thick anisotropic layer (and thus major plate-scale deformation), which Bastow et al. (2011b) cited as evidence that plate tectonics was in operation by Palaeoproterozoic times at the latest. Corroborating this conclusion, Snyder et al. (2013) presented evidence for gently dipping lithospheric layers deep beneath the SE Rae craton, consistent with under-thrusting beneath the Hearne and Meta Incognita Sugluk blocks (Fig. 3). The new tomographic images presented in this study also yield evidence for plate-scale under-thrusting beneath Meta Incognita (Fig. 17a, b, f ). In contrast there is little evidence in our new wavespeed maps for such large-scale tectonic processes across the Rae -Hearne domains to the west of Baffin Island (Fig. 17a, b, g).
Taken together, the results from the receiver function, SKS and tomography studies support a model for secular evolution in processes of crustal formation formation, with non-plate tectonic processes during the Palaeoarchaean evolving towards fully developed plate tectonics by the Palaeoproterozoic (Thompson et al. 2010).
Based on joint tomographic inversion of ambient-noise and teleseismic Rayleigh-wave dispersion measurements, Pawlak et al. (2012) documented systematic variations in anisotropic fabric moving from the upper crust, to the lower crust and uppermost mantle beneath Hudson Bay. In particular, while anisotropic fabrics in the upper crust and uppermost mantle are consistent with tectonic fabrics created during the THO, lower-crustal fabrics exhibit a more uniform north -south orientation. By analogy with models for crustal evolution in the modern Himalayan orogen, Pawlak et al. (2012) interpreted this pattern as evidence for frozen channel flow in the lower crust that formed during the culmination of the orogenic cycle. This model furnishes additional evidence that modernstyle plate tectonics was in operation c. 1.8 Ga during the THO.

Implications for the formation of the Laurentian keel
The tectosphere or lithospheric mantle beneath cratons commonly extends to depths of 250 km or more into the mantle. These 'keels' have remained stable, resistant to thermal and mechanical erosion during multiple Wilson cycles over billions of years. Archaean processes, such as the extraction of komatiitic magmas (e.g. Griffin et al. 2003), are often invoked to explain the intrinsic low density, but the Laurentian keel does not fit easily into the Archaean formation paradigm: it extends below both Archaean and Palaeoproterozoic terranes. Precisely how the keel beneath Hudson Bay formed is thus unclear. Is the present-day keel a relict of initial craton formation in the Archaean (e.g. Sleep 2003), or do later tectonic processes (specifically the THO) also play a role, perhaps resulting in a vertically stratified lithosphere? Stratification of cratonic lithosphere has been inferred previously from both geophysical (Fishwick & Reading 2008;Angus et al. 2009;Abt et al. 2010; and geochemical/thermobarometric studies (Griffin et al. 2003(Griffin et al. , 2004. The upper lithospheric layer has been interpreted as a highly meltdepleted peridotite chemical boundary layer; the lower, more fertile layer has been described as a thermal boundary layer, formed at a later stage in keel evolution (e.g. Lee et al. 2011). S to P converted phases isolated in receiver function studies, initially interpreted as the LAB (Rychert & Shearer 2009), are strong evidence for the latter hypothesis, and have been cited as evidence for midlithospheric discontinuities beneath the Canadian shield (Abt et al. 2010;. The new tomographic images presented in this study (Fig. 17) provide compelling evidence that Palaeoproterozoic mantle associated with the THO has been trapped between the Archaean cores of the Superior plate in the south and the Churchill collage (including Meta Incognita and the Sugluk block: Fig. 3) in the north. Depth resolution in such body wave models is never as good as lateral resolution (Fig. 16), but the inferred trapped lowwavespeed material between the Superior craton and the Churchill collage in Figure 17e is certainly a deep-seated lithospheric feature, not the smeared result of crustal structure.
The surface wave study of Darbyshire et al. (2013), which solved simultaneously for isotropic phase velocity heterogeneity and azimuthal anisotropy, also suggested a multi-stage assembly of the central Laurentian keel in the Archaean and Palaeoproterozoic. The Darbyshire et al. (2013) model (Fig. 12) shows that, in the uppermost mantle, probably associated with the development of the intracratonic basin, seismic wavespeeds immediately beneath Hudson Bay are lower than beneath the surrounding Archaean blocks, and anisotropic fast directions wrap around the Bay. In contrast, in the 70 -160 km depth range, the model shows two fast-wavespeed cores beneath the central Superior and Churchill cratons, separated by a near-vertical corridor of lower-wavespeed material. The geometry of this heterogeneous pattern prompted Darbyshire et al. (2013) to suggest that it arises from the terminal collision of the Superior and Churchill cratons, trapping more juvenile Proterozoic mantle between them during the THO. In this depth range, no large-scale coherency in anisotropic pattern is evident. However, in the deepest parts of the lithosphere the model reveals significantly more homogeneous velocity structure than the mid-lithosphere. Anisotropy within this layer is coherent, with a pattern strongly reminiscent of the inferred geometry of the THO suture. Given this distinctive deformation pattern, the basal layer is postulated to have accreted to the base of the Archaean cratonic lithosphere during or shortly after the terminal phases of the THO (Darbyshire et al. 2013).

Neotectonics, and implications for the development of the Hudson Bay basin
The Hudson Bay basin is the least studied of four major Phanerozoic intracratonic basins (including Williston, Illinois and Michigan) in North America and the mechanism by which it formed remains ambiguous. A number of hypotheses have been proposed over the years, each of which the HuBLE project has tested.
According to one hypothesis, subsidence occurred as a result of convective down-welling within the mantle (e.g. James 1992; Peltier et al. 1992). This explanation has also been championed to explain the long-wavelength negative gravity anomaly beneath Hudson Bay (Simons & Hager 1997;Mitrovica 1997;Tamisiea et al. 2007). The receiver function study of Thompson et al. (2011) addressed the issue of whether or not a mantle down-welling exists below the Bay by constraining the thermal structure of the mantle at transition zone depths. Elevated mantle temperatures are expected to reduce the depth interval between the olivineto-wadsleyite (the '410'), and ringwoodite-toperovskite + magnesiowüstite (the '660') phase transitions in the olivine system; reduced temperatures in the same depth range will increase transition zone thickness (e.g. Bina & Wood, 1987;Ita & Stixrude 1992;Helffrich 2000). The Thompson et al. (2011) study showed, however, that beneath the majority of cratonic North America, mantle transition zone seismic discontinuities are unperturbed with respect to the global mean, with the implication that there is little thermal variation at TZ depths beneath the Laurentian keel (Fig. 13). This suggests that any mantle down-welling associated with the long-wavelength negative gravity anomaly, and the Hudson Bay basin itself, must be confined to upper mantle depths. Small-scale convection around the Laurentian keel, if it operates at all, must also be confined to upper mantle depths of ,410 km (Thompson et al. 2011).
Two hypotheses for Hudson Bay basin formation concern the topography generated solely by crustal processes. In the first of these, basin subsidence was triggered by eclogite phase transformation within an orogenic crustal root (Fowler & Nisbet 1985;; according to the second, basin subsidence occurred in response to lithospheric extension that resulted in crustal thinning (Hanne et al. 2004). These hypotheses make different, testable, predictions about crustal thickness trends. The receiver function study of Thompson et al. (2010) was unable to constrain crustal structure beneath the Bay itself, but notably yielded no evidence for an eclogitized lowermost crust beneath its islands and coastal areas: such an assemblage around the Moho would result in gradational, not sharp, P to S conversions from the Moho, as are observed in Figure 5b & c (compare, for example, to results from the China craton presented by Zheng et al. 2008).
The study of ambient noise tomography by Pawlak et al. (2011) did, however, place new constraints on crustal structure beneath Hudson Bay. Tomographic maps and cross-sections obtained in the 5-40 s period range (Fig. 7) reveal markedly lower velocities at crustal depths beneath Hudson Bay than in the surrounding Archaean cratons (3.10 + 0.03 km s 21 within Hudson Bay v. 3.18 + 0.03 km s 21 in the Superior craton). The lowest mid-crustal velocities correspond to the region of maximum plate stretching near the centre of the basin, as inferred by Hanne et al. (2004). Pawlak et al. (2011) present the first compelling direct evidence (the upwarp of mantle material) for crustal thinning beneath the Bay, obviating the need for eclogitization of a remnant lower crustal root as a mechanism for basin formation.
Analysing intra-plate seismicity in the northern part of Hudson Bay (Fig. 8), Steffen et al. (2012) noted that waveforms from these small earthquakes are consistent with thrust-fault mechanisms, as inferred from previous studies for the region, as well as model predictions for glacial isostatic adjustment. Present-day crustal stresses are thus influenced by both glacial rebound and pre-existing faults. Hudson Bay itself is an ephemeral feature (c. 15 ka), caused by combined incomplete glacial rebound, mantle flow-driven dynamic topography and the thermochemical structure of the Laurentian keel. The seismological results from HuBLE, when synthesized in light of the sedimentary basin study of Hanne et al. (2004), and more recent backstripping studies of the Hudson Bay basin (Pinet et al. 2013), suggest strongly that plate stretching is the primary mechanism responsible for the formation of the Hudson Bay basin.

Conclusions
The Hudson Bay lithospheric experiment has placed fundamental new constraints on the crust and upper mantle seismic structure of Laurentia. Highlights of the experiment (so far) include: † The strong age-dependence on crustal architecture and bulk crustal properties (crustal thickness and V p /V s ratio) indicates a secular evolution in processes of crustal formation formation, with non-plate tectonic processes during the Palaeoarchaean evolving towards fully developed plate tectonics by the Palaeoproterozoic. † Seismic anisotropy, revealed by SKS splitting and surface wave analysis, indicates the preservation of a fossil lithospheric fabric, interpreted as evidence for modern-day-style plate tectonics operating by Palaeoproterozoic times. † Seismic tomographic images of mantle structure reveal an age dependence of seismic wavespeed: the Archaean cores of Laurentia exhibit faster seismic velocities than Proterozoic material within the heart of Laurentia, trapped since the THO. State-of-the-art anisotropic surface-wave inversions indicate depth-dependent anisotropy. Taken together with the variable wavespeeds across the shield, and evidence for mid-lithospheric discontinuities in S-P receiver function studies, this implies that cratonic keel formation is not confined to Archaean times. † Deep beneath the Laurentian lithosphere, a mantle transition zone study shows that the olivine phase transitions from olivine to wadsleyite (the '410') and from ringwoodite to perovskite + magnesiowüstite (the '660') are at remarkably normal levels across the Laurentian keel's c. 3500 km lateral extent. This implies that the keel has no significant thermal effect on the underlying mantle (≤50 K), and any small-scale convection or cold mantle downwelling associated with the large free-air gravity anomaly beneath the Laurentian shield must be confined to the upper mantle above the transition zone (i.e. shallower than 410 km). † Evidence for crustal thinning beneath the Bay, inferred from ambient noise tomography, indicates that the Hudson Bay basin probably owes its existence to shallow Phanerozoic tectonic processes (plate stretching), not deeper mantle down-wellings.