Regional Pliocene exhumation of the Lesser Himalaya in the Indus drainage

New bulk sediment Sr and Nd isotope data, coupled with U–Pb dating of detrital zircon grains from sediment cored by the International Ocean Discovery Program in the Arabian Sea, allow the reconstruction of erosion in the Indus catchment since ∼17 Ma. Increasing ε_Nd values from 17 to 9.5 Ma imply relatively more erosion from the Karakoram and Kohistan, likely linked to slip on the Karakoram Fault and compression in the southern and eastern Karakoram. After a period of relative stability from 9.5 to 5.7 Ma, there is a long-term decrease in ε_Nd values that corresponds with increasing relative abundance of >300 Ma zircon grains that are most common in Himalayan bedrocks. The continuous presence of abundant Himalayan zircons precludes large-scale drainage capture as the cause of decreasing ε_Nd values in the submarine fan. Although the initial increase in Lesser Himalaya-derived 1500–2300 Ma zircons after 8.3 Ma is consistent with earlier records from the foreland basin, the much greater rise after 1.9 Ma has not previously been recognized and suggests that widespread unroofing of the Crystalline Lesser Himalaya and to a lesser extent Nanga Parbat did not occur until after 1.9 Ma. Because regional erosion increased in the Pleistocene compared to the Pliocene, the relative increase in erosion from the Lesser Himalaya does not reflect slowing erosion in the Karakoram and Greater Himalaya. No simple links can be made between erosion and the development of the South Asian Monsoon, implying a largely tectonic control on Lesser Himalayan unroofing.

1 Introduction

The Western Himalaya region represents a classic example of an orogen formed by the collision of two continental plates. Sediment eroded from across this area records the development and unroofing of the mountains and has been accumulating in the Arabian Sea since the start of the collision, likely in the Eocene (DeCelles et al., 2014; Najman et al., 2010). Because erosion removes rocks from the mountains, the older history of the ranges is best reconstructed from the sedimentary record. Although some of this record is preserved onshore in the foreland basin (Najman, 2006), these sequences represent a relatively poorly dated and incomplete archive, with the erosion history of any particular section of the foreland basin being dominated by the immediately adjacent ranges rather than providing a complete regional record (Burbank et al., 1996). This makes it difficult to use the foreland sequences alone to address the ongoing debate regarding the competing roles of climatically modulated surface processes compared to solid Earth tectonic forces in controlling the structural evolution of the Himalaya (Beaumont et al., 2001; Robinson et al., 2006; Webb et al., 2011). Without a detailed erosional history, it is impossible to fully reconstruct how the mountain belt has evolved over long periods of geologic time and what role erosion has played in focusing exhumation and controlling the location of major structures (Beaumont et al., 2001; Wobus et al., 2003).

2 Marine erosion records

Scientific drilling conducted in 2015 by International Ocean Discovery Program (IODP) Expedition 355 now provides the opportunity to examine how Himalayan erosion has changed since ∼10.8 Ma (Pandey et al., 2015). Although drilling in the Laxmi Basin offshore western India was able to reach the Cretaceous basement at Site U1457 (Fig. 1) the recovered submarine fan sequences are mostly limited to the past 10.8 Myr as a result of a large mass transport complex (MTC), which removed much of the older record (Calvès et al., 2015). At Site U1456 a single sample of siliciclastic sand dated to ∼15.5 Ma was recovered below the carbonate breccias of the MTC, but otherwise the sediments span the same 10.8 Myr seen at Site U1457. Sediments from within the MTC are not included in this study. Nonetheless, the new core provides a much better-dated and more complete basin-wide history of erosion than has been previously available for the Western Himalaya covering the period since the middle Miocene. Here we present new provenance data from the IODP drill sites, as well as additional data from industrial borehole Indus Marine A1 drilled by Wintershall Holding (Germany) and located near the river mouth (Fig. 1), in order to reconstruct the erosion of the area in detail since ∼17 Ma, the date of the oldest material from Indus Marine A1.

Figure 1 (a) Shaded bathymetric and topographic map of the Arabian Sea area showing the location of the drilling sites within the Laxmi Basin, considered by this study. Map also shows the primary source terranes and the major tributary systems of the Indus River, as well as smaller Indian Peninsula rivers that may have provided material to the drill sites. (b) Inset map shows detail of the Laxmi Basin and location of the drill sites considered in this study. Numbered red circles indicate existing scientific boreholes from the Deep Sea Drilling Project (DSDP) and Ocean Drilling Program (ODP). Magnetic anomalies are from Miles and Roest (1993). White dashed lines show transform faults. KK: Karakoram; NP: Nanga Parbat. (c) Geological map of the Western Himalaya showing the major tectonic units that are eroded by the Indus River and its tributaries. Map is modified after Garzanti et al. (2005). Rivers are shown in thick black lines. ISZ: Indus suture zone; MCT: Main Central Thrust; MBT: Main Boundary Thrust; MFT: Main Frontal Thrust. Thick black line shows the boundary of the Indus drainage, while thinner lines demark the limits of the major Himalayan tributaries.

3 Background and prior studies

Age control at Sites U1456 and U1457 is provided by a combination of biostratigraphy and magnetic stratigraphy (Pandey et al., 2016), while at Indus Marine A1 age control is accomplished through biostratigraphy alone (Shuaib, 1982). We used bulk decarbonated sediment Nd and Sr isotope compositions coupled with single-grain detrital U–Pb zircon ages in order to constrain how the source of sediment changed over long periods of time. Zircon U–Pb dating is only possible when the siliciclastic sediment is sufficiently coarse grained, restricting the application of this method to the last ∼8.3 Myr. Nd isotopes have a history as robust provenance indicators of siliciclastic sediment and are especially effective in the Western Himalaya, where there is a wide range of isotopic values known from the different bedrock source terranes (Clift et al., 2002). This isotope system is particularly effective because it is not believed to be affected by sediment transport or chemical weathering processes (Goldstein et al., 1984), although it is moderately influenced by grain size in the Indus catchment (Jonell et al., 2018). Although it has been noted that HCl leaching during carbonate removal prior to analysis can affect the measured ¹⁴³Nd∕¹⁴⁴Nd values (Hein et al., 2017), in this study we treated all analyses with acetic acid so that no systematic bias was introduced. Sr isotopes are also employed because they provide an additional dimension for resolving sources, although the sediment has to be carefully decarbonated prior to analysis. Sr isotopic values are also influenced by chemical weathering and grain size variations (Derry and France-Lanord, 1996; Jonell et al., 2018). Earlier work from the Arabian Sea, including samples from Indus Marine A1, used only Nd isotopes to demonstrate a substantial change in source since ∼6 Ma, a trend that was attributed to large-scale drainage capture of the Eastern Himalaya draining tributaries of the Indus away from the Ganges Basin towards the east (Clift and Blusztajn, 2005; Zhuang et al., 2015).

More recently this late Miocene to Holocene, including present-day Nd isotopic trend was attributed in part to the onset of uplift of the Nanga Parbat Massif within the Western Syntaxis (Chirouze et al., 2015) since this source yields extremely radiogenic Nd (Whittington et al., 1999). Accelerated erosion elsewhere in the Himalaya was also inferred to have a played a part in causing the evolution in Nd isotopes. Comparison of cosmogenic and high-temperature thermo-chronometers (e.g., Ar∕Ar muscovite) implies that faster erosion has occurred during the last few million years in the Western Himalaya (Vance et al., 2003). We supplement bulk sediment Nd analyses with U–Pb dating of single zircon sand grains that are also recognized as effective provenance proxies within the Western Himalaya (Alizai et al., 2011; Gehrels et al., 2011) but which have never before been used at the basin-wide scale over million-year timescales due to lack of appropriate sections. Zircon grains are physically robust during transport and may be recycled many times from their original sources, but the U–Pb age is not easily reset because of the high temperatures required (∼950 ^∘C; Gehrels, 2014). Sources rarely have single distinctive ages, but rather a characteristic spectrum of zircons ages, and these are transferred to the sediments during erosion. Comparing sediment and bedrock age spectra allows the source of the sediment to be constrained.

4 Analytical techniques

4.1 Grain size analysis

Samples were prepared using standard procedures, with material being taken from each of the 11 samples that were the subject of zircon U–Pb dating. We put a small amount of sample into a cleaned 50 mL plastic centrifuge tube and added 5–7 mL of sodium phosphate solution. The tube was capped and vortexed to deflocculate clay-sized sediment and separate organic particles. The sample was poured through an 850 µm sieve and funneled into a 15 mL glass test tube. After centrifuging and removing the clear supernatant, 2–3 mL of sodium phosphate and 5 mL of 30 % H₂O₂ were added. Tubes were vortexed again and then put into a hot bath that was heated to 70 ^∘C. This step requires persistent monitoring to prevent the loss of reactant by spraying it with acetone until the reaction is stabilized. Reactants then sat overnight to completely oxidize organic matter. Reacted supernatant was removed and 5 mL of sodium phosphate was added. These treated samples were then rinsed with deionized water, transferred into clean 50 mL plastic centrifuge tubes, and topped with sodium phosphate into a sample solution of up to 40 mL. Samples were vortexed again prior to grain size analysis. Grain size analysis was conducted on a Beckmann Coulter LS13 320 laser diffraction particle size analyzer at Louisiana State University (LSU). The obscuration of all running samples in the aqueous liquid module (ALM) was between 8 % and 12 %. Results are provided in Table S1 in the Supplement.

4.2 Major element analyses

Bulk samples of the sands targeted for zircon U–Pb dating were analyzed for their major element contents by inductively coupled plasma emission spectrometry (ICP-ES) at Boston University, USA. Sediment samples were decarbonated with acetic acid, washed with distilled and deionized water with a purity of 9–12 MΩ cm, and hand-powdered before total fusion preparation. Glass beads for each sample were made in a muffle furnace under 1050 ^∘C by fusing 100±0.5 mg of sample mixed with 400±0.5 mg lithium metaborate (LiBO₂). The melted mixture was then dissolved in 5 % HNO₃, sonicated, manually shaken until no visible grains were observed, and further diluted for analysis (Dunlea et al., 2015). Precision for all elements was better than 1 % of the measured value, and accuracy was confirmed by repeated analyses of International Standard Reference Materials (Basalt, Hawaiian Volcano Observatory, BHVO-2) (Wilson, 1997). Results of the geochemical measurements are shown in Table S2.

4.3 Isotope analysis

Nd and Sr isotopes were measured from powdered whole sediment samples. After decarbonation with 10 % acetic acid and dissolution, Sr and Nd were concentrated using standard column extraction techniques, and isotopic compositions were determined by a Thermo “Neptune” multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS) at Woods Hole Oceanographic Institution. Sr results were corrected against NBS 987 standard ${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}\mathrm{Sr}=\mathrm{0.710240}$, and Nd data were corrected against JNdi-1 standard ${}^{\mathrm{143}}\mathrm{Nd}{/}^{\mathrm{144}}\mathrm{Nd}=\mathrm{0.512104}$. We calculate the parameter ε_Nd (DePaolo and Wasserburg, 1976) using a depleted-mantle model age and ¹⁴³Nd∕¹⁴⁴Nd value of 0.512630 for the Chondritic Uniform Reservoir (CHUR; Bouvier et al., 2008). Results of the geochemical measurements are shown in Table S3.

4.4 Zircon U–Pb dating

After standard mineral separation, zircon grains were sprinkle-mounted onto double-sided tape on 1^′′ acrylic discs and analyzed at random using depth-profiling LA (laser ablation)-ICP-MS U–Pb geochronology. For each sample at least 120 zircons were analyzed to obtain provenance datasets that resolve any component comprising >5 % of the total population (Vermeesch, 2004). The analyses were completed using a PhotonMachine Analyte G.2 Excimer laser (30 µm laser spot size) with a large-volume Helex sample cell and a Thermo Element2 ICP-MS at the UTChron facilities at the Jackson School of Geosciences at the University of Texas at Austin using procedures described in Hart et al. (2016). GJ1 was used as the primary reference standard (Jackson et al., 2004) and a secondary in-house zircon standard (Pak1 with a thermal ionization mass spectrometer (TIMS) ²⁰⁶Pb∕²³⁹U age of 43.0 Ma). The data from the analyses were then reduced using the Iolite data reduction software and VizualAge (Paton et al., 2011; Petrus and Kamber, 2012). For analyzed detrital zircons, the ²⁰⁶Pb∕²³⁸U age was used for grains younger than 850 Ma and the ²⁰⁷Pb∕²⁰⁶Pb age was used for grains older than 850 Ma (Gehrels et al., 2008). All ages reported use 2σ absolute propagated uncertainties, ²⁰⁷Pb∕²⁰⁶Pb ages are less than 30 % discordant, and ²⁰⁶Pb∕²³⁸U ages are less than 10 % discordant (Gehrels et al., 2011). The discordance reported is calculated with the ²⁰⁶Pb∕²³⁸U and ²⁰⁷Pb∕²³⁵U ages if <850 Myr and the ²⁰⁶Pb∕²³⁸U and ²⁰⁷Pb∕²⁰⁶Pb ages if >850 Myr. Results of the geochemical measurements are shown in Table S4.

5 Results

Sediment grain size can be assessed using the classification scheme of Folk (1974) (Fig. 2a). Sediments are dominantly silty sand and sandy silt, with one sample defined as a silt (U1456A-11H-6W). If we consider the range of grain sizes in any given sample (Fig. 2b), we see a generally good sorting (positive kurtosis) and a positive skew, i.e., a dominance of the finer grain sizes and a tail of coarser grains comprising a diminishing proportion of the sediment. Very little material of coarse sand size is seen in any of the material.

Figure 2 (a) Grain size range of all samples analyzed for U–Pb zircon dating from Laxmi Basin shown on the scheme of Folk (1974). (b) The detailed grain size spectra of the analyzed samples.

The general geochemical character of the sediments can be seen on CN-A-K ($\mathrm{CaO}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}-{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}-{\mathrm{K}}_{\mathrm{2}}\mathrm{O}$) ternary diagram (Fedo et al., 1995) (Fig. 3a). The data plot in an array close to a Chemical Index of Alteration (CIA) of ∼70 (Nesbitt et al., 1980). The geochemical analyses from the samples form a roughly linear array trending towards the illite end-member and suggestive of its progressive involvement as the primary mineral product of breakdown. Here we compare the Laxmi Basin samples with sediments from the Indus Delta (Clift et al., 2010) and Indus Canyon (Li et al., 2018), as well as sediments from the western Indian shelf and slope between Saurashtra Peninsula and Bhatikal (Kurian et al., 2013) (Fig. 1). The samples from the Indian shelf largely lie offshore from extensive exposures of the Deccan Trap flood basalt sequences. We note that the Laxmi Basin sands have very similar bulk compositions compared to the Quaternary Indus Canyon and Delta but plot below the array of the Indian shelf sediments. Only the sample taken near Bhatikal plots lower than the Laxmi Basin sediments. This plot confirms that the analyzed sands have little in common with material eroded from the Indian Peninsula and are consistent with an Indus River origin. Likewise, the Laxmi Basin sediments plot overlapping with the Quaternary Indus sediments on the discrimination diagram of Herron (1988) (Fig. 3b). These sediments form an array defined as wackes and litharenites, while the western Indian shelf sediments fall into the Fe shale and Fe sand fields.

Figure 3 (a) Geochemical signature of the analyzed samples illustrated by a CN-A-K ternary diagram (Fedo et al., 1995). CN denotes the mole weight of Na₂O and CaO* (CaO* represent the CaO associated with silicate, excluding all the carbonate). A and K indicate the content of Al₂O₃ and K₂O, respectively. Samples closer to A are rich in kaolinite, chlorite, and/or gibbsite (represented by Kao, Chl, and Gib). CIA values are also calculated and shown on the left side, with its values correlating with CN-A-K. Samples from the delta have the lowest values of CIA and indicate high contents of CaO and Na₂O and plagioclase. Abbreviations: sm (smectite), pl (plagioclase), ksp (K-feldspar), il (illite), m (muscovite). (b) Geochemical classification sediments from the Indus Delta (Clift et al., 2010), Indus Canyon (Li et al., 2018), and western Indian Peninsula rivers (Kurian et al., 2013) following the scheme of Herron (1988).

Nd and Sr isotope compositions of silty and muddy sediments from the boreholes are shown in Fig. 4 and are compared with source regions in the Greater Himalaya, Karakoram and Transhimalaya, as well as the Deccan Plateau of the neighboring Indian Peninsula. The Lesser Himalayan Series have very low ε_Nd and high ⁸⁷Sr∕⁸⁶Sr values and lie outside this plotted area, but still contribute material. The Lesser Himalayan influence on the fan sediment is not resolvable from the Greater or Tethyan Himalaya using the Nd and Sr isotopic systems. Many of the analyses overlap the range of bedrocks in the Karakoram, allowing for these ranges also to be an important source. However, on the basis of these isotopic data alone the relative influence of basement sources cannot be quantified. Interestingly, several IODP samples <10.8 Myr old, as well as those from the high-resolution study by Yu et al. (2019) spanning the past 0.6 Myr, often have higher ε_Nd values than samples from Indus Marine A1 or sedimentary rocks of the Kirthar and Sulaiman ranges of southern Pakistan (Fig. 5). Both of these latter sequences are generally considered to be ancient deposits of the lower reaches of the Indus River and are largely free from other influences (Zhuang et al., 2015). The IODP samples also have higher ε_Nd values and lower ⁸⁷Sr∕⁸⁶Sr values than many Quaternary sediments from the Indus Delta (Clift et al., 2010) (Fig. 4).

Figure 4 Diagram showing the Nd and Sr isotopic compositions of sediments recovered from the Laxmi Basin drill sites and Indus Marine A-1 compared to the major possible source terrains. Inset map shows the full isotopic range of the source ranges. Data from the Deccan Plateau compiled from Georoc (http://georoc.mpch-mainz.gwdg.de/georoc/, last access: 16 November 2015). Pink squares show the composition of selected Indian Peninsula rivers (Goswami et al., 2012). Compositions for the Indus Delta since 15 ka are from Clift et al. (2010). Data from the Kirthar and Sulaiman ranges are from Zhuang et al. (2015). Transhimalaya data are from Rolland et al. (2002), Singh and France-Lanord (2002), and Khan et al. (1997). Greater Himalayan data are from Ahmad et al. (2000), Deniel et al. (1987), Inger and Harris (1993) and Parrish and Hodges (1996). Karakoram data are from Crawford and Searle (1992) and Schärer et al. (1990).

Figure 5 Evolution of Indus related Nd isotope composition since 17 Ma. (a) New data from IODP Sites U1456 and U1457 and Indus Marine A-1 together with existing data from the latter, from IODP Site U1457 measured at high resolution by Yu et al. (2019), from ODP sites in the western Arabian Sea, and from the Kirthar and Sulaiman ranges. The yellow-shaded area shows the ±1ε_Nd uncertainty for those data not at Sites U1456 and U1457. Samples plotting to the right of this range could be influenced by flux from Indian Peninsula. (b) Data from all the sources in the previous diagram with peninsular-related samples excluded and a five-point running average plotted with the. ±1ε_Nd uncertainty. Note the decline in ε_Nd values after 6 Ma and especially after 3 Ma.

U–Pb zircon ages in all samples span a wide range, but all show a significant detrital age component <120 Myr, comparable to bedrock ages from the Karakoram and Kohistan, as shown in the kernel density estimate diagram of Vermeesch (2012) (Fig. 6). In general, the abundance of grains with U–Pb ages of 750–1250 and 1500–2300 Ma increases up-section. The 300–750 Ma component is present throughout, mostly at 20 %–30 % of the total, peaking briefly at 39 % at 5.9 Ma and falling from 30 % at 3.4 Ma to 18 % in the present day. The 1500–2300 Ma detrital age component shows an initial increase from 5 % to 17 % after 8.3 Ma and shows another jump from 20 % to 35 % of the total load after 1.9 Ma (Fig. 6).

Figure 6 Kernel density estimate (KDE) plots for the zircon U–Pb ages from Laxmi Basin compared with some of the major source terrains in the Western Himalaya, from the compilation of Alizai et al. (2011). Colored strips depict characteristic age peaks that can be linked to specific source areas within the Western Himalaya. Left column shows KDE plots to 3000 Ma and the right-hand column shows the data from 200 Ma.

6 Discussion

6.1 Shifting erosion patterns

The long-term temporal evolution in Nd isotope character can be used to reconstruct the erosion of source rocks and exhumation patterns. Long-term trends in Nd isotopic composition are highlighted by a five-point running average that defines a trend to slightly increasing ε_Nd values from 17 to ∼9.5 Ma, followed by a period of approximate stability until 5.7 Ma and declining values after that time. Phases of steep decline were identified, most notably after 5.7 Ma and after 3 Ma (Fig. 5). The tendency of some IODP samples to have higher ε_Nd values compared to previously analyzed Indus-sourced sediment could reflect a non-Indus contribution to the IODP sites from sediment eroded from the Deccan Plateau and delivered by the Tapti River (Fig. 1). Rivers draining the peninsula and located further south than the Tapti River are not considered to be likely suppliers to the area of the drilling sites because of the southward-directed longshore currents (Shetye et al., 1994), as well as the contrasting bulk sand samples' geochemistry character (Fig. 3). Alternatively, this difference might reflect relatively greater erosional flux from the Transhimalaya or Karakoram during some of the time since 10.8 Ma in a way not seen in earlier analyses (Zhuang et al., 2015). Significant flux from a paleo-Mahi River (Fig. 1) can be ruled out due to its low ε_Nd and high ⁸⁷Sr∕⁸⁶Sr values, although this river could have had higher ε_Nd values in the past when much of this catchment was covered by Deccan volcanic rocks that have now been removed by erosion. Because no independent record of Mahi River discharge exists, this possibility remains speculative. The Narmada River is isotopically indistinguishable from the post-15 ka Indus (Fig. 4).

Indus Marine A-1, located close to the river mouth, and the Kirthar and Sulaiman ranges that have preserved the lower reaches of the Indus River may be considered to represent a relative pure Indus signature (Zhuang et al., 2015). This allows for a better assessment of whether the more distal sediments in the Laxmi Basin were truly derived from this river or may be influenced by significant flux from the Indian Peninsula. Differences between the IODP sample compositions and the isotope record from closer to the Indus River mouth could reflect a number of processes. For example short-term changes in erosion patterns related to climate change and reproducibility of Nd isotope composition, generally considered to be ±1ε_Nd (Jonell et al., 2018), could also result in significant variability. We suggest that when analyses from the Laxmi Basin depart by >1ε_Nd from the trend of the Indus Marine A-1, Sulaiman and Kirthar, or western Indus Fan Ocean Drilling Program (ODP) samples, then they likely do not represent only flux from the Indus River. This is especially true when ε_Nd values increase, suggesting mixing with sediment from the neighboring isotopically positive Deccan Plateau, consistent with the interpretation of Yu et al. (2019) since 0.6 Ma at IODP Site U1457. Significant variations in Indus River ε_Nd values have been observed over recent glacial cycles and have been linked to changes in monsoon strength driving rapid change in erosion patterns in the mountains (Clift et al., 2010), suggesting that ε_Nd values more positive than this range may indicate some flux from the Indian Peninsula.

Sediments that have ε_Nd values more positive than those sediments deposited closer to the Indus River mouth require mixing with additional sources outside anything known from within the Indus catchment. Such analyses are considered to be at least partly derived from the Indian Peninsula (Fig. 5). Sediments with such positive ε_Nd values are only found in modern rivers in the Karakoram, not near the delta (Clift et al., 2002) and would require supply from a river with essentially no drainage of Himalayan source regions, which seems unlikely. Mixing with the sediment from the peninsula is more likely. Indeed, sedimentary petrography has identified peninsula-derived material in sands in the uppermost parts of the IODP section (Pandey et al., 2016), but the potential for erosion from these areas lower down in the section remains. Eolian sediment derived from Arabia is not considered to be important, even in those parts of the section where sedimentation was slow. When sedimentation rates were ca. 10 cm kyr⁻¹, as they are in the muddier parts of the stratigraphy, then using a dry density of ca. 1.6 g cm⁻³, equates to a dry mass accumulation rate (MAR) of 160 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{ky}}^{-\mathrm{1}}$ total sediment (Pandey et al., 2016). In contrast, the nearest sediment traps from the Arabian Sea showed lithogenic MARs of only 1.05–0.29 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{kyr}}^{-\mathrm{1}}$, with a clear reduction away from Arabia (Honjo et al., 1999). It thus seems unlikely that windblown material would have contributed much of the sediment in Laxmi Basin, even when MARs were low.

Based on our understanding of the primary sources of the Indus River the simplest explanation for the long-term increasing trend in ε_Nd values between 17 and 9.5 Ma would be a relative increase in erosion from either the Karakoram, Kohistan, or Transhimalaya compared to Himalayan sources, followed by a reversal since 7.0 Ma. A simple two-part mixing model using ${\mathit{\epsilon }}_{\mathrm{Nd}}=-\mathrm{14.7}$ as a representative value for the Greater Himalaya, based on the sources plotted in Fig. 4 and ${\mathit{\epsilon }}_{\mathrm{Nd}}=-\mathrm{8.0}$ for a combined Karakoram–Kohistan source based on bedrock analyses and the range of Indus Fan sediment compositions, yields a range of contributions from only 15.4 % in Karakoram at 1.13 Ma up to 96.3 % at 94.8 % at 6.74 Ma. This estimate is based on equal Nd contents of the sources, consistent with modern river sediment data (Alizai et al., 2011). The detrital U–Pb zircon ages also show temporal evolution, although these data are limited to the past 8.3 Myr with the exception of a single sample of 119 grains dated to ∼15.5 Ma, recovered from below the mass transport complex (Fig. 6). Because all of the sands contain grains <120 Myr, this requires derivation from the Indus River and not the Indian Peninsula (Fig. 6).

The implications of the zircon U–Pb dating can be better understood by subdividing the ages into provenance diagnostic groups. Zircon grains <25 Ma are found in igneous rocks at the core of the Nanga Parbat Massif (Zeitler et al., 1993). Such young grains are unknown elsewhere in the Indus catchment, although older grains dated to ∼1850 and 400–500 Ma are also known to be common in the Nanga Parbat Massif (Zeitler et al., 1989), reflecting the correlation of the massif with the Lesser Himalaya based on Nd isotopes (Whittington et al., 1999). Zircons dating to 40–70 Ma are associated with both Kohistan and the Karakoram, while 70–120 Ma grains are largely found only in the Karakoram (Searle, 1996). In both the Greater and Tethyan Himalaya, 300–750 and 750–1250 Ma grains are found (Alizai et al., 2011; DeCelles et al., 2000; Gehrels et al., 2011), while 1500–2300 Ma zircons are most frequent in Lesser Himalayan rocks (DeCelles et al., 2000). Zircon grains dating to 1500–2300 Ma are also found in other Himalayan units and comprise 38 % of the Greater and Tethyan Himalayan analyses compiled; however, 99 % of Lesser Himalayan grains fall in that age range. This observation indicates that Lesser Himalayan units are likely most critical in controlling the supply of 1500–2300 Ma grains, especially when 750–1250 Ma grains do not increase at the same time.

The major element analyses of the sands suggest that they are similar to Quaternary Indus Canyon and Delta sediments and unlike western Indian shelf sediments (Fig. 3). The Indian Peninsula is not affected by Himalayan mountain building, and magmatism and fission track data show progressive, modest cooling of western India since the Cretaceous that would not affect the U–Pb ages presented here (Gunnell et al., 2003). The youngest zircon U–Pb grain age that can reasonably be expected from an Indian Peninsular river is ∼65 Ma, derived from the Deccan Plateau (Schoene et al., 2015), but no such peak is noted in the data from the Laxmi Basin sands (Fig. 6). This is consistent with the sands not being derived from the Indian Peninsula. The changing abundance in the 1500–2300 Ma zircon group associated with the Lesser Himalayan does not precisely mirror changes in the Nd isotope character. For example, this group leaps from 19.6 % to 34.6 % of the total population between 1.9 and 1.5 Ma, while at the same time ε_Nd values fell only slightly after a large fall between 3.4 and 2.5 Ma. Nonetheless, it is this group in particular that shows a clear relative increase in abundance over the long term since 7.0 Ma and which is also associated with bedrock sources with very negative ε_Nd values. Zircons of this age comprise just 5 % at 8.3 Ma but 17 % by 5.9 Ma, even before the large increase to >34 % after 1.9 Ma. In contrast, the 300–750 Ma group falls slightly in abundance after 8.3 Ma. The 750–1250 Ma group shows an increase from 20 % at 5.9 Ma to 32 % by 3.0 Ma. This suggests that it is a source rich in 1500–2300 Ma grains, and to a lesser extent 750–1250 Ma, that changes most in its contribution to the fan since 8.3 Ma, especially after 1.9 Ma.

Provenance evolution can be further assessed by comparing the spectra of the detrital zircons with bedrock data using a multidimensional scalar diagram (Vermeesch et al., 2016). This is a form of principle component analysis that statistically compares the U–Pb age spectra. Samples that are similar to one another plot closely together on the diagram. Figure 7 shows that Himalaya sources plot on the left side of the diagram and Karakoram and Kohistan sources on the right. Nanga Parbat appears to plot in the upper center, but this is based only on bedrock data from the center of the massif, not including the older zircons that probably occur at the periphery of the massif. Not surprisingly the modern trunk Indus River plots closely to the Karakoram, as do the older samples from the Laxmi Basin. There is a clear long-term trend in the detrital samples towards more and more Himalayan compositions, consistent with the analysis of the KDE plots in Fig. 6. This figure indicates that the Greater and Tethyan Himalaya is likely the most important sources for the youngest samples (0.9 and 1.6 Ma), which are similar to the modern river mouth.

Figure 7 Multidimensional scalar (MDS) plot comparing the detrital samples from the Laxmi Basin (green dots) with potential source regions onshore (red dots) and both the modern river mouth and mainstream Indus downstream of Nanga Parbat (blue dots). Data sources are as for Fig. 6. Note progressive migration away from Karakoram–Kohistan end-members and towards the Greater and Tethyan Himalaya.

Sediment storage and recycling are unlikely to be affecting the erosional signal. Although storage in basins between source and sink does influence the Himalaya, this is largely of millennial duration and rarely >100 kyr (Blöthe and Korup, 2013), a value still too small to affect the reconstruction presented here. Recent work on intermontane basins in the Indus catchment estimates sediment buffering on timescales of 10³–10⁴ years (Jonell et al., 2017), while the erosional response to deglacial climate change allowed zircon travel times in the Indus to be constrained to 7–14 kyr (Clift and Giosan, 2014).

6.2 Unroofing the Lesser Himalaya

Nanga Parbat related 0–25 Ma grains are present throughout the sedimentary section in Laxmi Basin but are never very abundant, although they do increase from 1.8 % to 5.1 % of the total between 3.4 and 3.0 Ma (Fig. 8). Their modest overall contribution limits the role of Pliocene unroofing of Nanga Parbat in driving the Nd isotope evolution (Chirouze et al., 2015). However, the presence of 300–750 and 750–1250 Ma grains in all samples demonstrates the long-term flux from the Greater and Tethyan Himalaya to the Indus Fan and allows us to rule out large-scale drainage capture of eastern tributaries as the cause of falling ε_Nd values since 5.7 Ma (Clift and Blusztajn, 2005). We note that ε_Nd values started to fall after 5.7 Ma, close to the time when the first Inner Lesser Himalayan units, characterized by their very negative ε_Nd values ($<-\mathrm{22}$), were first exposed (Najman et al., 2009) but also at a time when the 750–1250 Ma zircon grains became more abundant. Although 1500–2300 Ma grains, characteristic of the Lesser Himalaya, started to become more abundant after 8.3 Ma, the sharpest increase in their population was after 1.9 Ma. The shift towards the Lesser Himalaya is shown by the fact that at 7.8 Ma the 750–1250 Ma (Greater and Tethyan Himalaya) and 1500–2300 Ma (Lesser Himalaya) age populations accounted for 26 % and 17 % of the total zircons, respectively, while by 0.9 Ma the proportions had changed to 32 % and 38 %. Such a change can only be achieved by increasing Lesser Himalayan erosion even more than erosion in the Greater and Tethyan Himalaya, which also rose. Sedimentation rates in the submarine fan increased from the late Miocene into the Pleistocene (Clift, 2006) so this change in relative influence cannot indicate reduced erosional flux from the Greater and Tethyan Himalaya.

Figure 8 Comparison of climate, erosion, and exhumation proxies in the Himalaya. (a) Smooth Nd isotope history for the Indus River from the Arabian Sea with grey background showing effective uncertainties from Fig. 5. (b) Breakdown of detrital zircon U–Pb age populations with each source diagnostic group shown as its proportion of the total. Percentage is recalculated to exclude nondiagnostic aged grains. (c) Carbon isotope character of pedogenic carbonate in Pakistan as an indicator of dominant vegetation in the Potwar Plateau of Pakistan (Quade et al., 1989) and NW India (Singh et al., 2011). (d) Relative exhumation rates of the Greater Himalaya tracked by bedrock Ar–Ar dating (Clift et al., 2008) and zircon fission track from foreland basin sediment (Chirouze et al., 2015). (e) Rates of sediment supply to the Arabian Sea calculated from regional seismic reflection surveys (Clift, 2006). KK: Karakoram; LH: Lesser Himalaya; LHCS: Lesser Himalayan crystalline series.

Furthermore, the abundance of zircons in the source rocks that contribute to the Indus mainstream is relatively low compared to the Eastern Himalaya draining tributaries. Zr can be used as a proxy for zircon content because of its high concentration in this mineral (Amidon et al., 2005). The Indus mainstream sample from Attock records a Zr concentration of ∼18 ppm compared to a range of 14 to 63 ppm for the Himalayan tributaries (Alizai et al., 2011). Using water discharge data as a proxy for sediment transport capacity, we estimate an average Zr concentration for the modern Punjabi tributaries of 40 ppm, more than double that in the mainstream (Alizai et al., 2011). Zr is even higher in the Laxmi Basin sediment (142–223 ppm) reflecting hydrodynamic sorting during transport.

This discrepancy in zircon source abundance is important because although the youngest Laxmi Basin sediment contains 79 % of grains older than 300 Ma, this does not mean that 79 % of the erosion is from the Himalaya, rather than the Karakoram, Kohistan, and the Transhimalaya. Indeed, petrographic data (Garzanti et al., 2005) and Nd isotope constraints (Clift et al., 2001) indicate that only ∼39 % and 41 %, respectively, of the total flux to the Arabian Sea is now from Himalayan sources, with the Karakoram accounting for much of the rest. The low zircon fertility along the mainstream upstream of Attock means that it is unlikely that uplift of Nanga Parbat could be driving the changes in zircon reported in this study despite the fact that there are Lesser Himalayan rocks within the massif. It is noteworthy that the same is not true of Nd concentrations because the Indus mainstream contains 27 ppm Nd at Attock, compared to a range of 20–29 ppm in the Himalayan tributaries (Alizai et al., 2011). Because the two areas have similar Nd concentrations the evolution in Nd isotope composition can be interpreted to indicate more erosion from both the Lesser Himalaya and Nanga Parbat since 5.7 Ma. The modern river data indicate that ∼10 % of the total modern zircon flux in the trunk river upstream of Attock is from Himalayan bedrocks (>300 Ma), i.e., not from the eastern tributaries. The >300 Ma zircons account for just 16 % of the zircon flux from the upper mainstream. In contrast, the mainstream supplies 61 % of the total sediment flux to the Arabian Sea based on Nd isotope constraints (Alizai et al., 2011). Overall, 96 % of the zircon flux from the eastern Punjabi tributaries is >300 Ma, but these streams account for only 39 % of the total Indus sediment load reaching the Arabian Sea (Alizai et al., 2011). If we estimate that the Punjabi tributaries are 2.2 times more enriched in zircon, but around the same in Nd, compared to the mainstream Indus River, then this implies that the evolution to more ε_Nd negative values seen offshore is 35 % driven by erosion from Nanga Parbat and 65 % driven by erosion in the Lesser Himalaya.

The lack of a tight correlation between Nd and U–Pb zircon data reflects the fact that these different analyses were not both performed on the same samples, although they both show the same long-term drift to more erosion of older Himalayan crust and less from the Karakoram and Kohistan. Short-lived changes in erosion pattern would account for mismatches in provenance between samples with similar but different depositional ages. We also recognize that Nd and Zr are not uniformly concentrated in the sources along the mainstream and Punjabi tributaries and that the fine-grained sediments analyzed for Nd may have a different provenance than the sandy zircon material; e.g., coarser material may be preferentially derived from the Karakoram, while finer material is more Himalayan (Jonell et al., 2018) or derived from the Indian Peninsula.

6.3 Controls on erosion

The importance of Karakoram and Kohistan for the erosional flux in zircon grains is clear until after 7.0 Ma and may be related to motion on the Karakoram Fault that started after 16 Ma (Phillips et al., 2013) and that drove rock uplift and unroofing along the length of that structure. The reduction in erosion is despite thermochronology data from the K2 mountain region that emphasize Pliocene cooling of the central Karakoram after slower exhumation earlier in the Cenozoic (Foster et al., 1994). The lack of any detectable signature in the fan related to material from K2 likely reflects its small area compared to the other potential sources to the Indus Fan and the fact that much of the sediment is transported north into the Tarim Basin (Clift et al., 2017). The fan record provides a basin-wide average image of erosion but is not good at identifying dramatic tectonic events in limited areas that do not produce much sediment, just as also noted in the Nanga Parbat Massif. Away from the Karakoram Fault itself the gneiss domes of the southern Karakoram show rapid uplift and cooling after 20–25 Ma (Rolland et al., 2001), potentially contributing to the trend to more positive ε_Nd values from 17 to 9.5 Ma. However, exhumation rates in the southern Karakoram continued to be rapid into the Pleistocene (Mahéo et al., 2004), while the sediment data require them to be outpaced by Himalayan erosion since 7.0 Ma. However, erosion in the eastern Karakoram has a markedly different history, with faster rates of exhumation in the middle and late Miocene followed by a slowing after 5 Ma (Wallis et al., 2016), consistent with the zircon and Nd isotope evolution presented here. Because the summer monsoon does not reach the Karakoram (Karim and Veizer, 2002), changes in the erosion of those mountains are unlikely to be linked to changes in the summer rains but could be related to the strength and location of the westerly jet which supplies the glaciers that control erosion in that region.

Neither Nd isotopes nor zircon U–Pb ages exhibit large changes in provenance between 15.5 and 8.3 Ma, implying that the Tethyan and maybe the Greater Himalaya was already unroofed by ∼15 Ma and continued to erode through that time interval. It is noteworthy that the initiation of greater relative contributions from the Lesser Himalaya after 7.0 Ma as tracked by zircon populations and after 5.9 Ma based on Nd isotopes commenced at a time of reduced monsoon rainfall across the foreland basin as tracked by the relative abundance of C₃ versus C₄ vegetation (Dettman et al., 2001; Quade et al., 1989; Singh et al., 2011) (Fig. 8). The Lesser Himalayan related 1500–2300 Ma grains increased from 5 % to 17 % of the total zircon population between 8.3 and 7.8 Ma, while the 750–1250 Ma group only rose from 23 % to 26 % and the 300–750 Ma group fell from 26 % to 20 %. The trend in erosion at this time is strongly towards the Lesser Himalaya. Subsequently, the 750–1250 Ma group shows a steadier rise from 5.9 to 3.0 Ma, while the 300–750 Ma group is more erratic and peaked at 5.9 Ma. The start of a progressive increase in zircon flux from the Himalaya after 8.3 Ma occurred during a time of overall slowing erosion and drying climate, as inferred from the regional sediment flux budget to the fan (Clift, 2006) (Fig. 7). Falling intensities of chemical weathering and slowing total sediment delivery rates from the late Miocene to the Pliocene (Clift et al., 2008) mean that the total weathering flux must be decreasing. Consequently, if Himalayan uplift has caused more silicate weathering, drawing down CO₂ and cooling the global environment (Raymo and Ruddiman, 1992), this is because of processes in the Eastern Himalaya and not those in the Indus catchment.

Lesser Himalayan erosion must be considered in two phases. The Outer Lesser Himalaya was likely eroding in NW India since 16 Ma (Colleps et al., 2018). The structurally overlying Inner Lesser Himalayan units lie directly under the Greater Himalaya sequences and contrast in Nd isotope character (ε_Nd $<-\mathrm{22}$) with the Outer Lesser Himalaya that overlaps with the Greater and Tethyan Himalaya in terms of Nd isotopes (ε_Nd $>-\mathrm{18}$; Fig. 4) (Najman et al., 2009). Nd isotope work on foreland basin sedimentary rocks has indicated that the unmetamorphosed Inner Lesser Himalaya was first exposed after 9 Ma, at least locally in the Beas River area (Najman et al., 2009) and in Pakistan (Meigs et al., 1995) (Fig. 1), earlier than the fall in ε_Nd values in the new marine data after 5.9 Ma. These foreland Nd data further indicate that the very ε_Nd negative rocks of the metamorphosed Inner Lesser Himalaya were exposed after 6 Ma (Fig. 8), which correlates more closely with the offshore record of a decrease in ε_Nd values at that time. However, the Beas area foreland records do not record the steep rise in Inner Lesser Himalayan erosion after 3.4 Ma observed in the Indus Fan (Najman et al., 2009). This implies that although the metamorphosed Inner Lesser Himalaya was exposed locally in the Beas area earlier in the Miocene, as shown by the change in isotope character in the fluvial sediments in that part of the foreland basin, their widespread unroofing across the entire Indus catchment is younger than previously recognized. This prediction is however consistent with recent palinspastic structural reconstructions of the Western Himalaya (Webb, 2013). Thermochronology from the bedrocks of the Lesser Himalaya in the Indus catchment indicate that exhumation was rapid from 10 to 6 Ma but that the rocks now exposed at the surface were still deeply buried (12–15 km) at that time (Caddick et al., 2007). Lower temperature thermo-chronometers including zircon and apatite fission track constraints from the Lesser Himalaya in the Sutlej Valley point to continued rapid cooling after 6 Ma (Vannay et al., 2004), but we can only constrain the exhumation of the rocks now exposed at the surface. They cannot inform us of when these units were first unroofed. This is the role of the erosional record we present here.

It is not apparent how the evolution in exhumation since the late Miocene might be linked to climatic evolution. More vigorous Lesser Himalayan erosion parallels a weakening monsoon and slower overall erosion (Fig. 8). When the monsoon is strong the occurrence of extreme summer rainfall increases (Turner and Annamalai, 2012) and when this occurs strong rainfall penetrates further north into Himalayan valleys and drives erosion of those rock units (Bookhagen et al., 2005). Consequently, late Miocene monsoon weakening would be expected to involve a southward retreat of strong rainfall from the edge of the Tibetan Plateau, resulting in more rainfall and erosion over the area where the Lesser Himalayan rocks are exposed.

Inner Lesser Himalayan thrust sheets were imbricated above a ramp in the basal decollement (Webb, 2013). This stacking resulted in surface uplift and generated a break in the topography against which the summer monsoon rains began to precipitate, focusing erosion and bringing these rocks to the surface. A positive feedback is likely because modeling has shown that focused denudation can encourage the formation of duplexes in thrust belts (Malavieille, 2010). Exposure of the Inner Lesser Himalaya on a regional scale happened later than has been previously recognized (Meigs et al., 1995; Najman et al., 2009). Zircon $(\mathrm{U}-\mathrm{Th})/\mathrm{He}$ thermochronology argues for a start of cooling in the Inner Lesser Himalaya after 11 Ma (Colleps et al., 2018), but this does not require widespread unroofing of this unit until the later times indicated by this study, i.e., 3.4–1.9 Ma. The final regional exposure does not correlate either with changes in the vegetation in the foreland basin and only roughly with the onset of Northern Hemisphere glaciation that is often believed to enhance erosion at 2–4 Ma (Zhang et al., 2001). Although it does correlate with a strengthening of the monsoon derived from nitrogen records from the same drill sites (Tripathi et al., 2017), it is unclear if this marine biogenic production record is closely related to Himalayan rainfall or not. Even if the Lesser Himalayan unroofing is not linked with monsoon intensity this does not preclude the older Greater Himalaya (Fig. 1c) from being more intimately linked to summer monsoon development. However, the evidence of this study suggests that the monsoon has been reacting more passively to structural changes in the Himalaya driven by solid Earth tectonic processes since the late Miocene and possibly since ∼15 Ma.