Stratigraphy of the Ghaggar alluvium

We established the subsurface stratigraphy of the Ghaggar alluvium along an ~300 km stretch of the basin based on field and age data from this work and earlier studies (Fig. 1B; Supplementary Fig. S1). Sedimentological analyses along the river valley helped us ascertain the river’s palaeo morpho-dynamics. The alluvium consists of three distinct sedimentary facies; grey micaceous coarse/medium grained sand (GS; grain size >250 μm), yellowish-brown fine-grained sand (YBS; grain size <125 μm) and brown colored mud (BM), deposited from bottom to top. The GS facies is composed of quartz, feldspar, heavy minerals (zircon, amphibole, kyanite, sillimanite, garnet, and pyroxene)21 and mica, which gives it a salt and pepper appearance. It contains poorly sorted coarse to medium, angular sand grains and no clay, which suggest their deposition as fluvial channel-fill/bar deposits under strong hydrodynamic conditions. The textural pattern of the GS facies resembles that of typical channel fills in neighboring Himalayan rivers16,22. The fining upward YBS facies contains clay, apart from fine sand and is overlain by the clay dominated BM facies.

The sharp changes in grain size, sorting and clay content from the GS facies to the overlying YBS and BM facies are indicative of the transformation of the fluvial system from a high-energy domain to a low-energy one. Similar observations have been made by many earlier studies in the Ghaggar alluvium16,21,23. However, what these studies lacked is the chronology of depositions of these three facies, which is presented in this work. We determined the timings of depositions of various facies by AMS C-14 dating of mollusk shells and optically-stimulated luminescence dating of sands (Supplementary Tables S1 and S2). Our results in conjunction with the existing information suggest that the suspended load, the BM facies, dominated the Ghaggar system since ~4.5ka, whereas sediments of the YBS facies got transported by the river during 20 and 3 ka (Fig. 1B)16,24,25. The deposition of the GS facies, on the other hand, can be traced back to the Pleistocene (Fig. 1B). Whereas most studies in the Ghaggar alluvium report cessation of deposition of the GS facies prior to the Holocene, a few report much younger occurrence of the GS facies, in the upper Haryana plains21,24. The lack of report of the younger phase of the GS facies in most of the earlier investigations could possibly have been due to sampling bias in a highly meandering fluvial system containing multiple buried channels. With these new data it can be clearly established that the GS facies occur in two time domains, during >80-20 ka and 9.0-4.5 ka with a break of about eleven millennia, along a long stretch of the Ghaggar palaeo-channel (Fig. 1B). The pre-Holocene deposits of the GS facies is generally accepted to be the relicts of a mega Himalayan paleo-river system15,16, however, the nature of the younger phase of the GS facies has not yet been explored. In this work we make an attempt to establish the origin of these Holocene sediments of the Ghaggar.

Provenance of the Ghaggar sediments

Since the coarse sands of the GS facies were deposited by a powerful channel of the river in the entire stretch of the alluvium (~300 km), it is reasonable to expect that the sediments were transported from a long distance (>600 km), most likely from the Himalayas. To determine the location of the sources, we utilized 40Ar/39Ar ages of coarse-grained muscovites present in these sands as provenance indicator. We dated three multigrain separates from the two depositional time intervals to capture the dominant age groups and their temporal distribution. Indistinguishable plateau and isochron ages, and intercepts showing atmospheric 36Ar/40Ar suggest that these muscovites were derived from igneous or metamorphic rocks which had reached argon closure temperature during 20.1-18.6 Ma (Supplementary Fig. S2). Presence of identical age mica clusters in both the sand horizons (Fig. 2A) suggests their derivation from similar sources. This age has a significant overlap with the mode of the 40Ar/39Ar age distribution for the muscovites in the Higher Himalayan Crystalline rocks (Fig. 2A), which establishes a most probable genetic link. The dominant source lithologies that could have provided the micas are the 20 ± 2 Ma leucogranites and associated migmatites located above the Main Central Thrust, in the glaciated regions of the Higher Himalaya26. The distribution of muscovite 40Ar/39Ar ages of the Lesser Himalayan rocks also shows a mode, albeit minor, at ~20 Ma (Fig. 2A), therefore, we cannot completely rule out the contribution of sediments from these rocks to the GS facies in the Ghaggar palaeo channels. The rocks of the Siwaliks are unlikely sources of the GS facies sediments, because the detrital muscovites in these rocks are much older than that found in the latter (Fig. 2A). In addition, removal of these mineral grains from the Siwalik rocks could only have contributed weathered mica grains or clays to the latter, which does not appear to be the case. An earlier study of 40Ar/39Ar ages of detrital mica from the Ghaggar alluvium16 had yielded similar results as ours.

Figure 2 Source fingerprinting with 40Ar/39Ar muscovite ages and Sr-Nd isotopic ratios of sediments. (A) Probability density plot of 40Ar/39Ar ages of muscovites from different Himalayan litho-tectonic units40 compared with our data from detrital muscovites from the Ghaggar alluvium - the orange band (18.6–20.1 Ma). (B) ε Nd vs. 87Sr/86Sr plot of sediments from the Ghaggar alluvium compared with binary mixing curves for the Higher Himalayan Crystalline (HHC) and the Lesser Himalayan Series (LHS) end-members46. LS: Local Sources in the Thar Desert27. SYG: Sutlej-Yamuna-Ganga. Full size image

To further constrain the provenance of the Ghaggar alluvium, we used Sr and Nd isotopic ratios of sediments for source fingerprinting (Supplementary Table S3). We found that the sediments in the GS facies horizons, irrespective of their depositional ages, have high 87Sr/86Sr (>0.75) and low ε Nd (<−17). These compositions overlap with that of sediments in the glacier-fed rivers like the Sutlej, Yamuna, and Ganga, which hint at their derivation from similar sources (Fig. 2B). A binary mixing model using possible Himalayan sources (Fig. 1A) suggests that the GS facies sediments had their sources in the Higher and Lesser-Himalayas (Fig. 2B). In contrast, the sediments of the BM and the YBS facies, having much different isotopic compositions (Fig. 2B), were likely derived from sources made up of already recycled sediments such as the Siwaliks, older alluvium and dunes of the Thar Desert (Fig. 1A)15,27. The dominance of smectite over illite in the clay mineral fraction of the YBS facies also supports a local origin for these sediments27, since smectite, a secondary clay, was believed to have been produced within the Ghaggar alluvium as a result of the second stage of chemical weathering process27. Isotopic data considered in conjunction with the 40Ar/39Ar ages of the detrital muscovite indicate that the GS facies sediments were directly deposited by a river which had its catchment in the Higher and Lesser Himalaya. An earlier provenance study in the Ghaggar alluvium based on detrital zircon U-Pb geochronology had also traced the sources of the GS facies to the Higher Himalayas16. However, the youngest GS facies horizon reported by these authors was older than 12.3 ka16, whereas we have encountered and studied much younger horizons of the facies (Fig. 1B; Fig. 3A; Supplementary Fig. S1).

Figure 3 Evolution of the Ghaggar from changes in sediment provenance and the Harappan settlement dynamics. (A) Stratigraphic changes in sediment Sr-Nd isotopic compositions in the Ghaggar alluvium during last 20 ka. Symbols and abbreviations are as in Fig. 1. (B) Evolution of the Harappan civilization in north-western India and eastern Pakistan as inferred from the settlement dynamics through ages (9.0-3.5 ka)7,10,47,48. Modern and inferred former courses (dashed lines) of the major Himalayan rivers are also shown. I: Indus; J: Jhelum; C: Chenab; B: Beas; S: Sutlej; SS: Saraswati (Ghaggar); Y: Yamuna; G: Ganga. Full size image

The evolving river system

A comprehensive reconstruction of the fluvial history of the Ghaggar, based on existing knowledge and results of this study (Fig. 1B; Fig. 3A), reveals that the river was perennial during the Marine Isotope Stages (MIS) 3 and 4 (80-20 ka), when it transported coarse-grained micaceous sand (GS facies) from the glaciated regions of the Higher-Himalaya onto the plains of north-western India. This fluvial phase ended during the peak aridity of the Last Glacial Maximum (LGM)25, which also had adversely affected the discharge in other major western-Himalayan rivers8. Subsequently, during 20-9 ka, the sediment load in the river got overwhelmed with material derived from the Siwaliks/local sources (YBS facies), which is clearly reflected in the shifts in 87Sr/86Sr and ε Nd of the Ghaggar alluvium (Fig. 3A). With the intensification of the ISM10 and melting of the Himalayan glaciers during the MIS-128 the GS facies made a reappearance at ~9 ka and continued its dominance until ~4.5 ka (Fig. 3A). Although, the Ghaggar received sediments from the Higher Himalaya during the river’s rejuvenated phase, there exists no evidence to suggest that any of its modern tributaries originated from glaciers. In this scenario, the possible pathways for delivery of such sediments into the Ghaggar could have been through the neighboring Sutlej and Yamuna, which currently flow through the Higher and Lesser Himalayan rocks. However, since it is believed that Yamuna had abandoned the Ghaggar channel during 49-10 ka15, the only remaining pathway for the Higher Himalayan sediments into the Ghaggar during the Holocene could have been the Sutlej. As a matter of fact, muscovites in the present-day Sutlej sand, which is derived from sources in both the Higher and Lesser Himalayas, have a major 40Ar/39Ar age mode at 16 ± 3 Ma16. Supporting evidence for such a hypothesis comes from the similarity of 87Sr/86Sr of in-situ mollusk shells from the grey-sand bodies (0.7184–0.7190; this work) with that of the water of the Sutlej (0.7166–0.7218)29,30. The U-Pb ages of detrital zircons from the Ghaggar alluvium also suggest that the majority of zircon populations in the younger deposits were predominantly originated from the Sutlej15,16. The Sutlej is known to have had substantial increase in its water and sediment flux in upper reaches during the Holocene31, which could have been delivered to the Ghaggar making the latter a perennial stream for the next few thousand years. At ~4.5 ka, 87Sr/86Sr and ε Nd of the alluvium revert to the pre-perennial values suggesting the end of supply of the Higher Himalayan sediments due to a break in the Sutlej-Ghaggar connection, which turned the Ghaggar into an ephemeral system. This final disruption of the river flow is also reflected in the sediment composition of its purported delta towards further south in the Great Rann of Kachchh32.

Dynamics of the Harappan settlements

A thorough scrutiny of the settlement dynamics of the Harappan Civilization reveals that the timing of the rejuvenated perennial phase of the Ghaggar (9-4.5 ka) coincides with that of the flourishing of the Pre-Harappan and Early Harappan cultures along its banks (Fig. 3B). Towards the end of the Mature Harappan phase (4.6-3.9 ka), there is a clear evidence of human migrations to the lower and upper reaches of the river, leaving the middle part sparsely populated (Fig. 3B), which could be attributed to the disorganization of the river as established in this work. The lower reaches of the river, in the Hakra sector, had possibly remained perennial, through a connection from the Sutlej, supporting mature and post-urban Harappan settlements (Fig. 3B). Our study brings to light the fact that the Harappans built their early settlements along a stronger phase of the river Ghaggar, during ~9 to 4.5 ka, which would later be known as the Saraswati. However, by the time the civilization matured, the river had already lost its glacial connection. These inferences confirm the observation of an earlier study, based exclusively on changes in the settlement patterns, that the Ghaggar first broke up at ~4.6 ka12. Interestingly, the timing of the ultimate disruption to the perennial phase of the Ghaggar roughly coincides with the beginning of the Meghalayan Stage (~4.2 ka)33. The urban Harappans abandoned their settlements in the Ghaggar valley within next few centuries and the civilisation declined by 3.9 ka. Although, the decline of the civilization in the Ghaggar-Saraswati valley postdates the exceptional changes to the flow of the river, a stronger perennial phase appears to have helped the early societies to sow the seeds of the earliest known civilization of the Indian subcontinent.

Methodology

In the arid to semi-arid climate of western India, preservation of organic carbon within the sandy sediments is rare. Therefore, for radiocarbon dating to constrain depositional ages we focussed on collection of mollusc shells buried in situ along with the sediments. Shells were cleaned with H 2 O 2 to remove organic matter. Subsequently, they were dipped in 0.1 N HCl for ~5 seconds and washed thoroughly with deionized water to remove (any) altered outermost layers. Sample powders for dating were micro drilled from the umbo regions, which are generally considered resistant to alteration. AMS C-14 dating was done at the Centro Nacional de Aceleradores, Sevilla, Spain34, and the ages were calibrated using INTCAL 1335. The typical reproducibility (1σ) of our radiocarbon ages is <1% (Supplementary Table S1).

Optically stimulated luminescence (OSL) of quartz grains in a sediment sample is proportional to its burial age. This dating method works with the premise that the quartz grains have negligible inherited OSL signal prior to deposition and that they subsequently acquire it from ambient radiation. The zeroing of OSL signal in these grains happens upon exposure to daylight prior to deposition, while on transport, or during any recycling/re-exposure event. For dating purpose, total OSL of quartz grains, separated from a sample in dark-room conditions, is measured and converted to units of radiation dose or total paleo-dose, which is then divided by the dose rate to derive the age since last exposure. The dose rate is determined from measured concentrations of U, Th and K and estimated contributions from cosmic rays36 in the immediate burial environment of the sample. Samples for of this study, collected from freshly cleaned up stratigraphic sections using specially designed aluminium pipes, were dated using the Single Aliquot Regeneration (SAR) protocol for quartz grains37. Concentrations of U, Th and K were measured on a High purity Germanium detector (HPGe). The OSL dating was done at the Physical Research Laboratory, India. The typical reproducibility (1σ) of our OSL ages is <10% (Supplementary Table S2).

The closure temperature of the 40Ar/39Ar isotope system for muscovites is ~350 °C38, which enables it to capture the timings of the latest tectonothermal events experienced by the parent rocks. Since the litho-tectonic units of the Himalaya had exhumed diachronously, they contain various age populations of muscovites representing each of the exhumation event39. Therefore, the 40Ar/39Ar ages of detrital muscovites, derived from these rocks, serve as powerful source indicators40. To determine 40Ar/39Ar ages of detrital muscovites in the Ghaggar alluvium we handpicked ~200 mg of grains, with grain size >500 μm, from grey sand horizons from three distinct stratigraphic horizons with depositional ages of 37, >16 and 6.3 ka. Smaller mica grains were avoided because they are more easily affected by secondary alteration as compared to larger grains41. Moreover, determination of 40Ar/39Ar ages of smaller grains is difficult and can become erroneous because of Argon loss due to recoil during irradiation42. Samples were packed in aluminium capsules and irradiated in a research reactor by fast neutrons. The Minnesota hornblende reference material (MMhb-1), of age 523.1 ± 2.6 Ma43, was used as the flux monitor and high-purity CaF 2 and K 2 SO 4 salts for interference corrections arising from the production of Ar from Ca and K isotopes. Argon was extracted by incremental heating between 450 °C and 1400 °C at steps of 50 °C and isotopic ratios were measured in a Thermo Fisher ARGUS-VI multi-collector mass spectrometer at the Indian Institute of Technology Bombay. Plateau and isochron ages were calculated and plotted using the software ISOPLOT 2.4944.