Salinity stratification controlled productivity variation over 300 ky in the Bay of Bengal.

The unique hydrographic setting of the Bay of Bengal (BoB) makes it an ideal tropical marine system to study the influence of regional and global forcings on productivity and [CO2aq] through the late quaternary. Enormous fresh water flux into the BoB and consequent salinity stratification significantly weaken the convective mixing and wind driven processes which are commonly responsible for transport of nutrients to the euphotic zone driving primary productivity. Here we present a high resolution organic carbon-CaCO3 MAR and δ13CTOC records for the last 300 ky from the BoB. The results show significant productivity variation at marine isotope sub-stages and millennial timescales. Colder sub-stages and stadials (Dansgard-Oeschger cycle) show a boost in productivity which may be attributed to thinning of low salinity cap, thereby facilitating efficient nutrient transport across the euphotic zone by the combination of wind driven processes (entrainment and upwelling), convective mixing and cold core eddies. The [CO2aq] was a net result of global pCO2 variation and regional processes. Our long term high-resolution data indicates a possibility of marked change in productivity/biogeochemistry of BOB in the future due to global warming, thus affecting the coastal economy.

Introduction

The Bay of Bengal, a tropical semi-enclosed basin in the northern Indian Ocean is the largest bay in the world, bordered by India, Sri Lanka, Bangladesh, Myanmar and the Andaman-Nicobar Islands. Ganga-Brahmaputra (G-B), Irrawaddy, Godavari, Mahanadi, Krishna and Kaveri rivers contribute 60% of the total freshwater received by the BoB of which the G-B river system contributes[1] 44%. Enormous riverine water flux (2.95 × 1012m3/yr)[1] and excess of precipitation over evaporation results in a stable water column (50–80 m) salinity stratification[2] in BoB, in contrast to other Indian Ocean regions. The surface salinity in BoB is lowest (~29 psu) above 20°N latitude and increases to 34 psu around 7°N latitude. The salinity gradient[2] from north to south decreases from 5.5 to 1. The salinity stratification results in a shallow mixed layer depth (MLD ~ 5 to 30 m) which shows marked regional and seasonal range[3]. The salinity stratification enhances the stability of the water column and prevents mixing with the underlying cooler waters leading to high sea surface temperatures (~28 °C) throughout the BoB[4]. The water column stability does not allow the prevailing wind (5–10 ms−1)[5] to disrupt the stratification[6] significantly except during cyclonic episodes. The salinity stratification affects the vertical distribution of heat in the near surface layers of BoB and can influence processes such as the active–break cycles of summer monsoons and development of regional tropical cyclones[7]. Owing to these enigmatic oceanographic characteristics, physical processes like convective mixing and wind driven processes (upwelling and nutrient entrainment) which are responsible for the transport of nutrients across the euphotic zone leading to primary productivity are significantly weaker in BoB compared to that of Arabian Sea[8]. As a result, maximum productivity in the open ocean BoB is associated with subsurface chlorophyll maxima (SCM ~40 to 90 mbsl)[2,5] in both the central and western bay throughout the year[2,9], whereas, surface productivity in BOB is limited to the coastal regions receiving nutrients along with perennial runoff from the rivers[9]. The depth of SCM is controlled by the vertical transport of nutrient by mesoscale cold core eddies (anticyclonic)[2,5,6,9,10]. The genesis of mesoscale eddies are attributed mainly to the interaction of Kelvin wave propagating along the coastal boundary with the northward flowing western boundary current and also by the breaking of Rossby waves propagated from the eastern boundary of BoB[10].

In view of this unique oceanographic setting and being one of the most vulnerable[11] regions in the world due to climate change, BoB is one of the hotspots for climate change research. Thus, it is pertinent to investigate how marine productivity and surface water dissolved CO2 ([CO2aq]) responded to highly variable past climatic conditions. Understanding the natural variability of the paleoceanic proxies at high temporal resolution will allow more accurate modeling of the consequences of future warming.

Study area

As part of this investigation, a giant Calypso piston corer was used on board ORV Marion Dufresne (MD-161) for the retrieval of the sediment core MD161-19 (core length: 39 m) off Mahanadi Basin (Lat.:18°59.1020; Long.:85°41′′.1669′′) in western BoB at a water depth of 1480 m (Fig. 1).

Figure 1

Location of cores MD161–19 and NGHP-19a (J 19) on the bathymetry map. The two core locations are indicated by thick blue arrows. The map is prepared using Generic Mapping tools software (GMT, version 4.5.1.15; www.gmt.soest.hawaii.edu/gmt4).

Results

The age-depth model (Supplementary Fig. 1 and Supplementary Table 1) for MD161–19 is based on the calibrated radiocarbon ages[12] and correlation of our δ18OG.ruber profile with the standard δ18O profile[13]. Eight marine isotope stages (MIS 1-8) representing glacial and interglacial changes are demarcated in Fig. 2. The linear sedimentation rate (LSR) varies from 1.5 to 38.6 cm/ky except in the last 1114 yr where the sedimentation rate reaches a peak value of 241 cm/ky. Total inorganic carbon (TIC) concentrations range from 0.01% to 5.6 wt% (Supplementary Table 2). The CaCO3 mass accumulation rate (MAR) reaches a peak value of 12.6 g/cm2/ky within the time window of 126885–145000 yr. Beyond this time window, the CaCO3 MAR ranges from 0.02 to 0.26 g/cm2/ky except during the period of 689–2000 yr where the MAR reaches a maximum of 0.66 g/cm2/ky. The Total organic carbon (TOC) value ranges from 0.15% to 2.9 wt% (Supplementary Table 2). The TOC MAR ranges from 0.01 to 0.5 g/cm2/ky barring the high flux of 2.0 g/cm2/ky within 689–2000 y. The δ13CTOC values range from −22.3 to −16.3‰ (VPDB).

Figure 2

Geochemical profile through core MD161–19. Blue sections indicate cold MIS and orange sections indicate warm MIS. (a) δ18O profile indicating marine isotope stage and sub-stages. Marine isotope stage (MIS) boundaries (1 to 8) and substages[13] are demarcated by dashed lines. (b) CaCO3 MAR. (c) TOC MAR and LSR. (d) δ13CTOC (‰VPDB). (e) pCO2 ppmv (composite Ice core) (f) Temporal variation [CO2aq].

Discussion

CaCO3 and TOC MARs (Fig. 2b,c) are considered here as proxies for past productivity variations in response to changes in regional forcing like surface water hydrography and nutrient supply in the BoB. The MAR deducts the siliciclastic dilution effect and indicates CaCO3 and TOC fluxes on the sea floor assuming no loss or gain of material in the sediment[14,15]. The colder sub-stages[13] (Fig. 2a) are associated with enhanced CaCO3 MAR, whereas, warmer sub-stages allied with diminished CaCO3 MAR. We attribute this marked temporal variations in CaCO3 MAR primarily to changes in calcite productivity[16,17] rather than differential preservation[18]. The core MD161–19 collected at a water depth of 1480 mbsl, is well above the reported depth of calcite saturation/lysocline (~3000 m) in the BoB[19]. Comparison of shallow (730–809 m) and deep (1727–2250 m) water CaCO3 particle flux data[20] from northern BoB (NBBT) show minor differences, suggesting lack of major calcite dissolution at water depths of ~2000 m. Based on foraminifera distribution in the surface sediments of BoB, ~2000 mbsl was suggested as the calcite lysocline depth[21]. Since, pre-industrial pCO2 was at least 100 ppmv less than the present day value, lysocline depth was definitely deeper than the present. Post depositional dissolution of calcite by anaerobic biogeochemical reactions like sulfate reduction is unlikely, since pore fluid chemistry data[12] of the current core shows a drop in pore-water calcium concentrations and increase in total alkalinity (TA) with depth below the sea floor. The increase in TA is attributed to organoclastic sulfate reduction and/or anaerobic oxidation of methane (AOM) which promotes calcite precipitation rather than dissolution[22]. However, anaerobic oxidation of pyrite by Fe3+ ion during deep burial diagenesis may result in a drop in pH leading to some calcite dissolution. Influence of such diagenetic processes cannot lead to the observed systematic temporal variation in calcite burial flux. The apparent increase in δ18OG.ruber values coupled with enhanced CaCO3 MAR and vice versa (Fig. 2a,b) and an overall positive statistical correlation between δ18OG.ruber and CaCO3 MAR (Supplementary Fig. 2a and 2b and Supplementary Table 3) suggests an underlying climatic control on productivity.

Total organic carbon and CaCO3 MARs show an overall similarity in temporal trends (except in MIS-2 and 4) which can be attributed to carbonate productivity variation coupled with enhancement of phytoplankton biomass mainly diatoms[23,24] in BoB[25]. Such a coeval pattern in productivity may be attributed to nutrient availability. The TOC MAR shows relatively stronger influence of sedimentation rate than CaCO3 MAR (Fig. 2b,c). This may be attributed to significantly lower %RSD of TOC wt% (~30.5%) compared to that of TIC wt % (105%). TOC MAR is also influenced by high sedimentation rate owing to both detrital dilution as well as enhanced preservation of organic carbon[15].

Here we establish a link between temporal variation in regional forcing and paleo-productivity in the BoB. Studies from BoB have shown a significant reduction in monsoonal intensity/fresh water flux during colder isotope sub-stages and opposite in warmer sub-stages[26]. This observation is also supported by multiple core data from BoB[27] suggesting enhanced riverine fresh water flux (due to increased monsoonal precipitation) in the northern BoB during early to mid-Holocene and diminished fresh water flux (reduced monsoonal precipitation) during the last glacial maxima (LGM). Overall diminished fresh water flux during colder and arid sub-stages caused thinning of low salinity cap leading to destabilization of the water column stratification which triggered enhanced nutrient entrainment by wind driven processes and convective mixing leading to enhanced productivity (enhanced CaCO3 MAR).

Exceptionally high CaCO3 productivity during colder sub-stage MIS 6a may be attributed to the sustained intensification of the physical forcings. Inferred[28] intensification of NE monsoon characterized by cold dry winds and drop in river discharge during the last glacial maxima (LGM) gives further credence to our hypothesis. Cold core eddy and wind driven shoaling of nitrate into the MLD and enhanced productivity has been reported during the winter monsoon in BOB[9]. On the other hand, the marked drop in CaCO3 MAR during warmer and humid sub-stages may be attributed to thickening of low salinity cap and stabilization of water column (shallow MLD). A similar relationship is also apparent at the stadial/interstadial time scale (Fig. 3). The interstadials (D-O events) commonly observed at high latitudes have also been identified in the BoB[26] and Arabian Sea[29]. The interstadials are supposed to have experienced relatively warm/ humid conditions and stronger monsoonal rainfall[26] in contrast to the cold/ arid and weak monsoonal conditions during the stadials. The sharp drop in CaCO3 MAR during the interstadials and enhancement during the stadials supports the influence of freshwater flux on productivity even at millennial scale. During the warmer events, enhanced glacial melt water[30] possibly contributed to the overall fresh water flux of the G-B river system for a short period of time.

Figure 3

Millennial scale changes in CaCO3 MAR. (a) δ18O ‰VPDB profile constructed from NGRIP (North Greenland Ice Core Project) data, Daansgard-Oeschger (D-O) events are numbered from 1 to 25[13]. The blue coloured zones indicate the colder events in between the D-O events. (b) Temporal variation in CaCO3 MAR.

The warmer/ humid and cold/arid time windows are also associated with relative increase and decrease in atmospheric pCO2 respectively (Fig. 2e). The near mirror image relation of δ13CTOC profile with the ice core pCO2 record (Fig. 2d,e and Supplementary Table 2) shows the influence of atmospheric pCO2 on the carbon isotope ratios of particulate organic matter. A positive statistical correlation between the temporal trends of δ13CTOC (−22.3 to −16.3‰ VPDB) and CaCO3 MAR (Supplementary Fig. 2a and 2b and supplementary text) indicates an overall underlying paleoclimatic control on both the parameters. The carbon isotopic composition of marine phytoplankton biomass during CO2 fixation via photosynthesis depends primarily on the ambient extracellular dissolved CO2 concentration ([CO2aq])[31-34] (Fig. 2f). However, the net cellular carbon isotope fractionation depends on phytoplankton species variation[35], solar irradiance and growth rate at seasonal scale[36]. Carbon isotope ratios of [CO2aq] normally range from ~ −9 to −11‰ VPDB[37]. Under [CO2aq] replete-conditions photosynthetic biomass show13C depletion, while, depletion of [CO2aq] leads to active-HCO3− (δ13C~1‰ VPDB) ion transport inside the cell via an active inorganic carbon concentrating mechanism (CCM)[33] and results in carbon isotopic enrichment of phytoplankton biomass[34]. The CCM pathway involves active uptake of HCO3− ions into the cell and subsequent conversion to CO2 (catalyzed by the carbonic anhydrase) and fixation by RubisCO. The calculated [CO2aq] values (~5 to 11.6 micromoles/kg) (Fig. 1f) using the equationare comparable to the estimated paleo marine [CO2aq] record[38]. The [CO2aq] profile (Fig. 1f) shows a gross similarity with the atmospheric pCO2 profile indicating dependence of [CO2aq] on atmospheric pCO2, on the other hand, observed non-synchronous variation in [CO2aq] and atmospheric pCO2 at finer time scale may be attributed to processes including CO2 draw down via productivity fluctuation[39,40] and/or influx of deep water CO2 through physical forcing. The calculated [CO2aq] values may also be affected by the presence of some terrestrial organic matter in the sediments. Variation in δ13C values (~0.4‰ VPDB) of ice core based atmospheric CO2 record[41] is minor compared to the fluctuation in δ13CTOC (~6‰ VPDB) reported in our work, however, this component may be considered for more precise calculation of [CO2 aq].

We conclude here that, in contrast to the present scenario, BoB experienced remarkable oscillations in productivity and surface water CO2 budget during the last 300 ky, controlled by variation in the intensity of global and regional physical processes. In view of the alarming influence of global warming on marine productivity[42], monsoonal variability[43] and ocean acidification[44], the high-resolution long term natural variability observed here will be useful in vulnerability modeling of BoB[11]. The paleo productivity data coupled with δ15N may also help in reconstructing paleoxygenation and denitrification[45] processes in BOB.

Method

Coring and sample preservation

A giant Calypso piston corer was used on board ORV Marion Dufresne (MD-161)[12] for the retrieval of the sediment core MD161–19 off Mahanadi Basin at a water depth of 1480 m (Fig. 1) at Lat: 18° 59.1092′′N Long: 85° 41.1669′′E. The sea bottom temperature at this location was 4.4 °C. The core was subsampled into 5 cm thick slabs. Aliquots for organic geochemistry was preserved at 2 °C to arrest further microbial activity. Samples were freeze dried at the earliest and stored at 2 °C (dark) for hydrocarbon extraction.

Separation of planktic foraminifera and oxygen isotope ratio measurement

200 dried and weighed aliquots of the samples of MD-161-19 were suspended in distilled water and gently sieved through a 63 µm mesh sieve. Approximately 15–20 clean tests of Globigerinoides ruber with the size ranging from 250 to 355 µm were picked from the oven dried > 63 µm fraction for carbon and oxygen stable isotope ratio measurement. Prior to the analyses, the foraminiferal tests were broken, cleaned in 5–10% H2O2 followed by ultra-sonication in distilled water and methanol to remove contaminants. Samples for oxygen and carbon stable isotope ratios were analyzed in a Kiel III carbonate preparation device interfaced with a Finnigan-MAT 252 isotope ratio mass spectrometer at the Department of Geological Sciences, University of Florida. This is a dual inlet system having a precision of 0.04‰ VPDB for δ13C and 0.08‰ VPDB for δ18O for calcite standard NBS-19. The results are presented in Supplementary Table 1.

TIC, TOC contents, and δ 13C measurement

600 samples were desalinated and powdered prior to compositional measurements. Total Inorganic Carbon (TIC) was determined by carbon coulometer (UIC-CM5130). The accuracy of TIC content of standard reference material (Ultrapure CaCO3 from Sigma-Aldrich) was 12.0 ± 0.25 wt%. Total carbon (TC) content was measured by the elemental analyzer (Thermo EA1112). Total organic Carbon (TOC) content was calculated by subtracting TIC from TC. 2, 4-DNP was used as a calibration standard for TC. Reproducibility for TC in NIST-SRM1944 sediment standard was found to be 4.4 ± 0.2 wt%. Carbon isotope measurement of total organic carbon (δ13CTOC) was carried out on decarbonated sediments. A Thermo-Finnegan Delta-V-Plus continuous flow isotope ratio mass spectrometer coupled to an elemental analyzer (Thermo EA1112) was used for C isotope ratio measurements. The external standard reproducibility calculated for δ13CTOC using IAEA-C3 cellulose standard was −24.7 ± 0.1‰ (VPDB). The results are presented in Supplementary Table 2.

Calculation of porosity and dry bulk density (DBD)

Dry bulk density for our core was calculated using the DBD-porosity relation for we have taken the help of porosity DBD relation obtained from core number NGHP-19[46]. MD161–19 and NGHP-19 are in close proximity and have similar lithology and porosity and profile for the top 40 m. The linear equation obtained between porosity and DBD data for NGHP-19 (DBD = [(porosity−98.86)/−35.75)] with an r2 of 0.99 is used to calculate the DBD for MD161–19. CaCO3 MAR was calculated as (DBD* Sed Rate* 8.33*TIC wt%/100) and TOC MAR was calculated as (DBD* Sed Rate* TOC wt%/100). For porosity measurements, a measured volume of sediment was dried at 105 °C. Moisture content in the sediment was calculated from the difference in wet and dry weight of sediments. Porosity was calculated as (volume of sediment pore-water/wet sediment volume weight) × 100. The results are presented in Supplementary Table 2.

Statistical analyses

Simple moving average filter: The moving average calculates the mean of the data in a particular period for a large dataset. It is simply used to reduce the random fluctuations generated in a big time series data. To calculate simple moving average (SMA), every data point is given equal weightage. The mathematical expression for estimating the SMA for a period of n in a time series data is aswhere, k is the position of the period and y is the variable to be filtered at the time t[47,48].

In Supplementary Fig. 2a, we have exhibited the raw data for CaCO3 MAR, δ18OG.ruber, and the δ13CTOC with their estimated SMA filtered data for the whole time series. It is to be carefully noted that we have lesser number of available data points for δ18OG.ruber with respect to the others. To calculate the Pearson’s correlation coefficient between CaCO3 MAR and δ18OG.ruber, we have truncated the CaCO3 MAR string as same as the size of the δ18OG.ruber data string. As a consequence of that, we have plotted the full CaCO3 MAR with the reduced one. In all the cases, the filtered lines can efficiently follow the trend of the raw data after removing the fluctuations.

Pearson’s correlation coefficient: To quantify the strength of association among two simultaneously evolving quantities, there are many statistical approaches like covariance, correlation, etc. We have opted the correlation coefficient as an index here[49,50]. The Pearson’s correlation coefficient (rxy) for n number of data points, can be written as

In the scatter plot (Supplementary Fig. 2b), we have plotted the δ18OG.ruber and the δ13CTOC raw data as well as the smoothened data with respect to CaCO3 MAR. We consider the axis of CaCO3 MAR in Log scale to cover the full range of CaCO3 MAR time series data points in a presentable form. The linear trends in the scattered points for both δ18OG.ruber and δ13CTOC with CaCO3 MAR encourages us to assess their mutual association through Pearson’s correlation coefficient (rxy). We have calculated rxy on the full data before and after using the SMA filter. The rxy values (Supplementary Fig. 2a) improves after removal of the local random fluctuations, but cannot affect the exact trend.

Supplementary information

Supplementary Data Table-1

Supplementary Data Table-2