Indian winter and summer monsoon strength over the 4.2 ka BP event in foraminifer isotope records from the Indus River delta in the Arabian Sea

2.1 Monsoon – land-based processes

Today, most of the annual precipitation over northwest South Asia stems from the ISM, and occurs mainly between June and September. The pressure gradient between the low-pressure Tibetan Plateau and high-pressure Indian Ocean is accompanied by the ITCZ (Intertropical Convergence Zone) reaching its northward maximum in summer, which draws in moisture over the subcontinent via southwesterly winds from the Indian Ocean (Gadgil, 2003). The summer rainfall gradient increases from the central Thar Desert (as little as 100 mm direct summer rainfall per year) to the Himalaya mountains in the north (>1000 mm) and the Aravalli range to the west (>500 mm) (Fig. 1b).

The IWM rain falls between December through March and is mainly the result of atmospheric Western Disturbances (Dimri and Dash, 2012; Yadav et al., 2012) originating over the Mediterranean and Black Sea (Hatwar et al., 2005) that allow for moisture incursion from the Arabian Sea (Rangachary and Bandyopadhyay, 1987). During the IWM, the pressure gradient is reversed from the summer condition, allowing the passage of Western Disturbances when the ITCZ moves southward. As winter transitions to spring, predominantly northeasterly winds shift to westerly winds (Sirocko, 1991) that result in peak winter rainfall over the plains of northwest India in February and March. Anomalously cool, evaporative conditions over the northern Arabian Sea (promoting deeper winter mixing) also correlate with increased winter precipitation in the western Himalayas (Dimri, 2006). The winter rainfall gradient increases from the southern Thar Desert (<10 mm yr⁻¹) up to the Himalayas in the northwest (>400 mm) (Fig. 1c). Overall, the IWM contributes between roughly 10 % and 50 % of the total annual rainfall of northwest South Asia today.

The Indus and the other rivers that make up Punjab are partly fed by winter snowmelt and ice melt from their upper mountain catchment areas. Melting peaks during the summer months around July–August (Yu et al., 2013), which coincides with the peak of ISM rainfall, and Indus River discharge reaches its maximum during August (Karim and Veizer, 2002). The proportion of winter to summer precipitation contributing to the Indus River is not entirely clear, although one study has estimated a 64 %–72 % contribution of winter precipitation from the deuterium excess of Indus River water (Karim and Veizer, 2002), whereas a previous study estimated a lower 15 %–44 % contribution of snowmelt to Indus tributaries (Ramasastri, 1999). Since the 1960s, the Indus River has seen more than a 50 % reduction in discharge because of the construction of barrages as well as the diversion of water for agricultural uses (Ahmad et al., 2001).

2.2 Hydrography – core site and ocean-based processes

Core 63KA was obtained by the PAKOMIN cruise in 1993 (von Rad et al., 1995). The laminated core from the northeastern Arabian Sea (24^∘37^′ N, 65^∘59^′ E) was taken at 316 m water depth on the continental shelf, ∼100 km west of the Indus River delta. The core has high sedimentation rates (equivalent to a temporal resolution of around 18 yr cm⁻¹ in the period of interest, 5.4–3.0 ka), and all foraminifer proxies were produced from the same laminated core with no bioturbation. An important aspect of core 63KA is that different components of the monsoon system are co-registered in the same sediment core, thereby permitting an explicit evaluation of the relative timing of different parts of the climate system (e.g., ISM and IWM).

Modern hydrographic conditions in the northeastern Arabian Sea are highly influenced by the seasonal monsoon. During summertime, highest sea surface temperatures (SSTs) are observed along with a shallow mixed layer depth <25 m (Schulz et al., 2002) (Fig. 2a). A low-salinity plume surrounds the Indus River delta and shoreline extending as far as the coring location (Fig. S1 in the Supplement). The reverse occurs in winter when the lowest SSTs are accompanied by surface water mixing to >125 m, resulting in warming of the deeper waters (Schulz et al., 2002). Northeasterly winds promote convection in the northeastern Arabian Sea by cooling and evaporation of surface water (Banse, 1984; Madhupratap et al., 1996), and during the transition from winter to spring, wind directions shift from northeasterly to westerly (Sirocko, 1991).

Figure 2(a) Seasonal surface water mixing depth based on station EPT-2 located nearby the coring site of 63KA (adapted from Schulz et al., 2002, who also used data from Hastenrath and Lamb, 1979). (b) Foraminifer depth ranges based on CTD profile. (c) Foraminifer abundances from EAST traps (overlapping peaks indicate data from multiple traps): G. ruber (orange), G. sacculifer (green), and N. dutertrei (blue) (adapted from Curry et al., 1992, using Zaric, 2005). (d) World Ocean Atlas (WOA) mean (1955–2012) temperature (red) and salinity (yellow) profiles at 24.875^∘ N, 65.875^∘ E, shown for summer (JAS) and winter (JFM) seasons (Locarnini et al., 2013; Zweng et al., 2013).

Download

The northern Arabian Sea is dominated by highly saline (up to 37 psu) surface waters known as Arabian Sea High-Salinity Water (ASHSW), which extend from the surface down to 100 m depth (Joseph and Freeland, 2005). The high salinity is explained by the high evaporative rates over this region. ASHSW forms in the winter but is prevented from reaching our coring site on the shelf by northerly subsurface currents until the summer (Kumar and Prasad, 1999). Along coastal areas, the ASHSW is starkly contrasted by the fresh water discharge of the Indus River, combined with direct precipitation. In contrast, surface waters in the Bay of Bengal on the eastern side of India have much lower surface water salinity because of overall higher precipitation and stronger stratification from weaker winds (Shenoi et al., 2002). The heightened evaporative conditions and highly saline surface waters of the northeastern Arabian Sea make it a sensitive study location to observe changes in discharge of the entire Indus River catchment area – ultimately tracking changes in monsoon strength. Unlike individual terrestrial records, which may be affected by local climatic processes, the marine record from core 63KA is more likely to integrate regional changes in the large-scale ocean–atmosphere system.

Planktonic foraminifera complete their life cycle within a few weeks (Bé and Hutson, 1977). Peak abundances indicate the time of year when each species tends to calcify, thereby recording the δ¹⁸O and temperature of the seawater in their CaCO₃ shells primarily during certain seasons. Foraminifer abundances in the eastern Arabian Sea have been studied by Curry et al. (1992) using sediment traps deployed at shallow (∼1400 m) and deep (∼2800 m) water depths (“T” in Fig. 1a). G. ruber and G. sacculifer have peak abundances during the summer months (June–September), whereas N. dutertrei lives mainly during the winter and has a secondary peak in summer months (Fig. 2c). Preferred depth ranges for each species reflect their ecological niches, including requirements for nutrients and tolerance for ranges of temperature and salinity (Bé and Hutson, 1977; Hemleben et al., 2012). G. ruber lives in the upper surface waters (0–10 m), G. sacculifer is found in slightly deeper surface waters (10–40 m), and N. dutertrei inhabits the base of the mixed layer near the thermocline (40–140 m) (estimates based on ranges from Farmer et al., 2007, and the local CTD – conductivity, temperature, and depth – profiles) (Fig. 2b).

3.1 Age model

The AMS (accelerator mass spectrometer) radiocarbon dates from Staubwasser et al. (2002, 2003) were obtained from 80 samples of mainly the foraminifer G. sacculifer and three samples of O. universa. In the interval of interest (5.4–3.0 ka), there are 15 radiocarbon dates with a 95 % confidence range of 30–130 years. The average sample resolution is 18 yr cm⁻¹. Bayesian age modeling software, BACON v2.3.3 (Blaauw and Christen, 2011), was used as an R package to update the age model of core 63KA. No major difference exists between the old and new age models, except for the period 13–11 ka (Fig. S5, Table S2). IntCal13 was used for radiocarbon calibration (Reimer et al., 2013) with marine reservoir ages provided by Staubwasser et al. (2002, 2003).

3.2 Stable isotope analysis

Oxygen and carbon isotopes were measured on three species of foraminifera selected from washed samples at 1 cm intervals throughout 132 cm of the core covering 5.4–3.0 ka: G. ruber (white, sensu stricto), G. sacculifer, and N. dutertrei. For G. ruber, 12±8 foraminifera were picked from the 400–500 µm size fraction with an average weight of 21.4±2.5 µg. The 400–500 µm size fraction was picked because too few specimens remained in the size fraction 315–400 µm used by Staubwasser et al. (2003). For G. sacculifer, 34±7 foraminifera were picked from the 315–400 µm size fraction with an average weight of 21.9±2.6 µg. For N. dutertrei, 34±4 foraminifera were picked from the 315–400 µm size fraction with an average weight of 25.9±2.2 µg. At some depth levels in the core there were insufficient foraminifera for measurement, along with outlier measurements in two cases, leaving 14 gaps in the G. ruber 400–500 µm record, 4 gaps in the G. sacculifer record, and no gaps for N. dutertrei. The published G. ruber is from the 315–400 µm size fraction and contains 17 gaps in the depth range examined (Staubwasser et al., 2003).

All foraminifera were weighed, crushed, and dried at 50 ^∘C. Samples were cleaned for 30 min with 3 % H₂O₂, followed by a few drops of acetone, ultrasonication, and drying overnight. Where sample weights exceeded 80 µg, oxygen and carbon isotopes were measured using a Micromass Multicarb Sample Preparation System attached to a VG SIRA mass spectrometer. In cases of smaller sample sizes, the Thermo Scientific Kiel device attached to a Thermo Scientific MAT253 mass spectrometer was used in dual-inlet mode. This method adds 100 % H₃PO₄ to the CaCO₃, water is removed cryogenically, and the dry CO₂ is analyzed isotopically by comparison with a laboratory reference gas. For both measurement methods, 10 reference carbonates and 2 control samples were included with every 30 samples. Results are reported relative to VPDB (Vienna Pee Dee Belemnite), and long-term reproducibility of laboratory standards (e.g., Carrara marble) is better than ±0.08 ‰ for δ¹⁸O and ±0.06 ‰ for δ¹³C. Reproducibility of foraminiferal measurements was estimated by five triplicate (three separately picked) measurements of G. ruber (400–500 µm) that yielded 1 standard deviation of ±0.12 ‰ (δ¹⁸O) and ±0.10 ‰ (δ¹³C). For G. sacculifer (315–400 µm) the standard deviation of eight triplicate measurements was ±0.07 ‰ (δ¹⁸O) and ±0.07 ‰ (δ¹³C), and for N. dutertrei (315–400 µm) the standard deviation of nine triplicate measurements was ±0.06 ‰ (δ¹⁸O) and ±0.07 ‰ (δ¹³C).

To calculate equilibrium values of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{calcite}\left(\mathrm{PDB}\right)}$, we used the CTD profile from station 11 (24.62^∘ N, 66.07^∘ E) taken in September 1993 during PAKOMIN Sonne cruise no. 90 (von Rad, 2013), which is nearly identical to the location of core 63KA (24.62^∘ N 65.98^∘ E). The $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{water}\left(\mathrm{SMOW}\right)}$ was calculated from salinity following Dahl and Oppo (2006), and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{calcite}\left(\mathrm{SMOW}\right)}$ was further calculated using the calcite–water equation of Kim and O'Neil (1997). We also used the equation of Shackleton (1974) as a comparative method for calculating $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{calcite}\left(\mathrm{PDB}\right)}$.

3.3 Statistical treatment

Statistical tests were applied to the raw data from the δ¹⁸O and δ¹³C time series, including the package SiZer (Chaudhuri and Marron, 1999; Sonderegger et al., 2009) in R software (2016), which calculates whether the derivative of a time series exhibits significant changes given a range of time spans. A Pearson's correlation test (confidence level 95 %) was done on paired samples from both size fractions of G. ruber. We also conducted Welch's t test to determine if the mean population of δ¹⁸O is significantly different before and after 4.1 ka.

Figure 3Core 63KA δ¹⁸O G. ruber from two size fractions shown in the context of the original record and also zoomed in over 5.4–3.0 ka: 400–500 µm (red) (this study), 315–400 µm (orange) (Staubwasser et al., 2003). δ¹⁸O of G. sacculifer 315–400 µm (green), and δ¹⁸O and δ¹³C of N. dutertrei 315–400 µm (blue) are shown over the interval 5.4–3.0 ka. Data are shown with a 210-year loess smoothing, and modern surface values ±1σ are plotted for comparison. Mean values for all species are denoted by the dotted line, and the pre- and post-4.1 ka mean values are indicated by an additional dotted line for N. dutertrei. Individual AMS radiocarbon dates are denoted by triangles near the timeline.

Download

As in the original data of Staubwasser et al. (2003), the oxygen isotope results show great variability, and distinguishing long-term trends in these data benefits from smoothing for visualization purposes. After completing all statistical tests and performing the differences on the raw data (132 depths), a loess (locally weighted) smoothing function was applied to the δ¹⁸O and δ¹³C data from 5.4 to 3.0 ka, using a 210-year moving window as described by Staubwasser et al. (2003). Loess smoothing uses weighted least squares, which places more importance on the data points closest to the center of the smoothing interval. The bandwidth of 210 years was considered a reasonable time window for capturing the overall trends in the dataset (other time windows are shown for comparison in Fig. S2).

Figure 4SiZer first derivative analysis (Chaudhuri and Marron, 1999; Sonderegger et al., 2009) applied to δ¹⁸O of (a) G. ruber 400–500 µm, (b) G. ruber 315–400 µm, (c) G. sacculifer 315–400 µm, and (d) N. dutertrei 315–400 µm. The red areas indicate statistically significant increases in δ¹⁸O, the blue represent decreases, and the purple no significant change. Black horizontal lines are the smoothing bandwidths (h=50, 80, and 200 years). The distance between the white lines denotes the change in smoothing bandwidth scaled to the x axis.

Download

5.1 Interpretation of foraminifer δ¹⁸O

The trends in the original δ¹⁸O record of G. ruber (315–400 µm) by Staubwasser et al. (2003) are reflected by our independent δ¹⁸O measurements of G. ruber in a larger size fraction (400–500 µm), although an important difference exists, suggesting a decrease in freshwater discharge as early as 4.8 ka. The larger size fraction is offset by approximately −0.2 ‰, which is similar to the size-related fractionation of −0.3 ‰ per +100 µm for G. ruber reported by Cayre and Bassinot (1998), and this could be attributed to size-related vital effects. Alternatively, part of the offset might be explained by interlaboratory calibration considering the data were produced using two different methods and mass spectrometers.

The observed 4.1 ka maximum in δ¹⁸O of G. ruber, living near the surface during summer months, could be attributed to either decreased SST or increased surface water salinity (Bemis et al., 1998). Staubwasser et al. (2003) acknowledged that a decrease in SST could cause the increase in δ¹⁸O in the G. ruber record but argued that this explanation is unlikely because a G. ruber δ¹⁸O record from core M5-422 in the northwestern Arabian Sea shows opposing trends over the same time period (Cullen et al., 2000), and a local alkenone SST proxy record shows relatively higher temperatures in the same period (Doose-Rolinski et al., 2001). If the ∼0.2 ‰ (relative to mean) increase in δ¹⁸O of G. ruber at 4.1 ka was caused by temperature change rather than salinity, a ∼1 ^∘C cooling of surface water would be required (Kim and O'Neil, 1997).

Following Staubwasser et al. (2003), we interpret the δ¹⁸O variations of G. ruber to be predominantly a salinity signal. Salinity at the core site is dependent on changes in Indus River discharge, local runoff, and direct precipitation. Although the ISM would be the main influence on direct precipitation and runoff at the coring location, changes in the IWM could also influence Indus River discharge.

The thermocline-dwelling foraminifera N. dutertrei show maximum abundances during winter and are interpreted to reflect winter mixing. During weak IWM conditions, colder unmixed water would result in higher δ¹⁸O values of N. dutertrei, whereas enhanced mixing and homogenization of the water column under strong IWM conditions would decrease δ¹⁸O. The minimum of δ¹⁸O in N. dutertrei occurs between 4.5 and 4.3 ka, pointing to a period of strengthened IWM. We interpret the stepped increase in δ¹⁸O of N. dutertrei at 4.1 ka to represent a decrease in IWM wind-driven mixing. Similarly, δ¹³C of N. dutertrei increases significantly after 4.1 ka (Fig. 3), which could indicate reduced upwelling of low δ¹³C intermediate water (Lynch-Stieglitz, 2006; Ravelo and Hillaire-Marcel, 2007; Sautter and Thunell, 1991); however, the interpretation of δ¹³C remains uncertain because of a poor understanding of the controls on the δ¹³C of planktonic foraminifera in this region. According to the δ¹⁸O signal of N. dutertrei, the temperature pattern in the thermocline implies surface cooling between 4.5 and 4.3 ka and surface warming after 4.1 ka interrupted only by a period of cooling between 3.7 and 3.3 ka, which is in broad agreement with records of alkenone sea surface temperature estimates from cores in the northeastern Arabian Sea (“E” in Fig. 1) (Doose-Rolinski et al., 2001; Staubwasser, 2012).

5.2 Interpretation of foraminifer Δδ¹⁸O

By using Δδ¹⁸O between foraminifer species, we can distinguish between additional processes affecting the surface waters and thermocline (Ravelo and Shackleton, 1995). This technique has been used previously to infer changes in the strength of the East Asian winter monsoon (EAWM) in the South China Sea (Tian et al., 2005), as well as mixed layer and thermocline depth in other studies (Billups et al., 1999; Cannariato and Ravelo, 1997; Norris, 1998). Here we use the difference in the δ¹⁸O of G. ruber and N. dutertrei ($\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$) to track changes in the surface-to-deep gradient. This gradient can be driven by either δ¹⁸O changes in the surface-dwelling (G. ruber) and/or the thermocline-dwelling species (N. dutertrei). During times of a strengthened winter monsoon, $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$ will decrease as surface waters are homogenized and the thermocline deepens. Similarly, $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$ will also decrease during times of a weakened summer monsoon, as decreased Indus River discharge will increase surface water salinity and δ¹⁸O of G. ruber will become more similar to N. dutertrei.

G. sacculifer is also a surface dweller but has a slightly deeper depth habitat than G. ruber. We thus expect G. ruber to be more influenced by surface salinity variations than G. sacculifer and suggest the δ¹⁸O difference between the two species ($\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$) reflects the influence of Indus River discharge on near-surface salinity. The smallest difference in $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$ occurs at 4.1 ka, which is interpreted as an increase in surface water salinity (Fig. 6).

The difference in δ¹⁸O between G. sacculifer and N. dutertrei ($\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-s}}$) also reflects surface mixing and thermocline depth, but G. sacculifer is less affected by surface salinity changes than G. ruber. Thus, the responses of $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$ and $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-s}}$ can be used to differentiate between surface water salinity changes and wind-driven mixing. Accordingly, simultaneously low $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-s}}$ and $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$ indicate a period of increased surface water mixing and increased IWM (such as the period between 4.5 and 4.3 ka), but times of relatively low $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-s}}$ but high $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$ and $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$ (around 5.0 ka) indicate periods of increased Indus discharge and strength of the ISM and IWM.

The following period of low $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-r}}$ from 4.1 to 3.9 ka is likely driven by increased salinity of surface water. This distinction becomes clearer when examining the $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$, where increased similarity from 4.8 to 3.9 ka (with a sharp increase at 4.1 ka) reflects the effect of increased sea surface salinity that reduces the δ¹⁸O difference between G. ruber and G. sacculifer. At the same time, weakened winter mixing increases $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{d-s}}$, which occurs from 4.2 to 3.9 ka. Importantly, the proxies also indicate that increased IWM mixing is generally positively correlated with increased Indus discharge and vice versa. The single time period when this does not hold true is 4.5–4.25 ka, when increased IWM mixing is coupled with decreased Indus discharge.

In summary, our multispecies approach using δ¹⁸O of G. ruber, G. sacculifer, and N. dutertrei allows us to differentiate between the strength of the IWM and freshwater discharge of the Indus River. We suggest that ISM strength decreased gradually from at least 4.8 ka, while the IWM strength peaked around 4.5–4.3 ka and then weakened afterwards. It is unlikely that the abrupt increase in G. ruber δ¹⁸O at 4.1 ka and low $\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$ could be caused solely by the decrease in IWM strength, even though IWM contributes to Indus River discharge. Weakening of the ISM must have played a substantial role in the 4.1 ka shift as well, indicated by the period 4.5–4.25 ka, when Indus discharge reflected a weak ISM ($\mathrm{\Delta }\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{s-r}}$) despite a phase of strengthened IWM.

5.3 Comparison to marine records

Other marine records from the Arabian Sea also suggest a gradual decrease in ISM strength from ∼5 ka (Gupta et al., 2003; Overpeck et al., 1996). Cullen et al. (2000) observed an abrupt peak in aeolian dolomite and calcite in marine sediments in the Gulf of Oman from 4.0 to 3.6 ka, and Ponton et al. (2012) also showed a shift to weaker ISM after 4.0 ka in the Bay of Bengal, based on δ¹³C of leaf waxes. Marine IWM reconstructions are not particularly coherent: although Doose-Rolinski et al. (2001) find a decrease in evaporation and weakening of the ISM between 4.6 and 3.7 ka, they argue this was accompanied by a relative increase in IWM strength. Giosan et al. (2018) inferred enhanced winter monsoon conditions from 4.5 to 3.0 ka based on a planktic paleo-DNA and percentage of Globigerina falconensis record close to our coring site (“C” in Fig. 1), which disagrees with our finding of decreased upper-ocean mixing after 4.3 ka. We suggest that the high stratigraphic (i.e., laminated) and chronological (i.e., 15 radiocarbon dates between 5.4 and 3.0 ka) resolution of core 63KA paired with a multispecies foraminifer δ¹⁸O record provides a robust history of the timing of changes in IWM and ISM strength, but additional studies are needed to resolve some of the discrepancies among the records.

Figure 7Comparison of the δ¹⁸O record of core 63KA with terrestrial records from the Indian subcontinent: (a, b) this study; (c) Berkelhammer et al. (2012); (d) Dixit et al. (2018); (e) Dixit et al. (2014); (f) this study; (g) Nakamura et al. (2016); (h) Kathayat et al. (2017). The mean value for each record indicated by the horizontal dashed lines is taken for all available data between 6.0 and 2.5 ka.

Download

5.4 Comparison to regional terrestrial records

The 63KA δ¹⁸O record obtained from three foraminifer species highlights several important ocean–atmosphere changes over the 5.4–3.0 ka time period. First, a sharp decrease occurred in both summer and winter precipitation at 4.1 ka, which is within a broader 300-year period of increased aridity spanning both rainfall seasons between 4.2 and 3.9 ka. In detail, we infer a relative decrease in Indus River discharge and weakened ISM between 4.8 and 3.9 ka, peaking at 4.1 ka, while a 200-year-long interval of strong IWM interrupted this period from 4.5 to 4.3 ka. Furthermore, the stepped change in δ¹⁸O of N. dutertrei suggests an enduring change in ocean–atmosphere conditions after 4.1 ka.

A relatively abrupt ∼4.2 ka BP climate event has been observed in several terrestrial records on the Indian subcontinent, most notably Mawmluh cave (∼4.1–3.9 ka) in northeastern India (Berkelhammer et al., 2012) and Kotla Dahar (∼4.1 ka) in northwestern India (Dixit et al., 2014) (Fig. 7). A less abrupt yet still arid period is documented in a peat profile (∼4.0–3.5 ka) from north central India (Phadtare, 2000), at Lonar Lake (∼4.6–3.9 ka) in central India (Menzel et al., 2014), and at Rara Lake (∼4.2–3.7 ka) in western Nepal (Nakamura et al., 2016). Finally, a recent study of oxygen and hydrogen isotopes in gypsum hydration water from Karsandi on the northern margin of the Thar Desert showed wet conditions between 5.1 and 4.4 ka, after which the playa lake dried out sometime between 4.4 and 3.2 ka (Dixit et al., 2018). Considering terrestrial records can record more local climatic conditions than marine records, it is remarkable that the records collectively agree on a period of regional aridity between 4.2 and 3.9 ka within the uncertainties of the age models that vary considerably among records.

However, not all records support this finding. For example, a reconstruction from Sahiya Cave in northwestern India shows an abrupt decrease in δ¹⁸O interpreted to reflect an increase in monsoon strength from ∼4.3–4.15 ka, followed by an arid trend after 4.15 ka (Kathayat et al., 2017). In addition, several other Thar Desert records do not identify a “4.2 ka BP event” sensu stricto but instead suggest that lakes dried out several centuries earlier (Deotare et al., 2004; Enzel et al., 1999; Singh et al., 1990) or later (Sinha et al., 2006) than 4.2 ka. This discrepancy may relate to nonlinear climate responses of lakes, which would not record a drought at 4.2 ka if they had already dried out earlier from the ongoing decrease in summer rainfall. In addition, there are also significant concerns about chronological uncertainties from the use of radiocarbon of bulk sediment for dating in some of these records. It is also possible that variations in the timing of climate change inferred from the terrestrial records may be real, reflecting a different sensitivity to ISM and IWM rain. As a marine record, core 63KA integrates large-scale ocean–atmosphere changes and therefore can help inform the interpretation of the more locally sensitive terrestrial records.

Figure 8Comparison of the δ¹⁸O record of core 63KA with more distant records: (a, b) this study; (c) Bar-Matthews et al. (2003); (d) Cheng et al. (2015); (e) Fleitmann et al. (2003). The mean value for each record indicated by the horizontal dashed lines is taken for all available data between 6.0 and 2.5 ka.

Download

More distantly, several terrestrial records in the Middle East also show a decrease in winter precipitation proxies around 4.2 ka: Jeita cave in Lebanon records a relatively dry period between 4.4 and 3.9 ka (Cheng et al., 2015) and Soreq cave in Israel shows a period of increased aridity starting at ∼4.3 ka (Bar-Matthews et al., 2003; Bar-Matthews and Ayalon, 2011) (Fig. 8). Lake Van in eastern Turkey also records reduced spring rainfall and enhanced aridity after ∼4.0 ka (Wick et al., 2003; Lemcke and Sturm, 1997). All of these records suggest a relatively arid period with reduced winter precipitation after ∼4.3 ka, as inferred from core 63KA. Qunf cave in Oman (Fleitmann et al., 2003), which is outside the range of IWM influence, instead shows a steady mid-Holocene weakening of the ISM that closely follows trends in summer solar insolation.

5.5 Cultural impacts

On the basis of our reconstruction of reduced IWM mixing after 4.3 ka, accompanied by decreased freshwater discharge of the Indus River, it is worth considering what impacts could be expected from a reduction in IWM and ISM precipitation. A weakened IWM overlying a reduced or more variable ISM would likely result in a distinct climate signal over the Indus River catchment, with broad implications for seasonal river flow and water availability throughout the year. The presence of the two rainfall systems creates a complex and diverse range of environments and ecologies across northwest South Asia (Petrie et al., 2017). In a situation when rainfall in both seasons is reduced over extended periods, step shifts in the natural environment may occur that are difficult to reverse (e.g., desertification, lake desiccation, regional vegetation changes, decline in overbank flooding, and shift in river avulsion patterns).

Societies reliant on IWM, ISM, or a combination of the two would have been vulnerable to years with monsoon failure, and a shift affecting both seasons will have challenged resilience and tested sustainability (Green et al., 2018; Petrie et al., 2017). Archaeological research into the transition from the urban Mature Harappan phase (∼4.6–3.9 ka) to the post-urban Late Harappan phase (∼3.9–3.6 ka) notes progressive de-urbanization through the abandonment of large Indus cities and a depopulation of the most western Indus regions, concurrent with a general trend towards an increase in concentrations of rural settlements in some areas of the eastern Indus extent (Green and Petrie, 2018; Petrie et al., 2017; Possehl, 1997) (Fig. S6). The relatively limited range of well-resolved available archaeo-botanical data suggests that there was a degree of diversity in crop choice and farming strategies in different parts of the Indus Civilization across this time span (Petrie et al., 2016; Petrie and Bates, 2017; Weber, 1999; Weber et al., 2010). Farmers in southerly regions appear to have focused on summer or winter crops, while the more northern regions of Pakistan Punjab and Indian Punjab and Haryana were capable of supporting combinations of winter and summer crops (Petrie and Bates, 2017). Although there is evidence for diverse cropping practices involving both summer and winter crops in the northern areas during the urban period, agricultural strategies appeared to favor the more intensive use of drought-resistant summer crops in the Late Harappan period (Madella and Fuller, 2006; Petrie and Bates, 2017; Pokharia et al., 2017; Weber, 2003; Wright, 2010). It has previously been suggested that a weakened ISM was a major factor in these shifts (e.g., Giosan et al., 2012; Madella and Fuller, 2006). Based on our reconstruction of a decreased IWM in northwest South Asia after 4.3 ka with a step shift at 4.1 ka, we suggest that both IWM and ISM climatic factors played a role in shaping the human landscape. This includes the redistribution of population to smaller settlements in eastern regions with more direct summer rain, as well as the shift to increased summer-crop-dominated cropping strategies.