The closure of ancient oceans created a dynamic setting suitable for craton formation via the thickening of continental material over a mantle downwelling. This process subjected the thickening lithosphere to extensive deformation, forming internal structure that can be preserved over the lifetime of the craton. Recent seismic imaging of cratonic lithosphere has led to observations of anomalous features colloquially known as midlithospheric discontinuities. These discontinuities are attributed to a range of sources, including the lithosphere-asthenosphere boundary, melt accumulation, and phase transitions. However, the internal structure imaged within these cratons might be reflective of their formation. In particular, the orientation and nature of the variable depths of the midlithospheric discontinuities suggest a more complicated origin such as that which could be introduced during the formation and thickening phase of cratonic lithosphere. Here, we present geodynamic models demonstrating the internal structures produced during the formation of cratonic lithosphere as well as new seismological observations of midlithospheric discontinuities in the West African craton, together with reassessment of midlithospheric discontinuities observed in the North American, South African, Fennoscandia, and Australian cratons. We suggest that the midlithospheric discontinuities observed in these cratons could be remnants of deformation structures produced during the formation of the cratons after ancient oceans closed.


The Wilson cycle model was introduced to describe the closing and opening of the Atlantic Ocean basin (Wilson, 1966). Wilson (1966) proposed that continents and their margins were continuously modified by collisional events that closed oceans and also by rifting events that then opened ocean basins. The closing of an ocean basin sutures continents together and produces new continental material through the accretion of arc material. During the suturing process, both the new and preexisting continental materials thicken and deform. In Wilson’s model, the newly created larger continent will then proceed to rift apart and break up (possibly in areas of preexisting weakness created during the collisional event) until a new ocean basin is created, and the cycle begins anew. We ask the question: Can the process of closing an oceanic basin form thick stable continental lithosphere that does not fall victim to the next stage of rifting in the cycle? Here, we present evidence that the closing of Paleoproterozoic oceans caused the formation of the cratons in Canada, Australia, Fennoscandia, West Africa, and South Africa, where there have been recent significant breakthroughs in seismic imaging.

Cratons are regions of long-lived, stable, thick lithosphere (Jordan, 1978; Pollack, 1986). Their stability comes primarily from their intrinsic strength determined by enhanced chemical buoyancy, viscosity, and finite strength (Jordan, 1978; Pollack, 1986; Lenardic and Moresi, 1999; Sleep, 2003). Stability might also be ensured by proximity to weaker material that would preferentially localize deformation, protecting cratonic lithosphere from destructive forces (Lenardic et al., 2003). In either case, craton stability is dependent on material properties to overcome the deforming forces exerted by the convecting mantle. Their long-term stability suggests that the inherent strength of cratons must be sufficient to offset deformation from initial stabilization forward; i.e., the craton must be able to sustain the ever-changing dynamic conditions within the mantle resulting from the cooling of Earth. In addition to intrinsic properties such as density, viscosity, and finite strength, thickness also may determine the stability and longevity of cratonic lithosphere. Therefore, thickness might not just be a characteristic of cratonic lithosphere, but a necessary requirement to help explain its stable nature (Lenardic and Moresi, 1999; Lenardic et al., 2003; Sleep, 2003; Cooper and Conrad, 2009).

Models of the formation of stable continental lithosphere often invoke the thickening of buoyant material over a mantle downwelling (e.g., Jordan, 1978, 1988; Bostock, 1998; Lee, 2006; Cooper et al., 2006; Rey and Houseman, 2006; Duclaux et al., 2007; Gray and Pysklywec, 2010). These studies suggest cratons are formed either through the stacking of buoyant oceanic lithosphere (Hart et al., 1997; Bostock, 1998; Musacchio et al., 2004; Cooper et al., 2006; Lee, 2006) or the progressive amalgamation of island arcs (Jordan, 1978, 1988; Rey and Houseman, 2006; Duclaux et al., 2007; Gray and Pysklywec, 2010). Regardless of the specific type of thickening or the type of material thickened, both mechanisms provide a means for stabilization as well as introduce complex internal structure to the newly thickened cratonic lithosphere. In addition, both scenarios are likely to occur during the closing of an ocean basin. The discrepancy between the two models primarily comes from the need to explain the chemistry of cratonic lithosphere and, in part, by the presence of seismically imaged internal structure within the craton (e.g., Bostock, 1998; Musacchio et al., 2004; Olsson et al., 2007; Wittlinger and Farra, 2007; Ford et al., 2010; Fischer et al., 2010; Miller and Eaton, 2010; Yuan and Romanowicz, 2010; Darbyshire et al., 2013; Snyder et al., 2013). For example, a midlithospheric seismic discontinuity observed in the North American craton has been interpreted as a remnant feature of craton formation (Abt et al., 2010; Fischer et al., 2010; Miller and Eaton, 2010; Yuan and Romanowicz, 2010; Snyder et al., 2013). The presence and nature of the internal structure of cratonic lithosphere should delineate between the two models and provide information on the formation and stabilization of cratons; i.e., the formation mechanism must produce structure that (1) can be preserved over billions of years and (2) can be imaged seismically, which requires sharp velocity contrasts or changes in anisotropy within the cold, thick craton keel.

The process of forming cratons through the thickening of buoyant material can be bracketed by two deformation mechanisms: thickening through thrust stacking and localized deformation or viscous thickening and distributed deformation. These two end members call upon differing physical mechanisms for thickening and stabilization. For example, if thick cratons form from thrust stacking of buoyant material, then that initial material must initially be highly viscous for deformation to primarily occur in localized shear zones (Cooper et al., 2006). Starting with a material with a high viscosity provides a means for stabilization and does not require subsequent modification of the material (Cooper et al., 2006). In addition, the localized shear zones have the potential to be preserved as ancient suturing or deformation features that could define the deep internal seismic structure imaged within some of Earth’s cratons (e.g., Bostock, 1998; Rychert and Shearer, 2009; Ford et al., 2010; Miller and Eaton, 2010; Darbyshire et al., 2013; Snyder et al., 2013). Conversely, if cratons form through the viscous thickening of buoyant material, this results in continuous deformation over a broad zone (Dewey and Burke, 1973; Jordan, 1978; Choukroune et al., 1995). The cratonic lithosphere then becomes viscously strong after a period of cooling, which then provides the necessary viscosity required for stabilization and longevity (Lenardic and Moresi, 1999). Formation of cratons through this type of thickening process would be less likely to produce deformational features that could be preserved over a long period of time, but it still could introduce compositional variations at depth that could be seismically observable.

While we do not rule out the possibility that cratons could be formed through an intermediary process between the two end-member cases, we have highlighted these cases to provide a dynamic framework to guide our interpretation of seismic observations. In other words, in order to use internal cratonic structure as an indicator of formation, stabilization, or modification, it is important to consider what structure we might expect given the dynamics of the proposed formation mechanisms. With that goal, we compare seismic observations of several cratons to each other and with geodynamic models that highlight the end-member processes of thickening of buoyant material over a mantle downwelling.


Seismology-based estimates of the depth to the base of the cratonic lithospheric appear to have large variations in thickness (e.g., Nataf and Ricard, 1996; Artemieva and Mooney, 2002; Rychert and Shearer, 2009; Eaton et al., 2009), but all find that fast shear and compressional wave velocities overlie a lower-velocity region, and that boundary resides between ∼100 and 300 km. On a global scale, the continental keels are imaged with surface wave tomography as fast velocity perturbations (>1%) in these depth ranges, although the depths to their bases are more difficult to determine due to the nature of the technique (e.g., Grand et al., 1997; Mégnin and Romanowicz, 2000; Houser et al., 2008; Ritsema et al., 2011). However, the locations of the fast velocities associated with the Precambrian cratons, in particular, those in Canada, Australia, Fennoscandia, West Africa, and South Africa, are quite clearly imaged by the global composite mantle tomography model SMEAN (Fig. 1; Becker and Boschi, 2002).

Recent regional surface wave tomography for these areas provides more detail on the cratonic structure (e.g., Bruneton et al., 2004; Fishwick et al., 2005, 2008; Chevrot and Zhao, 2007; Eaton and Darbyshire, 2010; Fishwick, 2010; Darbyshire et al., 2013) and finds that the fast velocity perturbations are greater (>3%) than those found in global models. Although the base of the craton remains challenging to define using surface wave–based techniques, these studies have inferred steps, gradients, and perhaps dipping layers at the edges and even within the craton. These regional surface wave imaging studies suggest that the deep keels appear to have a more complex internal structure than originally proposed or imaged in global studies.

Recent efforts in imaging the internal structure of the continents using receiver function techniques have detected boundary layers within cratons at variable depths between ∼60 and 160 km (e.g., Bostock, 1998; Olsson et al., 2007; Snyder, 2008; Hansen et al., 2009; Rychert and Shearer, 2009; Abt et al., 2010; Ford et al., 2010; Miller and Eaton, 2010; Snyder et al., 2013), which are identified by a negative-polarity pulse similar to the signal produced by the lithosphere-asthenosphere boundary. Deeper (>150 km) discontinuities are generally interpreted as the base of the lithosphere (for compilations of previous results, see Rychert et al., 2010; Fischer et al. 2010), although the shallower (60–160 km) features have also been suggested to be the lithosphere-asthenosphere boundary (cf. Rychert and Shearer, 2009) within the continental interior. Alternatively, this intermediate-depth feature evident in the cratons has been termed the midlithospheric discontinuity, but the origin of this feature has yet to be resolved. These structures appear to be present in many of the Paleoproterozoic- and Archean-aged cratonic keels (Fig. 1), but they are not consistent features that are evident as a distinct shear velocity decrease or a change in anisotropy from surface wave studies, making them a challenge to identify and interpret.

Large-scale scattered imaging studies of the continental interiors, such as those of Rychert et al. (2007), Abt et al. (2010), and Ford et al. (2010) for most of cratonic North America and Australia, illustrate that the Archean and Paleoproterozoic terranes have internal cratonic structure as identified by midlithospheric discontinuities. Other studies, also using receiver functions (Miller and Eaton, 2010; Snyder et al., 2013), have shown that the midlithospheric discontinuities beneath the Paleoproterozoic Superior craton appear to be dipping outward from the center of the craton, suggesting that they may be relict features resulting from the accretion of slabs and smaller cratonic terranes. These interpretations are also supported by more detailed and high-resolution array-based P receiver function (Bostock, 1998; Snyder, 2008) and reflection seismology imaging (Cook et al., 1997), which also detected dipping lithospheric structures within the Archean-age Slave Province. However, the apparent dip of these midlithospheric discontinuities in the younger Paleoproterozoic portion of the craton is away from the deepest part of the cratonic keel (∼260 km), making it difficult to reconcile with the model of ancient stacked or accreted oceanic slabs. This then leads to questions about whether these dipping structures could be another deformational feature, such as preserved sutures or shear zones and whether they are actually continuous structures.

There have been a series of studies using receiver functions to image cratonic structure in Fennoscandia and Africa that have also found internal cratonic structures (Olsson et al., 2007; Wittlinger and Farra, 2007; Kumar et al., 2007; Savage and Silver, 2008; Hansen et al., 2009) through identification of a negative conversion at ∼110–160 km depth (Fig. 1). The estimates of the base of the lithosphere are somewhat debated, but there are commonalities in all the receiver function images for these two regions, in which there are seismic discontinuities that produce the observed negative-polarity conversions internal to the cratons.


Previously, there had been no receiver function–based estimates for lithospheric structure in the West African craton, although there have been surface wave studies that target the region (e.g., Fishwick, 2010). Recent results from the PICASSO experiment in Morocco have found that stations on the craton (Fig. 1) have signatures of midlithospheric discontinuity structure (Miller and Becker, 2014; Miller et al., 2014). However, these stations are concentrated just on the northwest edge of the craton. Although there are very few broadband stations in West Africa, we were able to add additional observations across the region to help further understand the lithospheric structure of this large craton (Fig. 2A). For this study, we analyzed broadband seismic data for 15 broadband stations that recorded 36 earthquakes of magnitude 6 and greater at 55°–85° distance that occurred between 2010 and 2013 with high signal-to-noise ratios. We followed the S-receiver function methodology described in detail in Miller and Eaton (2010) and Levander and Miller (2012). The S-receiver functions were produced by deconvolving the SV component from the P component in the frequency domain following Langston (1977) and were depth converted using the one-dimensional (1-D) velocity model ak135 (Kennett et al., 1995), because it effectively represents the continental lithosphere. Examples of receiver function gathers from two stations (MBO in Senegal and IFE in Nigeria) are shown in Figure 2B.

We found that 10 of the 15 stations in Mali, Morocco, Ghana, and Nigeria indicate internal cratonic discontinuity structures or midlithospheric discontinuities (Table 1). We find that the stations off the edge of the craton have only a single clear lithosphere-asthenosphere boundary signal with no deeper (or intermediate depth) signal from the midlithosphere. The 10 stations that are clearly on the craton, as interpreted from fast shear wave velocities from SMEAN (Becker and Boschi, 2002), have a range of midlithospheric discontinuity signals between 75 and 160 km with an interpreted lithosphere-asthenosphere boundary signal extending to ∼300 km (Fig. 2A). Examples of receiver function gathers from stations IFE and MBO in Figure 2B show the existence and lack of midlithospheric discontinuity signals, respectively.

Due to the oblique raypaths that S waves from earthquakes between 55° and 85° have, the locations of conversion points (or piercing points) are at a distance away from the surface location of the station. Therefore, the interpreted midlithospheric discontinuity depths are plotted at their estimated piercing points based on back-azimuth bins, which demonstrate the variable depth range where the waves sample different parts of the cratonic lithosphere (Fig. 2A). Plotting the interpreted midlithospheric discontinuity depths in this way further shows that none of the stations is located very close to the continental shelf nor do those off the edge of the craton have clear internal lithospheric structure (MBO in Senegal and those in southernmost Ghana). In addition, all the stations in the interior of the craton, except TAM in Algeria, which is on the Hoggar swell (Liu and Gao, 2010), have midlithospheric discontinuity signals that are deeper than those placed on the edge of the craton (Fig. 2A). It is not surprising that TAM lacks any midlithospheric discontinuity signal because it is set on a region of Cenozoic volcanism and in the middle of a massif with high elevation and relatively slow shear wave velocities (Liu and Gao, 2010). However, we observe a general deepening of the midlithospheric discontinuities into the center of the craton, but due to the sparse stations, this can only be a very broad generalization.


To better understand how cratons formed and how internal structure can be created and preserved on a billion-year time scale, we use geodynamic models to highlight the processes that occur during craton formation. The closing of an ocean basin introduces a dynamic setting in which two continental masses are thickened over a mantle downwelling (initially existing in the form of a subducting slab). It is this process that we call upon to form thick, cratonic lithosphere, and as such, we can model this by exploring the thickening behavior of buoyant material over a mantle downwelling. An initial investigation addressed the dynamic feasibility of forming cratons via thickening processes (Cooper et al., 2006). We have expanded those models to track the predicted internal structure formed during thickening.

Model Description

Simulations were run using the numerical code Underworld (Moresi et al., 2007). The simulations were configured to isolate the key components driving the interaction between a chemically distinct layer and the convective mantle. The model domain consists of a 2 × 1 box with resolution of 640 × 128 elements (for initial model configuration, see Fig. 3A). Convection was driven from below with constant-temperature top and bottom boundaries. The side boundaries within the model domain are reflecting walls. We included a temperature-dependent viscosity within the mantle such that the viscosity contrast between the coldest and hottest mantle material spans five orders of magnitude. We started from a well-developed flow pattern and temperature field and emplaced a buoyant multilayered lithosphere within the upper thermal boundary layer of the system. The initial length of the emplaced layer spanned 75% of the surface of system, and the initial thickness was set to 7.5% of the system’s thickness. The layer was centered over the box and resided over a predeveloped mantle downwelling.

We used a viscoplastic rheology following that in Moresi et al. (2003). This rheology allows material to deform viscously as long as the local stresses are below a specified yield stress. Once that yield stress is achieved, then plastic deformation occurs, following the continuum representation of Byerlee’s law (Byerlee, 1968; Moresi and Solomatov, 1998; Moresi et al., 2003). Thus, the potential viscoplastic response of the chemically distinct layer will depend on the assigned viscosity and cohesive strength of the layer as well as convective stresses within the model. The buoyancy and thickness of the chemically distinct layer will determine whether the layer will thicken or recycle over a downwelling. Following the regimes outlined in Cooper et al. (2006), we chose a combination of parameters that highlight either a distributed or localized thickening response. All model parameter values are given within the captions in Figures 3 and 4. The simulations were run for several mantle overturns.

The main goal of these simulations was to highlight the internal structure that each end-member deformation style produces (localized vs. distributed thickening). As such, though the lithospheric layers within these models are visually different (demarcated by different colors), their material properties are the same. While varying the buoyancy, viscosity, and yield parameters of each layer will influence the thickening of the model lithosphere, we left this added complexity for future work. Rather, we chose uniform properties for all of the layers that we know will promote either localized or distributed deformation.

Localized Thickening as Applied to Craton Formation and Stabilization

The underthrusting of buoyant material proves to be a viable option for craton formation. Our simulations confirm that of Cooper et al. (2006) and show that the material will thicken along localized shear zones (Fig. 3). This behavior occurs if three conditions coexist: (1) Deformation must be localized (i.e., the layer must possess a viscosity high enough such that viscous deformation will not occur), (2) yielding must be possible, and (3) recycling into the mantle must be avoided. The first condition was constrained empirically by Cooper et al. (2006) and occurs if the layer’s viscosity is ∼1000 times that of the mantle. This is on order with the predicted viscosity contrast between the lithosphere and mantle if lithosphere is dehydrated of bound water (e.g., Pollack, 1986; Hirth and Kohlstadt, 1996). The last two requirements are simply stress and force balances set by the intrinsic properties of the material and the convecting mantle. In other words, a layer cannot yield unless the convective stresses exceed its yield strength, and a layer will not recycle if its positive chemical buoyancy exceeds the negative thermal buoyancy of the downwelling.

Within the simulations, yielding occurred on newly formed, localized shear zones and continued in those locations as long as the conditions for yielding existed. (Note, we defined the locations of shear zones within the simulations as regions that showed localized, linear regions of high strain rate.) If those conditions were no longer met, then the previously active shear zone was considered “healed” within the simulation. Thus, though it does not remain a zone of long-lived weakness within the simulation, it can still introduce long-lived structure within the lithosphere. As such, yielding often occurred in conjugate bands and evolved in the simulations as the magnitude and orientation of the driving stresses also changed (lower panels of Fig. 3B).

As stated in Cooper et al. (2006), the thickening process itself provides the stabilizing mechanism for a craton formed through the underthrusting of buoyant material. Given that the condition for yielding depends on lithospheric thickness, it is apparent how the thickening of a viscously strong material can create a nondeformable feature; the increased thickness increases the overall finite strength of the material such that it can exceed the local convective stresses, thus making it stable.

Localized Thickening: Expected Internal Structure and Potential to be Seismically Imaged

This deformation response introduces regions of localized shearing (highlighted regions in Fig. 3). These shear zones could be the source of the variable-depth midlithospheric discontinuities imaged in the cratons (e.g., Olsson et al., 2007; Whittlinger and Farra, 2007; Ford et al., 2010; Miller and Eaton, 2010; Abt et al., 2010). Ancient deformational features have also been seismically imaged in other regions and demonstrate that sutures and shear zones can be preserved over geologic time (Cook et al., 1997; Magnani et al., 2004; Karlstrom et al., 2005). The intense and localized shearing during this thickening process could introduce anisotropy into the cratonic lithosphere, providing a potential source for the observed midlithospheric discontinuities. This type of deformation may preserve relict features that possess a dip. The orientation of the dip of the relict shear zone depends on the direction of convergence as well as on the occurrence of conjugate shear zones. Thickening via localized deformation introduces conjugate shear zones (Fig. 3; see also Cooper et al., 2006). This type of deformation can introduce structures with dips that are in the opposing sense of direction of convergence (Fig. 3). Thus, caution should be used when interpreting these relicts of ancient deformation.

Distributed Thickening as Applied to Craton Formation and Stabilization

The process of forming a thick lithosphere can also be achieved through viscous deformation. During this thickening process in the simulations, little to no yielding is observed (Fig. 4), thus setting the first requirement of this type of behavior—the layer must not be too viscous such that it would allow localized deformation to dominate. As before, this localized/distributed cutoff occurs when the layer is ∼1000 times more viscous than the mantle (as determined empirically by Cooper et al., 2006). Finally, the layer must be able to thicken over the downwelling but not recycle wholesale back into the convecting mantle. This is determined by the force balance between the layer’s intrinsic buoyancy and negative thermal buoyancy of the downwelling as controlled by mantle dynamics. Craton formation via viscous thickening of buoyant material assumes that the entire chemically distinct lithosphere remains intact throughout the thickening process; i.e., the layer does not undergo any substantial dripping at its base.

However, we did observe small-scale dripping within those simulations that highlighted the viscous thickening deformation response (Figure 4). The occurrence of these drips is in large part a consequence of the lower values of viscosity required to allow the material to deform in a distributed, viscous manner. The development of these drips is potentially problematic because they counteract the thickening process. For example, Molnar et al. (1998) determined that a layer when thickened to twice its initial size could be thinned by potentially half in 10–20 m.y. through instability-driven drips. The size and presence of these drips can be suppressed by either a higher initial viscosity or temperature-dependent viscosity (Molnar et al., 1998; Conrad and Molnar, 1999; Dumoulin et al., 2005; Elkins-Tanton, 2007). A higher initial viscosity might be a feasible solution to avoid destructive thinning by drips; however, this runs counter to the lower viscosity needed to resist localized deformation. Given that cratons presently reside in a cooler thermal state (Jaupart et al., 1998), the increase in viscosity required to suppress dripping could be achieved by a temperature-dependent viscosity depending upon the cooling rate after the thickening process. The increase in viscosity must sufficiently suppress thinning via dripping before the layer thins below the critical thickness required for stability.

Distributed Thickening: Expected Internal Structure and Potential to be Seismically Imaged

This type of thickening is less likely to introduce deformation features that could be preserved for long periods of geologic time. The lower viscosity required to thicken the material in a distributed manner does not promote any localized yielding and also may promote viscous relaxation of any localized features. Depending on the composition of the protocratonic material, this type of deformation could still provide a source for the midlithospheric discontinuities, as it maintains compositional layering (Fig. 4). If there is a significant compositional change between the layers within the protocratonic lithosphere, then a seismic discontinuity could be induced. This discontinuity would follow the plane of compositional change. Distributed thickening produces relatively horizontal layering within the cratonic lithosphere (Fig. 4). As such, any compositionally driven seismic discontinuity would be expected to also be relatively horizontal. This could provide an explanation for regions that show relatively uniform depth of observed midlithospheric discontinuities, particularly in Australia, although the station density there is the lowest when compared to the other cratons (Ford et al., 2010).


Although station distribution is either quite sparse or concentrated within a subregion in the Archean and Proterozoic terranes of Earth’s cratons, the presence of internal structure within the cratons suggests that the midlithospheric discontinuities may be a common feature and possibly be relict internal structure from the formation of the oldest part of the continents. The range of depths of the midlithospheric discontinuities from craton to craton and even within a single craton (Fig. 1) strongly suggests that the discontinuities are related to the internal structure of the lithosphere rather due to any postformation modifications. Variations in the depth of midlithospheric discontinuities also limit the possibility that the discontinuities are marking a phase transition within the cratonic lithosphere. The present-day cool, thermal state of cratonic lithosphere (Jaupart et al., 1998) due to long-term cooling limits the presence of thermal perturbations that would influence the depth at which the phase transitions would occur. Thus, the thermal structure of cratonic lithosphere could dampen any variations in the depth of seismic discontinuities produced by phase transitions. However, we cannot rule out that this could be the case in regions that show more of a uniform midlithospheric discontinuity depth.

Finally, it also seems unlikely that the midlithospheric discontinuities represent the lithosphere-asthenosphere boundary in these cratonic regions, where other techniques such as global surface wave tomography indicate a much thicker lithosphere (Fig. 1). Rather, the simpler interpretation is that the midlithospheric discontinuities are relict features of the internal structure produced during the formation of a craton. Not all cratons may form with the same or similar combinations of buoyancy, viscosity, finite strength, and thicknesses, and, as such, different internal structures are likely to be present in different cratonic lithosphere. This would also then be reflected in any patterns that we observe in the midlithospheric discontinuities within cratonic lithosphere. Regions of cratonic lithosphere that show complexity and dip and depth variations in midlithospheric discontinuities may indicate that the craton was formed through processes capable of promoting localized deformation. Regions of cratonic lithosphere with uniform midlithospheric discontinuity depths might be more indicative of formation through viscous thickening.

The ability to delineate between the two interpretations allows us to use present-day observations to elucidate early Earth dynamics as well provide additional information about the type of material from which cratons were formed. In particular, cratonic regions with present-day variations in midlithospheric discontinuity depth and dipping midlithospheric discontinuities could be the end result of the closure of ancient oceans. The midlithospheric discontinuities in this case could be remnant features of the suturing and thickening of ancient arc material possessing significant inherent strength to be able to experience and preserve localized deformation. The complex midlithospheric discontinuities could also be recording the stacking of subducting slabs (e.g., Bostock, 1998; Cooper et al., 2006), which also could provide material strong enough to deform in a localized manner; however, even buoyant oceanic lithosphere may not provide sufficient buoyancy to achieve long-term stability (Lee, 2006; Cooper et al., 2006). Regions that show more uniform midlithospheric discontinuities may be cratonic lithosphere that was formed from a weaker initial material that underwent viscous thickening and obtained its stability during postformation cooling. Those regions formed via viscous thickening would be more susceptible to weakening and destabilization due to metasomatism or other postformation modification.

C.M.C. was supported by NSF grant EAR-1112820 and M.S.M was supported by NSF grant EAR- 1054638. The authors would like to thank the anonymous reviewers for the helpful suggestions and timely edits.

1GSA Data Repository Item 2014100, Table S1, station locations with MLD signals from receiver functions mapped in Figure 1, is available at www.geosociety.org/pubs/ft2014.htm, or on request from editing@geosociety.org, Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301-9140, USA.