The assembly of Rodinia involved the collision of eastern Laurentia with southwestern Amazonia at ca. 1 Ga. The tectonostratigraphic record of the central Andes records a gap of ∼300 m.y. between 1000 Ma and 700 Ma, i.e., from the beginning of the Neoproterozoic Era to the youngest part of the Cryogenian Period. This gap encompasses the time of final assembly and breakup of the Rodinia supercontinent in this region.

We present new petrographic and whole-rock geochemical data and U-Pb ages combined with Hf isotope data of detrital zircons from the volcanosedimentary Chilla beds exposed on the Altiplano southwest of La Paz, Bolivia. The presence of basalt to andesite lavas and tuffs of continental tholeiitic affinity provides evidence of a rift setting for the volcanics and, by implication, the associated sedimentary rocks. U-Pb ages of detrital zircons (n = 124) from immature, quartz-intermediate sandstones have a limited range between 1737 and 925 Ma. A youngest age cluster (n = 3) defines the maximum depositional age of 925 ± 12 Ma. This is considered to coincide with the age of deposition because Cryogenian and younger ages so typical of Phanerozoic units of this region are absent from the data.

The zircon age distribution shows maxima between 1300 and 1200 Ma (37% of all ages), the time of the Rondônia–San Ignacio and early Sunsás (Grenville) orogenies in southwestern Amazonia. A provenance mixing model considering the Chilla beds, Paleozoic Andean units, and data from eastern Laurentia Grenville sources shows that >90% of the clastic input was likely derived from Amazonia. This is also borne out by multidimensional scaling (MDS) analysis of the data.

We also applied MDS analysis to combinations of U-Pb age and Hf isotope data, namely εHf(t) and 176Hf/177Hf values, and demonstrate again a very close affinity of the Chilla beds detritus to Amazonian sources. We conclude that the Chilla beds represent the first and hitherto only evidence of Rodinia breakup in Tonian time in Andean South America.


The plate tectonic evolution of Earth at least since the beginning of the Proterozoic is marked by periods in which most of the continental fragments amalgamated episodically to form supercontinents (e.g., Nance et al., 2014; Evans et al., 2016). The processes and events recording the amalgamation and dispersal of the supercontinent Rodinia between 1200 and 700 Ma are relatively well established in broad terms (McMenamin and McMenamin, 1990; Hoffman, 1991; Pisarevsky et al., 2003; Cawood et al., 2016; Evans et al., 2016). However, precise data on the onset and progress of rifting and dispersal events are scarce (Evans et al., 2016).

The Sunsás orogen in southwestern Amazonia and the Grenville orogen of eastern Laurentia (present coordinates) formed by the collision of the two continents during the assembly of Rodinia (Fig. 1; Tohver et al., 2005, 2006; Li et al., 2013; Cawood et al., 2016; Evans et al., 2016). Onset of actual rifting leading to final separation of the two continents is recorded in eastern Laurentia at between 750 and 600 Ma (Fig. 1; Allen et al., 2010; Burton and Southworth, 2010; McClellan and Gazel, 2014).

At the southwestern margin of Amazonia, there is no known record of rifting related to Rodinia dispersal. After the end of the Sunsás orogeny at ca. 1000 Ma, there is a hiatus of ∼300 m.y. represented by an absence of known outcrops of any type of rock with an established Tonian age (1000–720 Ma) in the central Andes of southern Peru, Bolivia, northwestern Argentina, and northern Chile (e.g., Turner, 1970; Suárez Soruco, 1992, 2000; Ramos 2000, 2008; Hervé et al., 2007; Fig. 2). Here, the oldest dated Neoproterozoic and <900 Ma event is represented by a dacitic dike in the Arequipa massif of northern Chile, which gave a lower intercept age of 635 ± 5 Ma (Loewy et al., 2004).

The fact that Laurentia separated from Amazonia and the presence of rift-related successions on the eastern margin of Laurentia make it very likely that evidence of the dispersal of Rodinia may also be found at the southwestern margin of Amazonia.

In this contribution, we present results of a study of a hitherto enigmatic volcanosedimentary unit informally named the Chilla beds located on the Bolivian Altiplano west of the Bolivian capital La Paz. The Chilla beds comprise basaltic and andesitic lavas and tuffs in association with siliciclastic sandstones and conglomerates. From our new data on the petrography and geochemistry of the lavas, we conclude that they are the product of extension-related magmatism. New U-Pb age and Hf isotope data on the Chilla beds sandstones allow us to constrain the provenance of the clastic units to sources on Amazonia and to derive a Tonian age of the formation of the Chilla beds. We conclude that the Chilla beds represent the first, and presently only, evidence of extensional processes in early Neoproterozoic time at the southwestern margin of Amazonia. Our limited first data set provides important new insights further elucidating the dispersal history between Laurentia and Amazonia during Rodinia breakup.


The evolution of the central Andean region during the Neoproterozoic Era prior to the beginning of the Ediacaran Period at ca. 635 Ma is largely unknown. The Tonian and Cryogenian Periods encompass the time between final amalgamation of the Rodinia supercontinent between 1000 and 900 Ma (e.g., Cawood et al., 2016) and the initiation and activity of the Pampean magmatic arc and subduction system in central Argentina between 650 and 530 Ma (Fig. 1; Rapela et al., 1998; Escayola et al., 2007, 2011).

The transition to Cryogenian System rocks comprises two local occurrences of A-type granitoids in southern central Peru dated at (1) 751.7 ± 8.1 and 691 ± 13 and (2) 752 ± 21 Ma, respectively (Mišković et al., 2009). Locally limited deposition of the glacial siliciclastic rocks of the Chiquerío Formation occurred on the Arequipa massif in southern Peru at ca. 700 Ma (Chew et al., 2007a) and in northern Chile with the presumably glacial Limón Verde diamictites (Fig. 2; Morandé et al., 2012, 2018). The Limón Verde diamictites include granite clasts dated at 1040 Ma (Pankhurst et al., 2016).

During the amalgamation of Rodinia, southwestern Amazonia collided with eastern Laurentia (Tohver et al., 2005, 2006; Cawood et al., 2016). During the late Neoproterozoic Sunsás orogeny (Litherland et al., 1989; Cordani and Teixeira, 2007; Santos et al., 2008), at ca. 1000 Ma, southwestern Amazonia also experienced re-accretion of the native Paraguá terrane of eastern Bolivia and the collisional transfer of the Arequipa terrane now found in coastal southern Peru and northern Chile (Loewy et al., 2004; Boger et al., 2005). One of the few sedimentary units marking the Mesoproterozoic-Neoproterozoic transition in the central Andes is represented by metasedimentary rocks at Sierra de Moreno in northern Chile. This unit records a maximum depositional age of 1140 Ma (Fig. 2; Pankhurst et al., 2016). Its stratigraphic age and potential connection to Rodinia assembly and terrane transfers still needs to be established.

Circumstantial evidence derived largely from other continents suggests that the margin of southwestern Amazonia experienced a phase of protracted rifting during Tonian and Cryogenian time in the course of Rodinia breakup (Fig. 1; Cawood et al., 2016). Starting at ca. 650 Ma, the Pampean active margin developed in central and northwestern Argentina (Fig. 1). It was connected to the opening and closing of the marginal Puncoviscana basin in late Ediacaran and early Cambrian time (Rapela et al., 1998; Ježek et al., 1985; Keppie and Bahlburg, 1999; Zimmermann, 2005, 2018; Escayola et al., 2007, 2011). There is no direct evidence of the Pampean magmatic arc in either Bolivia or southern Peru. However, on the basis of detrital zircon age data, Chew et al. (2008), with support from Naidoo et al. (2016), argued for a mainly Ediacaran magmatic arc active between 650 and 550 Ma (Fig. 2), which is now buried beneath the eastern part of the central Andes and/or the western Amazon Basin.

Indirect evidence of Tonian and Cryogenian magmatic and metamorphic activity in the region of the central Andes is provided by detrital zircon U-Pb age data. Data sets from very late Neoproterozoic and Paleozoic siliciclastic sedimentary rocks contain between 18% and 35% ages between 1000 and 600 Ma (e.g., Chew et al., 2007b, 2008; Bahlburg et al., 2009, 2011; Reimann et al., 2010; Adams et al., 2011; Augustsson et al., 2015). These data bear witness to a plate-tectonic evolution not yet recognized or observed in outcrop.


The Chilla beds are present on the northern Bolivian Altiplano in the Serranía Chilla ridge to the southwest of the capital city La Paz (Figs. 3A, 3B). The highest peak of the ridge is formed by Cerro Chilla at an elevation of 4824 m, ∼800 m above the Altiplano surface. In the main outcrop area in the central-western part of the ridge, the Chilla beds are in tectonic contact with Cretaceous–Neogene sedimentary rocks of the Molina Formation. The bounding faults delimit a small lenticular popup structure along inward-dipping reverse faults (Fig. 3C; Zapata, 1992).

The Chilla beds consist of alternating volcanic and epiclastic sedimentary rocks. Lavas are low-grade metamorphic metabasalt occurring as massive and vesicular flows or as stacks of pillow lavas up to 6 m high (Figs. 4, 5B–5D). They are associated with epiclastic conglomerates (predominantly sedimentary and less-abundant granite, gneiss, and vein-quartz clasts) and graded sandstone interpreted as turbidites (Paton, 1990). Very fine-grained altered metatuffs were originally classified as phyllite (Figs. 4, 5; Paton, 1990; Zapata, 1992; Matos et al., 2000, 2002). The succession is very distinct from the surrounding and overlying weakly deformed and altered Devonian to Neogene rocks due to its marked and tight disharmonic folding, a cleavage well developed in the fine-grained rocks, and the low-grade metamorphic overprint (Paton, 1990; Matos et al., 2002). The basement to the Chilla beds is unknown.

Geochemical analysis of major and some trace elements suggested to Paton (1990) that the metabasalts have an affinity to calc-alkaline arc basalts. Permian whole-rock K-Ar age dates are between 294 and 278 Ma and were interpreted as defining the age of folding and metamorphism (Paton, 1990). Because the Chilla beds are regionally overlain by a conformable and much less-deformed succession of Devonian to Neogene sedimentary rocks (Paton, 1990), their stratigraphic age was recognized to be pre-Devonian. These field relations cast serious doubt on the interpretation of Permian K-Ar ages as the age of folding and metamorphism. Witschard (1992) speculated that the unit may represent an ophiolite complex in a major suture zone without providing any evidence in support. Outcrop observations together with the association of partly pillowed metabasalts with relatively quartz-rich sandstones and conglomerates led Matos et al. (2002) to assign this rock assemblage to an extensional, rift-like basin on continental crust. Díaz Martínez (1996) noted that none of the typical lithologies of an ophiolite are present in the Chilla beds, and correlated the unit with siliciclastic Neoproterozoic to Ordovician units in northwestern Argentina and the volcaniclastic Ollantaytambo Formation in southern Peru. The Ollantaytambo Formation is now recognized to be of Early Ordovician age (Carlotto et al., 1996; Bahlburg et al., 2006, 2011).


We apply a set of petrographical, geochemical and geochronological methods to samples of mafic lavas, tuffs, and epiclastic sandstones collected from the Chilla beds on the Bolivian Altiplano. We surveyed the outcrop area (Fig. 3) of the tightly folded rocks. We sampled the major lithologies for further analysis along a section where they appeared the least altered. We use geochemistry to elucidate the nature and tectonic setting of the magmatic rocks and the degree of alteration of the sandstones. Microscopy and electron microscopy permit the determination of the petrographical composition of the rocks. U-Pb geochronology of detrital zircon allows for establishing maximum depositional ages and the geochronological structure of the provenance signal stored in the clastic rocks. We use the Hf isotopic data here as a further means to constrain the provenance of the detritus contained in the Chilla beds.

Field Emission Gun Scanning Electron Microscopy (FEG-SEM) Coupled with Energy Dispersive Spectroscopy (EDS)

Field emission gun scanning electron microscopy (FEG-SEM) analyses were performed at the University of Stavanger (Norway) with a Zeiss Supra VP 35. Fifteen (15) polished thin sections from field samples were used to identify mainly strongly weathered minerals and imaged for textural characteristics. The samples were coated with palladium, and parameters of the microscope were set to an acceleration voltage of 15kV, 60 µm aperture, and a working distance between 10 and 11 mm. High current setting was used. For imaging of the samples, the secondary electron detector was used, while a backscattered electron detector was used to provide an image with different grayscale values based on the average atom number of each mineral. For energy dispersive X-ray spectroscopy (EDS), providing semiquantitative elemental analyses, an EDAX detector was used. To calibrate the EDS, an Island spar calcite crystal was used, along with standard dolomite, alkali feldspar, and plagioclase (provided by Astimex). Chemical data are not based on a standardized procedure.

U-Pb Geochronology

Four sandstone samples (Fig. 4) of ∼3 kg weight of different grain sizes were crushed, sieved, and processed. The basaltic tuffs consist mainly of clay minerals and a minor silt component; they did not yield any zircons. The grains <250 µm underwent routine magnetic and heavy-liquid separation using a Frantz magnetic separator and sodium polytungstate. The number of zircons obtained from each sample was relatively small. A combined total of 170 zircons from four samples were hand-picked onto epoxy mounts and polished for cathodoluminescence (CL) imaging. CL images were obtained using a scanning electron microscope JEOL 6610 to image the internal structure and zonation patterns of zircons (Fig. 6). Grains characterized by oscillatory or sector zoning are interpreted as of magmatic origin; unzoned grains or those with irregular to concentric and planar zoning are considered to be metamorphic in origin (Corfu et al., 2003; Wu and Zheng, 2004).

Zircons were analyzed by laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) at the Institute for Mineralogy at the University of Münster (Germany) using a ThermoFisher Element2 mass spectrometer coupled to a Photon Machines Analyte G2 Excimer laser. Whenever possible, cores and rims, if developed, were analyzed. Laser spots of 35 µm in diameter were used for grains with diameters >50 µm; to those <50 µm, a spot diameter of 25 µm was applied. Groups of 10 unknowns were bracketed with three analyses of the GJ-1 reference zircon (609 Ma; Jackson et al., 2004) for external calibration. To monitor precision and accuracy over the course of this study, the 91500 reference zircon (1065 Ma, 206Pb/238U age; Wiedenbeck et al., 1995) was measured as an unknown. Results of 91500 zircon analyses yielded a weighted mean average age of 1065.1 ± 5.5 Ma (n = 11; Fig. S1 in Table S11) that matches published values (Wiedenbeck et al., 1995).

Zircon analyses have an average relative uncertainty of <2.3% (2σ); analyses with relative errors >5% (2σ) and an elevated content of common lead (≥2%) were excluded. The latter criterion applied only to one analysis. Data reduction was performed with an in-house program using Microsoft Excel (Kooijman et al., 2012). Error propagation conforms to the standards outlined in Horstwood et al. (2016). We performed common 204Pb-based correction after Stacey and Kramers (1975). Zircon age data were analyzed using Isoplot/Ex 3.76 software (Ludwig, 2012), and age distributions are displayed as kernel density estimate and cumulative frequency plots using the provenance software of Vermeesch et al. (2016). Statistical evaluation of the data has been performed using the provenance software of Vermeesch et al. (2016) and the DZstat and DZmix tools by Sundell and Saylor (2017).

We uniformly applied a concordance filter of 90% and 101% to all our data (Table S1 [footnote 1]), with 124 U-Pb ages fulfilling this criterion. For zircons >1.6 Ga, the 207Pb/206Pb ages are used, and to those <1.6 Ga, preference is assigned to the 206Pb/238U ages. Here we follow Spencer et al. (2016), who evaluated the error dimensions of 38,000 published zircon ages and concluded that the crossover point from 207Pb/206Pb ages to 206Pb/238U ages should be placed at 1.5 Ga. We combine this with the approach of Gehrels et al. (2008, p. 8), who advised that the crossover point should “not artificially divide a cluster of analyses”. In view of the clustering of ages in our analyses, we placed the crossover point at 1.6 Ga.

Lu-Hf Isotope Analysis

We selected 53 of the dated grains for Hf isotope analysis, focusing on the main U-Pb age clusters. The Lu-Hf isotope analyses were performed by multicollector (MC) ICP-MS using a Thermo-Scientific Neptune at Goethe University Frankfurt (GUF, Germany) coupled to a New Wave Research UP-213 laser system with a teardrop-shaped, low-volume laser cell, following the method of Gerdes and Zeh (2006, 2009). The Lu-Hf laser spot of 40 µm diameter was drilled on top of or directly next to the U-Pb laser spot in the same zone of zonation. Multiple analyses of Lu- and Yb-doped JMC 475 Hf standard solutions show that results with a similar precision and accuracy can also be achieved if Yb/Hf and Lu/Hf is 5×–10× higher than in most magmatic zircons. All data were adjusted relative to the JMC 475 176Hf/177Hf ratio of 0.282160, and quoted uncertainties are quadratic additions of the within-run precision and the reproducibility of the 40 ppb JMC 475 solution (2 standard deviation = 0.0033%, n = 10 per day). Accuracy and daily reproducibility of the method was verified by repeated analyses of reference zircons GJ-1 and Temora, which yielded a 176Hf/177Hf of 0.282009 ± 0.000021 (2 SD, n = 30) and 0.282670 ± 0.000025 (n = 23), respectively. This is in perfect agreement the LA-MC-ICP-MS long-term average of GJ-1 (0.282010 ± 0.000025; n > 800) and Temora (0.282483 ± 0.000023, n > 250) reference zircons at GUF. The reproducibility of 0.007%–0.008%, which corresponds to 0.7–0.8 epsilon units, is similar or slightly better than the quoted uncertainties (mean = 0.95 epsilon units) of the detrital zircon grains in this study.

The initial 176Hf/177Hf values are expressed as εHf(t), which is calculated using a decay constant value of 1.867 × 10−11 yr−1 (Scherer et al., 2001; Söderlund et al., 2004), and for the chondritic reservoir, a 176Hf/177Hf and 176Lu/177Hf of 0.282785 and 0.0336, respectively (Bouvier et al., 2008). For the calculation of Hf two-stage model ages (TDM), a 176Lu/177Hfcrust,today value of 0.0093 (Amelin et al., 2000) was used for the average continental crust, and for depleted mantle (DM), 176Lu/177HfDM,today and 176Hf/177HfDM,today values of 0.03813 and 0.283224, respectively (Vervoort et al., 2000), were used. Hf-isotope data from the literature were recalculated using the parameters above for better comparison.

The εHf(t) notation of Hf isotope data results in positive or negative values, which refer to the difference to a value of 0 for the chondritic unfractionated reservoir (CHUR), reflecting the composition of the bulk Earth. There is an increase over time in the difference between the 176Hf/177Hf ratios of the DM reservoir and of the continental crust, and εHf = 0 does not correspond to a constant 176Hf/177Hf ratio. An εHf = 0 in a 3 Ga zircon corresponds to a crust that formed from the DM ∼0.3 b.y. earlier, corresponding to a difference ∼6 εHf units between zircon and the DM composition at that time. In contrast, an εHf = 0 in a 1 Ga zircon indicates growth in a crust that evolved for ∼0.65 b.y., corresponding to ∼12.5 εHf units difference. Analogous to the approach of Reimann et al. (2010), we consider residence times since separation from the depleted mantle of <0.3 b.y., equaling ∼6 εHf units, as juvenile. Residence times between ∼0.3 and 0.7 b.y. result in differences between ∼6 and ∼12.5 εHf units and are considered moderately juvenile. Larger differences from the depleted mantle are considered evolved.

Hafnium depleted-mantle model ages (HfTDM; Table S2 [footnote 1]), or average crustal residence times for the zircons, give a qualitative estimate of the time of separation of the source rocks, or their precursors, from a hypothetical depleted-mantle reservoir. While depleted-mantle model ages do not necessarily provide any real age information, they, along with εHf(t) values, are nonetheless useful in identifying older crustal versus juvenile mantle components (e.g., Bennett and DePaolo, 1987; Vervoort and Blichert-Toft, 1999; Vervoort et al., 1999). For the calculation of these zircon Hf model ages, it is generally assumed that the crust formed from a depleted mantle reservoir of known composition and had a known and uniform 176Lu/177Hf until the respective zircon, representing this crust, crystallized (Bouvier et al., 2008).


Ten basalt, one dacite, five tuff, and 12 sandstone samples were powdered in an agate mill for analysis of their whole-rock element geochemical compositions. Geochemical data were obtained by ICP-MS analysis at Acme laboratory (Vancouver, Canada; Table S3 [footnote 1]). Details for the analytical method and processing can be found at http://acmelab.com, but are compiled here. The milled sample was mixed with LiBO2/Li2B4O7 flux in crucibles and fused in a furnace. The cooled bead was dissolved in American Chemical Society–grade nitric acid and analyzed by ICP-MS. Loss on ignition (LOI) was determined by igniting a sample split then measuring the weight loss. A 1 g sample was weighed into a tarred crucible and ignited to 1000 °C for one hour, and then was cooled and weighed again. The loss in weight is the LOI of the sample. Total carbon and sulfur were determined by the LECO method. Here, induction flux is added to the prepared sample, then ignited in an induction furnace. A carrier gas sweeps up released carbon to be measured by adsorption in an infrared spectrometric cell. Results are total concentrations and attributed to the presence of carbon and sulfur in all components. An additional 14 elements were measured after dilution in aqua regia. The prepared sample was digested with a modified aqua regia solution of equal parts concentrated HCl, HNO3, and DI-H2O (DI stands for deionized water) for one hour in a heating block or hot-water bath. The sample volume was increased with dilute HCl solutions, and splits of 0.5 g were analyzed. None of the measured concentrations was far above the possible detection limit, but in standard range, and accuracy and precision were between 2%–3%.


Petrographic Compositions

Volcanic Rocks: Lavas

The lavas are composed of a fine groundmass in which medium-sized plagioclase needles (100–150 µm, mainly andesine according to EDS, with 5%–9% Ca) with a trachytic texture occur together with chlorite (Fig. 7A). In some samples, plagioclase crystal sizes are large (>200 µm), forming a porphyritic texture. Albite is strongly weathered. Hornblende occurs as smaller-grained blocky crystals (50–100 µm) characterized by Ca, Fe, and Mg (Fig. 7A). Some hornblende crystals contain traces of Na and K as well as Ti (according to EDS) and can be classified as kaersutite (Fig. 7A). Elongate orthopyroxene needles are classified as enstatite-ferrosilite, with higher amounts of Fe compared to Mg, hence ferrosilite (Fig. 7A). Large plagioclase crystals are mainly bytownite and andesine (6%–13% Ca; Fig. 7B). In the larger grains, exsolution can be observed toward K-feldspar (Fig. 7B). However, most of the plagioclase is strongly weathered and albitized, epidotized and/or transformed to calcite. Characteristic accumulations of illite are interpreted as former glass fragments. Accessory minerals are chalcopyrite (Fig. 7B), chlorite, magnetite, and small grains of K-feldspar (Fig. 7B). Epidotization occurs along fractures in plagioclase minerals.

Volcanic Rocks: Tuffs

The tuffs are strongly folded and foliated. The original texture is destroyed. Abundant dark minerals represent glass weathered to illite. The rocks contain silt-sized albite, small quartz grains (1–10 µm), K-feldspar, and notable amounts of phlogopite (according to EDS; Fig. 7C). Silt-sized albite displays slight rotation in some samples. Accessory minerals are apatite, rutile, iron oxide, and Ti-rich ilmenite. Some of the 10–30 µm rutiles (Fig. 7C) are strongly enriched in Nb (≤5.5%) and Ta (≤1.6%) but depleted in Cr (<0.2%) indicating, a derivation from felsic sources (Zack et al., 2004; Meinhold et al., 2008).


The sandstones are mostly poorly sorted and structureless. Some beds show coarse-tail or distribution grading but no structures from which current directions could have been derived. The sandstones commonly include granule-sized angular to well-rounded monocrystalline quartz grains, polycrystalline quartz, microcline, albite, sericitized plagioclase, and epidotized feldspar, broken metasedimentary clasts, and large muscovite grains. The matrix is ∼15%, classifying some of the rocks as wackes. Most of the matrix contains quartz, calcite, albite, illite, and some chlorite. One sample (CC20) contains a potentially secondary, minor calcite cement. Other samples are less rich in matrix and represent arenites (Fig. 7D). These rocks are better sorted and contain albite and quartz. Both K-feldspar and andesine are commonly rounded. Rutile is common; it has very low Cr contents (below detection limit) and is not enriched in Nb and Ta as observed in the tuffs. Slightly rounded fluorine apatite is rare.

The sample with the highest silica content (80.68 wt%, sample CC30; Table S3 [footnote 1]) represents the matrix of a conglomerate. It contains mainly large quartz and rounded to subrounded microcline clasts. Quartz shows resorption embayments, typical for felsic volcanic rocks. Exsolution of plagioclase to K-feldspar is present, similar to minerals in the basaltic lavas. Accessory minerals in the sandstones and conglomerates are magnetite, zircon (Fig. 7D), pyrite, and unspecified iron oxides.


Lavas and Tuffs

The lavas in the Chilla beds mostly have a subalkaline basalt to basaltic andesite composition with SiO2 between 47 and 51 wt%; however, one sample is a dacite with SiO2 of 65 wt% (Table S3 [footnote 1]). The SiO2 contents of the tuffs range between 53 and 60 wt%, which classifies them as basaltic to andesitic (Cox et al., 1979). Low K/Cs ratios between 1500 and 6000 indicate low degrees of weathering and alteration of the tuffs (Table S3). For comparison, the upper continental crust (UCC; McLennan, 2001) has a value of 6117 (Table S3).

The abundances of selected major and trace elements like Ti, Nb, Y, and Zr identify the rocks as subalkaline basalts and, to a lesser extent, as andesites (Fig. 8A; Pearce, 1996). According to their Zr and P2O5 abundances, their Nb/Yb and TiO2/Yb ratios, and AFM (A [Al2O3], F [FeO], and M [MgO]) systematics, the lavas and tuffs are classed as tholeiitic (Figs. 8B–8F; Irvine and Baragar, 1971; Winchester and Floyd, 1976; Cabanis and Lecolle, 1989; Pearce, 2008).

Chondrite-normalized element patterns of lavas and tuffs are rather flat and uniformly show a moderate enrichment of the light rare earth elements (LREEs), with LaN/YbN (N stands for chondrite normalized values) ratios between 2 and 6 for the lavas and between 3 and 10 for the tuffs, and a weak or absent negative Eu anomaly (Fig. 9A; Table S3 [footnote 1]) expected for plagioclase-rich rocks. The trace element patterns of lavas and tuffs lack the Nb and Ta anomalies typical of magmatic arc rocks (Fig. 9B).

Considering the immobile trace and rare earth elements in particular, the composition of the rocks is akin to enriched mid-ocean-ridge basalt (E-MORB) (Sun and McDonough, 1989). Relative to E-MORB, the rocks are slightly enriched in Th, which indicates some crustal assimilation (Fig. 8E; Pearce, 2008). In Nb/Yb versus TiO2/Yb space, the data array straddles the boundary between the MORB array linked to shallow melting and the oceanic island basalt array of deep melting. This is considered indicative of conditions linked to a thinning lithosphere, which may be connected either to plume-ridge interactions or volcanic rifted margins (Pearce, 2008). The concentrations of La, Yb, and Nb suggest that the Chilla volcanic rocks were derived from melting of a garnet peridotite (Liu et al., 2018). The geochemical element distributions of the Chilla lavas are also quite similar to those of tholeiitic within-plate basalts of Cenozoic age found along the Cameroon line (Pearce, 1982; Condie, 1982; Asaah et al., 2015; Figs. 8D, 9C), which dissects oceanic as well as continental parts of the western African plate.


The sandstones and wackes associated with the Chilla volcanics have SiO2 values mostly between 67 and 74 wt% (Table S3 [footnote 1]) and correspond to quartz-intermediate sandstones of Crook (1974). One sample with SiO2 = 81 wt% is quartz rich.

The chemical index of alteration (CIA; Nesbitt and Young, 1982; Bahlburg and Dobrzinski, 2011) is a valuable tool to gauge the effects of weathering and alteration on the composition of siliciclastic sedimentary rocks. In CN-A-K [molecular proportions of (CaO* + Na2O)–Al2O3–K2O] space, the samples of the Chilla sedimentary rocks display a relatively flat trend toward the A-K conode (Fig. 10A). This trend deviates from natural weathering trends (Nesbitt and Young, 1984; von Eynatten, 2004) and away from the post-Archean average Australian shale composite (PAAS; Taylor and McLennan, 1985) toward the K apex. It likely reflects the secondary addition of K through the conversion of aluminous clay minerals to illite or the transformation of plagioclase to K-feldspar through metasomatic processes (Fedo et al., 1995). This results in relatively low CIA values averaging 60 (Fig. 10A). If restored to the natural weathering trend, corrected values would lie between 60 and 80. High K/Cs ratios gauge the degree of weathering of clay minerals in siliciclastic sedimentary rocks (McLennan et al., 1990). Sandstone values between 15,000 and 50,000 indicate significant weathering and alteration at least regarding the clay mineral fraction (Table S3 [footnote 1]).

Using ratios of high field strength elements in the Zr/Sc versus Th/Sc diagram of McLennan et al. (1993) is a good tool to evaluate the degree of zircon concentration through sedimentary recycling by comparing the Zr abundance to tracers of mantle (Sc) and crustal (Th) components (Fig. 10B). The UCC (McLennan, 2001) has Zr/Sc and Th/Sc ratios of 13.97 and 0.79, respectively. The values of the Chilla sandstones cluster around averages of 29.74 and 1.28 and are thus both only mildly elevated relative to the UCC (Fig. 10B). They fall well within the fields of other upper crustally sourced sandstones of variable late Proterozoic and Paleozoic age in the central Andes (Fig. 10B). This coincides with the petrographic observations indicating slight reworking and the presence of a rather high amount of feldspar.

In view of the association of the Chilla beds sandstones with mafic lavas of continental tholeiite affinity, it is helpful to test the sedimentary rocks for respective signatures of ophiolite sources not necessarily apparent from microscopic inspection (Bock et al., 1998; Reimann Zumsprekel et al., 2015). The combination of a Cr/V ratio >10 with low Sc/Th (<0.2) and Y/Ni (<1) ratios would point to the presence of ophiolitic detritus (Floyd et al., 1990; McLennan et al., 1993). Similarly, a Cr/Th ratio >40 would indicate the presence of ultramafic detritus like chromite as an erosional product from ophiolites (Floyd et al., 1990). In comparison, the average value of the Cr/Th ratio is 3.9 in the studied rocks (Table S3 [footnote 1]). The average ratio in the UCC is 7.8 (McLennan, 2001). None of the Chilla beds samples is characterized by such source components. The data show that the Chilla sandstones do not contain detritus derived from ophiolites or ultramafic rocks.

The chondrite-normalized trace element patterns of the Chilla sandstones show a marked Nb and Ta anomaly typical of magmatic arc rocks and the UCC (Fig. 9B). This, together with the other presented data, indicates that the Chilla sandstones were derived from upper crustal sources outside and beyond the mafic volcanic part of the basin. The detritus does not show effects of significant recycling (Fig. 10B).

U-Pb Geochronology

Many of the zircons obtained from the Chilla beds sandstones were broken. The CL images revealed that fractures and other defects are common (Fig. 6). Many grains are metamict and display complex zoning patterns. Zircons are rounded to variable degrees. Idiomorphic grains, or broken pieces of such grains, are common. Metamorphic zircons represent <10% of the analyzed grains.

From 170 analyzed detrital zircons, we obtained 124 concordant U-Pb ages (Table S1 [footnote 1]). The ages range between 1750 and 800 Ma with 53% of all ages between 1300 and 1140 Ma and the main age mode at 1180 Ma (Figs. 11, 12A, 12B). The next-most abundant group occurs between 1120 and 1000 Ma (20% of all ages) with its main age mode at 1030 Ma (Fig. 12A). The youngest age of 796.6 ± 9 Ma was obtained from a slightly discordant single zircon grain (discordance = 6.6%, common Pb corrected). The concordia age of the youngest age cluster (n = 3) is 925 ± 12 Ma (Fig. 11). The oldest single zircon age is 1736.5 ± 37.7 Ma (Fig. 11). We obtained seven pairs of core and rim ages, which were similar within error for each pair (Table S1).

Hf Isotope Systematics

We analyzed a total of 50 dated zircons for their Hf isotope systematics (Table S2 [footnote 1]). The U-Pb ages of the analyzed zircons range from 1.7 Ga to 719 Ma. The grains’ 176Lu/177Hf ratios are <0.002, present-day 176Hf/177Hf ratios are between 0.281809 and 0.282199, and εHf(t) values are between −14 and +8. About 28% of the analyzed zircons have positive εHf(t) values.

There are no zircons with truly juvenile Hf isotope values corresponding to the depleted-mantle composition at the time t (Fig. 12). Five zircons with ages ranging from 1338 to 1487 Ma have the most positive εHf(t) values between +6.6 and +7.8. They form the major part of the group of oldest zircons in the analyzed samples. In view of the isotopic evolution having started at the time of separation of the parent melt from the protolith, these εHf(t) values may still be considered juvenile (Bahlburg et al., 2011).


Tectonic Setting of the Volcanic Rocks

The geochemical data on the mafic lavas and tuffs permit the inference of a broad geological context of the Chilla beds on the Altiplano west of La Paz, Bolivia, even though they form an isolated unit, in time and space, of the central Andes. The volcanic rocks have a uniform geochemical composition classifying them as tholeiitic basalts and basaltic andesites formed in a continental extensional setting (Fig. 8). In view of the absence of any coeval, true MORB-type rocks or ophiolites of similar age in the central Andes, we interpret the Chilla basalts to have formed as continental within-plate tholeiites indicative of an extensional setting (Fig. 8D). The only known oceanic basalts present in the central Andes are of Ordovician age and occur in northern Peru (Willner et al., 2014). Moreover, the lavas’ association with siliciclastic sedimentary rocks (Fig. 4) of upper crustal geochemical composition, bearing zircons indicating a derivation from continental crust, also attest against a mid-ocean-ridge environment. This agrees with the common presence of <35% of zircons of Neoproterozoic age in Paleozoic sedimentary units of the central Andes giving evidence of the wide distribution of upper crustal sources (Chew et al., 2007b; Bahlburg et al., 2009, 2011; Reimann et al., 2010; Adams et al., 2011; Augustsson et al., 2015) or even a Neoproterozoic magmatic arc (Chew et al., 2008).

There is a geochemically similar but stratigraphically and genetically completely unrelated unit cropping out in E-MORB–like mafic associations along the western margin of the Cuyania terrane in the western Argentinian Precordillera. The volcanic rocks of the Chilla beds are geochemically similar to these much younger, Ordovician rocks (Kay et al., 1984; Boedo et al., 2013). The origin of this Ordovician belt in Argentina is linked to early stages of oceanic rifting along the Gondwana (Precordillera) continental margin. The E-MORB–like character there reflects mixing of depleted and enriched mantle and continental lithospheric sources (Boedo et al., 2013) also recognized in the much older mafic lavas of the Chilla beds (Figs. 8E, 8F; Pearce, 2008).

Implications of the Sedimentological Features

The volcanosedimentary succession of the Chilla beds accumulated in a basin below storm wave base. The latter is indicated by the presence of turbidites and the absence of any wave-related sedimentary structures. The petrography and geochemistry of the sandstones characterize the rocks as relatively immature, little weathered, and composed of potentially first-cycle sediments of upper crustal composition, i.e., devoid of significant input from the associated lavas and tuffs. This gives us confidence that the petrographical composition of the clastic rocks reflects the largely unmodified provenance signal. The sedimentary rocks were likely derived from relatively nearby nonvolcanic margins of the depositional basin (Fig. 13).

Stratigraphic Age and Detrital Zircon Provenance

Maximum Depositional Age and Stratigraphic Age

Many of the dated zircons are idiomorphic and have very likely been transported only over short distances (Markwitz and Kirkland, 2018). They still could be polycyclic if they have been delivered to the Chilla site as parts of rock fragments derived from nearby sources supplying coarser, reworked material. Alternatively, one could speculate that source rocks of the indicated ages were exposed in Neoproterozoic time in the relative vicinity but are no longer exposed. This would coincide with the compositional immaturity and framework composition of the sedimentary rocks, the rather moderate alteration, and the geochemical composition close to that of typical UCC.

The age distribution of the Chilla beds stands out when compared to the large number of detrital zircon age distributions from late Neoproterozoic and Paleozoic sedimentary rocks in the central Andes (Fig. 14A), because ages >1.74 Ga are absent and it lacks late Neoproterozoic, including Brasiliano, ages in general and Pampean ages (650–530 Ma) in particular (Rapela et al., 1998; Escayola et al., 2007, 2011). It also lacks Famatinian ages ranging between 520 and 450–420 Ma (e.g., Bahlburg et al., 2009; Einhorn et al., 2015). Both Famatinian age groups are typical of Gondwana and very common in a large variety of late Neoproterozoic and Paleozoic sedimentary formations in the central Andes (e.g., Chew et al., 2007b, 2008; Bahlburg et al., 2009; Cardona et al., 2009; Augustsson et al., 2015; Einhorn et al., 2015; Naidoo et al., 2016; Pepper et al., 2016), but are notably absent from the Chilla beds (Fig. 14A).

In contrast, on the Arequipa massif north of Arequipa at Río Majes, there is a biostratigraphically dated shallow marine siliciclastic unit of Lower Devonian age belonging to the lower Yamayo Group (Torán-Aplao site; Boucot et al., 1980; Boekhout et al., 2013). The detrital age patterns obtained from two samples of this unit by Reimann et al. (2010) range from 1900 Ma to 800 Ma, with one sample including a single detrital zircon age of ca. 500 Ma (Fig. 14A). The Pampean and Famatinian ages usually abundant in coeval strata in the region (Fig. 14A; Chew et al., 2007b, 2008; Bahlburg et al., 2009; Augustsson et al., 2015) are missing. The age range at the Torán-Aplao site is somewhat similar to that of the Chilla beds; the age distribution, however, is highly dissimilar (Fig. 15), with a value of the statistical cross-correlation coefficient of 0.04 (Sundell and Saylor, 2017). The Devonian sediments of the lower Yamayo Group at the Torán-Aplao site have a unique detrital zircon age distribution most likely reflecting a provenance restricted to the Arequipa massif (Reimann et al., 2010). They appear to be unrelated to the Chilla beds based on the available age spectra.

In view of this Yamayo Group exception, we note the absence from the Chilla beds of Pampean and Famatinian ages so typical of central Andean sedimentary units (Fig. 14A). We are aware that, depending on the tectonic setting and characteristics of the drainage system, there may be a large time lag between the maximum depositional age and the stratigraphic age of deposition (Cawood et al., 2012; Krippner and Bahlburg, 2013). This is rather typical of passive margins where large hinterlands deliver a wide spectrum of source ages and an extended, rather episodic age distribution.

The new detrital zircon U-Pb age data presented here suggest a maximum depositional age between a concordia age of 925 ± 12 Ma (n = 3) and a slightly discordant single zircon age of 796.6 ± 9 Ma (Fig. 11) for the Chilla beds. Individual zircon ages such as our 796 Ma age cannot be relied on as a robust maximum depositional age indicator (Dickinson and Gehrels, 2009).

The age of 925 ± 12 Ma is, in turn, interpreted to represent the Tonian (1000–720 Ma, Cohen et al., 2013, updated) stratigraphic age of the succession. We consider this to be the most realistic interpretation of the stratigraphic age of the unit, which is in accordance with the field relations and the characteristics of the detrital zircon age data in terms of presence and absence of maxima in age abundances. However, further geochronological work is required to date this unit unequivocally.

Correspondingly, we consider the age of 925 ± 12 Ma robust and conclude that it defines the maximum depositional age in the Tonian system of the Neoproterozoic.

Derivation of the Zircon Grains

At first glance, it appears that one does not have to look far to locate potential source areas to which the obtained detrital zircon age record can be linked. The study site on the Bolivian Altiplano is located at the southwestern margin of the Amazonia craton (Fig. 14B). Considering all orogenies at the southwestern Amazonia margin since 2300 Ma (Rizzotto et al., 2014), the detrital zircon age spectrum of the Chilla beds sandstones ranging between 1736 and 800 Ma encompasses only the accretionary orogenic cycles from the Rio Negro–Juruena (1800–1600 Ma) and Rondônia–San Ignacio (1550–1300 Ma) to the Sunsás orogeny (1200–1000 Ma; Figs. 14A, 15, 16). The oldest cycle, the Rio Negro–Juruena, represents 5% of ages and includes the oldest peak at 1730 Ma, which may be related to the felsic Roosevelt volcanic rocks of 1740 Ma of the Juruena province (Santos et al., 2008). Most abundant are zircons (55% of ages) representing the Sunsás orogenic cycle. Here, the peak at 1180 Ma stands out, reflecting the Nova Brasilândia orogeny (Santos et al., 2008). Next in abundance are those ages originally derived from the Rondônia–San Ignacio belt (33% of ages).

The Sunsás orogeny has been connected to the collision of eastern Laurentia with southwestern Amazonia (e.g., Tohver et al., 2005; Cawood et al., 2016; Fig. 1). Consequently, the eastern Laurentian crust could be expected to have formed part of the source areas of the Chilla detritus. Eastern Laurentia at the time of the formation of the Chilla beds is characterized by the continental crust of the Grenville province (McLelland et al., 2010). Basement lithologies and detrital units derived from them have rather limited age distributions, mainly between 1.5 and 0.9 Ga with a pronounced maximum at 1.2 Ga and a minor one at 1.1 Ga (Figs. 16I, 16J; Quinn, 2012; Spencer et al., 2014, 2015; Howard et al., 2015; Petersson et al., 2015; Moecher et al., 2018). When inspected visually, the age patterns appear rather similar to the age distribution shown by the Chilla sandstones (Fig. 16).

The visual inspection of detrital zircon age spectra is a commonly applied method to compare different and sometimes large and variable data sets, even though this method is prone to introducing a subjective bias. Applying visual inspection as an exercise, we can say that the Chilla beds detrital zircon age spectrum is rather similar to the 2000–750 Ma segments of respective distributions of the central Andean units (Figs. 2, 16A–16G), notably (1) the data of Chew et al. (2008) which indicate the presence of a hitherto unknown late Neoproterozoic arc in the region of the eastern central Andes, (2) data from the glacial Chiquerío Formation of Cryogenian age in southern Peru (Chew et al., 2007a) which has a detrital zircon age spectrum and record somewhat similar to those of the Chilla beds, and (3) the detrital zircon age record of modern eastern Andean rivers (Pepper et al., 2016). Age data of the Sunsás orogen (Bettencourt et al., 1999; Boger et al., 2005; Santos et al., 2008) and its immediate precursors are shown in Figure 16H for comparison.

The age distributions of the Neoproterozoic metasedimentary rocks of the Sierra de Moreno in northern Chile (Pankhurst et al., 2016), the Ordovician formations of northern Bolivia (Reimann et al., 2010), and the eastern schists of the Ordovician Marañon complex of central Peru (Cardona et al., 2009), in turn, appear to be more dissimilar to the Chilla beds than the previously mentioned units.

Visual inspection of the age spectra like those shown in Figures 14A and 16 is subjective, offering a high potential for deriving misleading interpretations. This deficiency can be overcome by a statistical analysis of the similarities and differences of the considered age spectra and their age components (Gehrels, 2011; Saylor et al., 2013; Bahlburg and Berndt, 2016; Vermeesch et al., 2016; Sundell and Saylor, 2017; Huber et al., 2018).

The comparative statistical evaluation and analysis of dissimilarity measures of different zircon age distributions by way of multidimensional scaling (MDS) analysis aids in determining relative “sameness” between samples (Vermeesh et al., 2016; Sundell and Saylor, 2017; Wissink et al., 2018). There is disagreement, however, about which of the measures return satisfactory results (Wissink et al., 2018). Likeness and similarity, the Kuiper V test and Kolmogorov-Smirnov statistic, the cross-correlation coefficient, and the Sircombe-Hazelton dissimilarity measure are favored by different authors to return meaningful results on the degree of dissimilarity between samples (Vermeesch et al., 2016; Sundell and Saylor, 2017; Saylor et al., 2018; Wissink et al., 2018). Sundell and Saylor (2017) demonstrated the cross-correlation coefficient to be a reliable tool, whereas Vermeesch et al. (2016) concluded that the Sircombe-Hazelton dissimilarity measure and Kolmogorov-Smirnov distance return functional results. In another evaluation of the significance of the different criteria when applied to natural detrital zircon age distributions, Wissink et al. (2018) concluded that likeness and similarity return the most meaningful results and that the cross-correlation coefficient is not of similar value. A means to gauge the quality with which a MDS configuration preserves the inter-sample dissimilarities is the calculation of the stress (Kruskal, 1964) displayed commonly in Shepard plots including the stress value (e.g., Saylor et al., 2018; Wissink et al., 2018).

We apply MDS analysis using the DZmds software (Saylor et al., 2018) focusing on the criteria favored by different groups of authors, namely the cross-correlation coefficient, the likeness value, and the Kolmogorov-Smirnov distance. The cross-correlation coefficient and the likeness value analysis return “good” goodness-of-fit values (stress = 0.058 and 0.067, respectively; Fig. 17), whereas the Kolmogorov-Smirnov distance resulted in “poor” fits (stress = 4.61).

We apply this statistical analysis of the dissimilarity of zircon age spectra pursuing two goals. The first is to statistically compare the considered detrital zircon age spectra from 2000 to 750 Ma of a variety of Neoproterozoic and Paleozoic formations from the Andean region with the Chilla beds data (Fig. 16). The second is to discriminate between the two possible source regions of the Chilla beds detritus, namely the Grenville province of Laurentia and the Sunsás and Rondônia orogens of southwestern Amazonia (Figs. 16, 17).

The MDS plots in Figure 17 shows a clustering of formations in three groups. In both the cross-correlation and the likeness plots (Figs. 17A, 17B), the Grenville basement and Grenville detrital zircon data plot on the left side of the MDS maps. The center is taken up by a relatively closely spaced cluster comprising the Sunsás orogen, the inferred Neoproterozoic arc, the modern eastern central Andean rivers, and the Neoproterozoic Chiquerío Formation of coastal southern Peru. This cluster of samples is linked by nearest-neighbor relationships (Fig. 17A). In the center of the plots, the Chilla beds demonstrate dissimilarities to most other groups, with first-neighbor but still relatively distant relationships to the right side of the map with the Ordovician clastics of northern Bolivia and the eastern Marañon schists. The metasedimentary rocks of the Sierra de Moreno in northern Chile have a very limited age distribution and plot with a large dissimilarity from all the other groups.

The Sunsás data have an unequivocal Amazonian source from the Sunsás orogen. (Bettencourt et al., 1999; Boger et al., 2005; Santos et al., 2008). The low degrees of dissimilarity and first-neighbor relationships to the Chiquerío Formation, the inferred Neoproterozoic arc, and the modern eastern Andean rivers identify this group as clearly having an autochthonous South American provenance. This group of formations with a clear Amazonia provenance is significantly dissimilar to the Grenville basement. This group is also more dissimilar to the group forming the loose cluster on the right hand side of the diagrams, namely the Chilla beds, the eastern Marañon schists, and the Ordovician sedimentary rocks in northern Bolivia. However, in terms of first- or second-neighbor relationships and statistical dissimilarity, the Chilla beds appear to be very dissimilar to the Grenville sources and very similar to the autochthonous Amazonian sources (Fig. 17).

However, the detrital populations of siliciclastic sedimentary rocks are only rarely derived from a single source only. Mixing of sources is common, and unmixing of the source signal in detrital distributional data is a great challenge (Weltje, 2012). Quantifying source contributions to an observed composition is particularly important when sedimentary recycling is an important factor (Campbell et al., 2005; Satkoski et al., 2013; Sundell and Saylor, 2017). On the other hand, this also means that successful unmixing of such data strongly amplifies the interpretative power of a provenance analysis.

We tested the Chilla detrital zircon provenance results using the DZmix software of Sundell and Saylor (2017). This analysis tool combines an inverse Monte Carlo model with an optimized forward model. It returns relative source distributions based on the cross-correlation coefficient, the Kuiper test V value, and the Kolmogorov-Smirnov test D value (Sundell and Saylor, 2017). We tested the Chilla beds age spectrum by using the equivalent age ranges from the Sunsás data (Bettencourt et al., 1999; Boger et al., 2005; Santos et al., 2008), the two data sets from the Grenville basement (Quinn, 2012; Howard et al., 2015; Petersson et al., 2015; Moecher et al., 2018), and the detritals derived from the Grenville province (Fig. 18; Howard et al., 2015; Spencer et al., 2015; Moecher et al., 2018) as potential sources. The Monte Carlo model returned the following values of the relative contribution of the sources based on the cross-correlation coefficient: Sunsás, 0.98; Grenville basement, 0.01; and Grenville detritals, 0.01. Similarly, the Kuiper V values were 0.97, 0.02, and 0.01, respectively. If the Grenville basement and detrital data are combined into a single source, the respective values for Sunsás and Grenville are 0.99 and 0.01 respectively using the cross-correlation coefficient, and 0.99 and 0.01 respectively using the Kuiper V values. This, together with the previously discussed statistical assessments, makes it highly likely that almost all of the detrital zircons contained in the Chilla sandstones have been derived from primary or recycled Sunsás sources belonging to the Amazonian realm (Figs. 17, 18).

Implications of the Hf Isotope Systematics

In age versus εHf(t) space (Fig. 12), the εHf(t) data form a single array from relatively juvenile compositions between 1.6 and 1.4 Ga which evolved through time according to the crustal evolution trend (Vervoort and Patchett, 1996) to progressively more negative values, with the youngest zircons having values of −14 and −12 at ca. 0.9 Ga. This arrangement can be read in two ways. It may indicate that the zircons originally formed from a juvenile protolith between 1.6 and 1.4 Ga, with the oldest zircon projecting back to 2.0 Ga, and followed an uninterrupted crustal evolution path to ca. 0.9 Ga. On the other hand, between 1.4 and 0.95 Ga, the data form vertical arrays which are strongly indicative of crustal contamination of juvenile protoliths having formed at about the respective times. When viewed in combination with εHf(t) data from detrital zircons of late Proterozoic and early Paleozoic sedimentary rocks of the central Andean region, the arrangement in vertical arrays displays itself as a distinctive feature (Fig. 12). This arrangement reflects widespread crustal recycling in the different orogenic belts of southwestern Amazonia and the proto-Andean region (Bock et al., 2000; Zimmermann and Bahlburg, 2003; Willner et al., 2008; Reimann et al., 2010; Bahlburg et al., 2009, 2011, 2016; Boekhout et al., 2015), an interpretation also favored for the data obtained from the Chilla sandstones.

The analyzed zircons have HfTDM uniformly projecting back to crustal residence age estimates (model ages, TDM) between 1500 and 2300 Ma. In South American geology, this represents the times of the Marconi-Icantiunas (2.2–2.0 Ga), Ventuari-Tapajós (Transamazonian; 2.0–1.8 Ga), Rio Negro–Juruena (1.8–1.55 Ga), and Rondônia–San Ignacio (1.55–1.2 Ga) orogenies on the southwestern Amazonia craton (Fig. 14B; Cordani and Teixeira, 2007; Bettencourt et al., 2010). This also corresponds with Nd and Pb whole-rock isotope data (Tosdal, 1996; Pankhurst et al., 1998; Rapela et al., 1998; Bock et al., 2000; Cordani et al., 2000; Zimmermann and Bahlburg, 2003). However, relating model ages to actual orogenic events may be misleading because almost all εHf(t) values reflect some degree of crustal contamination, as is indicated by their offset from the depleted mantle curve and the arrangement in the mentioned vertical arrays. The data rather reflect juvenile magmatic activity ongoing from ca. 1.6 to 0.9 Ga, the melts of which experienced significant crustal contamination throughout.

The use of zircon Hf isotopic signatures in the analysis of provenance of central Andean lithologic units is strongly hampered by the limited availability of Hf isotope data. To the best of our knowledge, there are no Hf isotope data on any rocks from the Sunsás orogen. Referring to our statistical analysis of the U-Pb age zircon provenance discussed above, we decided to use the 1750–800 Ma modern river U-Pb and Hf isotope data of Pepper et al. (2016) as a proxy for the Sunsás provenance. We combined these data with the similar-age-range spectra from Ordovician sedimentary units from northern Bolivia and southern Peru (Reimann et al., 2010) and the Puncoviscana Formation (upper Neoproterozoic–lower Cambrian) in northwestern Argentina (Adams et al., 2011). These we compared to equivalent data from the eastern Laurentia Grenville basement (n = 261; Howard et al., 2015; Petersson et al., 2015; Hantsche, 2015) and detrital units derived from it (n = 605; Howard et al., 2015; Spencer et al., 2015).

In the analysis of the Hf isotopic signatures of the Chilla beds and their potential sources, we also applied the classical Kolmogorov-Smirnov MDS statistical analysis. MDS operates on single dissimilarity measures (Vermeesch et al., 2016). Because the Hf isotopic signature of a zircon is age dependent, we first combined the two by calculating the ratio between the U-Pb age and the initial Hf isotope value at the time t, age/εHf(t), to obtain the required single dissimilarity measure (Fig. 19A). We checked this approach by calculating the ratio of the U-Pb age to the initial 176Hf/177Hf value in order to minimize the influence of the U-Pb age on the Hf factor (Fig. 19B).

The result shows that the statistical dissimilarity in both models is largest between the Chilla beds and potential Grenville basement sources (Fig. 19). The dissimilarity is smaller to the detrital units derived from Grenville age sources, which may include input from other, as yet unspecified sources of Grenville age from outside the actual Grenville province. Sunsás sources are not likely candidates for these unspecified sources because the dissimilarity between the Sunsás and the Grenville detrital data is relatively large in both measures (Fig. 19). There is a relatively low statistical dissimilarity to the South American sources including the Sunsás domain, and the respective ages contained in the Puncoviscana Formation and the Ordovician sedimentary rocks of northern Bolivia and southern Peru. This relatively close connection is in keeping with the available geological evidence (Reimann et al., 2010; Reimann Zumsprekel et al., 2015; Augustsson et al., 2016; Pepper et al., 2016). The closest statistical neighbor to the Chilla beds is the Puncoviscana Formation, followed by the Ordovician sedimentary rocks of northern Bolivia and southern Peru (Fig. 19), all of which are significantly younger than the Chilla beds. This MDS analysis of the Hf isotope data, together with the U-Pb data, allows differentiation of an Amazonia provenance realm from the Grenville realm, with the Chilla beds data plotting with the Amazonia realm (Fig. 19). We conclude that the Hf isotope evidence also is in keeping with a derivation of the Chilla beds from Amazonia sources, including the Sunsás orogen and intermittent sedimentary storage units.

Regional Considerations

The basement to the Chilla beds is unknown. In fact, in the central Andes there are only very few, distant and isolated outcrops of the pre–latest Neoproterozoic basement (Ramos, 2008, 2010). However, this mostly unexposed basement has been probed in particular by younger, early Paleozoic igneous rocks that assimilated the older crust or incorporated it as xenoliths. Very diverse U-Pb zircon age spectra in Ordovician granitoids in northwestern Argentina, for example, give evidence of crustal inheritance and recycling ranging in age from 2680 to 450 Ma (Bahlburg and Berndt, 2016; Bahlburg et al., 2016). Also, the detrital age spectra of Ordovician sedimentary and metasedimentary rocks exposed in the relative vicinity in northern Bolivia and southern Peru all contain significant numbers of ages between 2000 and 900 Ma (Chew et al., 2008; Reimann et al., 2010; Bahlburg et al., 2011), indicating the former prominence of this older basement. We speculate therefore that the detrital age spectrum of the Chilla beds in fact indicates the presence and sampling of an underlying Paleoproterozoic to mainly middle to late Mesoproterozoic basement during Tonian time.

This Mesoproterozoic crust records the collision of eastern Laurentia and southwestern Amazonia during the amalgamation of Rodinia (Fig. 1; Loewy et al., 2004; Tohver et al., 2005; Casquet et al., 2010; Cawood et al., 2016). Breakup of Rodinia commenced in the Tonian (Cawood et al., 2016) and led to continental separation of Laurentia from southwestern Amazonia by ca. 760 Ma (Cawood et al., 2016; Fig. 1). We interpret the rift-related tholeiitic basalts of the Chilla beds as having erupted during this period of progressive continental breakup in an extensional continental basin. In this basin, immature quartz-poor to quartz-intermediate sandstones of regional, southwestern Amazonian provenance were deposited as turbidites. We infer from our analysis that this detritus has been derived almost exclusively from Amazonian sources. We consider the Chilla beds therefore to presently represent the only direct evidence of Rodinia breakup at the evolving southwestern Amazonian passive margin. This passive margin was transformed into an active one after 650 Ma when the Puncoviscana marginal basin and associated Pampean magmatic arc formed in what is now northwestern Argentina (Rapela et al., 1998, 2016; Keppie and Bahlburg, 1999; Zimmermann, 2005; Escayola et al., 2011; Fig. 1).


The collision of eastern Laurentia and southwestern Amazonia in the course of the formation of the supercontinent Rodinia transformed the Grenville-age Sunsás orogen at the southwestern Amazonia margin into an interior orogen. Breakup of Rodinia started in this region after 0.9 Ga. West of the Sunsás orogen, in what is now the Altiplano west of La Paz, Bolivia, rifting led to the deposition of the volcanosedimentary Chilla beds in a continental extensional basin on Mesoproterozoic crust. Interpretation of the detrital zircon U-Pb age data indicate that deposition of the Chilla beds occurred in the Tonian period.

By 760 Ma, toward the end of the Tonian, the passive margin of Gondwana in the region of southwestern Amazonia had evolved very likely outboard of the present position of the Chilla beds (Fig. 1; Cawood et al., 2016) because no true passive margin sequence or oceanic crust of Neoproterozoic age is known from the region. The detrital zircon age distribution and statistical characteristics of the Chilla beds sandstones point to predominant Amazonian sources linked directly or indirectly to the Sunsás and precursor orogens.


We thank Peter Cawood for supplying us with a digital version of figure 2 of Cawood et al. (2016), which we adapted for our Figure 1. We thank P. Castillo (Westfälische Wilhelms-Universität Münster, Germany), for her constructive comments on an earlier version of this manuscript. We acknowledge support from the Open Access Publication Fund of the University of Münster. We thank Z. Sickmann and an anonymous reviewer for their insightful comments, and associate editor T. LaMaskin for editorial guidance. This study has been funded partly by German Research Foundation grant BA 1011/43-1 to HB. This is a contribution to IGCP (International Geoscience Program) Project 648, “Supercontinent Cycles and Global Geodynamics”.

1Supplemental Tables. Table S1: Chilla beds, geochronological data. Table S2: Chilla beds, Hf isotope data. Table S3: Chilla beds, whole-rock geochemical data. Please visit https://doi.org/10.1130/GES02151.S1 or access the full-text article on www.gsapubs.org to view the Supplemental Tables.
Science Editor: Shanaka de Silva
Associate Editor: Todd LaMaskin
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.