Visualizing the sedimentary response through the orogenic cycle: A multidimensional scaling approach

Sedimentary rocks contain denuded and transported remnants of precursor igneous, metamorphic, and sedimentary material. Characterizing the age composition of detrital minerals within sediments has become a widespread tool applied to a range of questions within earth sciences (Fedo, 2003). In particular, zircon is a refractory phase highly amenable to high-precision U-Pb geochronology. Dates from individual zircon grains typically reflect the magmatic crystallization age of a rock in its provenance, although dates indicating other zircon genesis events (e.g., metamorphism) also occur. The age spectra of detrital zircon in sedimentary deposits contain a snapshot of the geological evolution of the basin system at a given moment in time, frequently recording components from multiple crystalline sources, themselves derived from magmatic episodes in the hinterland (Thomas, 2011). Detrital ages may characterize a range of processes, including mixing of drainage systems and capturing of different basement blocks as uplift and erosion progress. Additionally, due to the resistant nature of zircon, multiple recycling events of older sedimentary deposits frequently leave an enduring legacy in the detrital record (Rainbird et al., 1992). Additionally, when considered over time, changing provenance trends reflect the local and potentially more distal geodynamic environment, which will have characteristic responses in the detrital record, capturing such information as terrane accretion, basement uplift, basin intercommunication, and reworking. The combination of the ubiquity of detrital zircon grains and the mechanism to analyze sufficient material, with precision, to characterize components from sediment has made detrital zircon geochronology the tool of choice in many provenance studies (Gehrels et al., 2011). The evaluation method of this age spectrum information is very important, and most studies rely on visual comparison of age information displayed via a range of graphical means, including U-Pb concordia plots, probability density diagrams, histograms, or pie charts. Some more quantitative approaches, including principal component analysis (Sircombe, 2000), age spectrum deconvolution (Sambridge and Compston, 1994), kernel functional estimates (Sircombe and Hazelton, 2004), and statistical descriptors and tests (Kirkland et al., 2012) have also been applied. However powerful these methods can be, additional significant geological information can be extracted through quantitatively evaluating the changing provenance signature through time. Multidimensional scaling (Vermeesch, 2013) is one such tool that has specific potential to further characterize the sedimentological processes associated with the orogenic cycle. To this end, we apply multidimensional scaling to three case study examples, two of which have a well-understood geodynamic framework (Albany-Fraser orogen and Grenville orogen) with sedimentary basins rooted on Archean cratonic crust that has undergone several cycles of orogenesis. In the final case study from the North Atlantic region, we explore the likely geodynamic environment as elucidated by detrital provenance evolution.

The Albany-Fraser orogen is located along the southern and eastern margins of the Archean Yilgarn craton (Fig. 1A; Spaggiari et al., 2011). The orogen formed during reworking of the Yilgarn craton, but it also importantly underwent phases of variable juvenile mantle addition, from at least 1800 Ma to 1140 Ma. Within the Albany-Fraser orogeny, at least two regionally extensive basin systems are preserved—the ca. 1815–1600 Ma Barren Basin and the ca. 1455–1305 Ma Arid Basin (Spaggiari et al., 2015; Waddell et al., 2015). Spaggiari et al. (2015) recently presented extensive secondary ion mass spectrometry (SIMS) U-Pb zircon geochronology from sediments of the Barren and Arid Basins. These data reveal that the Barren Basin was dominantly filled with Neoarchean zircon detritus with age and Hf isotopic signature identical to granitic rocks of the Archean Yilgarn craton (Kirkland et al., 2011, 2015a; Wyche et al., 2012). An additional local source of Paleoproterozoic to Mesoproterozoic detritus was derived from felsic magmatism of the 1815–1800 Ma Salmon Gums event, 1780–1760 Ma Ngadju event, and 1710–1650 Ma Biranup orogeny, with which it shares strong age and isotopic affinity (Kirkland et al., 2011; Spaggiari et al., 2011; Spaggiari and Tyler, 2014). Sedimentary deposition of the Barren Basin in the Biranup zone was initiated at a similar time to intrusion of granitic magmas at 1686 ± 8 Ma, possibly related to extensional processes within the Yilgarn craton edge (Spaggiari and Tyler, 2014; Spaggiari et al., 2015). The Biranup zone also includes Archean tectonic fragments, interpreted to reflect rifted lozenges of the Yilgarn craton now isolated with Proterozoic crust (Kirkland et al., 2011a; Spaggiari et al., 2011). The younger Arid Basin reflects a new cycle of deposition on the margin of the Yilgarn craton. The Arid Basin contains at least some similar detrital zircon cargo as the Barren Basin, but it also contains a ca. 1455–1375 Ma detrital age component that is unrecognized within the older Barren Basin. This age component is not found within any autochthonous unit within the orogen or elsewhere in the Yilgarn craton. However, similar ages and juvenile isotopic Hf signatures are recognized in the Madura Province to the east of the Albany-Fraser orogen (Spaggiari et al., 2012; Spaggiari and Smithies, 2015). Spaggiari et al. (2015) suggested that closure of a marginal basin, via east-dipping subduction, accreted a Madura Province arc terrane with oceanic affinity (Loongana arc) at ca. 1330 Ma (Smithies et al., 2015). This outboard terrane supplied the exotic 1455–1375 Ma juvenile detritus into the Arid Basin (Spaggiari et al., 2015; Kirkland et al., 2015b). Thus, in a basin dynamic model, the Barren Basin reflects a prolonged 1815–1600 Ma extensional phase where the basin was fed by zircon detritus derived from the Yilgarn craton and coeval extension-related magmatism within it. In contrast to this, the younger Arid Basin reflects the development of a marginal ocean basin on the edge of the Yilgarn craton that in part was fed by zircon detritus derived off an exotic arc system accreted onto the cratonic margin.

The eastern margin of Laurentia is composed predominantly of a series of accretionary orogens that developed during the Paleoproterozoic to Mesoproterozoic Eras. These include the orogenic intervals ca. 1850–1700 Ma (i.e., Central Plains, Yavapai, Penokean, Makkovik, Ketilidian orogens), ca. 1700–1600 Ma (i.e., Mazatzal and Labradorian orogens), and ca. 1500—1340 Ma (Granite-Rhyolite igneous province), which culminated during continental collision at ca. 1090–980 Ma between Laurentia and Amazonia, referred to as the Grenville orogeny (Fig. 1B; Rivers et al., 2012). The cratonic sedimentary basins associated with the Grenville orogeny record sedimentation during the assembly and breakup of the supercontinent Rodinia. The sedimentary successions deposited in these basins are classified as: (1) inboard sedimentary successions reflecting the active convergent margin predating the collisional Grenville orogeny, (2) intracratonic sedimentary successions deposited synchronously with the collisional assembly of Rodinia, and (3) basins formed following the assembly and during the breakup of Rodinia.

The sedimentary successions that predate Rodinia, found inboard of the convergent margin, for which detrital zircon age spectra are available include the Siamarnekh Formation (Wheeler, 1964; Spencer et al., 2015), Wakeham Group (Madore et al., 1998; van Breemen and Corriveau, 2005), Seal Lake Group (Reardon et al., 2009), Siamarnekh Formation (Spencer et al., 2015), Appalachian Inliers (Carrigan et al., 2003; Gates et al., 2004; Ownby et al., 2004), and Composite Arc Belt (Sager-Kinsman and Parrish, 1993; Friedman and Martignole, 1995; Wodicka et al., 1996; Corrigan and van Breemen, 1997). These successions show multimodal age spectra that include zircon grains derived from both the crystalline basement of the cratonic interior (>1500 Ma) but also younger populations associated with magmatism related to subduction zone(s) preceding the Grenville orogeny (ca. 1500 Ma to ca. 1200 Ma).

The Grenville orogeny (sensu stricto) spanned 1085 Ma to 985 Ma (Gower and Krogh, 2002). Preserved sedimentation that possibly occurred during this interval is very sparse, in part as a response to the compressional tectonic regime at this time. Units with a signature indicative of deposition during this period come from within the cratonic foreland and also within the orogenic interior. The Middle Run Formation reflects a cratonic foreland setting (Santos et al., 2002), whereas the Battle Harbour Psammite is an orogenic interior deposit (Kamo et al., 2011; Spencer et al., 2015). The Middle Run Formation is only exposed in drill core found in Ohio, within ∼50 km of the Grenville orogenic front (Santos et al., 2002), and it is interpreted to have been deposited in the foreland of the collisional orogeny (Cawood et al., 2007a). The Battle Harbour Psammite has a maximum depositional age of ca. 1045 Ma (Spencer et al., 2015) and was intruded by a 1024 ± 3 Ma pegmatite (Kamo et al., 2011). These units contain dominant detrital zircon age peaks between ca. 1300 Ma and ca. 1100 Ma, likely derived from the magmatic arc(s) that preceded the Grenville orogeny. Concomitant with the Grenville orogeny, sedimentation occurred within the Midcontinental rift of the Midwest region of the United States. Sedimentation was constrained to after major volcanic episodes at ca. 1085 Ma, although the youngest zircon ages in the rift extend to ca. 950 Ma (Craddock et al., 2013). The detrital zircon age spectra of the Midcontinent rift sedimentary successions are dominated by local Paleoproterozoic and Archean basement and underlying Mesoproterozoic volcanic rocks in basal units and transitions to greater proportions of Mesoproterozoic ages similar to that of the Middle Run Formation and Battle Harbour Psammite toward the top of the sequence. Sedimentary units of a similar age are found in West Texas along what is assumed to be the western extension of the Grenville orogeny (Spencer et al., 2014a), though the timing of magmatism and sedimentation differs from that of the Grenville orogeny (sensu stricto) (Mosher et al., 2008) and is not included in this study.

The sedimentary successions that formed after the assembly of Rodinia and those following the breakup of the supercontinent are traced along the length of the eastern Laurentian margin (Cawood et al., 2007a). These include successions deposited within failed intracratonic rift basins in the Eastern Blue Ridge of the Appalachians (e.g., Dahlonega gold belt, Sauratown Mountains window, and Smith River allochthon; Bream et al., 2004; Carter et al., 2006). This was followed by successions deposited during the initial stages of rifting and those recording the rift to drift transition, e.g., the Rome Formation of the Southern Appalachians (Thomas, 2004), the Unicoi, Erwin, Hardystron Formations and the Poughquag Quartzite of the central and north-central Appalachians (McLennan et al., 2001; Eriksson et al., 2004), and the Bradore Formation of western Newfoundland (Spencer et al., 2015; Cawood and Nemchin, 2001). The deposition of units along the passive margin of the Iapetus Ocean was followed shortly thereafter by the Taconic orogeny in the Middle Ordovician (McLennan et al., 2001).

In the circum–North Atlantic region, within and on the foreland of the Ordovician–Silurian Caledonian orogen, thick successions of Mesoproterozoic to Neoproterozoic sedimentary rocks are preserved (Fig. 1C; e.g., Cawood et al., 2007b). These deposits track the development of a series of successor basins rooted on the distal margin of the denuding Grenville orogen. Many of these deposits have been subjected to a range of Neoproterozoic compressional orogenic episodes. Latest Mesoproterozoic to Neoproterozoic sedimentary strata within the North Atlantic region have been dispersed through Caledonian orogenic events and are now found on and between East Greenland, Scotland, Ireland, Norway, the North Sea, and Svalbard. These sedimentary rocks chart three grand first-order cycles of sedimentation. The first cycle constrains early Neoproterozoic to late Mesoproterozoic Lithotectonic Group 1 units. These sediments are dominated by late Paleoproterozoic and Mesoproterozoic detrital zircon populations and an absence of Archean material. Granitic intrusive units within these units yield a range of zircon crystallization ages between 980 and 915 Ma (Strachan et al., 1995; Leslie and Nutman, 2003; Johansson et al., 2005; Kirkland et al., 2006, 2008; Myhre et al., 2009; Cutts et al., 2009). Examples of this lithotectonic group include the Krummedal of Greenland, Svaerholt succession in Norway, and terranes on Svalbard (Kirkland et al., 2007; Cawood et al., 2010). In the mid- to late Neoproterozoic, Lithotectonic Group 2 units were deposited and also characterized by minimal Archean detritus but dominant late Paleoproterozoic and Mesoproterozoic detrital zircons (e.g., Cawood et al., 2004; Dhuime et al., 2007; Kirkland et al., 2007). This lithotectonic group was originally suggested to record deposition in extension-related, intra-Rodinia basins (Cawood et al., 2004, 2007a), but more recently, it has been regarded as successor basins on Rodinia’s exterior margin with proximal oceanic crust and accreted outboard units (Kirkland et al., 2007; Cawood et al., 2010). The Moines, Grampians, and the Torridon Group of Scotland, the lower Eleonore Bay Supergroup of East Greenland, the Sørøy succession of Norway, and potentially other Scandinavian Caledonide tectonostratigraphic levels all can be ascribed to this package (Cawood et al., 2007b; Kirkland et al., 2007; Strachan et al., 2010; McAteer et al., 2010; Be’eri-Shlevin et al., 2011). These rocks are characterized by late Archean, late Paleoproterozoic, and Mesoproterozoic zircon detritus consistent with the removal of the Grenville orogenic barrier by this stage in the evolution of the North Atlantic region. Sedimentary sequences of late Neoproterozoic to early Paleozoic age define Lithotectonic Group 3, which reflect deposition during the main phase of Rodinia breakup. This package includes the upper Dalradian sequence in Scotland, Neoproterozoic to Cambrian siliciclastic sequences on the Caledonian foreland in Newfoundland, NW Scotland, and East Greenland (Cawood et al., 2007b), and the upper Eleonore Bay Supergroup in East Greenland (e.g., Cawood et al., 2003).

Latest Mesoproterozoic to mid-Neoproterozoic (1030–710 Ma) sedimentation and orogenic activity in the North Atlantic region have been considered in terms of the Valhalla orogeny (Cawood et al., 2010). The Valhalla orogen represents an accretionary orogen that developed along the margin of Laurentia on its free oceanic margin. Sedimentary successions within the orogen were stabilized during two episodes of tectonothermal activity in which Renlandian (980–910 Ma) activity was associated with Lithotectonic Group 1 and Knoydartian (830–710 Ma) orogenic events were associated with Lithotectonic Group 2. The Valhalla orogen is regarded as a distinct exterior system isolated from the interior Grenville orogen formed between collision of Laurentia, Baltica, and Amazonia.

The similarity of detrital zircon age spectra is frequently assessed by visual comparison. Largely, the object of this has been to ascertain the presence or absence of a particular zircon age component. There have been several attempts to quantify similarity between age spectra beyond binary qualification—that is, the presence or absence of an age population. Gehrels (2000) and Gehrels et al. (2002) proposed a set of metrics to measure both overlap and similarity between two samples, in essence a quantified visual comparison. Berry et al. (2001) utilized the Kolmogorov-Smirnov test, which provides a single-dimension statistic that calculates the maximum difference between two cumulative probability functions. At issue with the now popular probability density plot/distribution is the fact that the “probability” is not determined by data density but rather is a function of the precision of the analysis. This reduces the perceived importance of imprecise data despite these components potentially having a greater data density that arguably may have greater significance for provenance investigations (Vermeesch, 2012). We suggest that in order to avoid subjective comparisons between various age probability plots, intercomparison is best afforded through objective statistical methodologies.

In this study, we use a statistical tool referred to as multidimensional scaling (MDS; Vermeesch, 2013). This tool is based on a dissimilarity matrix for a suite of samples derived by using the D value from the Kolmogorov-Smirnov (K-S) test to create a two-dimensional “map” of points, where similar samples lie close together, and dissimilar samples lie far apart. This method is particularly useful when visualizing the dissimilarity of a large set of samples. The K-S test converts a detrital zircon probability spectrum to a cumulative density arrangement, which is the sum of probabilities with increasing age. The importance of the cumulative density arrangement is that it provides a means with which to compare strongly nonparametric distributions. The K-S test utilizes this projection of data and then statistically compares the maximum difference between two cumulative density functions, which is the maximum vertical difference between two curves, i.e., the D statistic.

The next step of the process uses a matrix of the D statistics (the dissimilarity matrix) to visualize dissimilarity in Euclidean space. As an example, we use the distances between cities from around the globe to represent dissimilarity (Table 1; Fig. 2A). MDS is used to arrange the data (e.g., location of cities) in a way in which the dissimilarity (or distance in this case) is represented in two dimensions. In the simple city example (e.g., a geographic map), there are essentially only two dimensions (latitude and longitude), but one can see that when there is a multitude of dimensions, such as different detrital age components, this approach will significantly simplify the complexity while still retaining the important age similarity/dissimilarity relationships between sedimentary samples. When MDS is applied to the dissimilarity matrix of city locations, their geographic locations are rearranged in such a way that mimics their true geographic locations (Fig. 2B). It is important to note that the orientation of the axes in MDS is arbitrary, as the distance between data remains constant if axes are rotated.

While this illustration is helpful to visualize the function of MDS, the real power of this technique comes when dealing with data containing multiple components (e.g., various zircon age populations). We demonstrate this using synthetic data of three age “sources” with normally distributed age populations (1000 ± 50 Ma, 2000 ± 50 Ma, 3000 ± 50 Ma) within four hypothetical samples, each with equal proportions of two of the three age populations (Fig. 3A). These samples, along with the hypothetical end members, are then plotted in Euclidean space using MDS (Fig. 3B) and the dissimilarity matrix of D statistics from the cumulative density functions (Fig. 3C).

Analyses from individual units were combined and treated as a single metric. For each case study, we also calculated synthetic unimodal age spectra with assumed normal distributions to provide an aid in the visual comparison of the MDS maps. This provides a vector on the figure that can be related to the increasing contribution of a given age component. The construction of the nonmetric MDS maps was performed via a modified MATLAB script (Vermeesch, 2013). MDS maps are presented in Figure 4. The dissimilarity matrices (K-S test D values), QQ plots, Shepard plots, and measure of stress (i.e., goodness of fit) are provided in the online supplementary tables and figures.¹ Only nonmetric MDS was used in this study because it is most appropriate for large data sets (see Vermeesch, 2013).

Holmes (1926) proposed that orogenic zones are affected by a cyclic repetition of events, part of a process called the orogenic cycle. Although his views have been significantly modified, the idea of some consistent sequence of events during orogenesis is appealing and expected to lead to characteristic responses in the first-order basin provenance signature. Convergent plate boundaries are characterized by the development of orogenic wedges and plateaus, which transition in later stages to gravitational collapse of previously thickened crust, producing intermontane and eventually oceanic backarc basins. Foreland and extensional sedimentary basins in the plate boundary region are filled by the erosional products of orogenic crust. These features portray the crustal orogenic cycle and are indicators of the thermal and mechanical evolution within the plate-boundary region. Although controlled by complex interactions between plate-tectonic, gravitational potential energy, and buoyancy forces, through this orogenic cycle, detrital mineral provenance trends will chart its progress. For example, unroofing processes after the orogenic thickening phase should readily be resolved through greater deep basement detritus incorporated into sedimentary packages. In contrast, one of the clearest indications of compressional processes is the arrival, and legacy within the sedimentary record, of exotic elements derived from outboard units. Thus, specific patterns on MDS plots may be expected to reflect the basin’s evolution and a signal from the original geodynamic setting. MDS plots have already been used to compare detrital zircon age spectra in sediments with similar depositional settings (e.g., loess deposits—Vermeesch, 2013; intrarift basin fluvial deposits—Spencer et al., 2014b). However, when viewed as an evolving system, provenance patterns visualized on MDS plots can be used to track the spatial and temporal location of a region through the orogenic cycle. We demonstrate this with the three case study examples as follows.

When comparing the sedimentary successions associated with the Albany–Fraser orogeny with the synthetic age spectra representative of basement rocks in the region (Fig. 4A), the MDS map reveals a transition from those rift-related units deposited in the on-craton Barren Basin to those units with a detrital signal derived from an arc-continent collision. This pattern is consistent with the “front end” of the idealized orogenic cycle as initially proposed by Wilson (1966), in that there is transition from rifting to subduction and ultimate accretion. The “back end” of the orogenic cycle, a process from subduction to collision to rifting, is seen in the foreland basin systems of the Grenville orogeny sensu stricto, in which there is a transition from subduction and continental collision to rifting and the formation of a passive margin (Fig. 4B; see McLelland et al., 1996; van Breemen and Corriveau, 2005). We contend that within these orogenic settings, the evolution of sedimentation can be characterized by increasing dissimilarity with time and therefore can provide an important metric through which the geodynamic regimes of ancient orogens can be deconvolved.

In stark contrast to case studies 1 and 2, the sedimentary successions along the periphery of the North Atlantic region reveal a MDS signature characterized by initial clustering progressing toward more variable provenance through time, which is reflected in a spreading MDS pattern. Although the siliciclastic-dominated sequences found along the Mesoproterozoic to Neoproterozoic margins of Laurentia and Baltica (Fig. 4C) are commonly invoked as representing stages in tectonism associated with the assembly and rifting of Rodinia (Cawood et al., 2007a, 2010, 2015; Kirkland et al., 2008; Bingen et al., 2011; Lorenz et al., 2012; Gee et al., 2014, 2015; Spencer et al., 2015), there is no clear pattern apparent in the dissimilarity of the detrital zircon age spectra through time. This is likely due to early deposits (Lithotectonic Groups sub-1 and 1) being characterized by a restricted provenance strongly linked to proximal basement rocks, whereas later deposits (Lithotectonic Groups 2 and 3) relate to much greater interbasin communication driven through sedimentary recycling processes. In this context, it is important to note that the Mesoproterozoic to Neoproterozoic history of the North Atlantic region is currently without consensus as to the tectonic setting and extent of orogenesis (e.g., Bingen et al., 2008; Cawood et al., 2010; Lorenz et al., 2012; Slagstad et al., 2013a, 2013b; Möller et al., 2013, 2015; Roberts and Slagstad, 2015; Gee et al., 2015). Some favor interpretations that include sedimentary components of the Mesoproterozoic to Neoproterozoic North Atlantic deposited in an interior position on a northward arm of the Grenville-Sveconorwegian orogen stretching into the High Arctic (Gee et al., 2015), while others favor either an exterior position on the distal margin of Rodinia (Cawood et al., 2015; Kirkland et al., 2007), or in even more outboard locations (Corfu et al., 2007). Although currently proposed tectonic models for the North Atlantic region are in flux, it is nonetheless clear that transcurrent tectonics among relatively small crustal fragments played an important role in the various episodes of Mesoproterozoic to Neoproterozoic orogenesis, given the increasing spreading pattern in the MDS map through time. Additionally, it is apparent that determining a unique Baltica or Laurentian signature for many units based on age correlations alone is futile given the strong similarity between specific lithotectonic groups across much of the North Atlantic region (Fig. 4C). Such detrital similarity reflects deposition during the time frame of Rodinia, the paleogeography and geodynamics of which strongly controlled many aspects of Earth’s geological record at this time (Van Kranendonk and Kirkland, 2013; Spencer et al., 2013). Hence, the question of outboard versus endemic derivation of many North Atlantic region units reflects the positioning of basins after rifting of Iapetus, as the Rodinia margin was not broken into discrete crustal fragments until after this point in time. We propose that it is due to these complexities, associated with loosely spaced cratonic margins, that no coherent pattern of increasing dissimilarity is seen through time in the North Atlantic region.

Detrital zircon populations through time record important information on basin dynamics that are not readily explored by simple comparison of probability or kernel density plots. However, MDS provides a useful tool in this regard because it allows visualization of the changing detrital zircon provenance through time. The case studies presented herein demonstrate characteristic differences in the evolutionary response of the detrital zircon record to the geodynamic environment. Specifically, unambiguous age progressive trends are noted in the Albany-Fraser and Grenville (sensu stricto) orogens due to the availability of progressively younger detritus through the orogenic cycle. In contrast, the North Atlantic region displays clustered to anticlustered distributions through time, reflecting greater basin promiscuity through time in the distal depositional environments of a complex orogenic system.

We thank Brendan Murphy, Pieter Vermeesch, and an anonymous reviewer, whose reviews greatly improved the clarity of this paper.