Orogenesis and the exhumation of mountain belts, combined with erosion and distribution of their sediments, are important Earth surface processes that respond to forcing mechanisms such as climate change, tectonics, and internal fluvial system dynamics. There continues to be much debate on the relative importance of these various processes in shaping Earth's landscapes (e.g., Molnar and England 1990;Raymo and Ruddiman 1992;Burbank et al. 2003;Clift 2006;Willenbring and Jerolmack 2016;Hessler et al. 2018), the rates over which they operate (e.g., Molnar 2004;Ganti et al. 2016;Perron 2017), and the potential for storage of these signals in the stratigraphic record (e.g., Jerolmack and Paola 2010;Romans et al. 2016). The Himalaya in particular has served as a natural laboratory for emerging understanding of coupled tectonicsedimentary systems (Clift et al. 2001;Najman 2006;Inam et al. 2007;Yin 2010), with the deep sea Bengal-Nicobar Fan (BNF) potentially preserving environmental signals propagated across the source-to-sink system (Hessler and Fildani 2019). Furthermore, hydroclimate and temperature proxies stored in the sediments of large delta-fan systems can inform the magnitudes and rates of carbon storage and release that modulated climate change throughout the late Cenozoic (Galy and Eglinton 2011;Hein et al. 2017;Hein et al. 2020). In this paper we examine the stratigraphic record of the Indo-Burman Ranges (IBR) to constrain the timing of deposition, burial, and exhumation of shallow and marginal marine sediments of the ancestral Brahmaputra delta that fed the BNF, and illuminate eastern India's location as central to evolving patterns of sediment dispersal and tectonic deformation through the Neogene. Ultimately, as our understanding of coupled surface and tectonic processes in this region deepens, so too does our understanding of how the planet responds to perturbations associated with climate change, coastal vulnerability, and geohazards that impact the lives of millions in this densely populated region.
Sediments of the Ganges-Brahmaputra Delta (GBD) and BNF record the Cenozoic to present uplift and denudation of the Himalayas. Detailed knowledge of the facies distribution, stratigraphic architecture, and surface processes in the GBD is emerging (Wilson and Goodbred 2015;Reitz et al. 2015;Sincavage et al. 2018;Pickering et al. 2019), and recent drilling across the BNF (IODP 354;362) provides a wealth of new insights into the propagation of environmental signals from the Himalayan source area into the Bay of Bengal sink (France-Lanord et al. 2016a;McNeil et al. 2017;Blum et al. 2018;Hessler and Fildani 2019). While these studies greatly advance our knowledge of sediment dispersal and their relationship to tectonics in the region, they are restricted to spatial and temporal limitations of borehole data largely from Pliocene-Pleistocene (BNF) and Holocene (GBD) deposits. For example, the IODP BNF record spans mostly Pleistocene sediment, piercing Miocene sediment in only one core (U1451), underscoring the need for onshore studies of the up-dip Mio-Pliocene stratigraphic record of this Himalayan depositional system.
Ongoing deformation along the active Himalayan front and flanking Indo-Burman ranges (IBR) exposes a rich stratigraphic archive of the Oligocene-Pliocene delta-fan system that formed by progradation of the ancestral Brahmaputra Delta (Evans 1932;Johnson and Alam 1991;Reimann 1993;Alam et al. 2003;Cina et al. 2009;Lang and Huntington 2014;Bracciali et al. 2015Bracciali et al. , 2016;;Vadlamani et al. 2015;Lang et al. 2016;Govin et al. 2018aGovin et al. , 2018b)). Furthermore, this deformation likely played a key role in a major drainage basin reorganization whereby the ancestral Brahmaputra River migrated westward around the rising Shillong Massif sometime in the late Miocene to Pliocene (Johnson and Alam 1991;Govin et al. 2018a). Despite this deep record of the Neogene eastern Himalayan peripheral basin and its central location in this important source-to-sink system, outcrops of the IBR have received little attention until recently due to their limited accessibility.
In this paper, we present new lithofacies definitions, stratigraphic mapping, and basin thermal modeling of the Miocene-Pliocene IBR stratigraphy from outcrops located in an ~4,000 km 2 region of the IBR in central Mizoram, India (Fig. 1). We interpret the lithofacies architecture of the Neogene fluvial Tipam Group and shallow-marine Surma Group of the IBR in the context of the modern fluvial and shallow-marine environments of the GBD and constrain maximum depositional ages and basin history with detrital-apatite fission-track (dAFT) analyses and (U-Th)/He (dAHe) thermochronology and thermal modeling. Our analyses illuminate the shifting depositional environments and rates of sedimentation and exhumation of the terrestrial and shallow-marine part of the late Miocene-Pliocene ancestral Brahmaputra Delta. Furthermore, thermal modeling results reveal the onset of exhumation in the developing IBR fold-thrust belt, which led to a major reorganization of sediment delivery to the BNF in this Himalayan source-to-sink system (cf. Najman et al. 2012Najman et al. , 2016;;Govin 2018a).
Sediments in the Burma accretionary wedge near latitude 23.58 N record an evolution from Late Cretaceous to Paleocene deep marine conditions of the Disang Formation (Nandy 1999) to syn-Himalayan (Oligocene-present) shallow-marine and intertidal facies of the Barail (Oligo-Miocene) and Surma groups (Miocene), Bhuban and Boka Bil formations (Alam et al. 2003), and fluvial deposits of the Tipam and Dupi Tila groups (Miocene-Pleistocene; Alam et al. 2003;Betka et al. 2018a) (Fig. 2). A similar sequence of facies is being deposited now in the modern delta, from marine sediments in the Bengal Fan (France-Lanord et al. 2016a;McNeill et al. 2017), to shallow marine to shelf and fluvial environments farther north (Barua 1997;Michels et al. 1998;Davies et al. 2003;Kuehl et al. 2005). Importantly, the ages of the syn-Himalayan deposits in the IBR are not well known, as the existing age data for the units (as well as the unit names themselves, e.g., ''Tipam Sandstone'' vs. ''Tipam Group'') are generally extrapolated over the . 100,000 km 2 region (Alam et al. 2003) and the lithostratigraphic units are time-transgressive in this highly progradational system (Fig. 2). Because knowledge of this system has expanded since the original stratigraphic nomenclature was coined, we develop here a new nomenclature for IBR facies and associated subenvironments that are correlative with observable depositional styles on the modern GBD.
The Brahmaputra River avulsed westward to its present course in Bangladesh by the end of the Miocene or the early Pliocene to feed sediment into the modern GBD and Bengal Fan (Uddin and Lundberg 1999;Najman et al. 2012Najman et al. , 2016;;France-Lanord et al. 2016b;McNeill et al. 2017). Miocene-present uplift of the Shillong massif (Biswas et al. 2007;Clark and Bilham 2008;Govin et al. 2018a) and subsidence of the Sylhet basin on the Indian foreland (Najman et al. 2016) initiated this major shift in depositional patterns across the Bengal Basin (Johnson andAlam 1991, Uddin andLundberg 1999).
Results from IODP Expedition 354 provide insights into distal sediment dispersal patterns in the Bay of Bengal from the Neogene to Quaternary. BNF deposits overall are characterized by pelagic and hemipelagic units with relatively slow (1-2 cm/ky) accumulation rates interspersed with turbidite units with accumulation rates an order of magnitude higher (10-100 cm/ky). These facies changes are interpreted to represent channel migration across the fan through the Pliocene and Pleistocene (France-Lanord et al. 2016a;Bergmann et al. 2019). The earliest evidence of turbidite deposition is found on the eastern BNF in the late Oligoceneearly Miocene (France-Lanord et al. 2016a, 2016b), initiating a shift from predominantly pelagic Paleogene carbonates to Neogene turbidites sourced from both the ancestral Ganges and Brahmaputra rivers (France-Lanord et al. 2016b). An order-of-magnitude increase in turbidite deposition rates in the BNF (from ~1.5 cm/ky to 12 cm/ky) is noted starting at about 13.5 Ma and lasting until about 6.8 Ma, when similar increased rates of turbidite deposition shift further west (France-Lanord et al. 2016a, 2016b). Geochemical and detrital-zircon U-Pb evidence suggests temporal variability in the amount of mixing of source material for these turbidites between the two rivers (France-Lanord et al. 2016b). Autogenic avulsions on the delta plain (particularly during sea-level highstands) and the shelf margin (during sea-level lowstands) are suggested as the drivers of the amount of mixing (Blum et al. 2018). In spite of relatively rapid lateral shifts in depocenter, the overall rate (10-20 cm/ky) of fan growth is fairly uniform when integrated over long (. 1 Ma) timescales (France-Lanord et al. 2016b).
To characterize the depositional sub-environments of the paleodelta in northeast India, we undertook a targeted field mapping campaign of: 1) shallow-marine to intertidal deposits of the Surma Group (Bhuban and Boka Bil Formations), and 2) fluvial deposits of the Tipam and Dupi Tila groups. Published lithofacies assemblages from the modern GBD (Michels et al. 1998;Kuehl et al. 2005;Pickering et al. 2014Pickering et al. , 2019;;Goodbred et al. 2015;Wilson and Goodbred 2015;Reitz et al. 2015;Patrick et al. 2015;Sincavage et al. 2018Sincavage et al. , 2019) ) were used as a reference for interpreting field descriptions of IBR paleodelta outcrops in terms of equivalent modern environments (Fig. 2). SRTM 1 arcsecond (30 m) data (NASA, Jet Propulsion Laboratory) generated a hillshade digital elevation model used as a basemap for geologic contacts.
Due to thick vegetation throughout the region, mapping was focused along road-cut exposures along a series of transects crossing fault-cored anticlines of the fold-thrust belt (e.g., Betka et al. 2018aBetka et al. , 2018b)). This vegetative cover also limited the ability to access lateral continuity of facies in most instances. Architectural elements such as smaller (meter-scale) channel bodies are observable at the outcrop scale, whereas larger (tens to hundreds of meters in scale) stratal geometries can be inferred from panoramic imagery of some of the better exposures (Fig. 3). Grain-size variability was assessed qualitatively at the outcrop scale on a bed-by-bed basis, as were the nature of contacts between outcrop-scale units. Preservation of fossil material is poor in these units, so biofacies descriptions are not included, although the occurrence of fossils and trace fossils were noted where observed. Mineralogical descriptions are limited to field observations and thin-section descriptions from the quarry south of Aizawl.
Stratigraphic sections and lithofacies observations were made along two road cuts, north and south, respectively, from the city of Aizawl, Mizoram (Fig. 2). Each section documents the shift from shallow-marine to fluvial conditions across the paleodelta. South of Aizawl, a continuous (~2 km) road-cut exposure along the World Bank Road between the villages of Thiak and Sumsuih (Fig. 2) exposed ~500 m of section that was measured for stratigraphic thickness and facies variability (Fig. 3, Fig. S1, see Supplemental Material). North of Aizawl, a road cut from the Lengpui airport to the south-southeast (Fig. 2) with discontinuous exposure over approximately 450 m of section provided an additional dataset from which to examine the shallow-marine to fluvial transition (Fig. 4). Stratigraphic thicknesses along World Bank Road were measured using a Jacob's staff and Brunton compass, corrected for apparent thickness based on bed attitudes and angle between road orientation and strike of bedding. Field mapping was complemented with a composite image of the World Bank Road section from 102 digital photos stitched together using Gigapan image-processing software (www.gigapan.com/gigapans/206848) to highlight the stratigraphic architecture of the deposits (Fig. 3). Exposures along the airport-road section were logged in detail, and stratigraphic thicknesses were estimated using the elevation of mapped contacts based on subhorizontal structural dips observed throughout the valley (Fig. 4).
We collected four samples for detrital-apatite fission-track (dAFT) analysis and detrital-apatite (U-Th)/He (dAHe) dating (Fig. 2, Table 1) for comparison with previously published detrital-zircon fission-track (DZFT) ages on the Tipam Group (Betka et al. 2018a). Samples 16AIZ1 and 16AIZ2, near Aizawl, sample the lower and middle Surma Group, respectively, and samples 16KZL5 and 16KZL4 were taken from the same stratigraphic intervals from the eastern part of Mizoram near the village of Khawzawl (Fig 2). dAFT analyses were conducted at the University of Arizona Fission Track Laboratory. Samples were irradiated at the Oregon State University Triga Reactor and analyzed with the external-detector method. Fission tracks were counted and Dpar values were measured using an Olympus BX51 microscope with automated stage system driven and calibrated using FTStage software. See supplemental files for a detailed description of dAFT methods. dAHe dating was performed at the University of Arizona Radiogenic Helium Dating Laboratory. Analyses followed the procedures outlined in Reiners et al. (2004) and Reiners (2005). See supplemental files for a detailed description of dAHe methods.
dAFT and dAHe single-grain ages were analyzed to determine maximum depositional age (MDA) and/or cooling age after burial resetting for the lower and middle Surma Group. The software RadialPlotter (Vermeesch 2009) was used for peak fitting and/or to determine minimum age using the model of Galbraith and Laslett (1993) from radial plots for dAFT and dAHe data. Maximum depositional ages of ~8 Ma for fluvial (Tipam Group) deposits in the study area were estimated from previously published detrital-zircon fission-track (dZFT) data (Betka et al. 2018a). HeFTy Time-Temperature Path Modeling dAFT and dAHe time-temperature (tT) modeling was conducted using the software HeFTy, Version 1.9.3 (Ketcham 2005). For apatite FT, the 5.5M etchant AFT annealing model of Ketcham et al. (2007) with length reduction in standard of 0.911 (personal value of SNT calculated using induced and spontaneous tracks in a Durango standard) was used to predict (forward model) fission-track data from randomly generated tT paths. For dAHe, the RDAAM model of Flowers et al. (2009) was used. Samples were modeled individually to determine the tT pathways that best match the measured data for each site. Each modeling run was constrained using the dAFT age and the six highest precision dAHe grains, as HeFTy allows a maximum of seven inputs per sample. Track lengths were measured for samples 16AIZ1 (n ¼ 69), 16AIZ2 (n ¼ 13), 16KZL4 (n ¼ 10), and 16KZL5 (n ¼ 7). Track lengths were included in tT modeling runs for sample 16AIZ1 but were omitted from tT modeling runs for samples 16AIZ2, 16KZL4, and 16 KZL5 because HeFTy was unable to find acceptable fits with track length data. We infer that this is because the track lengths reflect source-terrane thermal histories that are unique for each grain in our unreset or partially reset dAFT samples.
For each randomly generated tT path, the predicted dAFT and dAHe data were matched against measured ages, diameter values of fission-track etch pits (Dpar), and single-grain effective uranium (eU Flowers et al. 2009) values until 100 good fits were found. For sample 16AIZ1 with track length data, the model was run until 100 acceptable fits were found. The goodness-of-fit value is a merit function that combines tests for the fit of model-predicted dAFT age against the measured data using a Kuiper's test to indicate the probability of failing the null hypothesis that the model and data are different (Willett 1997;Ketcham et al. 2000). We used the default merit values of 0.5 for ''good'' fit and 0.05 for ''acceptable'' fit. When fitting multiple statistics HeFTy uses a mean merit value of 0.5 to define good fit, with an additional criterion imposed with the minimum merit value for any one statistic being 1/(Nþ1), where N is the number of statistics being evaluated (Ketcham et al. 2009).
Three tT boxes (and a present-day surface temperature constraint of 20 6 58 C) were used to constrain each randomly generated HeFTy tT path. A depositional age constraint box is based on the sample MDAs (determined below) and an assumed near surface temperature of 25 6 58 C. An earlier constraint box between the deposition age and maximum age of 50 Ma, and temperature range of 208 C to 2008 C, was used to account for a full variability in detrital-apatite grain ages incorporated into the sediment. A third constraint box was used between the depositional age and present day with temperature range of 20-1208 C to allow HeFTy to explore a full range of possible burial tT histories following deposition.
Inverse thermal modeling was also performed using the QTQt software package (Gallagher 2012) to verify HeFTy results and because QTQt allows a greater number of inputs per sample. With QTQt, the AFT data and all ten of the single-grain AHe results were used to constrain model runs for each site. Track length data were included for samples 16AIZ1 and 16AIZ2. Otherwise, model runs were set up identically to the HeFTy models, using the Ketchum et al. (2007) annealing model for AFT and the Flowers et al. (2009) RDAAM diffusion model for apatite He. QTQt models were also constrained based on the sample MDAs and an assumed near-surface temperature of 25 6 58 C. The current surface or subsurface temperature (final point in each tT path) was similarly set to 20 6 58 C. The QTQt Markov chain Monte Carlo (MCMC) algorithm sampling parameters and thermal-history path relative probabilities are presented in supplemental materials (Fig. S2).
We described the facies of the Tipam and Surma groups in terms of modern analog environments of the GBD (Fig. 2). We identified and mapped seven facies units, including four that represent shallow marine to intertidal deposition (M1-M4) and three that represent deposition in fluvial environments (F1-F3). We present a generalized overview of facies in the first subsection, with more detailed descriptions and associated departures from these generalized descriptions from the well-exposed sections along the World Bank Road (Fig. 3) and the airport road (Fig. 4) in the subsections that follow. 1) Analyses by external detector method using 0.5 for the 4p/2p geometry correction factor.
2) Ages calculated using dosimeter glass: IRMM540R with f 540R ¼ 368.1 6 14.9 (apatite).
3) Pv 2 is the probability of obtaining a v 2 value for v degrees of freedom where v ¼ no. of crystals -1.
Shallow-Marine Lithofacies in the Surma Group (M1-M4)
The M1 facies contains mm-cm-scale planar-laminated gray to tan silty shale, interpreted to represent deposition on the bottomset of the subaqueous clinoform (Figs. 2A,B,5A). Bottomset beds of the modern GBD subaqueous delta, found in water depths from ~80 m to the shelf edge, exhibit similar patterns of grain size and lamination, and there is little evidence of distorted bedding or slumping in this facies, suggesting deposition at some distance from the slope of the clinoform (Kuehl et al. 2005). Contextually, the M1 facies is often found grading upwards in the M2 facies (Fig. 4). We interpret this facies association as the progradation of the subaqueous delta across the paleo-shelf (Fig. 2A,B).
The M2 facies consists of cyclical beds of alternating fine sands and silts with muds. Coarser units exhibit occasional unidirectional-current indicators (Fig. 5B). Outcrop scale (15-20 cm, Fig. 5B) layering scales with annual (monsoon-intermonsoon) cycles as observed in the Pride Formation (Upper Mississippian) of southern West Virginia, and superimposed mm-scale laminae share characteristics with spring-neap cycles identified in Pride Formation rhythmites (Eriksson and Romans 2015). Alternating coarse-fine layers are consistent with cycles of varying sediment input, suggesting that the M2 facies was deposited in an environment similar to the modern foreset beds of subaqueous delta of the GBD (Kuehl et al. 2005). Preliminary thin-section analysis from samples at a quarry south of Aizawl indicate ~95% quartz composition. We classify these beds as rhythmites, and interpret them to represent the prograding foreset beds of the subaqueous delta (Figs. 2A,B,5B,S3).
The M3 facies contains cm-scale beds of light gray sands, silts, and muds with flaser, wavy, and lenticular bedding (Fig. 5C). The presence of isolated Thalassinoides and Skolithos ichnofossils indicate deposition on the shallow shelf and subtidal to intertidal region (Figs. 2A, B, 3, 5C), but the general lack of bioturbation suggests rapid deposition and reworking of sediments on the shelf. These characteristics are consistent with topset bed deposits of the modern GBD, which are frequently disrupted by cyclones and shallow currents (Kuehl et al. 2005). A general increase in sand content from the M2 to the M3 facies is consistent with observations from the modern GBD subaqueous delta, where sand:silt:mud ratios increase from 1:7:2 on the foreset beds to 4:3:3 on the topset beds (Kuehl et al. 2005) (Fig. 5B,C). Flaser, wavy, and lenticular bedding found in some units is interpreted to represent tidal influence in the shallowest deposits in the M3 facies.
The M4 facies consists of massive, thick-bedded (. 1 m) sandstone with mixed trough crossbeds and bidirectional current indicators (Figs. 2A,B,5D). The sandy cross-bedded tabular sandstone beds are generally limited to tens of meters in lateral extent, such that the M4 facies is not always found near the transition from shallow-marine to fluvial deposits. For example, a prominent M4 exposure is visible on the airport-road section directly underlying fluvial deposits of the Tipam Group (Fig. 4) but the World Bank Road section (~25 km to the south) does not contain any evidence of these tabular sand bodies (Fig. 3). We interpret the crossbedded sand units of the M4 to represent large tidal channels in an intertidal environment inboard of the shallow shelf deposits of the M3 facies.
The F1 facies contains tan to brown fine-to medium-grained sandstone beds with meter-scale thicknesses and 0.5-m-scale trough cross bedding (Figs. 2A,B,3,5E). The trough cross-bedded sandstones are often contained within small (meter scale) scour-and-fill structures (Fig. 5E). Current indicators are generally unidirectional, and the red and brown colors (suggesting conditions of oxidation) of this facies is in striking contrast to the blue-gray colors of the M1-M4 marine facies. The F1 often exhibits a gradational contact with the M3 facies below (Fig. 3). We interpret this unit to represent distributary channels of the lower delta plain (Fig. 2A,B).
The F2 facies consists of fine-medium tan sandstones with large-scale (up to several meters) trough and tabular cross beds and erosional scours (Fig. 5F). The scale of the internal features in the F2 facies scale with similar features found on cut banks on the modern Brahmaputra River (Fig. 5F). Well sorted fine to medium sands are also consistent with the deposits found in the Holocene-modern braidbelt of the Brahmaputra River in Bangladesh (Pickering et al. 2014;Sincavage et al. 2018). This facies shares characteristics with the Tipam Group and is interpreted as the main braidbelt of the late Miocene-Pliocene ancestral Brahmaputra River (Fig. 2A,B).
The F3 facies consists of thin-bedded (tens of centimeters) multi-colored silts and fine sands with small-scale (, 0.5 m) cross beds (Figs. 2A,B,5G). This facies is poorly exposed due to its susceptibility to rapid weathering, and as such has limited lateral extent in the study area. From the limited outcrops available, the scale of the fluvial features and grain sizes observable in this facies are comparable to smaller distributary channels and splays off the main trunk channels of the lower delta plain as observed near Kuakata Peninsula in southern Bangladesh (Allison et al. 2003).
The exposed stratigraphy along the World Bank Road near Sumsuih (Figs. 2,3) illustrates the progradation of the ancestral Brahmaputra River delta front in central Mizoram (Fig. 3; complete stratigraphic section in Fig. S1) expressed in sub-environments of the M3 facies. Decimeter-scale alternating beds of blue-gray silts and fine sands with generally conformable contacts dominate the base of the section. Fine units contain evidence of soft-sediment deformation, whereas coarse units contain bioturbation, long-wavelength (~30 cm) ripples, and channel-form geometries (~1/2 m thick). The lowest ~100-m-thick unit contains two coarsening-and thickening-upwards packages of ~25 and ~60 m thickness, respectively (Figs. 3,S1). The gradational transition from bioturbated silts and sands to thicker and coarser units exhibiting softsediment deformation and slumping is interpreted to represent progradation of the advancing prodelta.
The overlying ~100 m thick section is dominated by gray-tan very fine to fine sands with generally coarsening and thickening upwards to beds up to 1 m thick and occasional evidence of erosive basal contacts (Figs. 3,S1). The prevalence of wavy and lenticular bedding increases upsection, as well as low-angle cross bedding and ripples. Thin (10-15 cm) silty interbeds are found between the thicker and coarser sand beds near the top of this section. We interpret this section to represent sandy distributary-mouth bars (Figs. 3,S1).
The next section upwards is a relatively thick (~160 m) succession of alternating thin (10-15 cm) blue-gray silts and fine sands of the M3 facies, with an overall coarsening and thickening up-section (Figs. 3,S1). Poorly indurated basal interbeds exhibit wavy and flaser bedding interspersed with bioturbated beds. Thicker (2-3 m) tabular, structureless sands (Figs. 3,S1) punctuate this section. The thicker and coarser tabular beds near the top of the section exhibit low-angle cross bedding with some soft-sediment deformation near the top of the section. This unit represents interdistributary-bay deposits punctuated by splays (Figs. 3,S1).
The uppermost ~50 m of the section shows the transition from shallowmarine facies (M3) into the predominantly fluvial facies (F1) of the Tipam Group (Figs. 3,S1). Meter-scale gray-tan fine-medium sand beds dominate this section, with thin (mm to cm scale) mud laminae. Sand beds alternate from structureless to exhibiting faint low-angle cross bedding, and generally coarsen and thicken upwards with oxidation present in some of the upper sand layers. A 5-m-thick clay unit in the village of
) caps the top of the section, interpreted as distributary channels of the lower delta plain.
The airport-road section near Lengpui does not contain continuous exposures like the World Bank Road to the south, but there are enough high-quality outcrops to document the shallow-marine to fluvial transition (Figs. 2,4). The gentle (10-158) NNE dip and excavated topography from the river provide several vantage points to view the succession of facies (Fig. 4). Exposures of the lowest stratigraphic unit (M1 facies) crop out along the road from the hairpin turn south of the river crossing near Lengte Road (Fig. 4). The lowest sections of the M1 here contain thin-bedded (1-2 cm) light gray silty fissile shales, grading upward into 3-5 cm sandier units. Near the new road cut just north of the Lengte Road intersection (Fig. 4), sand bodies transition from tabular geometries to scoured bases, with increased thicknesses (up to 8-10 cm thick) and more periodic character associated with the M2 facies. This section of the road illustrates the transition from bottomset (M1) to foreset (M2) beds of the prograding ancestral delta (Fig. 2A,B).
Farther up the road towards the sharp west bend (Fig. 4), grain size increases to silt and fine sand in decimeter-scale sand packages of the M3 facies. Current indicators begin to shift from unidirectional to more bidirectional, punctuated by 0.5-1 m tabular sand bodies starting near the bend in the road. Exposure quality is variable here, but with careful observation the transition from foreset beds (M2) to topsets (M3) and finally intertidal environments (M4) is evident (Fig. 2A,B).
The best exposures along the road are at a prominent cliff face on the north side of an east-west drainage between two hairpin turns in the road 1 km south of the airport (Figs. 4,5D). At road level, blue-gray cm-scale silts and fine sands with periodic character and bidirectional current indicators (interpreted as M4 facies) are overlain by a large (5-10 m thick) sand body with large-scale bidirectional current indicators (F1 facies). A dirt path provides access to exposures of the F1 facies at the top of the cliff, where tan-orange 1-m-thick channel bodies are located. This location provides excellent exposures of an abrupt transition from intertidal (M4) deposition to the delta-plain distributary channels of the Tipam Group (F1). As the road continues to meander towards the Lengpui airport, exposure quality decreases, but occasional outcrops show that the road traverses across the M4-F1 facies contact. The village adjacent to the airport is predominantly underlain by the F1 facies of the Tipam Group, exhibiting low-angle cross bedding. While quality exposures for a detailed paleocurrent analysis are limited in the region, the data available here indicate a predominantly NW direction of transport.
The dAFT and dAHe ages were used to constrain the maximum depositional age (MDA) and thermal history of the Surma Group. dAFT results from all of the samples show a distribution of single-grain ages between ~60 Ma (16KZL4) and ~10 Ma (Fig. 6). The wide (~40-50 Myr) distribution of single-grain ages in each sample indicates that none of the samples have experienced enough burial reheating to fully anneal fission tracks, constraining burial temperatures to , 1208 C and indicating that the age distributions can help to constrain MDAs. Samples 16AIZ2 and 16KZL4 from the middle Surma Group (upper Bhuban Fm.) give minimum peak ages 19.6 6 4.6 and 25.7 6 1.1 Ma, respectively. Samples 16AIZ1 and 16KZL5 from the lower Surma Group (lower Bhuban Fm.) give minimum peak ages 17.97 6 0.91 and 18.5 6 1.3 Ma, respectively (Fig. 6). The minimum peak ages from dAFT are younger from the lower Surma Group, limiting the deposition of the Surma Group to , 17 Ma, assuming an ~1 Ma AFT lag time for exhumation from the AFT closure temperature (~1208 C) to the surface) typical for the Himalayan source rocks of these sediments (Lang et al. 2016). This result is consistent with previously published biostratigraphy studies that place the Surma Group in the middle to late Miocene (Alam et al. 2003;Gani et al. 2003;Tiwari et al. 2013;Lalnuntluanga et al. 2014).
Detrital AHe results differ for samples collected near the cities of Aizawl (Fig. 2; 16AIZ1 and 16AIZ2) and Khawzawl (16KLZ4 and 16KLZ5). Near Aizawl, both samples 16AIZ1 and 16AIZ2 show single grain ages between ~3 and 9 Ma. Near Khawzawl, single-grain ages from samples 16KZL4 and 16KZL5 range from ~9 to 60 Ma (Figs. 7,8). In samples 16AIZ1 and 16AIZ2, effective uranium concentrations (eU ¼ [U ppm] þ 0.245*[Th ppm], Fig. 7) show a generally positive correlation with age which suggests partial resetting in the AHe partial retention zone following deposition (. 608 C;Reiners and Brandon 2006). Assuming a geothermal gradient from 16 to 238 C/km for the upper 3 km of the basin (Zahid and Uddin 2005) and ~208 C surface temperatures yields minimum burial depths from ~1.7 to 2.5 km for samples 16AIZ1 and 16AIZ2. The thickness of the eroded Surma Group intertidal (M4) and shelf (M1-3) facies overlying the samples has been estimated to be ~2 km based on geologic cross sections (Betka et al. 2018a(Betka et al. , 2018b). An additional 1.5 km of now eroded fluvial (F1-3) deposits may have also contributed to burial, consistent with our results. In contrast, the Khawzawl samples (16KZL4 and 16KZL5) yield uncorrelated eU vs. AHe age relationships consistent with unreset detrital AHe grain ages. Combined with the differing dAHe age distributions between Aizawl and Khawzawl samples, the eU vs. age patterns suggest burial resetting for the Aizawl samples near the center of the basin, and detrital ages (limited burial to , 608 C) for the Khawzawl samples near its eastern margin.
Near Khawzawl, we use the secondary peak at 9.9 6 0.1 Ma of middle Surma Group sample 16KZL4 to constrain the MDA (Fig. 7). This sample does have one young (4.95 6 0.21 Ma) AHe grain age, but it is a small grain (~44 lm radius) and contains a very low (14 ppm) eU. The increased susceptibility of low eU apatite to He diffusion could cause a young grain to experience partial resetting with minimal burial (Flowers et al. 2009). The minimum AHe grain age for the lower Surma near Khawzawl (16KZL5) is 11.9 6 0.15 Ma. Accounting for an ~1 Myr lag time in the rapidly eroding Himalayan source of these apatite grains during the Miocene before establishment of the Shillong Plateau rain shadow in Bhutan (e.g., Coutand et al. 2014;Lang et al. 2016) gives maximum depositional ages of ~9 and ~11 Ma for the middle and lower Surma, respectively. These results further support a mid-late Miocene maximum age for the Surma Group consistent with the dAFT data in this study, but slightly younger than the conventionally reported middle Miocene age for the lower Surma Group (Gani and Alam 2003;Alam et al. 2003).
Near Aizawl (16AIZ1 and 16AIZ2), the cooling ages after burial resetting were estimated using the minimum peak fitting routine in Radial Plotter (Galbraith and Laslett 1993). This assumes that the youngest grain ages are fully reset, and thus record the time of cooling onset. The youngest peak for sample 16AIZ1 is 3.14 6 0.13 Ma, and for sample 16KZL5 it is 3.36 6 0.14 Ma (Figs. 7,8). The resetting of 16AIZ1 and 16AIZ2 ages near Aizawl indicates that they were buried by a thicker (. 2 km) stratigraphic section near the center of the basin than equivalent unreset deposits located near Khawzawl (16KZL4 and 16KZL5,Figs. 7,8).
The postdepositional thermal history of each sample was modeled using both HeFTy and QTQt modeling programs (Ketchum 2005;Ketchum et al. 2000Ketchum et al. , 2009;;Gallagher 2012). Figure 9A shows the temperature-time history modeled with QTQt, and Figure 9B shows the results from HeFTy. By comparison, both QTQt and HeFTy predict similar tT histories for each sample. For Aizawl sample 16AIZ1, the best-fitting tT paths indicate reheating beginning ~11 with maximum reheating occurring ~8.5 Ma (QTQt) or 7 Ma (HeFTy), followed by steady-state cooling to the present day (Fig. 9A,B). For 16AIZ2, both HeFTy and QTQt indicate reheating beginning ~9 Ma followed by steady state cooling from ~6.5 Ma to present for the middle Surma at this location (Fig. 9A,B). QTQt results suggest maximum reheating to ~658 C for the middle Surma (16AIZ2) and ~858 C for the lower Surma (16AIZ1). Similarly, the HeFTy best-fit result indicates maximum reheating at ~6 Ma to ~758 C for the middle Surma (16AIZ2), and ~858 C at 8 Ma for the lower Surma (16AIZ1). The best-fit path determined by HeFTy for the middle Surma (16AIZ2, black curve) is just in the margin of error of the QTQt model; however, the weighted mean of good and acceptable fitting paths from HeFTy (yellow curve, Fig. 9B) suggests reheating to ~858 C, ~108 warmer than the QTQt model, and all the goodfit paths indicate rather more variability in the exact timing of maximum reheating. For the lower Surma (16AIZ1), both the HeFTy best-fit (black curve) and weighted-mean (yellow curve) paths predict maximum reheating of ~858 C, which is within the margin of error of the QTQt model but on the warmer side of the mean path (blue curve). QTQt model results for both 16AIZ1 and 16AIZ2 predict maximum reheating 0.5-1 Myr earlier than the respective HeFTy results. This more rapid burial curve reflects track-length data included in the QTQt model runs but not in HeFTy; thus the QTQt reheating path is better constrained than HeFTy.
To estimate burial depths from the modeled best-fitting reheating curves, we consider modern geothermal gradients reported from gas wells in the outer part of the IBR that range from 16-238 C/km at 1-3 km depths (Zahid and Uddin 2005). Taking the mean and standard deviation of these presentday gradients for the upper 3 km of the basin yields 20.3 6 2.68 C/km. Assuming burial reheating along this average geothermal gradient indicates burial to 2:20 þ0:33 À0:25 km depth for the upper Surma Group (16AIZ2). Calculation of rates is more uncertain given the variability in time of maximum reheating of the HeFTy good fit paths, but ~2.5 Myr of reheating from ~9-6.5 Ma implies a burial (accumulation) rate of 0:90 þ0:13 À0:10 mm/yr for the middle Surma (16AIZ2, Fig. 9). For the lower Surma (16AIZ1), burial depths reached 3:20 þ0:48 À0:37 km, and an approximate burial (accumu-lation) rate of 1:10 þ0:16 À0:12 mm/yr between ~11.5 and 8.5 Ma. Note that a cooler paleo-geothermal gradient would increase the depth range, for example, a 158 C/km gradient yields depths of ~3 and 4.3 km for samples 16AIZ2 and 16AIZ1 from the middle and lower Surma Group, respectively.
In contrast to the Aizawl samples, the single-grain AHe and AFT results from the Khawzawl samples (16KZL4 and 16KZL5) indicate that, with the exception of one young 5 Ma AHe grain age, they did not experience enough burial reheating to reset AHe and AFT thermochronometers (previous section, Figs. 7,8). The lack of postdepositional reheating is confirmed in the best-fitting QTQt and HeFTy thermal histories for these samples. With the applied depositional constraints, results indicate that samples 16KZL4 and 16KZL5 were not heated to postdepositional temperatures within the AHe partial-retention zone (. 608 C; Reiners and Brandon 2006; best-fit model paths show reheating between ~40 and 508 C). In 16KZL4, inclusion of the one reset 5 Ma AHe grain age in the HeFTy modeling results in both the best-fit (black curve) and weightedmean (yellow curve) model paths requiring minor reheating to ~558 C from 9-2 or 9-4 Ma, respectively (Fig. 9B). The QTQt model for 16KZL4 (including the same reset 5 Ma AHe grain age) also favors a tT path with slight reheating to ~60 þ25 À30 8 C between ~10-8 Ma. However, the lack of reset grains means that the exact shape of the postdepositional best-fit tT path in both HeFTy and QTQt is very poorly constrained by the data; thus we did not estimate burial rates for sample 16KZL4. Sample 16KZL5 does not have any reset AFT or AHe grains. Accordingly, the QTQt model results indicate that it resided at or near the surface temperatures since ~11 Ma (Fig. 9A). HeFTy did not find any acceptable fits for sample 16KLZ5.
Analogue depositional environments for all of the seven new facies mapped in Miocene-Pliocene outcrops of the IBR are found on the modern GBD (Fig. 2). The scale, lithologic character, stratigraphic architecture, and facies associations of modern and ancient deposits are strikingly similar, particularly when comparing the F2 facies with the modern Brahmaputra braidbelt (Fig. 5F). Large (meter scale), high-angle sweeping cross beds in fine to medium sands are characteristic of both the F2 facies as well as Holocene sediments of the Brahmaputra mainstem (Pickering et al. 2014;Sincavage et al. 2018). Deposits of this nature have also been described at the western edge of the Shillong Massif, interpreted as remnants of early post-avulsion Brahmaputra River braidbelt deposits (Pickering 2016).
The M2 facies contains stacked sets of alternating light and dark bands that are ~15 cm thick. Michels et al. (1998) document annual sedimentation rates on the GBD of ~10 cm/yr, comparable in scale to the observed banding in the M2 facies (Fig. 5B). Annual cyclical sedimentation rates of ~10 cm/yr were also noted in analogous central Appalachian Basin in the eastern United States in the Pride Formation, interpreted as tidally influenced distal prodelta deposits (Miller and Erikkson 1997;Erikkson and Romans 2015). While these rates are substantially higher than those estimated from our thermal modeling (which does not account for compaction), sediment accumulation rates are known to be highly variable across the modern subaqueous delta, from tens of cm/yr along the delta foresets to less than 1 cm/yr on the topsets and bottomsets (Kuehl et al. 2005). The sedimentation rates are also likely to be highly cyclic with sea-level changes, while the thermal modeling will reflect the long-term net sediment accumulation rates. Thin-section images of the M2 facies (Figs. 5B,S3) reveal finer (mmscale) cyclicity, a possible indicator of spring-neap tidal cycles and even dominant vs. subordinate currents in daily cycles (Erikkson and Romans 2015). Of five thin-section images, three show clear evidence of cyclicity, one shows faint evidence, and one does not contain evidence of finerscale cyclicity (Fig. S3). Combined, these observations from analogous deltas coupled with our documented facies associations provide reasonable justification for our interpretation of the M1 and M2 facies representing the prograding subaqueous clinoform of the ancestral Brahmaputra delta.
Exposures of the two ~500 m sections on World Bank Road and the Lengpui Airport road show evidence of progradation of a large delta system, from M1 and M2 subaqueous-delta facies, to M3 and M4 intertidal facies, capped by F1 fluvial facies (Figs. 3,4). Prior to progradation of the paleodelta associated with Surma and Tipam sediments, distal deposits of the Eocene-Miocene Barail system were deposited in the region. The Barail stratigraphy records the progradation of a delta smaller than that represented by deposits in the Surma and Tipam groups that was focused farther north in Assam (Alam et al. 2003;Biswas and Mukhopadhyay 2011;Bezbaruah and Muzamil 2013). Based on the burial depths determined from the thermal modeling, we deduce that the late Miocene-Pliocene sedimentary pile (Surma and Tipam groups) thinned toward the east, off-lapping the Oligo-Miocene Barail Group as the ancestral Brahmaputra Delta prograded southward (Fig. 10). The MDA constraints indicate that the lower to middle Surma Group was deposited from ~11-9 Ma, consistent with a previously published MDA of ~8 Ma for the overlying Tipam Group (Betka et al. 2018a). In the center of the basin (near Aizawl) the middle to lower Surma Group was buried to depths of ~2.2-3.2 km (using the mean depths predicted from the thermal models, and assuming the modern GBD near-surface geothermal gradient) at a rate of ~1 mm/yr, without accounting for compaction. In contrast, at the eastern margin of the basin (near Khawzawl), the middle to lower Surma Group was buried to a depth of only ~1.7 km based on thermal modeling of 16KZL4 with a single reset AHe grain (Figs. 7,8) at a lower sedimentation rate of ~0.3 mm/yr, without compaction.
Compaction rates in deltaic settings are spatially variable across delta subenvironments and notoriously difficult to measure (Meckel et al. 2007), but sedimentation and subsidence rates of millimeters to tens of millimeters per year are not uncommon on large modern deltas, including the GBD (Törnqvist et al. 2008;Rogers et al. 2013;Reitz et al. 2015, Jankowski et al. 2017;Rogers and Overeem 2017;Grall et al. 2018;Zoccarato et al. 2018;Becker et al. 2020). Fine-grained sediments of the Surma Group could be expected to experience 40-80% compaction (Allen 2000), such that actual sediment accumulation rates may have ranged from ~1.5-5.3 mm/yr near the center of the paleodelta (Aizawl) and ~0.6-1.7 mm/yr along its eastern margin (Khawzawl).
Post-burial cooling likely reflects tectonic uplift and exhumation of the Surma Group with the advance of the IBR fold-thrust belt (cf. Betka et al. 2018a). Both 16AIZ1 and 16AIZ2 from Aizawl show reset cooling ages of ~3 Ma (Figs. 7,8). Thermal-modeling results suggest steady-state cooling from ~7-8.5 Ma for 16AIZ1 and ~6.5 Ma for 16AIZ2. These results indicate initiation of uplift associated with the advancement of the IBR fold-and-thrust belt sometime between ~8-6.5 Ma (Fig. 10), in good agreement with the findings of Betka et al. (2018a) which limit the maximum age of deformation to the ~8 Ma MDA for the Tipam Group. Assuming that ~2.2 km of sediment were removed over ~6.5 Ma of uplift yields an average exhumation rate of 0.34 mm/yr. The deformation also resulted in uplift and exhumation of the Barail Group to the east of Khawzal near the Myanmar border (Figs. 2,10).
While the kinematic and mechanical evolution of the IBR fold-thrust belt are beyond the scope of this paper, we speculate that this latest Miocene to Pliocene exhumation event may record a kinematic response of the internal part of the thrust wedge to the newly applied load of the Shillong massif and its nascent Sylhet Basin. We infer that the ~8-6.5 Ma onset of cooling and exhumation of the Surma Group in the central part of the basin marks both the westward propagation of the thrust wedge and the end of deposition of ancestral Brahmaputra deposits with the avulsion of the river to the north and west of uplifting Shillong Massif (Najman et al. 2016;Govin et al. 2018a). We also infer that post ~8 Ma rapid (! 15 mm/yr) westward propagation of the thrust front (Betka et al. 2018a) reflects a critical wedge response to the westward-shifting depocenter. For example, both analog (Wu and McClay 2011) and numerical (Fillon et al. 2013) models of thrust wedges indicate that high sedimentation rates near the foreland promote both out-of-sequence deformation in the internal part of the wedge and favor rapid foreland propagation of the frontal thrust.
Our results, combined with observations from IODP 354, distinguish eastern India as an important location not only for the routing of Miocene-Pliocene sediments to the BNF, but also in initiating a major avulsion of the ancestral Brahmaputra River. Mid to late Miocene (~13.5-6.8 Ma) intensification of turbidite deposition to the eastern BNF (France-Lanord et al. 2016b) brackets our determined MDAs for shallow-marine deposition of the Surma Group (~9-11 Ma). The progradation of shallow marine (M1-M4) and fluvial (F1-F3) facies to produce sediment burial depths deep enough for grain resetting of our westernmost samples, combined with high rates (10-20 cm/ky) of turbidite deposition at IODP Site U1451, suggest a depocenter for the ancestral delta near Aizawl, India (Fig. 10). Focused sedimentation here led to increased burial depths (2.2-3.2 km) and resetting of AHe ages, whereas sediments on the eastern margin were not buried as deeply and thus were not reset (Figs. 9, 10).
Thermal-modeling results indicate that cooling and the onset of exhumation of the Surma Group at ~8-6.5 Ma was coincident with a decrease in deposition rate (to ~1-2 cm/ky) and initiation of increased carbonate deposition in the eastern BNF (France-Lanord et al. 2016a, 2016b). This likely represents a forced westward shift of the depocenter initiated by the advancing deformation front of the IBR FTB (Fig. 10), and eventual avulsion of the river around the rising Shillong Massif, leading to increased turbidite deposition on the western BNF through the Pliocene and Pleistocene (France-Lanord et al. 2016a;Najman et al. 2016;Govin et al. 2018a). Therefore, allogenic (i.e., tectonic) processes forced the Brahmaputra River into closer proximity to the Ganges River, which set the stage for additional allogenic (i.e., climatic and/or eustatic) and autogenic constraints on source mixing and sediment dispersal patterns throughout the Plio-Pleistocene (Blum et al. 2018). Our results connect marginal marine to shelf deposits of the Neogene paleo-Brahmaputra Delta to the BNF, underscoring the importance of deep-sea fans as receptors of tectonic and climate signals from the source area (Hessler and Fildani 2019). By constraining the timing of initiation of the ancestral Brahmaputra Delta, we provide pathways to link changes in climate forcing signals to delta growth and terrestrial carbon sequestration and establish the GBD as a potentially critical climate-change modulator throughout the late Cenozoic (Galy and Eglinton 2011;Hein et al. 2017Hein et al. , 2020)).
Outcrops of the IBR record a rich history of sediment dispersal, burial, and exhumation associated with the progradation of the Miocene ancestral Brahmaputra Delta and the advancing deformation front of the IBR fold belt. Detailed stratigraphic mapping reveals a progression of facies representing a shift in depositional environment from the subaqueous delta front, to intertidal, and finally fluvial deposition on the paleodelta plain. dAFT and dAHe analyses constrain the MDA of the shallow-marine facies of the Surma Group to ~9-11 Ma, consistent with previously published depositional ages for the fluvial Tipam Formation at ~8 Ma, and bracketed by periods (~13.5-6.8 Ma) of increased turbidite deposition rates (10-20 cm/ky) on the eastern BNF. Thermal modeling suggests deeper (2.2-3.2 km) burial of shallow-marine sediments near the center of the basin near Aizawl than on its eastern flank (near Khawzawl). Thermal histories further indicate initiation of exhumation at ~8-6.5 Ma, marking the end of deposition in the basin and the onset of both the avulsion of the Brahmaputra River around Shillong Massif and continued westward propagation of the thrust front. Timing of this deformation coincides with a decrease in turbidite deposition rates (to ~1-2 cm/ky) and increased carbonate deposition on the eastern BNF. Our results constrain the timing and rates of tectonic processes that set the stage for climatic and/or eustatic and autogenic forcing of sediment mixing and dispersal to the BNF through the Plio-Pleistocene and establish the GBD as a potentially critical modulator of late Cenozoic climate change.
Supplemental materials are available from the SEPM Data Archive: https:// www.sepm.org/supplemental-materials.