Changes in productivity and intermediate circulation in the northern Indian Ocean since the last deglaciation: new insights from benthic foraminiferal Cd ∕ Ca records and benthic assemblage analyses

We have measured $\mathrm{Cd}/\mathrm{Ca}$ ratios of several benthic foraminiferal species and studied benthic foraminiferal assemblages on two cores from the northern Indian Ocean (Arabian Sea and northern Bay of Bengal, BoB), in order to reconstruct variations in intermediate-water circulation and paleo-nutrient content since the last deglaciation. Intermediate water Cd_w records estimated from the benthic $\mathrm{Cd}/\mathrm{Ca}$ reflect past changes in surface productivity and/or intermediate–bottom-water ventilation. The benthic foraminiferal assemblages are consistent with the geochemical data. These results suggest that during the last deglaciation, Cd_w variability was primarily driven by changes in intermediate-water properties, indicating an enhanced ventilation of intermediate–bottom water masses during both Heinrich Stadial 1 and the Younger Dryas (HS1 and YD, respectively). During the Holocene, however, surface primary productivity appears to have influenced Cd_w more than intermediate water mass properties. This is evident during the early Holocene (from 10 to 6 cal ka) when benthic foraminiferal assemblages indicate that surface primary productivity was low, resulting in low intermediate-water Cd_w at both sites. Then, from ∼ 5.2 to 2.4 cal ka, surface productivity increased markedly, causing a significant increase in the intermediate-water Cd_w in the southeastern Arabian Sea and the northeastern BoB. The comparison of intermediate-water Cd_w records with previous reconstructions of past Indian monsoon evolution during the Holocene suggests a direct control of intermediate-water Cd_w by monsoon-induced changes in upper-water stratification and surface primary productivity.

1 Introduction

During the last deglaciation, a two-step rapid increase in atmospheric CO₂ occurred during the 17–13.8 and 12.3–11.2 cal ka time intervals (e.g., Monnin et al., 2001). Several studies suggest that variations in the Southern Ocean circulation contributed to these increases in atmospheric CO₂ by transferring deep-ocean carbon to the upper ocean and atmosphere, through enhanced upwelling and increased northward penetration of the Antarctic Intermediate Water (AAIW) in all ocean basins (e.g., Marchitto et al., 2007; Anderson et al., 2009; Skinner et al., 2014). Different proxies have been used to reconstruct past changes in intermediate circulation, such as radiocarbon activity (Δ¹⁴C) (e.g., Marchitto et al., 2007; Bryan et al., 2010), benthic δ¹³C (e.g., Pahnke and Zahn, 2005; Jung et al., 2009; Ma et al., 2019), foraminiferal ε_Nd (e.g., Pahnke et al., 2008; Xie et al., 2012; Yu et al., 2018) and benthic foraminifera $\mathrm{Sr}/\mathrm{Ca}$ (Ma et al., 2020). These studies have focused on the close relationship between enhanced ventilation in the Southern Ocean and rising atmospheric CO₂ during the last deglaciation period. Furthermore, it has been shown that glacial–interglacial transfer of CO₂ between the oceans and the atmosphere could also be linked to changes in the efficiency of the oceanic biological pump (Pichevin et al., 2009; Ziegler et al., 2013; Bauska et al., 2016; Hertzberg et al., 2016; Jaccard et al., 2016; Yu et al., 2019), which may contribute to up to half of the observed CO₂ flux (Kohfeld et al., 2005).

The oceanic biological pump and nutrient upwelling are at least partly controlled by intermediate–deep-water circulation, contributing to the observed CO₂ changes (e.g., Toggweiler, 1999; Marchitto and Broecker, 2006). To track past changes in the nutrient concentration of intermediate water masses, benthic foraminifera $\mathrm{Cd}/\mathrm{Ca}$ has been used in many recent studies (e.g., Came et al., 2008; Poggemann et al., 2017; Valley et al., 2017; Umling et al., 2018); indeed, the benthic foraminifera $\mathrm{Cd}/\mathrm{Ca}$ is a robust proxy for seawater cadmium concentrations (Cd_w) (Boyle, 1988, 1992), which show a positive linear correlation with labile nutrients (phosphate and nitrate) in the modern ocean (e.g., Boyle et al., 1976; Boyle, 1988; Elderfield and Rickaby, 2000). The benthic foraminifera incorporate Cd as a function of Cd_w with a species-dependent partition coefficient (e.g., Tachikawa and Elderfield, 2002). Thus, the Cd measured in the fossil tests reflects the paleo-nutrient concentrations of the surrounding water masses and can be used to investigate past changes in intermediate-to-deep-ocean properties (e.g., Boyle and Keigwin, 1982; Oppo and Fairbanks, 1987; Came et al., 2008; Poggemann et al., 2017; Valley et al., 2017; Umling et al., 2018).

Complementary to the geochemical proxies, the type of benthic foraminifers and their abundance, both of which are related to organic flux and ecosystem oxygenation, make benthic foraminifer assemblages a powerful proxy for estimating past variations in bottom-water conditions (e.g., Corliss et al., 1986; Schmiedl et al., 1998; Almogi-Labin et al., 2000) in conjunction with organic matter fluxes to the seafloor (e.g., Altenbach et al., 1999; Van der Zwaan et al., 1999; Fontanier et al., 2002; Caulle et al., 2015). Benthic foraminifera have been successfully used as indicators of surface productivity, especially in high-carbon-flux regions (Schnitker, 1994). By comparing past benthic foraminiferal assemblages to modern ones, changes in food supply and oxygen concentrations of the bottom water can be reconstructed (e.g., Corliss, 1979; Peterson, 1984; Murgese and De Deckker, 2005). Recently, the combining of benthic foraminiferal assemblages and geochemical proxies has received increasing attention and has been used to reconstruct the evolution of surface productivity and upwelling intensity in the Indian Ocean (e.g., Hermelin 1991, 1992; Hermelin and Shimmield, 1995; Den Dulk et al., 1998; Murgese and De Deckker, 2005).

The Arabian Sea is one of the most productive regions of the ocean today (Banse, 1987; Marra and Barber, 2005). Surface productivity is dominated by the monsoon system, which has a strong impact on the distribution and dynamics of stratification and vertical mixing (Lévy et al., 2007). Numerous studies have focused on the reconstruction of the paleo-productivity of the Arabian Sea in relation to past changes in monsoon intensity (e.g., Prell and Kutzbach, 1987; Naidu and Malmgren, 1996; Gupta et al., 2003; Singh et al., 2006, 2011; Bassinot et al., 2011; Saraswat et al., 2014). By contrast, little is known about the paleo-productivity of the Bay of Bengal (BoB), especially its links to changes in monsoon precipitation (Phillips et al., 2014; Zhou et al., 2020). Consequently, studying paleo-productivity and past nutrient concentration of intermediate water masses in the northeastern Indian Ocean will also allow us to completely understand the influence of monsoon climate changes in tropical ocean ecology at different timescales. Moreover, as the benthic foraminiferal $\mathrm{Cd}/\mathrm{Ca}$ is a promising proxy to reconstruct the intermediate–deep-water nutrient content (e.g., Boyle and Keigwin, 1982; Tachikawa and Elderfield, 2002; Came et al., 2008; Poggemann et al., 2017; Valley et al., 2017), most of the studies referred to above have reconstructed deep–intermediate water masses in the past (e.g., Came et al., 2008; Bryan and Marchitto, 2010; Poggemann et al., 2017; Valley et al., 2017), and only few works investigate the relationship between the intermediate water mass nutrient and surface productivity (Bostock et al., 2010; Olsen et al., 2016). Furthermore, the evolution of the nutrient content of intermediate water masses since the last deglaciation has never been reconstructed in the Indian Ocean, where only two low-resolution $\mathrm{Cd}/\mathrm{Ca}$ records are available for deep water depths (Boyle et al., 1995), and, to our knowledge, none are available for intermediate water depths.

In this study, we provide, for the first time, two benthic foraminifera $\mathrm{Cd}/\mathrm{Ca}$ records at intermediate water depths in the northern Indian Ocean (Arabian Sea and northern Bay of Bengal). These data make it possible to estimate past changes in the nutrient content since the last deglaciation. We have also investigated benthic foraminiferal assemblages obtained from core MD77-191 (southeastern Arabian Sea) to help us reconstruct the conditions at the seafloor. Combined with planktonic foraminiferal δ¹⁸O, benthic δ¹³C and $\mathrm{Cd}/\mathrm{Ca}$ records obtained from the same core as well as with results already published in the Bay of Bengal (Ma et al., 2019, 2020), this study aims to document past variations in intermediate- and deep-water conditions and to decipher their links with surface paleo-productivity and intermediate-water ventilation.

2 Material and modern hydrological setting

We analyzed sediment core MD77-191 (07^∘30^′ N, 76^∘43^′ E, 1254 m) located in the Arabian Sea (off the southern tip of India), and core MD77-176 (14^∘30^′ 5^′′ N, 93^∘07^′6^′′ E, 1375 m) retrieved in the northeastern Bay of Bengal (BoB). These cores were collected in 1977 during the OSIRIS III cruise of the French R/V Marion Dufresne (Fig. 1).

Figure 1 (a) Oceanographic setting and locations of core MD77-191 in the Arabian Sea (red star), core MD77-176 in the Bay of Bengal (red star) and reference site SK237 GC04 (red circle, Naik et al., 2017). The black arrows represent the general surface circulation direction in the northern Indian Ocean during the summer southwest monsoon (Schott and McCreary, 2001). (b, c) Phosphate distribution along depth–latitude sections during the southwest monsoon and northeast monsoon periods, for the Arabian Sea and the Bay of Bengal, respectively. Data (in µmol kg⁻¹, colored scale; Schlitzer, 2000) were contoured and plotted using the Ocean Data View (ODV) software (Schlitzer, 2015). On these two figures are shown the distribution and circulation of water masses in the Arabian Sea and Bay of Bengal (black arrows). ASHS: Arabian Sea High Salinity Water; EIOW: Eastern Indian Ocean Water; BoBSW: Bay of Bengal surface waters; AAIW: Antarctic Intermediate Water; RSOW: Red Sea Overflow Water; IDW: Indian Deep Water.

The age model of core MD77-191 was established by using accelerator mass spectrometry (AMS) ¹⁴C dates obtained on nine monospecific samples of planktonic foraminifera Globigerinoides bulloides (Bassinot et al., 2011), one sample of pteropods (Mléneck, 1997) and three samples of the planktonic foraminifera Globigerinoides ruber (Ma et al., 2020). The average sedimentation rate of core MD77-191 is about 53 cm kyr⁻¹ and up to 90 cm kyr⁻¹ during the Holocene, providing a high-resolution, continuous record since 17 cal kyr BP.

The age model of core MD77-176 was previously established by using 31 planktonic foraminifer (G. ruber) AMS ¹⁴C dates combined with the core MD77-176 oxygen isotope record obtained on planktonic foraminifer G. ruber, which were correlated to the GISP2 Greenland ice core record (Marzin et al., 2013). Core MD77-176 displays high accumulation rates (average ∼25 cm kyr⁻¹ and up to 40 cm kyr⁻¹ during the Holocene).

In the modern ocean, the surface waters of the Arabian Sea and BoB are characterized by seasonally reversing currents that are driven by the monsoon winds (Fig. 1a). The surface water masses shallower than 150 m in the Arabian Sea are mainly Arabian Sea High Salinity Water (ASHS, 36.5 psu) (Talley et al., 2011). In the BoB, the surface waters above 100 m are designated Bay of Bengal surface waters (BoBSW), which have a low salinity (31 psu) due to large river inputs (Talley et al., 2011). Today, the northward extension of AAIW in the Indian Ocean rarely reaches beyond 10^∘ S (Lynch-Stieglitz et al., 1994). The sites of cores MD77-191 and MD77-176 are mainly bathed, therefore, by the North Indian Intermediate Water (Olson et al., 1993; Reid, 2003) with a potential contribution from the Red Sea Outflow Water (RSOW) for the site MD77-191 (Beal et al., 2000).

Due to the land–sea configuration in the north by Asia, the deep waters of the northern Indian Ocean originate from the south, including the Circumpolar Deep Water (CDW) and North Atlantic Deep Water (NADW) (You, 2000; Tomczak and Godfrey, 2003; Talley et al., 2011). Thus, between 1500 and 3800 m, the dominant deep water in the northern Indian Ocean is Indian Deep Water (IDW), originating from the CDW admixed with NADW (You, 2000; Tomczak and Godfrey, 2003; Talley et al., 2011). Then, on their pathway, the bottom water upwells when it expands northward in the northern Indian Ocean, returning to shallower depths (You, 2000; Fig. 1c). Therefore, variations in deep water masses can also influence the intermediate-depth waters in the northern Indian Ocean.

As far as surface waters are concerned, during the summer monsoon, the clockwise circulation in the Arabian Sea drives high-salinity waters from the northern to the southeastern Arabian Sea. By contrast, during the winter monsoon, the northeastern winds bring low-salinity water (BoBSW) from the BoB. The northern Indian Ocean, especially the Arabian Sea, is characterized by highly variable seasonal productivity (Shankar et al., 2002). Southwest winds during the summer season induce a strong Ekman pumping resulting in very active upwelling along the western coasts of the Arabian Sea and thus promoting strong surface productivity (Shankar et al., 2002; Fig. S1 in the Supplement). By contrast, the surface productivity in the BoB is generally weak compared with the Arabian Sea (e.g., Prasanna Kumar et al., 2001; Thushara and Vinayachandran, 2016; O'Malley, 2017; Fig. S1). In the BoB, large river inputs of freshwater and direct monsoon precipitation lead to more stable stratification in the upper ocean (Vinayachandran et al., 2002), and hence the vertical mixing of nutrients from the subsurface to the euphotic zone is generally limited (Gomes et al., 2000). However, the primary productivity of the western BoB shows a slight increase during the winter monsoon, as indicated by the distribution of chlorophyll in the surface water (Thushara and Vinayachandran, 2016; O'Malley, 2017; Fig. S1).

Modern data indicate that the southern-sourced intermediate water (AAIW) in the Indian Ocean has a phosphate concentration of about 2–2.5 µmol kg⁻¹ (Fig. 1b and c). In the northern intermediate Indian Ocean, the phosphate concentration is significantly higher, ranging from 2.75 to 3 µmol kg⁻¹ in the Arabian Sea during the summer monsoon, and from 2.5 to 2.75 µmol kg⁻¹ in the BoB during the winter monsoon (Fig. 1b and c). The higher phosphate in the northern Indian Ocean can been linked to increased primary productivity (Banse, 1987; Marra and Barber, 2005).

3 Methods

3.1 Cd $/$ Ca analysis

In order to improve understanding of possible inter-species differences and microhabitat effects on the benthic $\mathrm{Cd}/\mathrm{Ca}$ records, we analyzed $\mathrm{Cd}/\mathrm{Ca}$ in three calcite (Cibicidoides pachyderma, Uvigerina peregrina and Globobulimina spp.) and one aragonite (Hoeglundina elegans) benthic foraminiferal species from core MD77-191. C. pachyderma is a shallow infaunal species U. peregrina and Globobulimina spp. are endobenthic species with intermediate and deep microhabitats, respectively (Fontanier et al., 2002). In core MD77-176, due to the limitation of calcitic species, we only measured $\mathrm{Cd}/\mathrm{Ca}$ ratios in H. elegans shells. Moreover, $\mathrm{Mn}/\mathrm{Ca}$, $\mathrm{Fe}/\mathrm{Ca}$ and $\mathrm{Al}/\mathrm{Ca}$ ratios were also measured in all benthic foraminiferal samples to check the robustness of $\mathrm{Cd}/\mathrm{Ca}$ results and the potential influence of contamination (i.e., oxides and sedimentary clay, Barker et al., 2003).

Each sample contained between 10 and 15 individuals picked from the 250–315 µm size fraction. Samples were gently crushed and cleaned to remove clays, organic matter and elemental oxides by using reductive and oxidative cleaning following previously published methods (Boyle and Keigwin, 1982; Barker et al., 2003). Each sample was dissolved in 0.075 N HNO₃ and analyzed using a single collector sector field high-resolution inductively coupled plasma mass spectrometer (HR-ICP-MS) Thermo Element XR hosted at the GEOPS Laboratory (University Paris-Saclay, France).

The detailed instrumental settings and mother standard solutions are described in Ma et al. (2020). A blank consisting of the same 0.1 N HNO₃ used to dilute the standards and samples was also analyzed. We removed the blank intensity values from all the raw intensities (including standards), and raw data were linearly drift-corrected by interspersing a drift standard every four samples. Standard curves were used to calculate elemental concentrations, coefficients of determination (r²) always being >0.9999 for all elemental ratios. The mean reproducibility and accuracy are 3.6 % and 7.5 %, respectively.

3.2 Faunal analysis

Benthic foraminiferal assemblages from core MD77-176 have already been published in Ma et al. (2019). For core MD77-191, a total of 72 samples were collected for benthic foraminiferal assemblage determinations. In each sample, benthic foraminifera (>150 µm) were extracted, counted and identified to species level following the taxonomical descriptions of various authors (e.g., Loeblich and Tappan, 1988; Jones, 1994; Holbourn et al., 2013). For core MD77-191, there is no material left in this old, small-diameter core, and so we used samples obtained earlier for stable isotope studies. Since the bulk weights of these samples were not recorded prior to sieving, we could not perform the calculation of the absolute abundance of foraminifera or accumulation rates. Thus, we only converted the individual counts to percentages with respect to the total benthic foraminifera present in each sample. In order to describe major faunal variations, we performed principal component analysis (PCA) on the variance–covariance matrix using the PAST software (Version 3.0, Hammer et al., 2001). Species present with a percentage >1 % in at least 1 sample were used for statistical analysis and diversity calculation.

4 Results

4.1 Elemental ratio results

To check the influence of oxide contaminants on the elemental ratios, $\mathrm{Mn}/\mathrm{Ca}$ was systematically measured. The $\mathrm{Mn}/\mathrm{Ca}$ of H. elegans from cores MD77-191 and MD77-176 ranges between 6.5–10 and 1–30 µmol mol⁻¹, respectively. Such ranges are much lower than the 100 µmol mol⁻¹ limit proposed by Boyle (1983). The $\mathrm{Mn}/\mathrm{Ca}$ values obtained on the three calcite benthic foraminifera species from core MD77-191 – C. pachyderma (5–18 µmol mol⁻¹), U. peregrina (3–23 µmol mol⁻¹) and Globobulimina spp. (4–69 µmol mol⁻¹) – are also all below 100 µmol mol⁻¹ (Boyle, 1983). The $\mathrm{Fe}/\mathrm{Ca}$ ratios are also lower than 1 mmol mol⁻¹ in all samples from cores MD77-191 and MD77-176, in agreement with the limit proposed by Barker et al. (2003). In addition, Barker et al. (2003) concluded that no significant pollution by clay minerals would be expected when $\mathrm{Al}/\mathrm{Ca}$ is <0.5 mmol mol⁻¹. In all our samples, $\mathrm{Al}/\mathrm{Ca}$ is below 0.5 mmol mol⁻¹, indicating that the sample cleaning procedure was efficient.

All of the above results indicate that our samples were not affected by contamination.

Cd $/$ Ca

The $\mathrm{Cd}/\mathrm{Ca}$ records of C. pachyderma, U. peregrina and Globobulimina spp. from core MD77-191 range between 0.07–0.2, 0.07–0.14 and 0.03–0.09 µmol mol⁻¹, respectively (Fig. 2d; Table S1 in the Supplement).

Figure 2 (a) GISP2 Greenland ice core δ¹⁸O signal (Stuiver and Grootes, 2000). (b, c) Globigerinoides ruber δ¹⁸O records of cores MD77-191and MD77-176, respectively (Marzin et al., 2013; Ma et al., 2020). (d) $\mathrm{Cd}/\mathrm{Ca}$ records of the benthic foraminifera Hoeglundina elegans (black), Cibicidoides pachyderma (green), Uvigerina peregrina (blue) and Globobulimina spp. (orange) obtained from core MD77-191. (e) $\mathrm{Cd}/\mathrm{Ca}$ records of the benthic foraminifera H. elegans from core MD77-176. EHCO stands for Early Holocene Climate Optimum, YD for Younger Dryas, B-A for Bølling–Allerød and HS1 for Heinrich stadial 1.

The $\mathrm{Cd}/\mathrm{Ca}$ records for the calcite benthic species C. pachyderma and U. peregrina have very low time resolutions during the last deglaciation. However, some common patterns can be observed. The $\mathrm{Cd}/\mathrm{Ca}$ records of C. pachyderma and U. peregrina show lower values during the Heinrich stadial 1 (HS1, 17–15.2 cal ka) and the Younger Dryas (YD, 13–11 cal ka) cold periods, with average values of ∼0.08 µmol mol⁻¹ for C. pachyderma and ∼0.09 µmol mol⁻¹ for U. peregrina. By contrast, these two species display higher $\mathrm{Cd}/\mathrm{Ca}$ ratios (∼0.12 µmol mol⁻¹) during the Bølling–Allerød warm period (B-A, 15–13.3 cal ka) compared with the HS1 and YD. Then, lower values (∼0.1 µmol mol⁻¹ for C. pachyderma; 0.11 µmol mol⁻¹ for U. peregrina) are observed during the early Holocene (10–5 cal ka) compared to larger variations occurring in the late Holocene (5.2–2.4 cal ka). The $\mathrm{Cd}/\mathrm{Ca}$ record of deep infaunal Globobulimina spp., obtained at a lower time resolution, shows different variations compared with the two other taxa without any clear trend during the Holocene.

The H. elegans $\mathrm{Cd}/\mathrm{Ca}$ values of core MD77-191 range from 0.05 to 0.31 µmol mol⁻¹ since 17 cal kyr BP (Fig. 2d; Table S1). Depleted values at about 0.07 µmol mol⁻¹ are recorded from the last deglaciation to the early Holocene (17–5 cal ka time interval). During the HS1 and the YD time intervals, a significant decrease of ∼0.06 µmol mol⁻¹ occurred (even when taking into consideration the analytical error bar of ±0.02, 2σ), and a slight increase (0.09 µmol mol⁻¹) is observed between 15 and 13.3 cal ka (B-A period). A rapid increase in the $\mathrm{Cd}/\mathrm{Ca}$ values beginning at 5.2 cal kyr BP reaches a maximum (0.31 µmol mol⁻¹) during the late Holocene.

For core MD77-176, the H. elegans $\mathrm{Cd}/\mathrm{Ca}$ records range between 0.06 and 0.17 µmol mol⁻¹ over the past 18 cal kyr BP (Fig. 2e; Table S1), without no clear trends and average benthic $\mathrm{Cd}/\mathrm{Ca}$ values of ∼0.09 µmol mol⁻¹ during the different periods (HS1, YD and Holocene). However, the benthic $\mathrm{Cd}/\mathrm{Ca}$ record during the Holocene seems to exhibit a slight increase both in value and range of variations after 6 cal kyr BP.

4.2 Foraminifera assemblages of core MD77-191

Benthic foraminiferal species richness ranges between 16 and 36, and the total abundance fluctuates between 82 and 642 specimens (Table S2). Hyaline species are the dominant constituents (>80 %) and mainly consist of Bulimina aculeata, H. elegans, C. pachyderma, Uvigerina spp., Gyroidina broeckhiana, Globocassidulina subglobosa, Sphaeroidina bulloides, Gyroidinoides spp., Lenticulina spp., Melonis barleeanum and Globobulimina spp. (including Praeglobobulimina spp.) (in decreasing order of relative average abundance). Agglutinated taxa reach on average about 1.6 % and consist of Textularia sp., Martinottiella communis and Eggerella bradyi. The average percentage of porcelaneous species, characterized by Pyrgo elongata, Pyrgo murrhina, Pyrgo depressa, Pyrgoella irregularis, Quinqueloculina spp., Sigmoilopsis schlumbergeri and Spiroloculina spp., is about 5.1 %.

Furthermore, we merged species that share an ecological similarity, such as Globobulimina affinis, Globobulimina pacifica and Praeglobobulimina spp. into Globobulimina spp. A total of 74 samples and 55 groups and/or species were adopted to perform principal component analysis (PCA) in order to identify major faunal trends. The PCA analysis suggests that the benthic foraminifera could be grouped into three assemblages, with PC1 (positive and negative loadings) and PC2 (positive loadings) representing 42 % and 19 % of the total variance, respectively (Table 1). Moreover, compared with the total variance of PC1 and PC2, PC3 is the largest one and only explains 8 % of the total variance for the rest of the principal component scores (PCs). The species composition consists of H. elegans, Globobulimina spp. (positive loadings), Uvigerina peregrina and C. pachyderma (negative loadings) (Table 1). It seems that the main composition of assemblages (PC3) is quite similar to PC1 and does not show more information about the bottom conditions. Therefore, we only focus on PC1 and PC2 in the paper for the interpretation and do not present other PCs in the discussion.

Table 1 Species composition of benthic foraminiferal assemblages from core MD77-191.

Assemblage 1 can be defined as the combination of Bulimina aculeata and C. pachyderma, together with Pullenia bulloides and Ehrenbergina trigona (Figs. 3 and S2) and display high positive PC1 loadings. This assemblage dominated the foraminiferal record during the late Holocene (between 6 and 1.4 cal ka).

Figure 3 Downcore variations in PC scores and the percentages of major species. (a) Sphaeroidina bulloides and (b) Gyroidinoides orbicularis are dominated assemblage 3, (c) Cibicidoides wuellerstorfi and (d) Hoeglundina elegans are the main associated species of assemblage 2, (e) Cibicidoides pachyderma and (f) Bulimina aculeata are major species in assemblage 1. The color-shaded intervals and abbreviations are the same as in Fig. 2.

By contrast, assemblage 2, dominated by H. elegans and Bulimina manginata, exhibits high negative PC1 loadings, and corresponds to the record during the early Holocene (Figs. 3 and S2). Other quantitatively important contributors are C. wuellerstorfi and Globocassidulina subglobosa (Fig. S2).

Then, assemblage 3, dominated by Sphaeroidina bulloides and Gyroidinoides orbicularis, corresponds to the positive loadings of PC2, which is more important during the last deglaciation (Figs. 3 and S2). The associated species of assemblage 3 are Bulimina mexicana and Gyroidinoides soldanii (Fig. S2).

However, as the main composition of PC2 negative loadings is dominated by the same benthic species in assemblages 1 and 2, it is difficult to glean any additional information from this analysis. Thus, to clarify the discussion, we prefer to use three assemblages in the following rather than the two PCs.

5 Discussion

5.1 Past intermediate-water Cd_w concentrations from the northern Indian Ocean

In the modern ocean, benthic foraminifera $\mathrm{Cd}/\mathrm{Ca}$ shows a positive correlation with Cd_w and dissolved nutrients (phosphate and nitrate) (Boyle et al., 1976; Hester and Boyle, 1982). As aragonitic benthic foraminifera H. elegans faithfully records the bottom-water Cd concentrations (Cd_w), $\mathrm{Cd}/\mathrm{Ca}$ ratios can be converted to seawater Cd_w with the appropriate relationship (Eq. 1), where the partition coefficient D_p≈1 for all water depths (Boyle et al., 1995; Bryan and Marchitto, 2010).

$$\begin{array}{}\text{(1)}& D_{\mathrm{p}}={\displaystyle \frac{{\left(\mathrm{Cd}/\mathrm{Ca}\right)}_{\mathrm{foram}}}{{\left(\mathrm{Cd}/\mathrm{Ca}\right)}_{\mathrm{water}}}}\end{array}$$

In contrast, the partition coefficient for calcite species changes with water depth. For water depths between 1150 and 3000 m, D_p was calculated based on the equation of Boyle (1992; Eq. 2). The seawater Ca concentration is assumed to be at a constant mean value of 0.01 mol kg⁻¹ (Boyle, 1992).

$$\begin{array}{}\text{(2)}& D_{\mathrm{p}}=\mathrm{1.3}+\left(\text{depth}-\mathrm{1150}\right)\times \left(\mathrm{1.6}/\mathrm{1850}\right)\end{array}$$

The intermediate Cd_w results based on the H. elegans $\mathrm{Cd}/\mathrm{Ca}$ values of core MD77-191 range from 0.5 to 3.1 nmol kg⁻¹ since 17 cal kyr BP (Fig. 4a), with a core top value of 0.80 nmol kg⁻¹ in agreement with the estimated intermediate water depth modern Cd_w (∼0.83 nmol kg⁻¹) in the northern Indian Ocean (Boyle et al., 1995). The intermediate Cd_w was also calculated from calcite benthic species C. pachyderma, U. peregrina and Globobulimina spp. from core MD77-191, with values ranging between 0.53–1.48, 0.52–1.04 and 0.26–0.65 µmol mol⁻¹, respectively (Fig. 4a). The Cd_w values of C. pachyderma and U. peregrina are within the same range. However, the H. elegans Cd_w values are higher than those from the two calcite species, especially during the late Holocene. Moreover, the core top data of C. pachyderma and U. peregrina are also lower (∼0.7 and 0.69 nmol kg⁻¹, respectively) than the modern estimated Cd_w data (∼0.83 nmol kg⁻¹) in the northern Indian Ocean (Boyle et al., 1995) (Fig. 4a). These depleted Cd_w values may be related to the benthic foraminiferal microhabitat effect; indeed, U. peregrina is known to be strictly a shallow infaunal species, as well as C. pachyderma (Fontanier et al., 2002), differing from strictly epifaunal taxa, such as Cibicidoides wuellerstorfi (Mackensen et al., 1993).

Figure 4 (a) Cd_w records calculated based on the $\mathrm{Cd}/\mathrm{Ca}$ of benthic foraminifera Hoeglundina elegans (black), Cibicidoides pachyderma (green), Uvigerina peregrina (blue) and Globobulimina spp. (orange) obtained from core MD77-191. (b) Cd_w record from core MD77-176 reconstructed using H. elegans $\mathrm{Cd}/\mathrm{Ca}$; the red line is the smoothed curves using a two-point moving average. The red stars represent the modern Cd_w (∼0.83 nmol kg⁻¹) in the northern Indian Ocean (Boyle et al., 1995). The color-shaded intervals and abbreviations are the same as in Fig. 2.

Moreover, the deep infaunal Globobulimina spp. Cd_w displays relatively much lower values and does not exhibit strong variations compared to the other species investigated in this study, displaying a general increasing trend from the last deglaciation to the Holocene. As Globobulimina spp. correspond to deep benthic infaunal species, this result may indicate a stable nutrient content of pore water as compared to other benthic taxa associated with bottom water (Fig. 4a). Thus, when tracking past changes in the bottom-water Cd_w concentrations, the use of a strictly epifaunal species living at the water–sediment interface such as H. elegans appears to be more robust than using endofaunal species that live in contact with pore water.

Relative variations in the Cd_w obtained from C. pachyderma and U. peregrina are in good agreement with the records obtained on H. elegans. Variations in H. elegans Cd_w during the last deglaciation indicate a decrease of about ∼0.6 nmol kg⁻¹ in the HS1 and YD periods, with a slight increase (0.9 nmol kg⁻¹) during the warm B-A. Cd_w results from core MD77-191 indicate a shift from the last deglaciation (∼0.7 nmol kg⁻¹) to the late Holocene (∼1.59 nmol kg⁻¹). During the Holocene, the Cd_w records display relatively low values of around 0.9 nmol kg⁻¹ in the 10–6 cal ka time interval and show a major shift at around 6.4 cal kyr BP with values rising up to 3.1 nmol kg⁻¹.

For core MD77-176, the intermediate-water Cd_w calculated from the H. elegans $\mathrm{Cd}/\mathrm{Ca}$ records ranges between 0.6 and 1.7 nmol kg⁻¹ over the past 18 cal kyr BP (Fig. 4b). Compared with intermediate Cd_w from MD77-191, the Cd_w record of core MD77-176 does not display any clear trend from the last deglaciation to the Holocene. However, a slight increase is observed since 6 cal kyr BP, in agreement with the MD77-191 intermediate Cd_w records. In addition, even though the MD77-176 record has a lower time resolution, it displays a shorter maximum (1.3 nmol kg⁻¹) during the 13.4–11 cal ka time interval.

To summarize, among the three calcite benthic taxa and the aragonitic benthic species H. elegans, the $\mathrm{Cd}/\mathrm{Ca}$ records of H. elegans appear to be the most suitable for tracking past Cd_w changes at intermediate water depth through time. Thus, in the following discussion, we will only focus on the intermediate Cd_w calculated from the H. elegans $\mathrm{Cd}/\mathrm{Ca}$ from both studied cores.

5.2 Comparison between geochemical records and benthic foraminiferal assemblages

Comparing the geochemical records to the benthic assemblages, we can observe similar patterns. For core MD77-191 from the southeastern Arabian Sea, three benthic assemblages were identified since the last deglaciation. S. bulloides and Gyroidinoides orbicularis are major components of assemblage 3 (during the last deglaciation), together with B. mexicana and Gyroidinoides soldanii (Figs. 3 and S2). S. bulloides and B. mexicana are found in intermediate to high organic carbon flux rate regions (e.g., Schmiedl et al., 2000; Eberwein and Mackensen, 2006, 2008), while G. orbicularis and G. soldanii are associated with well-oxygenated and oligotrophic environments (Peterson, 1984; Burmistrova and Belyaeva, 2006; De and Gupta, 2010). Thus, assemblage 3 reflects mesotrophic environments and/or well-ventilated conditions during the last deglaciation. Although millennial-scale changes in the benthic foraminiferal assemblages during the last deglaciation could not be observed, assemblage 3 seems at least partly consistent with previous studies in the northern Indian Ocean based on multiple geochemical proxies (e.g., benthic δ¹³C, intermediate water [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] and ε_Nd records); these studies have revealed the presence of better-ventilated waters, which might correspond to AAIW, during the HS1 and YD (e.g., Yu et al., 2018; Ma et al., 2019, 2020).

Benthic foraminiferal assemblage 2 predominates during the early Holocene and is characterized by H. elegans and B. manginata as major contributors (Figs. 3 and S2). The other important contributors are C. wuellerstorfi and G. subglobosa. B. manginata is found in high organic carbon flux rate conditions (De Rijk et al., 2000; Eberwein and Mackensen, 2006, 2008). However, previous studies on H. elegans, C. wuellerstorfi and G. subglobosa indicate that these species correspond to high levels of dissolved oxygen and oligotrophic settings (e.g., Altenbach et al., 1999; Fontanier et al., 2002; Murgese and De Deckker, 2005, 2007; De and Gupta, 2010). Periods dominated by these taxa probably indicate high oxygen levels and an oligotrophic environment. This is consistent with previous studies in the area, based on benthic foraminiferal δ¹³C and Δ¹⁴C age difference (e.g., Naqvi et al., 1994; Bryan et al., 2010) (Fig. S3). Indeed, the glacial to Holocene benthic δ¹³C shifts (0.35 ‰–0.4 ‰, vs. Pee Dee Belemnite, PDB) at intermediate–deep water depth in the northern Indian Ocean are interpreted as reflecting an increased contribution of better-ventilated deep water NADW in IDW, during the Holocene (e.g., Naqvi et al., 1994; Ma et al., 2019) (Fig. S3). Furthermore, the increased benthic–planktonic (B-P) age offsets and depleted ε_Nd records obtained from the same core site could also reflect the enhanced influence of NADW in IDW during the Holocene, which is characterized by well-ventilated conditions and depleted nutrient concentrations (modern Cd_w, ∼0.2 nmol kg⁻¹) (Poggemann et al., 2017; Yu et al., 2018; Ma et al., 2019). The impact of this change in the IDW composition can be recorded at intermediate water depth since the deep water masses are transformed to an upward flow during their pathway, thus being a potential contribution to intermediate-depth water masses (Naqvi et al., 1994; You, 2000). Although the intermediate benthic δ¹³C record from core MD77-191 is missing for the Last Glacial Maximum (LGM), the average value for the Holocene (∼0.31 ‰, vs. PDB) is consistent with previous studies carried out in the northern Indian Ocean; combined with the opposite trend between δ¹⁸O_ivc records and intermediate-water temperature from MD77-191 (Ma et al., 2020), all these records suggest well-ventilated conditions (Fig. S3). To summarize, the predominance of benthic foraminifera assemblage 2 in the early Holocene seems to reflect better-ventilated water masses, related to an enhanced contribution of NADW in IDW at the core site, as already observed in previous studies (Poggemann et al., 2017; Yu et al., 2018; Ma et al., 2019, 2020).

By contrast, B. aculeata and C. pachyderma are major components of assemblage 1 (during the late Holocene), together with P. bulloides and E. trigona (Figs. 3 and S2). Living B. aculeata have a widespread distribution, with a preference for water depths ranging from 1500 to 2500 m, and are typically associated with high organic carbon fluxes (Mackensen et al., 1995; Almogi-Labin et al., 2000; Caulle et al., 2015). P. bulloides is a shallow infaunal species, which prefers mesotrophic environments and shows adaptability with respect to oxygen concentration in the Arabian Sea (Gupta and Thomas, 1999; Caulle et al., 2015). E. trigona is commonly recorded in low-oxygen habitats (Caulle et al., 2015). We thus interpret assemblage 1 as indicating relatively low-oxygen and meso- to eutrophic bottom-water conditions during the late Holocene (6–1.4 cal ka). However, the lower oxygen concentrations reflected by benthic fauna 1 seem to be the opposite of what would be expected under an enhanced influence of better-ventilated NADW in IDW during the Holocene in the northern Indian Ocean. Thus, another process has to be explored to combine our observations. To do that, we can use the relative abundance of Globigerina bulloides, a proxy for upwelling activity, that increased in the late Holocene in core MD77-191, suggesting an increased productivity in the southeastern Arabian Sea (Bassinot et al., 2011) (Fig. 5). This record is synchronous with the benthic foraminiferal assemblage 1 (during the late Holocene). Thus, increased surface productivity during the late Holocene could have introduced more organic matter into the intermediate water, leading to depleted oxygen conditions. By contrast, benthic assemblages 2 and 3 (during the last deglaciation and early Holocene; 17–6 cal ka) are associated with low G. bulloides abundances, suggesting lower productivity in the southeastern Arabian Sea during this period (Bassinot et al., 2011) and thus indicating that intermediate water masses were characterized by higher bottom-water oxygen conditions and a lower flux of organic matter. Therefore, all of these elements suggest that changes in primary productivity seem to be an important factor impacting the distribution of benthic assemblages at the core MD77-191 site, especially during the Holocene.

Figure 5 (a) Organic carbon weight percentage (%C_org) and (b) G. bulloides percentage from core SK237 GC04 (1245 m, Arabian Sea, Naik et al., 2017). (c) Relative abundance of G. bulloides (Mléneck, 1997; Bassinot et al., 2011), (d) PC 1 scores and (e) Cd_w records from core MD77-191 (Arabian Sea). The color-shaded intervals and abbreviations are the same as in Fig. 2.

In order to examine the relationships between intermediate Cd_w and these different processes (surface productivity and/or water mass ventilation) in the eastern Arabian Sea, we can compare the MD77-191 Cd_w values with the relative abundance of G. bulloides and benthic foraminiferal assemblage analyses from the same core MD77-191, together with the records for C_org and the G. bulloides percentage obtained from core SK237 GC04 (1245 m, southeastern Arabian Sea, Naik et al., 2017) (Fig. 5). Indeed, the total organic carbon (C_org) could also be used as a qualitative indicator of past productivity and/or bottom-water ventilation changes (Naidu et al., 1992; Canfield, 1994; Calvert et al., 1995; Naik et al., 2017). Despite a lower resolution for MD77-191 H. elegans Cd_w records, when compared to the C_org and the G. bulloides percentage from core SK237 GC04, all of them seem to exhibit similar trends at a long timescale even though some small-scale discrepancies can be observed at a millennial timescales (Fig. 5). From the last deglaciation to the late Holocene, the Cd_w record displays a significant shift from ∼0.7 nmol kg⁻¹ to about twice the values of ∼1.59 nmol kg⁻¹. The intermediate Cd_w values are thus extremely high during the late Holocene and synchronous with the higher values of C_org and G. bulloides percentage records. These observed similar trends suggest that the increased surface productivity at the core site during the late Holocene is associated with higher intermediate Cd_w values. Moreover, previous studies have suggested that increased Cd_w values (>1 nmol kg⁻¹) could correspond to elevated surface productivity (Bostock et al., 2010; Olsen et al., 2016). However, at a millennial timescale, we also observed several decreases in intermediate Cd_w values (∼0.81 nmol kg⁻¹) during the late Holocene, reaching nearly similar values during the last deglaciation (Fig. 5). Thus, the variations in the Cd_w values cannot be fully associated with variations in the surface productivity.

As mentioned before, during the Holocene, an increased influence of NADW in IDW was observed in the northern Indian Ocean (Yu et al., 2018; Ma et al., 2019, 2020). NADW is characterized by a depleted nutrient content (modern Cd_w, ∼0.2 nmol kg⁻¹; Poggemann et al., 2017), and its contribution to IDW may affect the intermediate Cd_w by deep water masses upwelling when flowing northward. However, during the late Holocene, benthic foraminiferal assemblage 1 is associated with lower oxygen concentrations, which seem to be inconsistent with an enhanced influence of better-ventilated NADW in IDW in the northern Indian Ocean. Therefore, this apparent discrepancy seems to indicate that deep–intermediate water mass variations are not an important control during the Holocene in this area, although we could not fully exclude the influence of NADW in IDW at a millennial timescale. Moreover, there is no clear evidence for such a millennial-scale variability in the IDW and/or NADW circulation in the studied area. Thus, we suggest the intermediate Cd_w at the core MD77-191 site may be mainly influenced by surface productivity, especially during the Holocene.

In the Bay of Bengal, the benthic assemblages of core MD77-176 suggest that the intermediate water masses were characterized by oligotrophic to mesotrophic conditions and/or well-ventilated environments during the Holocene (Ma et al., 2019), associated with much lower surface productivity (Fig. S4). This observation is in agreement with low primary productivity during the Holocene reconstructed by the relative abundance of coccolith species Florisphaera profunda from the same core MD77-176 in the northeastern BoB (Zhou et al., 2020). In the modern ocean, Prasanna Kumar et al. (2001) indicate that primary productivity in the BoB is much lower than in the Arabian Sea, the lower surface productivity resulting from the large freshwater input from rivers and direct rainfall resulting from enhanced Indian summer monsoon precipitation (e.g., Vinayachandran et al., 2002; Madhupratap et al., 2003; Gauns et al., 2005). Moreover, when we compare the average Cd_w value of core MD77-176 from the BoB (∼0.9 nmol kg⁻¹) with results from core MD77-191 in the Arabian Sea (∼1.2 nmol kg⁻¹), lower values, especially during the late Holocene, are in agreement with the benthic assemblages.

To sum up, variations in the benthic assemblages seem to be associated with changes in the deep water mass ventilation and/or organic matter flux, linked to surface productivity. The benthic foraminiferal fauna are consistent with the Cd_w record of core MD77-191 particularly during the late Holocene (6–1.4 cal ka). Thus, our results seem to show that the Cd_w record is mainly controlled by changes occurring at the surface, especially during the Holocene. However, at millennial timescales, such during the HS1 and YD, the percentages of planktonic species G. bulloides from cores MD77-191 and SK237 GC04 all indicate modest paleo-productivity, the opposite of what is suggested by the results of core MD77-191 Cd_w and the C_org record obtained from core SK237 GC04. This interval is also marked by enriched G. ruber δ¹⁸O values, indicating a weaker monsoon and reduced freshwater inputs (Naik et al., 2017). This apparent discrepancy may be related to changes in the intermediate water mass sources and/or ventilation during the last deglaciation.

So, in the next sections, we discuss (i) processes controlling surface productivity and (ii) changes in the intermediate-water circulation, both of them being potential drivers of the observed variations.

5.3 Relationships between primary productivity and monsoon intensity

During the Holocene, the intermediate-water Cd_w records obtained from cores MD77-191 and MD77-176 seem to display depleted values in the early Holocene, followed by an abrupt increasing trend at the middle Holocene and then reaching higher values on average (despite a short-timescale variability) during the late Holocene.

Of the two cores, core MD77-176, located in the northeastern BoB, shows the lowest intermediate Cd_w (down to ∼0.83 nmol kg⁻¹) during the 10–6 cal ka time interval. Observations described above suggest that this low in Cd_w resulted from low primary productivity and thus reduced fluxes of organic matter to the intermediate depths. We attribute this evolution to monsoon variation. Indeed, the early Holocene Climate Optimum (10–6 cal ka) is characterized by enhanced monsoon precipitation (Marzin et al., 2013; Contreras-Rosales et al., 2014) (Fig. 6d–f) that resulted in increased freshwater discharge from the Ganges–Brahmaputra river system and from the Irrawaddy River. However, the distribution of chlorophyll in surface water of the western BoB suggests a low annual productivity, indicating that the BoB is not significantly influenced by the riverine nutrient input (Zhou et al., 2020). Thus, it is likely that this increase in freshwater drove pronounced ocean stratification in the northeast BoB, which could impede the nutrient transfer from the intermediate and/or deep layer to the euphotic upper seawater column and then induce low productivity. A similar low in Cd_w values is observed in the reconstructed intermediate-water Cd_w record from core MD77-191 during the early Holocene, with values descending to ∼0.92 nmol kg⁻¹ in the 10–6 cal ka time interval. These low values of intermediate Cd_w are coeval with low surface productivity as recorded by the G. bulloides percentage and low values in C_org content from SK237 GC04 in the Arabian Sea (Fig. 5). These variations are also recorded in changes in benthic assemblages, with the occurrence of assemblage 2 associated with high oxygen levels and an oligotrophic environment (Fig. 3). Off the southern tip of India, we cannot reject the possibility that increased monsoon precipitation and enhanced freshwater runoffs in the BoB during the early Holocene, inducing a stronger stratification, could explain part of the decrease in surface primary productivity. Yet, at this site, another explanation prevails which is related to the decrease in summer monsoon wind intensity that drives local Ekman pumping. As shown by Bassinot et al. (2011), the productivity variations at the southern tip of India are inversely related to the evolution of upwelling activity along the Oman Margin, to the west of the Arabian Sea. Based on a data–model comparison, Bassinot et al. (2011) showed that this anti-correlation can be attributed to the northward shift in the intertropical convergence zone (ITCZ) when boreal summer insolation reached a maximum in the early Holocene (Fig. 6a); this ITCZ location results in enhanced summer monsoon wind intensity and an increase in the associated Ekman pumping in the west of the Arabian Sea and along the Oman margin, while it weakens at the southern tip of India. This process may thus induce a decrease in surface productivity in the southeastern Arabian Sea.

Figure 6 (a) The solar insolation at 10^∘ N in summer (Laskar et al., 2004). (b, c) Intermediate Cd_w calculated from H. elegans obtained from MD77-176 and MD77-191, respectively. (d) Lonar Lake δ¹³C_wax record (Sarkar et al., 2015). (e) δD_Alk−ic record from core SO188-342KL (Contreras-Rosales et al., 2014). (f) Seawater δ¹⁸O anomaly obtained from MD77-176 (Marzin et al., 2013). The color-shaded intervals and abbreviations are the same as in Fig. 2.

In addition, Naik et al. (2017) pointed out the co-existence of low productivity during the early Holocene in the BoB and to the south of India, in agreement with our data that clearly show the impact of such a reduction in surface primary productivity on the intermediate-water Cd_w. These authors suggested a direct relationship between intense monsoon rainfall and reduced surface productivity. However, the northeastern BoB received a much larger amount of river input than the southern tip of India during the early Holocene (Marzin et al., 2013). Thus, it seems reasonable to propose that the northeastern BoB is more affected by the salinity-related stratification effect, while the southern tip of India is more affected by the decrease in wind intensity (Bassinot et al., 2011) with enhanced stratification being potentially made stronger by an additional freshwater effect, although weaker than in the BoB. Ultimately, both climatic features (summer wind intensity and precipitation) are directly under the control of monsoon evolution resulting from the orbital forcing of low-latitude boreal summer insolation.

By contrast, higher intermediate Cd_w values from core MD77-191 associated with higher G. bulloides relative abundances and C_org from core SK237 GC04 during the 5.2–2.4 cal ka time interval could indicate enhanced productivity during the mid to late Holocene (Naik et al., 2017) (Fig. 5). To a lesser extent, this is also observed in the records from the northern BoB for the same time period. These changes are consistent with the weakened summer monsoon intensity, with less rainfall during the late Holocene, as observed in the BoB using core MD77-176 seawater δ¹⁸O and core SO188-342KL δD_Alk−ic records (Marzin et al., 2013; Contreras-Rosales et al., 2014; Fig. 6e–f). In addition, this is also strongly supported by the δ¹³C_wax records from the Lonar Lake on the Indian continent (Sarkar et al., 2015; Fig. 6d) and a progressive increase in monsoon summer winds to the south of India (Bassinot et al., 2011). These observations could also strongly support the hypothesis that the major control on surface productivity is linked to monsoon evolution in the BoB and at the southern tip of the Arabian Sea during the Holocene (Bassinot et al., 2011; Naik et al., 2017; Zhou et al., 2020).

5.4 Millennial-scale changes in intermediate-water circulation during the deglaciation

During the last deglaciation, short events have been recorded at the site of core MD77-191 during the 16–15.2 (HS1) and 12.6–11 cal ka (YD) time intervals (Fig. 5). The low Cd_w values in the MD77-191 record are coeval with reductions in C_org in core SK237 GC04 during the HS1 and YD periods (Fig. 5). According to previous studies, extremely high Cd_w values (>1 nmol kg⁻¹) were reported to have been associated with enhanced surface productivity (Bostock et al., 2010; Olsen et al., 2016). However, the range of values of intermediate Cd_w (0.58–0.85 nmol kg⁻¹, HS1; 0.5–0.8 nmol kg⁻¹, YD) from core MD77-191 during the last deglaciation is much lower compared with the Holocene Cd_w values (>1 nmol kg⁻¹) and thus may be associated with other processes such as a better ventilation, changes in the water mass source and/or depleted surface productivity (Fig. 6). Significant decreases in G. bulloides relative abundance of SK237 GC04 (Naik et al., 2017) and MD77-191 records were observed from the HS1 to B-A (Bassinot et al., 2011), and thereafter slight increases occurred in the YD (Fig. 5). These high values at both core sites during the HS1 and YD may indicate an enhanced surface productivity during these intervals (Fig. 5). This should have led to increased intermediate Cd_w and organic matter preservation under low oxygen concentration conditions during the HS1 and YD. However, despite a low resolution for the MD77-191 Cd_w record during the last deglaciation, we do not observe high values of intermediate Cd_w during the HS1 and YD (∼0.7 nmol kg⁻¹) compared with the late Holocene (∼1.59 nmol kg⁻¹), especially at 16.5–16 cal ka. Although we cannot fully discard the influence of surface productivity on the intermediate Cd_w in these time intervals, this apparent discrepancy seems to provide another piece of evidence for the influence of changes in water masses and/or ventilation during the HS1 and YD, in line with previous studies and proxies in the northern Indian Ocean (Bryan et al., 2010; Yu et al., 2018; Ma et al., 2019, 2020).

Moreover, an increase in benthic δ¹³C values is observed during the HS1 and YD in the northern Indian Ocean (e.g., Duplessy et al., 1984; Curry et al., 1988; Naqvi et al., 1994; Jung et al., 2009; Ma et al., 2019) (Fig. S3). The increase in the different benthic δ¹³C records during the HS1 and YD in the western Arabian Sea, Pacific Ocean and BoB is interpreted as reflecting the northward expansion of AAIW (Pahnke and Zahn, 2005; Jung et al., 2009; Ma et al., 2019) (Fig. S3). The decreased B-P age obtained from marine sediment cores from the Arabian Sea and the Bay of Bengal during the same intervals could confirm enhanced vertical mixing in the Southern Ocean (Bryan et al., 2010; Ma et al., 2019). The transition in the ε_Nd and Δ¹⁴C records during the deglaciation also indicates a strong northward penetration of AAIW within the North Atlantic and Bay of Bengal (e.g., Cao et al., 2007; Pahnke et al., 2008; Pena et al., 2013; Yu et al., 2018). In addition, during the HS1 and YD, a decrease in the [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] record from core MD77-191 also suggests the release of CO₂ from the deep ocean in the deglacial period through the expansion of AAIW (Ma et al., 2020). These time intervals are associated with better ventilation in the Southern Ocean (e.g., Anderson et al., 2009; Skinner et al., 2010), which led to enhanced vertical ventilation resulting in increased production of intermediate water masses (AAIW) (Anderson et al., 2009).

As mentioned before, previous studies have suggested an enhanced northward flow of southern sourced intermediate water mass AAIW, observed also in the Atlantic, Pacific and Indian oceans during the last deglaciation (e.g., Pahnke et al., 2008; Bryan et al., 2010; Poggemann et al., 2017; Yu et al., 2018; Ma et al., 2019, 2020), indicating that the source of intermediate water masses may be partly the same in these oceans. Thus, as the benthic δ¹³C values collected from the north Indian Ocean could better constrain the influence of AAIW in the two studied cores (Naqvi et al., 1994; Jung et al., 2009; Ma et al., 2019, 2020) (Fig. S3), we can also compare the range values of AAIW Cd_w from both studied cores with data from the Atlantic and Pacific oceans at intermediate water depth during the HS1 and YD (Cd_w, 0.3–0.9 nmol kg⁻¹; Umling et al., 2018; Valley et al., 2017). Thereafter, we could get the ranges of Cd_w–δ¹³C values of AAIW during these intervals, based on the benthic δ¹³C records in the Indian Ocean (Naqvi et al., 1994; Jung et al., 2009; Ma et al., 2019, 2020) as well as benthic Cd_w values from the Pacific and Atlantic oceans (Valley et al., 2017; Umling et al., 2018) at intermediate water depths (Fig. 7). Unfortunately, the resolution of both intermediate Cd_w and benthic δ¹³C from core MD77-176 (northeastern BoB) are very low for the HS1 and YD events, making it difficult to extract reliable information. Thus, we have decided to focus on the results from core MD77-191 (0.5–0.85 nmol kg⁻¹) during these two time intervals; these results are in good agreement with the collected dataset (Fig. 7). Thus, the benthic Cd_w results provide new evidence for tracking the northern flow of AAIW in the northern Indian Ocean, which increased during the HS1 and the YD.

Figure 7 Intermediate Cd_w versus benthic δ¹³C obtained from core MD77-191 located off the southern tip of India. The yellow shaded area represents the ranges of Cd_w–δ¹³C values of AAIW during the HS1 and YD, which were reconstructed in the Indian Ocean (benthic δ¹³C, Naqvi et al., 1994; Jung et al., 2009; Ma et al., 2019, 2020) and the Pacific and Atlantic oceans (benthic Cd_w, Valley et al., 2017; Umling et al., 2018) at intermediate water depths. The abbreviations are the same as in Fig. 2.

Taken together, Cd_w, the B-P age offset, benthic δ¹³C, ε_Nd and Δ¹⁴C records reported from the northern Indian Ocean all suggest strong upwelling and enhanced northern flow of AAIW from the Southern Ocean during the HS1 and the YD. Thus, the variations in these records can provide strong evidence for the hypothesis that Southern Ocean upwelling played a vital role in the increase in atmospheric CO₂ in the deglacial period (Anderson et al., 2009; Skinner et al., 2010, 2014). However, Kohfeld et al. (2005) suggested that although physical processes (such as ventilation) are involved in the glacial–interglacial atmospheric CO₂ change, the biological pump may also contribute nearly half of the observed changes in CO₂ during the glacial–interglacial transitions. As shown above, the HS1 event is characterized by reduced surface productivity, as revealed by the lower percentage values of G. bulloides in core MD77-191 (Bassinot et al., 2011) and by several studies of cores located in the eastern and western Arabian Sea within the Oxygen Minimum Zone (e.g., Schulz et al., 1998; Altabet et al., 2002; Ivanochko et al., 2005; Singh et al., 2006, 2011; Naik et al., 2017). This reduced productivity at a millennial timescale suggests that the entire biological factory was related to the reduced monsoon intensity during the North Atlantic Heinrich events (e.g., Singh et al., 2011; Naik et al., 2017). Thus, a weaker biological production could also have contributed to the two-step increase in atmospheric CO₂ during the last deglaciation, at least for the HS1 period.

6 Conclusions

Changes in benthic foraminiferal $\mathrm{Cd}/\mathrm{Ca}$ and assemblages were reconstructed on core MD77-191 (1254 m water depth) located off the southern tip of India, as well as on core MD77-176 (1375 m water depth) from the northern BoB, in order to reveal the evolution of intermediate-water circulation and paleo-nutrient changes in the northern Indian Ocean since the last deglaciation. We reconstructed seawater Cd_w concentration by converting H. elegans $\mathrm{Cd}/\mathrm{Ca}$. Benthic $\mathrm{Cd}/\mathrm{Ca}$ ratios are mainly influenced by changes in surface productivity and intermediate–bottom-water ventilation.

Results indicate that assemblages 2 and 3, reflecting high bottom-water oxygen conditions and a low flux of organic matter, dominated between 17 and 6 cal ka, corresponding to a poor-productivity time period. The typical late Holocene assemblage indicates a relatively low oxygen level and meso- to eutrophic deep-water conditions, associated with high surface productivity. The early Holocene (10–6 cal ka) corresponds to a low in productivity associated with depleted Cd_w in intermediate water. These observations seem to result from enhanced monsoon precipitation and increased river inputs from the Himalayan rivers, which led to more marked stratification in the BoB and a reduction in primary and export productivity. At the southern tip of India, the decrease in vertical mixing is also associated with a reduction in summer wind forcing resulting from the northward displacement of the ITCZ during summer (Bassinot et al., 2011). During the late Holocene (5.2–2.4 cal ka), the increased intermediate Cd_w concentrations of cores MD77-191 and MD77-176 indicate enhanced surface productivity in the southeastern Arabian Sea and in the northeastern BoB, corresponding to weakened monsoon intensity and rainfall, in agreement with other local records and reconstructions of the paleo-monsoon strength. Thus, our results clearly show the strong control of intermediate-water Cd_w during the Holocene by orbitally driven changes in summer monsoon productivity.

As far as millennial-scale variability is concerned, during the last deglaciation, decreased intermediate Cd_w concentrations during the HS1 and the YD are coeval with increased benthic δ¹³C, depletion in [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] and decreased B-P age offsets. These observations indicate that the low Cd_w values in intermediate water mainly resulted from the increased northward flow of AAIW during the HS1 and YD intervals. These signals also provide strong evidence for the important role of enhanced Southern Ocean ventilation in the CO₂ increase during the last deglaciation. The declined intermediate Cd_w obtained from the southeastern Arabian Sea (Core MD77-191), combined with the published eastern and western Arabian Sea paleo-productivity results, provides evidence for the important influence of decreased monsoon intensity at a millennial timescale during cold events in the North Atlantic region, associated with the increase in atmospheric CO₂ during the last deglaciation.