Nature of the Moho transition in NW Canada from combined near-vertical and wide-angle seismic-reflection studies

The crust-mantle boundary, or the Moho, is one of the most distinct manifestations of a differentiated Earth. Major changes in petrology, mineralogy, chemistry, seismic wave velocity, density, and rheology occur across it (Jarchow and Thompson, 1989). Originally, the Moho was defined on the basis of seismic head waves (Mohorovičić, 1910). It was described as the depth at which the compressional P-wave velocity increases rapidly or discontinuously to 7.6–8.6 km/s (Steinhart, 1967). This suggested a petrologic interpretation of the Moho as a first-order (i.e., zero thickness) discontinuity in rock composition from predominantly mafic lower-crustal rocks to predominantly ultramafic upper-mantle rocks. Ophiolite studies indicate that the lithologic sequence of the oceanic Moho includes pelagic sediments, pillow basalts, massive and layered gabbro, and residual ultramafic tectonites of the upper mantle (Casey et al., 1981; Karson et al., 1984; Boudier and Nicolas, 1995; Jousselin and Nicolas, 2000; Dilek et al., 2008). These rock cumulates form laterally discontinuous lenses with thicknesses ranging from 1 cm to several tens of meters and aspect ratios ranging from 10:1 to 100:1. Prominent reflections from within the oceanic Moho transition zone obtained by seismic-reflection methods also support a geologically complex boundary (e.g., Hasselgren and Clowes, 1995; Nedimović et al., 2005).

Unlike abundant ophiolites, which allow direct observation of the oceanic Moho, exposed sections of the continental Moho are rare (Fountain and Salisbury, 1981). Nevertheless, important constraints on the nature of the continental Moho are provided by the small continental data set of xenoliths (McGetchin and Silver, 1972; Debari et al., 1987; Kopylova et al., 1998). Xenolith samples indicate a lower crust composed of mafic granulite, gabbro, and amphibolite and an upper mantle composed of several varieties of peridotite with some eclogite.

Despite all the geologic studies of the Moho by means of exposed sections (ophiolites) or xenolith samples, seismic data still remain the main source for our understanding of the structure of the Moho and its development. As techniques applied to lithospheric studies have improved considerably in terms of seismic data acquisition technologies (e.g., Meissner, 1986; Mooney, 1987; Mjelde et al., 1997; Panea et al., 2005) and computational capabilities (e.g., Fuchs and Müller, 1971; Spence et al., 1984; Levander and Holliger, 1992; Zelt and Smith, 1992; Carbonell et al., 2002), it has become apparent that the Moho is typically not a uniform, sharp, first-order discontinuity, but rather a complex and variable transition zone. From the seismic perspective, a frequently cited model for the Moho is characterized by a variable-thickness transition zone composed of anastomosing layers progressing with depth from mafic to ultramafic rocks (Jarchow and Thompson, 1989; Braile and Chiang, 1986). Alternatively, the Moho transition may represent a change in scale of vertical layering and horizontal extent of inhomogeneities (Tittgemeyer et al., 1999; Hurich, 2003; Carpentier and Roy-Chowdhury, 2007).

Although the Moho is well imaged by near-vertical incidence (NVI) and refraction/wide-angle reflection (R/WAR) data, the exact nature and characteristics of the transition are not well understood. The wave field (i.e., waveform, duration, and amplitude) of a seismic signal (reflected or refracted) from the lower crust and upper mantle may contain significant information about the structure of the Moho transition zone, its fabric, and scale of heterogeneities. Only a few studies have specifically addressed this issue. These studies can be grouped into two categories, qualitative and quantitative. Examples of qualitative analyses of reflection patterns near the Moho are those of Gibbs (1986) on Consortium for Continental Reflection Profiling (COCORP) seismic-reflection profiles, Meissner (2000) using Deutsches Kontinentales Reflexionsseismisches Programm (DEKORP) data, and Hammer and Clowes (1997), Cook (2002), and Eaton (2006) on Lithoprobe reflection transects that sample regions of diverse age and tectonic history. From these studies, reflection patterns appear to be divided into three primary classes: (1) a distinct, multicyclic band of reflectors that are laterally continuous or piece-wise continuous over tens of kilometers, (2) fading out of crustal reflectivity with depth, and (3) no clear reflections in the vicinity of the Moho. These studies did not involve wide-angle data. Although the comparisons have not clearly linked a specific tectonic process with a distinct reflection pattern, sharp Moho reflections have been associated with lower-crustal deformational processes including compression, extension, and intrusive features.

Quantitative studies use synthetic modeling to construct the seismic response of the continental Moho using seismic velocity data (laboratory seismic-refraction results) together with petrologic information derived from exposures of the crust-mantle boundary when available. In earlier studies, the Moho was modeled as a finite-thickness laminated zone of alternating high- and low-velocity layers (Clowes and Kanasewich, 1970; Hale and Thompson, 1982; Fountain, 1986; Braile and Chiang, 1986). More recently, the Moho has been considered as a laterally and vertically heterogeneous transition where the heterogeneities are represented by randomly distributed ellipsoids (Carbonell et al., 2002; Tittgemeyer et al., 2003). Another representation of the Moho follows a statistical approach, in which the transition zone is characterized by a stochastic heterogeneity distribution (Hurich, 2003; Nielsen and Thybo, 2006; Carpentier and Roy-Chowdhury, 2007). Only a few of these studies included consideration of wide-angle data (Carbonell et al., 2002; Tittgemeyer et al., 2003; Nielsen and Thybo, 2006).

In this paper, we analyzed the near-vertical incidence and wide-angle seismic data acquired along line 1 of Lithoprobe's Slave–Northern Cordillera Lithospheric Evolution (SNORCLE) transect in the Paleoproterozoic–Archean domains of the Northwest Territories, Canada (Cook et al., 1999; Fernández Viejo and Clowes, 2003). We investigated the dynamic characteristics of Moho reflections by calculating the synthetic seismic signature for both near-vertical and wide-angle data of a number of crust-mantle transition models using one- and two-dimensional wave propagation algorithms. The models range from a simple first-order discontinuity to more complex, laterally and vertically heterogeneous structures. Unraveling the detailed structure of the crust-mantle boundary and the way in which it forms can help us understand some major tectonic processes, such as crustal growth, accretion, and delamination (Morozov et al., 2001; Carbonell et al., 2002).

The Wopmay orogen in northwest Canada (Fig. 1) is a Paleoproterozoic assembly of domains that developed during the Calderian orogeny between 1.92 Ga and 1.84 Ga (Hoffman, 1989). The orogen is divided into four major north-south–trending domains: the Coronation Supergroup, the Great Bear magmatic arc, the Hottah terrane, and the Fort Simpson terrane. Rifting to the west of the Slave Province at ca. 1.90 Ga initiated deposition of the Coronation Supergroup, a west-facing shelf rise and foredeep succeeding prism, which shortly thereafter was intruded by a suite of plutons and translated eastward by the Hottah collision to form a foreland fold-and-thrust belt (Hoffman, 1988).

The Hottah terrane formed as a magmatic arc and collided with the western Slave craton ca. 1.89–1.88 Ga (Hoffman and Bowring, 1984; Hildebrand and Bowring, 1999). The collision caused compression, shortening, and eastward thrusting of the Coronation Supergroup onto the Slave craton (Hoffman and Bowring, 1984). The lack of coeval arc magmas on the western Slave craton indicates that the accretion of the Hottah terrane onto the craton was the result of west-dipping subduction of an ocean basin. The Hottah terrane consists of Paleoproterozoic sedimentary and intermediate volcanic rocks metamorphosed to amphibolite grade and intruded by calc-alkaline granitic plutons during the period 1.940–1.902 Ga (Hildebrand et al., 1987). To the west, another oceanic plate was converging against the western Hottah terrane, generating an east-dipping subduction zone below the Hottah terrane and Coronation Supergroup. It led to the formation of the Great Bear magmatic arc during the period between 1.875–1.843 Ga (Hildebrand et al., 1987; Gandhi et al., 2001).

Between the Coronation Supergroup and the Hottah terrane, there lies the 100-km-wide Great Bear magmatic arc, which consists of volcanosedimentary sequences and plutonic rocks (Hildebrand et al., 1987). Aeromagnetic data indicate that rocks of the Great Bear arc extend from north to south beneath the younger Proterozoic and Paleozoic cover for a total length of ∼900 km (Coles et al., 1976). At its southern end, the magmatic arc is ∼2.0–4.5 km thick, as interpreted from seismic-reflection data along SNORCLE line 1 (Cook et al., 1999). These data show a strong reflection from the top (∼0.2 s) of the arc, which has been drilled, and another strong reflection at ∼1.5 s from its interpreted base. Whether the thin basin-like feature represents the complete Great Bear section at this latitude or whether there are additional Great Bear rocks below is not known. In either case, development of the arc does not appear to have disrupted the older structure beneath it, at least within the resolution of the reflection survey. One possible explanation for the intact structure underneath is that the arc magmas migrated through the crust along relatively narrow preexisting zones of weakness without disrupting the primary crustal structure (Cook et al., 1999).

The entire Great Bear–Hottah–Slave assembly is deformed by a set of brittle conjugate faults, which have been interpreted as the result of the terminal collision of the Fort Simpson exotic terrane with the western margin of the Hottah terrane (Hoffman and Bowring, 1984). The collision of the Fort Simpson terrane with the Hottah terrane occurred after the formation of Fort Simpson magmatic rocks, ca. 1.845 Ga (Villeneuve et al., 1991), and pre–1.71 Ga mafic intrusions of sedimentary rocks (Abbott et al., 1997). Formation of the Fort Simpson basin, identified as such on SNORCLE line 1 by seismic-reflection and seismic-refraction data (Cook et al., 1999; Welford et al., 2001), was the result of lithospheric extension, which followed the collision of the westernmost terranes with the older Paleoproterozoic domains. From the reflection data, the basin is interpreted as a monocline with a base that dips westward from a depth of ∼3 km at the eastern end to a depth of ∼24 km at the western end, revealing at least 20 km of subsurface relief (Cook et al., 1999).

The enlarged continental assembly was affected by a number of crustal extension phases that continued to Neoproterozoic time (Thompson et al., 1987). The Proterozoic rocks of the Wopmay orogen along SNORCLE line 1 are overlain by Phanerozoic sedimentary rocks that thicken from their feather edge near the eastern boundary of the Great Bear magmatic arc to ∼1000 m over the Nahanni domain (Fig. 1; Cook et al., 1999). The Nahanni domain, west of the Fort Simpson terrane, is the westernmost component of the Wopmay orogen in the study area. It has no surface exposure and is only distinguished by its magnetic signature.

The main purpose of the SNORCLE transect studies was to investigate the various processes involved in the westward growth of North America from the Archean to the present using a multidisciplinary approach (Cook and Erdmer, 2005). Vibroseis near-vertical incidence reflection data were acquired along SNORCLE seismic line 1, one of four such lines in the transect. Line 1 is 725 km long and crosses the southwest Slave craton and Wopmay orogens (Fig. 1). In total, 404 geophones per record with a station spacing of 60 m were used. Shot points were every 90 m with sweep frequencies ranging between 10 and 80 Hz. Complete field acquisition parameters are given in Cook et al. (1999). For such data, the nominal vertical resolution is given by λ/4, where λ is the wavelength, and the horizontal resolution is given by the radius r of the first Fresnel zone ⁠, where z is the depth (Yilmaz, 2001). Thus, at a Moho depth of ∼32 km where the seismic velocity is ∼6.8 km/s (Fernández Viejo and Clowes, 2003) and the dominant frequency determined from Moho reflections is ∼21 Hz, the vertical resolution is ∼80 m and the horizontal resolution is ∼2.2 km.

Refraction/wide-angle reflection data were recorded from 13 explosive shots along the same SNORCLE line 1. The total number of seismographs was 594, with station spacing of 1.2 km. Further details of the field acquisition can be found in Fernández Viejo and Clowes (2003). Interpretation of the R/WAR data followed standard procedures (Zelt and Smith, 1992) and provided information on the two-dimensional (2-D) seismic velocity structure. The estimated uncertainties in velocities based on procedures discussed by Fernández Viejo and Clowes (2003) are given in their Table 3. For the WAR data, traveltime modeling indicated vertical resolution at lower-crustal depths on the order of 1–2 km and horizontal resolution between 60 and 80 km (Fernández Viejo and Clowes, 2003).

The reflection data were reprocessed following standard procedures with parameters (Table 1) that are similar to, but not identical with, those used by Cook et al. (1999). The main difference in our procedure is with handling noise. Incoherent noise present in the data corrupts the quality of the signal and may lead to misinterpretation. Thus, separation of signal and noise is an important step in data processing, particularly in the lower crust where the signal-to-noise ratio (SNR) is low. Traditional seismic processing takes advantage of the redundant (multifold) data to improve SNR. In this study, we adopted an additional and new technique, based on the curvelet transform (Hennenfent and Herrmann, 2006; Neelamani et al., 2008), to suppress incoherent noise and improve the resolution of the seismic data at deeper levels, including the Moho (e.g., Fig. 2).

The curvelet transform belongs to a family of multiscale and multidirectional transforms. A curvelet is localized in frequency and pseudo-localized in space (rapid spatial decay). In the physical domain, it looks like a small plane wave that is oscillatory in one direction and smooth in the perpendicular direction. The signal and noise have minimal overlap in the curvelet domain (Neelamani et al., 2008), making the curvelet transform an ideal choice for our purpose. The underlying assumption of constructing curvelets is that any object with wavefront-like structure (e.g., seismic images) can be represented by a relatively few significant transform coefficients (Candès et al., 2006).

The stacked and denoised section of the seismic image beneath the Great Bear magmatic arc shows a distinct Moho transition at ∼11.0 s from ∼15 to 60 km and a slightly deeper Moho at 11.6 s from 0 to 15 km distance (Fig. 2). The transition consists of a multicyclic band of reflections separating a highly reflective crust from a relatively transparent upper mantle. Ideally, we would like to simulate this Moho transition response by calculation of synthetic seismograms from model representations of the Moho. To simulate the stacked section accurately, calculation of a large number of shot gathers along a high-resolution 2-D velocity model is required. This process is computationally very expensive, and a high-resolution 2-D velocity model is not available. In addition, some of the processing steps preceding the stack may introduce some artifacts to the data.

An alternative approach is to identify Moho reflections on individual shot gathers and then attempt to simulate these through synthetic seismogram modeling. Previously, this has not been attempted because individual shot gathers are characterized by low SNR ratios, and deeper reflections, particularly at the base of the crust, are usually obscured by random noise. We overcame this problem by utilizing the aforementioned curvelet denoising technique to attenuate incoherent noise on the shot gathers. Our results show remarkable reflections at ∼11.0 s that we interpret as the Moho. The reflections are characterized by a narrow band of 5–9 cycles that are laterally variable and piecewise continuous for several kilometers or more (Fig. 3).

Processing of the refraction data involved editing/killing noisy traces and band-pass filtering (3–15 Hz). To offset distances of ∼160 km, the crustal refracted phase, Pg, is the first arrival (Fig. 4). Prominent secondary wide-angle reflections from the crust-mantle boundary (PmP) beneath the Great Bear magmatic arc are observed on three shots (1104, 1105, and 1108). They are characterized by a prominent reverberatory pattern with duration of ∼0.5 s between 80 and 220 km offset (Fig. 4). At offsets greater than 160 km, the refraction phase Pn traveling within the uppermost mantle becomes the first, but rather weak, arrival. Because of the ambiguity of the Pn phase, further amplitude/frequency analysis of this phase was not pursued. A normal-moveout–corrected and stacked section of shot gathers on which the secondary phase PmP was observed showed clear, horizontal reflections at ∼11 s that extend laterally over ∼200 km below the Hottah terrane and Great Bear magmatic arc (Fernández Viejo et al., 1999). These horizontal reflections correspond remarkably well with the horizontal Moho reflections observed on the NVI data.

For amplitude analysis, we calculated the PmP amplitudes within a 0.5 s window length after each traveltime pick for this phase (Fig. 5A). Only amplitudes associated with pick uncertainties <0.15 s were considered (Table 2). All trace amplitudes were initially scaled such that “true relative amplitude” between data of various instrument types was achieved; see Ellis et al. (2002) for details. The amplitude versus offset plots for the three shots show high variability and no obvious distance at which a common amplitude peak is observed.

Both near-vertical and wide-angle reflection seismic data show very distinct Moho reflections from beneath the Great Bear magmatic arc. The dynamic characteristics of these reflections contain valuable information that can help us shed some light on the structure of the crust-mantle transition (Moho). In order to investigate wavefield characteristics (traveltime, duration, amplitude, frequency, etc.), we calculated the synthetic seismic signature of different crust-mantle transition models using 1-D and 2-D wave propagation algorithms. We examined a suite of models ranging from simple to more complex transitions suggested by the many studies mentioned earlier. These include: first-order discontinuities, gradient transitions, layered transitions, and laterally heterogeneous transitions. The key parameters of these models, which control the seismic signature of the synthetic seismograms, are the total thickness of the transition zone, the velocity contrast between the transition zone and the surrounding material, and the correlation length of the lateral heterogeneities. By modifying these parameters, synthetic seismograms that best simulate the real data can be obtained. To examine the seismic signature of the first-order discontinuity, gradient, and layered transitions, we used a 1-D reflectivity modeling algorithm (Fuchs and Müller, 1971). For the laterally heterogeneous models, a 2-D viscoelastic finite difference algorithm (Robertsson et al., 1994) was required. The 1-D and 2-D velocity models that were incorporated into the 1-D reflectivity and 2-D viscoelastic modeling, respectively, were based on the interpretation of the R/WAR study along SNORCLE line 1 (Fernández-Viejo and Clowes, 2003). The crustal P-wave velocities increase from 5.9 km/s at the surface to 6.8 km/s at the crust-mantle boundary (∼32 km depth). The velocity of the uppermost mantle is 8.1 km/s and linearly increases to 8.4 km/s at 45 km depth. Corresponding S-wave velocities were calculated using the ratios Vp/Vs = 1.78 and Vp/Vs = 1.82 for the crust and upper mantle, respectively (Fernández-Viejo et al., 2005), whereas density values were based on laboratory measurements by Christensen and Mooney (1995) and Christensen (1996).

The P- and S-wave quality factors, Q_p and Q_s, respectively, were more difficult to assign because no attenuation studies have been carried out in the region. Values of Q_p = 500 and Q_s = 200 for the crust and Q_p = 1000, Q_s = 400 for the upper mantle were used for generating the synthetic seismograms displayed. These are “generic” values typically used within the algorithms by other scientists who have used it for crustal studies. The Q_p value for the crust is consistent with values of 400 ± 200 determined for the continental crust in England (Scheirer and Hobbs, 1990). In general terms, Q_s is about one half the value of Q_p (Fowler, 2005). Although Q_s must be specified within the algorithm, this study does not include shear-wave arrivals, so the value specified is not significant. Values of the quality factor for the uppermost mantle are generally considered higher. The generic value of Q_p = 1000 is consistent with that shown by Fowler (2005) for a general Earth model. It should be noted that quality factors of 1000 or more are indicative of very limited attenuation. Moreover, changes in the values of Q by factors of one half or two do not materially affect the dynamic characteristics of the seismograms, as established by additional calculations using different Q values.

For the wide-angle data, a zero-phase Ricker wavelet was used as the source signal with a bandwidth of 4–16 Hz to simulate an explosive source. For the Vibroseis data, which is characterized by higher frequencies, a zero-phase Ricker wavelet with a 20 Hz central frequency of a two-octave wavelet (10–40 Hz) was used as the input source.

As represented by WAR and NVI data sets, the characteristic features of the observed wave field of the Moho below the Great Bear magmatic arc are:

• (1) prominent PmP phase between 80 and 220 km offsets, followed by ∼0.5-s-long coda, and the absence of this phase at near-vertical incidence angles (Fig. 4);

• (2) amplitude behavior of the PmP phase with offset that is highly erratic and no common amplitude peak observed among the different shots (Fig. 5A);

• (3) variable dominant PmP frequency, measured on individual traces, at different offsets between 10 and 13 Hz (however, no clear relation between amplitude spectra and offset is observed; Fig. 5B);

• (4) a very weak amplitude Pn phase at offsets >160 km offset; and

• (5) a Moho signature on the NVI data characterized by a very distinct, mutlicyclic band of reflectors that are laterally variable and piecewise continuous over several kilometers (Fig. 3). Amplitude spectra measured over 25 traces show a dominant frequency centered at ∼21 Hz.

We calculated the synthetic seismograms of 1-D and 2-D velocity models to determine the extent to which these models can replicate the wave-field characteristics noted here. Both early studies of the properties of the Moho (e.g., Clowes and Kanasewich, 1970) and more recent ones (e.g., Carbonell et al., 2002) have demonstrated that representation of the crust-mantle transition by a first-order discontinuity, or a linear increase in velocity from lower-crustal to upper-mantle values over a limited depth interval, does not generate synthetic seismograms that replicate observed wave-field characteristics. Our studies of such Moho models, not discussed herein, further confirmed these results. Accordingly, we initiate our discussion with layered transition models.

The crust-mantle transition is represented by alternating high- and low-velocity layers with a velocity contrast of 0.6 km/s (Fig. 6). Alternative models with different/variable velocity contrasts (0.2–0.8 km/s) were also examined. Following Carbonell et al. (2002), these layered models are characterized by two parameters, the total thickness of the transition zone (H) and the thickness of the internal layers (h). Initial modeling suggests that H controls the length of PmP coda, whereas the number of wiggles is controlled by h, as well as the velocity contrast, similar to the results obtained by Carbonell et al. (2002). We kept the total thickness of the transition zone fixed at H = 3 km, which produces an ∼0.5 s coda length, comparable to the observed data, and we varied the internal thickness (h) between 100 and 500 m. Layered models produce prominent, reverberatory PmP phases beyond 80 km offset but rather weak Pn phases at around 140 km offset (Fig. 6). In a test of two models with the same H = 3 km, but one featuring h = 300 m and the other with h = 500 m, the latter generates a phase characterized by a more reverberatory pattern (Fig. 7). For very thin internal layers (h ≤ 100 m), the Moho signature is very similar to that generated by a first-order discontinuity model. From models with varying velocity contrasts, we note that smaller values produce fewer numbers of wiggles.

Amplitude-versus-offset (AVO) curves are relatively simple compared to the observed PmP amplitudes, regardless of the internal thickness of the layers. They are characterized by a clear peak amplitude located at ∼120 km offset that corresponds, as expected, to the critical distance (Fig. 8). Beyond 120 km, there is a nearly steady decrease in amplitudes with offset.

The amplitude spectrum of the PmP phase was calculated on single traces and over a 0.5 s window (here and in all subsequent models). The spectrum features a varying dominant frequency between 7 and 12 Hz at different offsets. There is no correlation between the dominant frequency and varying offsets (Fig. 9A). Also, the internal thickness of the layers has a limited effect on the shape of the spectrum (Fig. 9B).

The synthetic NVI seismograms feature a distinct Moho characterized by sharp, mutlicyclic bands of reflectors (4–9 bands) depending on the thickness (h) of the internal layers (Fig. 10). Higher velocity contrasts between the internal layers result in higher amplitude reflections and vice versa. The dominant frequency of these reflections, measured from averaging the amplitude spectrum over 25 traces and for different layer thicknesses, ranges between 19 and 23 Hz, compared with the 20 Hz central frequency of the source wavelet. This is due to the fine-tuning effect of the varying thicknesses of the internal layers (Fig. 11). We notice that thicker layers produce more spiked-shape spectra, whereas layers <100 m thick produce a smoother spectrum similar to that of the input source.

Previous results and those from our study indicate that the crust-mantle transition is not a simple, uniform discontinuity in which velocities increase suddenly or linearly from crustal to mantle values, but rather a complex and variable transition zone. One possible scenario that more closely replicates the observed data is a transition zone that consists of alternating high- and low-velocity layers, as described in the preceding based on 1-D modeling. However, a 1-D model does not generate the lateral variability that is observed, particularly on the NVI reflection gathers and stacked sections (Figs. 2 and 3). In order to account for all the characteristics of the wave field, a 2-D wave propagation algorithm is needed, one that can handle both vertical and lateral velocity variations.

To study the effects of both lateral and vertical heterogeneity at the crust-mantle boundary, we investigated a range of velocity models that represent two end members. In one model, we use laterally discontinuous layers with variable thicknesses and velocities. This represents an anastomosing layered structure with juxtaposed rocks from the lower crust and upper mantle (e.g., Jarchow and Thompson, 1989). In the second approach, the velocity model is assumed to be made up of two components (Mereu and Ojo, 1981): (1) a deterministic background component in which the velocity is increasing with depth and (2) a random velocity component consisting of a field of small elliptical reflectors embedded in the background velocity field. Random numbers are generated, which in turn are used to distribute the ellipses throughout the transition zone and to assign velocity contrast values to each ellipse (Mereu, 2003) with 2% standard deviation. The ellipses are characterized by their correlation length, which has a horizontal (a_x) and vertical (a_z) component. For the WAR data, we calculated the synthetic seismograms for three shots that simulate shots 1104, 1105, and 1108. In the case of NVI data, several synthetic seismograms simulating normal shot gathers are generated. For all cases, the velocity structure on either side of the region of the crust-mantle transition is derived from the refraction analysis of line 1 data (Fernández Viejo and Clowes, 2003).

The Moho transition consists of a complexly layered velocity structure over a 3-km-thick zone. Individual layers have thicknesses varying between 100 and 500 m, and their horizontal extent ranges between 5 and 20 km, with varying P-wave velocities of 6.9, 7.5, and 8.1 km/s (with 2% standard deviation), characteristic of lower-crust and upper-mantle values (Fig. 12A).

Synthetic seismograms generated by the viscoelastic finite-difference algorithm display a prominent PmP phase followed by an ∼0.5-s-long coda, in agreement with the observed data (Figs. 12B and 4). Through a variety of tests, layers as short as 3 km long still produce a PmP phase with similar dynamic characteristics. On the other hand, an increase in the total thickness of the transition zone to 6 km results in a much longer coda (∼1.2 s), which is not observed. In general terms, the synthetic seismograms of Figure 12 compare favorably with the observed data of Figure 4, including variability of the PmP phase and a weak Pn phase beginning at ∼170 km offset. The PmP amplitude analyses show complex behavior with offset in all three shots (Fig. 13). This resembles the amplitude analysis results for the observed data (Fig. 5). Amplitudes were corrected for geometrical spreading by multiplying all traces by r^½ (r = source-receiver distance), assuming a point source in 2-D as equivalent to a line source in three dimensions (3-D).

On the synthetic NVI data, the Moho signature features a strong, multicyclic band of reflectors. However, only models that have layers with shorter horizontal extent (3–7 km) display laterally variable and piecewise continuous reflections, similar to the observed data (Fig. 14).

The Moho is represented by a 3-km-thick transition zone characterized by randomly distributed ellipses with correlation lengths of a_x = 1000 m and a_z = 100 m and velocities ranging between 7.1 and 7.8 km/s (Fig. 15A). This velocity structure produces a strong and reverberatory PmP phase characterized by ∼0.5-s-long coda, comparable to the observed data (Fig. 15B). Increasing the dimensions of these ellipses (a_x = 3000 m and a_z = 300 m) or replacing their varying velocities with one fixed value at 7.2 km/s does not change the dynamic characteristic of the PmP phase substantially. On the other hand, when a_x < 700 m, the synthetic seismograms feature a simple PmP phase lacking a visible coda, similar to that generated by a first-order discontinuity. A velocity contrast as small as 0.1 km/s still produces a reverberatory PmP phase. Distributing the ellipses over a 6-km-thick zone produces a distinct PmP phase with a very long coda (∼1.5 s), whereas a 1-km-thick transition zone produces a PmP phase without a coda. Table 3 summarizes all the different combinations of parameters that were tested. All of the scenarios generate a weak Pn phase observed at offsets >160 km. As in the previous model, the PmP amplitude curves of all the different shots and for different correlation lengths vary considerably with offset (Fig. 16).

On the NVI data, the Moho resulting from this lamellae structure is imaged as a narrow (0.4–0.7 s two-way traveltime [TWTT]) band of 4–8 cycles of reflections (Fig. 17), which is comparable to the observed data (Fig. 3). These reflections are laterally variable and piecewise continuous over ∼10 km distance. Our results show that increasing the total thickness of the transition zone from 3 to 6 km would result in a much longer (∼0.9 s TWTT) band of 12–13 cycles of reflections. Different correlation lengths change the number of cycles, and even with values as small as a_x = 300 m and a_z = 30 m, distinct bands of reflections can still be detected, attesting to the higher resolution of the NVI data. However, a minimum velocity contrast of 0.2 km/s is needed between the ellipses and the surrounding material in order to produce visible reflections.

From the results of synthetic modeling and based on the duration of the PmP coda, a 3-km-thick Moho is required to generate 0.5-s-long coda. Thicknesses greater than 5 km or less than 2 km produce either a longer coda (∼1.0 s) or a PmP phase lacking any visible coda, respectively. Furthermore, observed NVI data at Moho levels feature a distinct band of reflections over 0.45 s traveltime (0.9 s TWTT). Considering average seismic velocity at this level to be ∼7.4 km/s, this band translates to an ∼3.3-km-thick zone, which is in good agreement with the results from the synthetic WAR data.

Only a laterally and vertically heterogeneous crust-mantle boundary can simulate all the dynamic features observed on the wide-angle and near-vertical data (Figs. 12 and 15). This transition consists of discontinuous layers/ellipses characterized by variable correlation lengths and velocities. Synthetic models for both WAR and NVI data allow us to place lower and upper limits on the values for such correlation lengths and velocities. A minimum thickness of 100 m and a 0.2 km/s velocity contrast are needed for the layers/ellipses in order to reproduce the reverberatory pattern of the PmP phase. Models with layers/ellipses thicker than 600 m do not simulate the multicyclic band of reflections observed on the NVI data. On the other hand, layers/ellipses that are at least 1 km but no more than 6–7 km long are required to generate a distinct PmP coda and to achieve the lateral variability seen on the NVI data.

Calculated amplitude versus offset curves for the PmP phase from the layered 1-D velocity model are characterized by one prominent peak amplitude at around 120 km. In contrast, amplitude curves calculated from the heterogeneous models have more complicated shapes with multiple high amplitudes, in agreement with the calculated amplitudes from the observed data (Figs. 8, 13, and 16). The deviation between these curves and the ones calculated for simple 1-D velocity models is most probably the result of constructive and destructive interference effects due to the lateral heterogeneity as well as the vertical complexity of the crust-mantle boundary.

Based on the spectral analysis, the dominant frequency of the PmP phase varies randomly between 7 and 12 Hz at different offsets. Thus, no relation can be inferred between frequency and offset. In addition, the wide-angle data are less sensitive to the fine structure of the crust-mantle boundary (e.g., internal thickness of the layers) than the near-vertical data (Figs. 9B and 11). Unfortunately, spectrum analysis of the reflection waveforms is insufficient to distinguish between different heterogeneous models.

The Wopmay orogen is an amalgamation of Paleoproterozoic domains that developed and accreted onto the Slave craton between 2.1 Ga and 1.84 Ga by complex geologic processes (Hoffman, 1989). Thus, in this region, reflection patterns near the crust-mantle boundary are expected to be similarly complex (e.g., Clowes et al., 1987). However, observations along SNORCLE line 1 reveal simple Moho reflections characterized by subhorizontal surfaces and relatively uniform arrival time (∼11.0 s) over the entire profile, even beneath regions with different ages such as the Archean Slave Province and the Paleoproterozoic Wopmay orogen (Cook et al., 1999). This uniformity implies that postaccretion, regionally extensive processes have effectively modified rocks and reset the reflection Moho to a uniform depth (∼33 km). Therefore, a complete interpretation for the development of the Moho must consider all possible thermal and/or mechanical processes in order to account for its formation and evolution from initial complexities during accretion to increasing coherency over time (Cook, 2002).

Several hypotheses have been proposed for the formation and development of the continental Moho that can be classified under either igneous or nonigneous origins (Eaton, 2006). Igneous origins include thermal front (e.g., partial melt and magmatic intrusion) and magmatic underplating, whereas nonigneous origins include metamorphic or metasomatic front (e.g., mafic granulite to eclogite) and regional décollement (Eaton, 2006; Cook et al., 2010).

The spatial correlation between the Moho and crustal reflections provides important constraints on the origin of the crust-mantle boundary. For instance, lower-crustal reflections that are listric into the Moho suggest that the crust-mantle boundary is most probably a regional décollement (Cook and Vasudevan, 2003). The transition zone acts as a rheologically weak layer and is underlain by a stronger layer. On the other hand, truncation of lower-crustal reflections at Moho levels could be caused by thermal or metamorphic fronts (e.g., Carbonell et al., 2002; Cook et al., 2010). Interestingly, seismic data from the lower crust of the Great Bear magmatic arc, part of SNORCLE line 1, exhibit both types of spatial correlation between the Moho and lower-crustal reflections near 10–11 s (Fig. 2). Thus, in our synthetic models, the lamellae structure with randomly distributed ellipses represents a thermal/metamorphic front, whereas the laterally discontinuous layering represents a regional décollement.

Many deep seismic profiles reveal lower-crustal reflections that are listric into the Moho (Cook et al., 1992, 1999; Calvert et al., 1995; Kukkonen et al., 2008). This suggests that the transition zone near the crust-mantle boundary behaves as a mechanically weak layer (more plastic) into which upper-crustal detachments sole (Eaton, 2006). This process could be associated with either regional extension and/or contraction (Cook, 2002).

Within the lithosphere, the strength of minerals decreases considerably from olivine, amphibole, and garnet to feldspar and quartz, the latter being the weakest mineral (Austrheim et al., 1997). The lower crust of the Hottah terrane consists primarily of amphibolite and paragranulite (Fernández-Viejo et al., 2005). Consequently the rheology is presumably controlled by amphibole and granulite, whereas olivine controls the rheology of the mantle. Therefore, at the crust-mantle boundary, one expects a distinct strength contrast that may tend to localize detachment surfaces (Cook, 2002; Cook et al., 2010). Based on published geotherms and lithospheric composition, Meissner and Mooney (1998) estimated the depth of low-viscosity (weak) zones within the crust and upper mantle along which decoupling may occur. They found these weak zones at three depths: (1) the base of the felsic upper crust; (2) just above the Moho; and (3) some tens of kilometers below the Moho.

In the lower crust beneath the Great Bear magmatic arc, a series of discrete, northeast-dipping reflections (Fig. 2) is listric into the Moho. Moreover, improved seismic images of the same profile (fig. 3d inCook and Vasudevan, 2003) show distinct upper-mantle structures, which were interpreted as a synform fabric overlain with angular discordance by the flat Moho. Such observations support the interpretation of the Moho as a structural boundary along which lower-crustal and upper-mantle rocks are displaced, sheared, and blended. However, as Cook and Vasudevan (2003) explained, the Moho here cannot be regarded only as a regional décollement for the following two reasons: (1) the crust is much more reflective than the upper mantle, and with a simple detachment, the rocks above the detachment surface should be equivalent to, but dislocated from, rocks below; and (2) a detachment surface does not explain how rocks below the Moho have significantly higher seismic velocities (∼8.1 km/s; Fernández Viejo and Clowes, 2003) than those above it (∼6.8 km/s).

From the preceding discussion, consideration of the Moho as a regional décollement represents only one stage of its development. A logical sequence of stages or events is probably necessary to explain the varying characteristics of the crust-mantle boundary observed on both NVI and WAR data.

Mafic igneous rocks may be injected into the lower crust and upper mantle in the form of sills during the amalgamation of island arcs at collisional orogens (e.g., Wopmay orogen). These sills are characterized by relatively high density and high seismic velocity (>7.2 km/s). The intermediate-composition crust above acts as a density filter for these intrusives, and hence they pond within the lower crust and along the crust-mantle boundary, producing seismically laminated lower crust and a sharply defined Moho (Furlong and Fountain, 1986; Nelson, 1991).

Partial melting can modify the chemical and physical properties of the intermediate-composition rocks at elevated lower-crustal temperatures and pressures by preferentially removing lower-density minerals such as quartz and feldspar. This will leave mafic residuals with higher density and higher seismic velocity (>7.5–7.8 km/s) appropriate for the upper mantle (Hynes and Snyder, 1995; Rudnick and Fountain, 1995). The solidus temperature for these restites is considerably elevated, and any further partial melting requires substantially higher temperatures. As a result, the Moho could be “frozen” unless significantly higher temperatures occur (Cook 2002). Because partial melting might not remove all of the lighter fraction from the lower crust, structural geometry near the Moho can be preserved (Cook et al., 2010).

As mentioned previously, termination or truncation of some of the dipping structures in the lower crust of the Great Bear magmatic arc by the Moho, in addition to the discordance between the underlying synform fabric and the Moho (Cook and Vasudevan, 2003), strongly suggests that the Moho postdates the synform and lower-crustal structures. Thus, it may have been superimposed on them by thermal processes. The approach and subsequent subduction of oceanic lithosphere beneath the Hottah terrane was probably the source of thermal activity responsible for partial melting and magmatic intrusion in the study area, which in turn was associated with the deposition of the Great Bear magmatic rocks between ca. 1.88–1.84 Ga (Cook et al., 1999; Oueity and Clowes, 2010).

Since the work of Green and Ringwood (1967) and Ito and Kennedy (1971), it has been recognized that under relatively high pressure-temperature (P-T) conditions, mafic rocks (e.g., mafic granulite, gabbros) might transform into eclogite at depths of ∼35 km. Hence, the Moho could be interpreted as a metamorphic front associated with a phase change. Eclogitization results in an increase in the density and seismic velocity of the lower crust such that it becomes indiscernible from the uppermost mantle (e.g., Austrheim, 1987; Mengel and Kern, 1992; Fischer, 2002; Kukkonen et al., 2008). Consequently, lower-crustal eclogitization has been invoked to explain the observation of a relatively shallow Moho at several places (e.g., the Scottish Caledonides [Hynes and Snyder, 1995] and the Appalachians off northeastern Newfoundland [Chian et al., 1998]).

Along the southwestern coast of Norway, eclogites are exposed over a large area in the Western Gneiss region (WGR). They form pods and lenses ranging in size from decimeters to several hundreds of meters in length embedded in amphibolites (Austrheim, 1987). This is very similar to the lamellae structure we tested in our synthetic models. Although in the Western Gneiss region, most eclogites are mylonitic, examples of partial eclogitization are abundant. Therefore, in our model, the wide range of velocities/densities characterizing the randomly distributed ellipses could be explained by different degrees of eclogitization, depending on the rock composition and availability of water.

Observed wide-angle and near-vertical seismic-reflection data and 1-D and 2-D synthetic seismogram modeling for the crust-mantle boundary offer different but complementary information about its nature and formation. A unified interpretation, in which both wide-angle and near-vertical data sets derived from the same Moho location, can provide a more complete explanation concerning the structure of the Moho, one with higher resolution and perhaps less ambiguity.

Our synthetic models lead to the conclusion that beneath the Great Bear magmatic arc, an ∼3-km-thick, laterally and vertically heterogeneous crust-mantle boundary properly simulates the dynamic characteristics of the observed wave fields recorded on the near-vertical and wide-angle reflection data. The heterogeneity can be achieved by either laterally discontinuous layering or a lamellae structure with randomly distributed ellipses. These two models may represent different re-equilibration mechanisms, one characterized by thermal/metamorphic fronts and another that features a regional décollement. Due to insufficient details on the seismic-reflection image, it is difficult to establish the exact geometric relationships between the Moho and lower-crustal layers and consequently the re-equilibration process. However, assuming that our images of the Moho beneath the Great Bear magmatic arc represent “snapshots” of the crust-mantle boundary “taken” at different stages of its evolution (Keller et al., 2005), we believe that the Moho has formed as a combination of the aforementioned processes.

Here, we propose a scenario for the evolution of the Moho (Fig. 18) that builds on, but differs somewhat from, a previous study (Cook and Vasudevan, 2003) as follows: (1) Prior to 1.9 Ga, the Moho represented the base of the Hottah terrane crust, which was generated as a magmatic arc on cryptic 2.4–2.0 Ga crust (Hoffman and Bowring, 1984). The depth to the Moho is unknown but was perhaps at ∼35 km. (2) The collision of the Hottah terrane onto the Slave craton caused substantial compression and led to the formation of the listric structures that flatten into the Moho. This implies that the Moho at this stage acted as a ductile decoupling zone or a detachment surface. (3) Between 1.88 and 1.84 Ga, the Great Bear magmatic arc was formed as a result of east-dipping subduction of oceanic lithosphere beneath the Hottah terrane and the Slave craton generated by incipient of the Fort Simpson terrane. The thermal activity associated with the previous subduction produced partial melt and magmatic intrusions at the base of the crust. A horizontal Moho was developed as a thermal/metamorphic boundary.

We wish to thank Phil Hammer for stimulating discussions and constructive comments. We thank the editor for comments that helped clarify some parts of the text. The data were acquired as part of Lithoprobe, which was funded through a Research Networks grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Geological Survey of Canada. This study was supported by an NSERC Discovery grant (7707-2009) to Clowes.