Glacial–interglacial changes and Holocene variations in Arabian Sea denitrification

At present, the Arabian Sea has a permanent oxygen minimum zone (OMZ) at water depths between about 100 and 1200 m. Active denitrification in the upper part of the OMZ is recorded by enhanced δ¹⁵N values in the sediments. Sediment cores show a δ¹⁵N increase during the middle and late Holocene, which is contrary to the trend in the other two regions of water column denitrification in the eastern tropical North and South Pacific. We calculated composite sea surface temperature (SST) and δ¹⁵N ratios in time slices of 1000 years of the last 25 kyr to better understand the reasons for the establishment of the Arabian Sea OMZ and its response to changes in the Asian monsoon system. Low δ¹⁵N values of 4–7 ‰ during the last glacial maximum (LGM) and stadials (Younger Dryas and Heinrich events) suggest that denitrification was inactive or weak during Pleistocene cold phases, while warm interstadials (ISs) had elevated δ¹⁵N. Fast changes in upwelling intensities and OMZ ventilation from the Antarctic were responsible for these strong millennial-scale variations during the glacial. During the entire Holocene δ¹⁵N values > 6 ‰ indicate a relatively stable OMZ with enhanced denitrification. The OMZ develops parallel to the strengthening of the SW monsoon and monsoonal upwelling after the LGM. Despite the relatively stable climatic conditions of the Holocene, the δ¹⁵N records show regionally different trends in the Arabian Sea. In the upwelling areas in the western part of the basin, δ¹⁵N values are lower during the mid-Holocene (4.2–8.2 ka BP) compared to the late Holocene (< 4.2 ka BP) due to stronger ventilation of the OMZ during the period of the most intense southwest monsoonal upwelling. In contrast, δ¹⁵N values in the northern and eastern Arabian Sea rose during the last 8 kyr. The displacement of the core of the OMZ from the region of maximum productivity in the western Arabian Sea to its present position in the northeast was established during the middle and late Holocene. This was probably caused by (i) reduced ventilation due to a longer residence time of OMZ waters and (ii) augmented by rising oxygen consumption due to enhanced northeast-monsoon-driven biological productivity. This concurs with the results of the Kiel Climate Model, which show an increase in OMZ volume during the last 9 kyr related to the increasing age of the OMZ water mass.

1 Introduction

The marine nitrogen (N) cycle is highly dynamic due to the many chemical compounds of reactive N and their rapid transformation processes (Casciotti, 2016). Its feedback mechanisms are able to respond to external perturbations, possibly stabilizing the marine inventory of fixed N (Deutsch et al., 2004; Gruber, 2008; Sigman et al., 2009). The range of both oceanic N sources and sinks is still uncertain due to the poor data coverage of rate measurements and the large uncertainties of the water mass ages. The estimates of total N sources and sinks vary by factors of up to 4 and it has been debated whether the recent marine N cycle is in balance (Brandes and Devol, 2002; Codispoti, 2007; Codispoti et al., 2001; Gruber, 2008; Gruber and Sarmiento, 1997). Models have constrained the major oceanic N sinks (total water column and benthic denitrification) to 120–240 Tg N year⁻¹ and brought them close to equilibrium with estimates of diazotrophic dinitrogen (N₂) fixation, the main oceanic N source (Deutsch et al., 2004; DeVries et al., 2013). New measurements have at the same time led to higher global estimates of N₂ fixation (Grosskopf et al., 2012).

A period of fundamental change in oceanic N cycling (among other element cycles) occurred during the transition from the last glacial to the Holocene due to adjustment to changes in wind forcing, ocean circulation, sea level, and nutrient and trace metal supply from land (Deutsch et al., 2004; Eugster et al., 2013). The present equilibrium was probably attained only a few thousand years ago (Deutsch et al., 2004; Eugster et al., 2013). Understanding the response of the N cycle to this complex reorganization is important to facilitate our present understanding of N cycling on global and regional scales (Gruber and Galloway, 2008).

Sedimentary δ¹⁵N values integrate signals derived from N sources and the fractionation processes occurring during N cycling so that δ¹⁵N records have to be carefully deciphered (Altabet, 2006; Brahney et al., 2014; Nagel et al., 2013). Locally, eolian and riverine N supply can impact δ¹⁵N values in sediments (Kendall et al., 2007; Voss et al., 2006), but generally sedimentary δ¹⁵N values reflect the role of denitrification vs. N₂ fixation in ocean basins. Denitrification in the water column OMZ reduces nitrate in several steps to N₂. These reactions strongly discriminate against the heavy ¹⁵N isotopes so that the residual nitrate is isotopically enriched to δ¹⁵N values above the oceanic average of 5 ‰ (Brandes et al., 1998; Cline and Kaplan, 1975). Convective mixing and especially upwelling force nitrate-deficient water masses to the surface so that the enriched δ¹⁵N signal of nitrate is effectively transported into the euphotic zone. After assimilation into biomass by phytoplankton, ¹⁵N-enriched particulate matter sinks through the water column to the seafloor where the signal of denitrification and OMZ intensity is preserved in the sediments (Altabet et al., 1995; Gaye-Haake et al., 2005; Naqvi et al., 1998; Suthhof et al., 2001). The nitrogen deficit produced by denitrification can be counteracted by N₂ fixation from the atmosphere, which introduces nitrogen with a δ¹⁵N only slightly lower than the atmospheric value of 0 ‰, as the process is associated with little isotopic fractionation (Carpenter et al., 1997).

Sedimentary δ¹⁵N records show that during the glacial denitrification was less intense than today (Galbraith et al., 2013). Models suggest – albeit with many uncertainties and unknowns – that both denitrification and N₂ fixation were lower during the glacial (Deutsch et al., 2004; Eugster et al., 2013; Galbraith et al., 2013; Schmittner and Somes, 2016; Somes et al., 2017). However, due to a stronger reduction of denitrification than of N₂ fixation, total export production was higher and increased the glacial oceanic N inventory by 10–50 % over that of the Holocene, also enhancing the carbon storage in the ocean (Deutsch et al., 2004; Eugster et al., 2013; Schmittner and Somes, 2016; Somes et al., 2017). Distinct changes in sedimentary δ¹⁵N values during deglaciation are interpreted to reflect the major changes in the N inventory (e.g., Galbraith et al., 2013). The decreasing iron supply after the LGM (19–26.5 ka BP; Clark et al., 2009) may have significantly reduced N₂ fixation, leading to a rise of δ¹⁵N (Eugster et al., 2013). Enhanced upwelling at about 15 ka BP led to abrupt onsets or increases in denitrification in the eastern tropical North and South Pacific as well as in the Arabian Sea (Altabet et al., 1995; Ganeshram et al., 2002, 2000, 1995; Suthhof et al., 2001). The corresponding signal of enhanced δ¹⁵N values was dispersed and registered in many parts of the global ocean from the glacial to early Holocene and was followed by a smooth decrease in δ¹⁵N from enhanced N₂ fixation stimulated by the delayed increase in benthic denitrification caused by sea level rise (Deutsch et al., 2004; Galbraith et al., 2013; Ren et al., 2012). This sequence of events is very well recorded in cores from the east Pacific upwelling areas, but differs from the temporal pattern seen in sedimentary records from the Arabian Sea that show stable or increasing δ¹⁵N values in the Holocene (e.g., Galbraith et al., 2013).

Table 1 Station number, locations, water depth (m), data sources (references), and variables used: SST A (alkenone sea surface temperatures), δ¹⁵N ratios of total N, and SST Mg ∕ Ca (Mg ∕ Ca sea surface temperatures).

In order to (i) discern why N cycling in the Arabian Sea differs from the global trend and to (ii) better understand the response of the OMZ to changes in the monsoon system we present a summary of δ¹⁵N records from the Arabian Sea including two new records from the Oman upwelling area (Table 1; Fig. 1a; Sect. S1 in the Supplement). The records are from different areas and trace the regional history of mid-water oxygenation over the last 25 kyr. To relate the records of mid-water oxygenation to the history of southwest (SW) monsoon upwelling and northeast (NE) monsoon winter cooling, we compiled SST records from the literature and generated a new temperature reconstruction for the Oman upwelling area (Table 1; Fig. 1b, Sect. S2). Based on these integrated δ¹⁵N and SST records for different regions of the Arabian Sea we examine contrasts between glacial and Holocene conditions over the entire basin and contrasting regional evolution within the basin during the Holocene. Finally, we discuss our conclusions with the results of the global climate and ocean biogeochemistry model (KCM–PISCES) for the Holocene Arabian Sea.

2 Study area

The Arabian Sea hosts one of the most pronounced mid-water OMZ of the world's ocean and is a major oceanic N sink due to denitrification and anammox (Bulow et al., 2010; Codispoti et al., 2001; Gaye et al., 2013a; Jensen et al., 2011; Ward et al., 2009). Suboxic conditions between the thermocline and 1200 m are maintained by the balanced interaction of oxygen demand (organic matter degradation) and oxygen supply (ventilation; e.g., Olson et al., 1993; Sarma, 2002). The degradation of organic matter sinking out of the surface mixed layer consumes oxygen in the upper sub-thermocline water column. Primary productivity and particle flux in the Arabian Sea are highly seasonal and more than 50 % of annual particle fluxes occur during the summer season (Haake et al., 1993; Nair et al., 1989; Rixen et al., 1996), when strong SW monsoon winds induce upwelling of cold, nutrient-rich water masses along the coasts of Somalia and Oman (Fig. 2a). Upwelling ceases as changing wind patterns reverse surface circulation from clockwise during the SW monsoon (June–September) to anticlockwise during the NE monsoon (December–March; Schott et al., 2002). During the NE monsoon the temperature minimum and the productivity maximum occur in the northeastern Arabian Sea off Pakistan (Fig. 2b) caused by deep convection due to winter cooling (Rixen et al., 2005; Wiggert et al., 2005).

Figure 1 Stations of sediment cores for δ¹⁵N (a) and for SST (b) reconstructions with colors indicating surface sediment samples (purple) and cores from the Oman upwelling (dark blue), Somali upwelling (light blue), and the western (green), northern (yellow), eastern (red), and southeastern (orange) Arabian Sea.

Figure 2 SST in ^∘C from Jul–Sep (SW monsoon) (a) and Jan–Mar (NE monsoon) (b) from the World Ocean Atlas (Locarnini et al., 2013). Solid arrows indicate major wind directions and broken arrows indicate surface general currents; WICC is the West Indian Coastal Current.

Upper water masses (< 1500 m of water depth) in the Indian Ocean have a net westward circulation, while deep waters follow an eastward route as part of the Great Ocean Conveyer Belt connecting the Pacific and Atlantic Ocean (Broecker, 1991). However, this general direction comprises rather complicated routes and pathways (Durgadoo et al., 2017). Four water masses contribute to the subsurface waters of the OMZ: Arabian Sea High-Salinity Water (ASHSW), Persian Gulf Water (PGW), Red Sea Water (RSW), and Indian Ocean Central Water (ICW). The ASHSW forms during the NE monsoon due to enhanced evaporative cooling driven by dry air from the Himalayas (Prasad and Ikeda, 2002). During September–October the core of the ASHSW at a salinity > 36.5 psu and σ_θ=23.9 kg m⁻³ is found in the upper 100 m in the northern and eastern Arabian Sea (Kumar and Prasad, 1999). PGW forms a salinity maximum of 35.1–37.9 psu at σ_θ=26.6 kg m⁻³ within the core of the denitrification zone (100 to 350 m) and undergoes strong isopycnal mixing (Prasad et al., 2001) so that its proportion of the salinity maximum on the 26.6 σ_θ surface is less than 40 % (Morrison et al., 1998). During the SW monsoon PGW core salinity is additionally reduced due to the northward flowing Somali Current providing the less saline ICW (Prasad et al., 2001; You, 1998). The RSW forms an intermediate salinity maximum between 600 and 900 m and is characterized by salinities between 35.1 and 35.6 psu and σ_θ of 27.2 kg m⁻³ (Kumar and Prasad, 1999). PGW and RSW are saturated with oxygen from atmospheric contact shortly before their passage into the Arabian Sea through the Strait of Hormuz (50 m sill depth) and the Strait of Bab-el-Mandeb (137 m sill depth), respectively (Rohling and Zachariasse, 1996; Sarma, 2002). ICW combines Subantarctic Mode Water (SAMW) and Indonesian Intermediate Water (IIW) and enters the Arabian Sea from the southwest as part of the Somali Current (Resplandy et al., 2012; You, 1998). Intermediate water from the southern sources is originally rich in oxygen, but becomes increasingly oxygen depleted and nutrient rich on its path to the Arabian Sea owing to oxygen loss during the mineralization of sinking organic matter. Progressive oxygen loss is reflected by the observed pattern with higher oxygen concentrations in the NW basin than in the NE basin of the Arabian Sea (Morrison et al., 1999; Pichevin et al., 2007; Rixen et al., 2014). The deep water below about 1500 m is fed by circumpolar deep water (CDW) with σ_θ of 27.8 kg m⁻³ and a salinity of 34.8 psu (Bindoff and McDougall, 2000; Schott and McCreary, 2001).

The intensity of the OMZ and denitrification is seasonally variable. Oxygen concentrations in its core and its volume vary in response to the seasonality of ventilation, probably related to the seasonality of isopycnal mixing (Banse et al., 2014; Rixen et al., 2014). Models produce similar patterns with major ventilation from the south during the SW monsoon, while reversing currents, progressive oxygen consumption, and isopycnal mixing reduce oxygen concentrations in the entire basin during the winter monsoon (Resplandy et al., 2012; Rixen et al., 2014). Reoxygenation during the SW monsoon occurs via the invigorated Somali Current in the western Arabian Sea and a northward-flowing undercurrent below 150 m of water depth along the SW coast of India (Resplandy et al., 2012), which was found up to 20^∘ N (Shetye et al., 1990). This undercurrent reverses its direction semiannually, opposing the direction of the West Indian Coastal Current (WICC; Fig. 2) at the surface (Shetye et al., 1990). At present, the strongest denitrification prevails in the NE Arabian Sea although productivity and particle fluxes are highest in the western part of the basin (Gaye-Haake et al., 2005; Nair et al., 1989). Denitrification that reduces nitrate to nitrite and gaseous dinitrogen is triggered when oxygen concentrations fall below 4–5 µM O₂ (Cline and Richards, 1972; Devol, 1978). In general, oxygen-deficient conditions enable denitrification below 100 m of water depth in the Arabian Sea and active denitrification indicated by the accumulation of nitrite (Naqvi et al., 2008) was found between about 100 and 400 m of water depth (Gaye et al., 2013a; Martin and Casciotti, 2017). The intrusion of PGW that flows in a southward direction along the coast of Oman can occasionally supply oxygen and suppress denitrification (Morrison et al., 1998), as was the case during the late SW monsoon in 2007 between 250 and 400 m of water depth (Gaye et al., 2013a).

Paleoceanographic studies from the Arabian Sea report the existence of a pronounced OMZ and elevated denitrification during IS and interglacial stages, whereas the Arabian Sea OMZ was better ventilated and denitrification was suppressed during the LGM and stadials (Altabet et al., 1995; Higginson et al., 2004; Möbius et al., 2011; Pichevin et al., 2007). Many studies used productivity proxies and SST reconstructions, often in combination with denitrification proxies such as sedimentary δ¹⁵N values, to show that warm periods (IS and interglacials) were characterized by a stronger SW monsoon inducing upwelling and higher productivity than cold periods so that denitrification was switched on in the OMZ (Altabet et al., 1999, 2002; Pichevin et al., 2007; Reichart et al., 1997; Schulte et al., 1999b; Suthhof et al., 2001). After the transition from glacial to interglacial conditions with the warm and cold excursions during the Bølling–Allerød and Younger Dryas (YD; 11.7–12.9 ka BP), respectively, the Holocene was evidently a more stable period of permanent upwelling and denitrification (Böll et al., 2015; Overpeck et al., 1996; Pichevin et al., 2007; Tierney et al., 2016). There are, however, indications that millennial-scale climate oscillations similar to the North Atlantic cold periods detected by Bond et al. (1997) also occurred in the monsoon realm, but with reduced amplitude (Azharuddin et al., 2017). These Holocene cold periods were found to be characterized by reduced precipitation on land (Menzel et al., 2014) and reduced monsoonal upwelling in the Arabian Sea (Gupta et al., 2003). The volume and intensity of the mid-water OMZ appear to have oscillated related to SW monsoon strength, intensity of winter cooling by the NE monsoon, and changes in OMZ ventilation (Böll et al., 2015; Das et al., 2017; Pichevin et al., 2007). Thus, understanding Holocene OMZ dynamics is indispensable to evaluate the recently observed OMZ intensification in the Arabian Sea (Banse et al., 2014; Rixen et al., 2014).

Indian Ocean water masses and circulation were quite different at glacial conditions as proxy studies of benthic foraminifera indicate. Deep water was evidently less oxygenated than at present (Duplessy, 1982; Kallel et al., 1988; Schmiedl and Leuschner, 2005; Waelbroeck et al., 2006). This was reproduced by models that showed a generally more sluggish bottom water ventilation from the Antarctic (Rickaby and Elderfield, 2005) with reduced oxygen contents due to the increase in sea ice cover (Buchanan et al., 2016; Somes et al., 2017). A better ventilation of the upper water column in the glacial ocean was explained by the better oxygen solubility in the colder water (Somes et al., 2017). Studies from the southern Arabian Sea further suggest that there was much stronger formation of AAIW during the Heinrich events. This glacial AAIW (GAAIW) would ventilate intermediate waters (∼ 1600 m) in the Arabian Sea where AAIW is not detectable today (Jung et al., 2009).

3 Material and methods

3.1 Sample collection

The two new core records included in this summary of SST and δ¹⁵N records from the Arabian Sea are the gravity core SL163 merged with the multicore MC681 taken at the same location (21^∘55.97^′ N, 59^∘48.15^′ E; 650 m of water depth) and multicore MC680 (22^∘37.16^′ N, 59^∘41.50^′ E; 789 m of water depth). Cores were retrieved in 2007 during Meteor cruise M74/1b from the continental margin off northern Oman. MC680 and MC681 were sampled in 1 cm intervals and the first 400 cm of core SL163 were sampled in continuous 3 cm intervals. We analyzed alkenones in all sediment samples of SL163; δ¹⁵N was analyzed in all samples of MC681 and SL163 and in every second sample of MC680. All sediment samples were freeze-dried and homogenized prior to chemical treatment and analyses.

3.2 Analyses of the new cores SL163–MC681 and MC680

3.2.1 Stable nitrogen isotopes

The ratio of the two stable isotopes of N (¹⁵N ∕ ¹⁴N) is expressed as δ¹⁵N, which is the per mil deviation from the N-isotope composition of atmospheric N₂ (δ¹⁵N = 0 ‰):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}=\left[\right(R_{\mathrm{Sample}}-R_{\mathrm{Standard}})/R_{\mathrm{Standard}}]\times \mathrm{1000},\end{array}$$

where R_Sample is the ¹⁵N ∕ ¹⁴N ratio of the sample and R_Standard is the ¹⁵N ∕ ¹⁴N ratio of atmospheric N₂; δ¹⁵N values were determined using a Finnigan MAT 252 isotope ratio mass spectrometer after high-temperature flash combustion at 1100 ^∘C in a Carlo Erba NA 2500 elemental analyzer. Pure tank N₂ was calibrated against the International Atomic Energy Agency reference standards IAEA-N-1 and IAEA-N-2, which were, in addition to an internal sediment standard, also used as working standards. Replicate measurements of a reference standard resulted in an analytical precision better than 0.2 ‰. The mean standard deviation based on duplicate measurements of samples is 0.07 ‰.

3.2.2 Alkenones

Sample preparation and the detailed analytical procedure for alkenone identification are described in Böll et al. (2014). Purified lipid extracts of between 1.5 and 5 g of freeze-dried and homogenized sediment samples were analyzed for alkenone concentrations using an Agilent 6850 gas chromatograph (GC) equipped with a split-splitless inlet system, a silica column (30 m × 0.25 µm film thickness × 0.32 mm ID; HP-1; Agilent), and a flame ionization detector (310 ^∘C). Alkenone unsaturation ratios were translated into sea surface temperature using the core top calibration for the Indian Ocean of Sonzogni et al. (1997b):

$$\begin{array}{ll}\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{SST}=({\mathrm{U}}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}-\mathrm{0.043})/\mathrm{0.033}{\displaystyle }& {\displaystyle }\text{with}{\mathrm{U}}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}={\mathrm{C}}_{\mathrm{37}:\mathrm{2}}/({\mathrm{C}}_{\mathrm{37}:\mathrm{2}}+{\mathrm{C}}_{\mathrm{37}:\mathrm{3}}).\end{array}$$

All lipid extracts were analyzed twice, resulting in a mean standard deviation of 0.2 ^∘C. The mean standard deviation of estimated SST based on replicate extraction and measurement of a working sediment standard is 0.5 ^∘C.

3.2.3 Sediment core age models

The age model for SL163 was published by Munz et al. (2017) and is based on 15 accelerator mass spectrometry (AMS) ¹⁴C datings. In this study the upper 50 cm were taken from multicore MC681. Four additional datings and visual matching of variations in element concentrations were used to correlate SL163 and MC681 in their overlapping parts (Table S3 and Fig. S3 in the Supplement).

The age model of core MC680 is based on four accelerator mass spectrometry (AMS) ¹⁴C datings from different core depths measured at Beta Analytics in Miami, FL (see Table S4 and Fig. S4). Calibration and reservoir age correction were done in the same way as for SL163 (Munz et al., 2017). Both cores have a conspicuous sedimentation hiatus around 5700 years BP. In core SL163 the hiatus was positioned at 57 cm; in MC680 the hiatus was at 37 cm based on a change in facies from olive foraminiferal nanofossil ooze below to olive brown organic-rich nanofossil silty clay above.

3.3 Integration and averaging of SST and δ¹⁵N reconstructions

The temperature and δ¹⁵N curves of most cores used here were taken from the literature (Table 1) except those for the new records of cores SL163–MC681 and MC680 presented for the first time in this paper (see above). All original data and all calculations are presented as Tables S1 and S2 in the Supplement.

In our compilation, temperature reconstructions from the eastern and southeastern Arabian Sea are based on Mg ∕ Ca ratios of planktonic foraminifera (see methods in, e.g., Govil and Naidu, 2010; Mahesh and Banakar, 2014; Saraswat et al., 2013; Tiwari et al., 2015) except core MD90963, which has alkenone temperatures (Rostek et al., 1997). All other records are alkenone temperatures calculated with the core top calibration of Sonzogni et al. (1997b). Using the published age models, we averaged the temperature records available from the northern, western, eastern, and southeastern Arabian Sea and for the Oman and Somali upwelling systems (Fig. 1b). Composites are based on two to five different core records. The data were binned in time slices of 1000 years for each individual sediment core. Next, all time slices of an age interval in a defined study area were averaged. The standard deviations of the calculated average SST curves rarely exceed the analytical precision of 0.5 ^∘C of alkenone-based temperature reconstruction.

Temperature reconstructions based on different methods may differ as proxies may be seasonally biased or impacted by dissolution or diagenesis (Huguet et al., 2006; Regenberg et al., 2014). TEX₈₆- and alkenone-based SST reconstructions of cores NIOP905 and SO42-KL74 were shown to differ in magnitude and phase as TEX₈₆ temperature reconstructions seem to have a SW monsoon bias (Huguet et al., 2006). Munz et al. (2015) showed that winter temperatures derived from planktic foraminiferal assemblages had stronger amplitudes than the alkenone-based annual average SST reconstructions in SO130-275KL of Böll et al. (2014). Two of the records used in this study have both, alkenone and Mg ∕ Ca temperature. In core P178-15P from the Gulf of Aden, alkenone and Mg ∕ Ca temperatures have uniform trends and are significantly correlated (P<0.05) with a slope of 1, and Mg ∕ Ca values are on average 0.5 ^∘C lower than alkenone temperatures (Tierney et al., 2016). A comparison of alkenone temperatures of Huguet et al. (2006) and Mg ∕ Ca temperatures of Anand et al. (2008) from core NIOP905 also shows that Mg ∕ Ca temperatures are lower. In contrast, the alkenone temperatures of core MD90963 (Rostek et al., 1997) are about 1 ^∘C lower than the two Mg ∕ Ca temperature records of nearby cores SK129-CR04 (Mahesh and Banakar, 2014) and SK157-4 (Saraswat et al., 2005) in the southeastern Arabian Sea. As both Mg ∕ Ca- and alkenone-based temperature reconstructions are calibrated with annual average surface layer temperatures (Regenberg et al., 2014; Sonzogni et al., 1997a, b) and as we can identify no trend in our comparison of results of the two methods in Arabian Sea sediment cores, we have compiled SST reconstructions of both methods.

The average δ¹⁵N values were calculated per time slice in a similar way as SST curves and averaged for the same areas (Fig. 1a). Before averaging the results of all curves of the selected areas, δ¹⁵N values were normalized to the average δ¹⁵N value of the respective core (Table S1). Some records were too short to use their average δ¹⁵N values as they did not cover the main δ¹⁵N shift from the glacial to Holocene. In these cases the normalization was done with the average value of a nearby core with δ¹⁵N in the same range. This procedure was carried out as δ¹⁵N in sediments is impacted by several factors in addition to the δ¹⁵N of nitrate upwelled or mixed from subsurface waters. Nitrogen fixation or an allochthonous supply of nitrogen from rivers or the atmosphere can reduce δ¹⁵N in particulate matter (Agnihotri et al., 2011, 2008, 2009; Lückge et al., 2012; Montoya and Voss, 2006). Upwelling from different water depths and the incomplete utilization of nitrate supplied by upwelling may furthermore lead to a gradient with increasing δ¹⁵N values offshore of the upwelling areas (Naqvi et al., 2003). Diagenesis increases δ¹⁵N values in the deep Arabian Sea by up to 3 ‰ in distal sediments (Gaye-Haake et al., 2005; Möbius et al., 2011). The normalization procedure makes the relative changes in δ¹⁵N comparable within each area despite differences in the diagenetic imprint and in δ¹⁵N sources so that relative changes may be interpreted with respect to the relative intensity of denitrification. Average δ¹⁵N curves of normalized values have a standard deviation of up to 1.5 ‰ with most values far below 1 ‰. The standard deviation is generally largest during deglaciation when δ¹⁵N changed rapidly. The curves represent averages of four to seven individual records except for the Somali upwelling system for which only two records were found. For the construction of the present δ¹⁵N chart, the results from surface samples published by Gaye-Haake et al. (2005) were included (Fig. 1a).

3.4 Climate and biogeochemistry model

We use results from an experiment with a global coupled atmosphere–ocean–sea ice model (the Kiel Climate Model, KCM; Park and Latif, 2009; Park et al., 2009) consisting of ECHAM5 (Roeckner et al., 2003) and NEMO (Madec et al., 2008) to force a global model of the marine biogeochemistry (PISCES; Aumont et al., 2003) in offline mode. KCM has been used for time-slice simulations of the preindustrial and the mid-Holocene climate (Schneider et al., 2010; Khon et al., 2010, 2012; Salau et al., 2012; Jin et al., 2014). Here, orbital parameters (eccentricity, obliquity, and precession) accelerated 10 times were varied transiently as forcing according to the equations of Berger (1978). The greenhouse gas concentrations follow the standard Paleoclimate Modelling Intercomparison Project Phase III (PMIP3) protocol (Braconnot et al., 2012) based on Indermühle et al. (1999). Changes in the ice sheets are neglected.

The ocean component (OPA9) uses a tripolar grid with a 2^∘ zonal resolution, and a meridional resolution varying from 0.5^∘ at the Equator to 2^∘ × cos (latitude) polewards of 20^∘. The water column is divided into 31 layers, with 20 layers in the upper 500 m (known as the ORCA2 configuration). ECHAM5, the atmospheric component of KCM, is run in T31 resolution with 19 layers corresponding to a grid cell size of about 400 km × 400 km. PISCES (Aumont et al., 2003) simulates the marine biogeochemistry including processes that determine dissolved oxygen concentrations based on the oceanic circulation as provided by NEMO (Madec et al., 2008) and an NPZD-type (nutrient–phytoplankton–zooplankton–detritus) description of the marine ecosystem. NEMO–PISCES in the ORCA2 configurations has been used to study monsoon and biological production interconnections in a recent study by Le Mézo et al. (2017).

Here we restrict the description of PISCES to the processes relevant for the oxygen concentration. Sources of oxygen are gas exchange with the atmosphere at the surface and biological production in the euphotic zone. The production of two phytoplankton groups (representing nanophytoplankton and diatoms) is simulated based on temperature, the availability of light, and the nutrients P, N (both as nitrate and ammonium), Si (for diatoms), and the micronutrient Fe. There are three nonliving components of organic carbon in PISCES: semi-labile dissolved organic carbon (DOC) and large and small particulate organic carbon (POC), which are fueled by mortality, aggregation, fecal pellet production, and grazing. In deviation to the standard PISCES setup, small POC sinks to the seafloor with a constant settling velocity of 2 m day⁻¹, while large POC settling is simulated depending on the calcite and opal ballast effect following Gehlen et al. (2006). Oxygen loss occurs through the respiration of organic matter in the entire water column. The respiration rate depends on temperature with a Q₁₀ of 1.8 and on the oxygen level, with a reduced rate for O₂ concentrations below 6 µmol L⁻¹. We also use an idealized age tracer that is set to zero in the surface layer and increases with time elsewhere. Advection and mixing are also applied to the age tracer. This tracer gives an indication of the subsurface circulation strength. Here we use it to analyze the change in the average age of the water masses in Arabian Sea OMZ over time.

Figure 3 Millennial regional averages of SST (^∘C) and ±1σ standard deviation in the northern, western, eastern, and southeastern Arabian Sea and the Oman and Somali upwelling areas of the last 25 kyr. Regions are indicated in Fig. 1. Times of interstadial 2 (IS2), Heinrich Event 1 (H1), and Younger Dryas (YD) are indicated by grey bars. The last glacial maximum (LGM), early Holocene (EH), mid-Holocene (MH), and late Holocene (LH) are also marked. Lines mark the beginning of MH and LH.

We do not attempt to provide a full model analysis of the Arabian Sea OMZ in this paper, but will mainly use it as an additional tool to estimate the most likely causes for the sediment core analyses of OMZ intensity changes. Some basic features of the simulated OMZ can be found in the Supplement (Figs. S5–S7).

Figure 4 Annual SST distribution in ^∘C from the World Ocean Atlas (Locarnini et al., 2013) (a); alkenone- and Mg ∕ Ca-derived SST reconstruction for the time slice from 17–18 ka BP (b) from cores shown in Fig. 1.

4 Results

4.1 Temperature reconstruction

All temperature reconstructions indicate lower SST during the Pleistocene compared to the Holocene (Fig. 3). Warming started at about 16–17 ka BP during the period defined as Termination 1 (Stern and Lisiecki, 2014) except for the southeastern region where SSTs rise at about 22 ka BP in response to rising summer insolation over the Northern Hemisphere (Berger and Loutre, 1991). The largest SST increases from the glacial to the Holocene can be observed in the northern (4 ^∘C) and the eastern (3 ^∘C) Arabian Sea. The increase is about 2.5 ^∘C in the Oman upwelling area and less than 2 ^∘C in the open western Arabian Sea, the Somali upwelling, and the southeastern Arabian Sea south of India. Some small-scale temperature variabilities exceeding the analytical error of 0.5 ^∘C are visible. There is an increase of different amplitude during the warm IS 2 (∼ 23.4 ka BP) and a small temperature drop during the YD in the available records of higher resolution (Fig. S8). In the average curves (Fig. 3) the YD is visible only in the east and west.

Figure 5 The δ¹⁵N values in ‰ in recent surface sediments (a); δ¹⁵N in sediments at 17–18 ka BP (b) from surface and core locations shown in Fig. 1.

In order to compare the modern and glacial SST distributions (Fig. 4a, b) we plotted the SST map from the World Ocean Atlas (Fig. 4a; Locarnini et al., 2013) and the time slice at 17–18 ka BP from the core records (Fig. 4b). This time slice is neither an IS nor a Heinrich event and has the best data density of the glacial (see Sect. S2). There is a change in the SST pattern in the basin between glacial conditions and the Holocene. During the last glacial, the SST minimum was situated in the northern Arabian Sea and Oman upwelling and there was generally a north–south temperature increase (Figs. 3; 4b). During the Holocene the SST pattern deviates from this north–south increase (Figs. 3; 4a): (i) SSTs in the Oman and Somali upwelling areas are lower than northern Arabian Sea temperatures, and (ii) SSTs in the eastern and southeastern Arabian Sea are high and in the same range. Small drops in SST occur in some of the curves at about 9 and 4–5 ka BP (Fig. 3).

4.2 Patterns of δ¹⁵N

The absolute δ¹⁵N values in surface sediments in the present Arabian Sea are elevated with values between 6 and > 12 ‰ compared with those of the last glacial with values between 3.5 and 7 ‰ (Fig. 5a, b). Holocene δ¹⁵N values are highest in the central part of the basin and in the Oman upwelling area and lower in most shelf and slope sediments outside upwelling areas (Figs. 5a and 6). Glacial shelf and slope sediments have δ¹⁵N values below 6 ‰. Similar to the present situation δ¹⁵N increased towards the center of the basin. However, there are no glacial δ¹⁵N records from the deepest part of the central Arabian Sea (Figs. 1a; 5b).

The δ¹⁵N values increase between 16 and 14 ka BP in all sectors except in the eastern Arabian Sea, where the increase occurs at about 8 ka BP (Fig. 6). The normalized highest relative increase in δ¹⁵N values by about 3.5 ‰ is observed in the northern Arabian Sea. All other normalized δ¹⁵N curves increase by ≤ 2 ‰. Most integrated curves show a relative minimum during the YD when δ¹⁵N almost returned to the low glacial values. The general pattern of the δ¹⁵N curves is thus similar to the GISP ice core δ¹⁸O record during the glacial and deglaciations (Fig. 6). The three δ¹⁵N curves of high resolution (Fig. S9) follow the GISP core with distinct minima during Heinrich Event 1 (H1) and the YD and a maximum during IS1, whereas IS2 maxima are found neither in the Somali upwelling area nor in the eastern Arabian Sea (Fig. 6). The Holocene δ¹⁵N patterns differ across the basin. An early Holocene (> 8.2 ka BP; Walker et al., 2012) maximum is observed in the open western Arabian Sea including the upwelling areas, whereas a late Holocene (< 4.2 ka BP; Walker et al., 2012) maximum is visible in the northern and eastern part of the Arabian Sea. An early and late Holocene δ¹⁵N peak occurs in the Oman and Somali upwelling areas.

5 Discussion

5.1 Temperature differences between glacial and the Holocene

The temperature rise from the LGM to the Holocene in the northern and eastern coastal regions of the Arabian Sea of 3–4 ^∘C is 1–2 ^∘C larger than modeled for the tropical ocean (Annan and Hargreaves, 2013; Hopcroft and Valdes, 2015; Jansen et al., 2007). This may be induced by the much lower glacial land temperatures of central Asia (Annan and Hargreaves, 2013), which weakened the SW and strengthened the NE monsoon compared to the Holocene (Duplessy, 1982). Changes in annual average temperatures in the northern Arabian Sea were shown to be determined mainly by the intensity of winter cooling and the resulting deeper thermohaline mixing; they thus added to the cooling induced by lower insolation during the glacial (Böll et al., 2014, 2015; Reichart et al., 2004).

Figure 6 The δ¹⁸O values in ‰ from the GISP2 ice core (Grootes and Stuiver, 1997) and sea level (m above NN) reconstruction from the Red Sea (Siddall et al., 2003) (a) compared with millennial regional averages of normalized and averaged δ¹⁵N (‰) values from the northern (b), eastern (c), and western (d) Arabian Sea, the Oman (e), and Somali (f) upwelling area. Total organic carbon mass accumulation in the Arabian Sea (TOC MAR in Tg C year⁻¹; data from Cartapanis et al., 2016, and Tim Rixen, unpublished) and insolation at 30^∘ N in W m⁻² (Berger and Loutre, 1991) (g) during the last 25 kyr. Error bars denote ±1σ. Grey bars and abbreviations as in Fig. 3; interstadial 1 (IS1).

The observed regional differences in temperature rise from the LGM to the Holocene (Fig. 3) led to a change in the general SST pattern (Fig. 4a, b). The SW monsoon SST pattern in the modern Arabian Sea with its NW–SE-oriented gradient (Fig. 4a) is strongly modulated by upwelling off Oman and Somalia and the inflow of warm and low-saline surface water from the Bay of Bengal via the WICC (Vijith et al., 2016) (Fig. 2a). The WICC inflow is fed by the North Equatorial Current and starts in the post-SW-monsoon period, probably forced by local winds around the southern tip of India (Suresh et al., 2016). It is related to prevailing sea level height difference between the Arabian Sea and Bay of Bengal, which is due to the enhanced precipitation and river discharge to the bay (Shankar and Shetye, 2001). A reason for the more latitudinal gradient of glacial isotherms (Fig. 4b) was a strengthened NE and a weakened SW monsoon so that winter cooling in the northern Arabian Sea was much stronger (Reichart et al., 2004) and upwelling was reduced or even inactive during the glacial; the cold water source in the western Arabian Sea was strongly reduced (Böll et al., 2015; Duplessy, 1982). In addition, salinity reconstructions indicate that there was less advection of low-salinity, warm surface waters by the WICC into the eastern Arabian Sea from the Bay of Bengal, probably due to the low glacial precipitation and river runoff (Mahesh and Banakar, 2014).

Glacial SSTs off Somalia are in a similar range as in the western and eastern Arabian Sea and almost 2 ^∘C higher than off Oman. This suggests that upwelling off Somalia was weaker than off Oman or even shut down. At present, SW monsoon upwelling is driven by the positive wind stress curl induced by the Findlater Jet – a low-level, cross-equatorial jet stream recurring during each SW monsoon over eastern Africa and the western Indian Ocean (Brock et al., 1991; Findlater, 1977). The strength of the Findlater Jet is directly coupled with the moisture transport to the Indian monsoon region (Fallah et al., 2016). Precipitation on land, such as that from the All-India Rainfall Series, is thus used as a measure for SW monsoon strength (Mooley and Parthasarathy, 1984). However, even at present there is no straightforward coupling between high rainfall on land and low SST (Levine and Turner, 2012) and thus also no direct correlation with productivity. Further, sediment proxies indicate and modeling studies suggest that the position of the Findlater Jet changed with monsoon intensity and this could decouple SST, productivity, and monsoon strength (Anderson and Prell, 1992; Le Mézo et al., 2017; Sirocko et al., 1991). The Himalayan ice shield during glacial conditions led not only to a reduced temperature gradient between land and sea so that the Findlater Jet weakened, but also to an eastward shift of the continental low-pressure cell so that the Jet had a more latitudinal orientation (Le Mézo et al., 2017; Sirocko et al., 1991). Our data suggest that its glacial position was not favorable for upwelling off Somalia.

Figure 7 Millennial regional-averaged SST from the northern Arabian Sea minus Oman-upwelling-averaged SST in ^∘C (black line) and northern Arabian Sea minus Somali-upwelling-averaged SST (red line; regions are indicated in Fig. 1) compared with reconstructions of the mean effective moisture in southern and central Asia (blue line) and the Indian monsoon region (green line) from continental archives (Herzschuh, 2006). Grey bars and abbreviations as in Fig. 3.

The SST differences between the northern Arabian Sea on the one hand and the more southern Oman and Somali upwelling areas on the other hand can be used as upwelling indices as they show deviations from the insolation-driven, southward temperature increase (Böll et al., 2015). Enhanced upwelling is indicated by a positive or rising index as it shows lower temperatures in the more southern upwelling areas compared to the northern Arabian Sea. In Fig. 7 we compare it with the index of effective moisture calculated from a large number of lake, peat, loess, and river records from the Asian continent by Herzschuh (2006). Peaks in the upwelling indices at 22–23 ka BP suggest that upwelling prevailed during IS2. During this short warm interval the upwelling was enhanced off Oman and became active for about 2 millennia off Somalia. The Somali upwelling notably started at about 16 ka BP, which could reflect a shift of the Findlater Jet to a position more parallel to the western margin of the Arabian Sea and thus favorable for upwelling in both upwelling areas. The temperature minimum during the YD, which is quite pronounced in SST records of high resolution (Böll et al., 2015; Saraswat et al., 2013; Schulte and Muller, 2001; Tierney et al., 2016; see Fig. S8), is reflected as a minimum in the upwelling indices at 11–13 ka BP (Fig. 7). Both the Oman and Somali upwelling indices increased during further warming after the YD at 11 ka BP. The moisture index drops after the early Holocene to the present in parallel with temperature records from the continent (Herzschuh, 2006; Marcott et al., 2013; Peterse et al., 2014), but upwelling indices remain on about the same level or even increase during the mid-Holocene (Fig. 7). We surmise that this is a signal of a mid-Holocene shift in the position of the Findlater Jet. During the warmest period in the early Holocene the Findlater Jet reached its position closest to the coast so that the Oman upwelling may have been restricted to a much smaller area (Le Mézo et al., 2017). In the mid-Holocene the Findlater Jet shifted offshore and upwelling remained high, covering a larger area so that SST minima prevailed and productivity was enhanced (Le Mézo et al., 2017).

Holocene temperature minima do not coincide in the basin and may have different causes (Fig. 3). SST minima near coastal stations in the upwelling centers during the early Holocene climatic optimum could indicate enhanced upwelling (Böll et al., 2015). The second Holocene temperature minimum at around 4–5 ka BP in the eastern part of the basin coincides with a severe drought on the Indian peninsula (Prasad et al., 2014) and a colder climate in other terrestrial climate records from central Asia (Berkelhammer et al., 2012; Hong et al., 2003; Ponton et al., 2012). The SST minimum may thus be due to cooling by a strengthened NE monsoon due to colder winters rather than to warming-related enhanced upwelling.

5.2 Nitrogen cycling in the glacial

At present, nitrate reduction between 100 and 400 m of water depth leaves residual nitrate with δ¹⁵N values up to > 20 ‰, and upwelling can transport enriched nitrate from 250 to 300 m of water depth into surface waters in the western Arabian Sea upwelling areas (Gaye et al., 2013a, b; Yoshinari et al., 1997). Therefore, near shore sediments from the upwelling area off Oman have δ¹⁵N elevated to > 10 ‰ (Fig. 5a), and δ¹⁵N values in all other recent sediments collected at water depths of < 1000 m, i.e., at depths where the diagenetic effect on sedimentary δ¹⁵N is small or negligible (Altabet and Francois, 1994; Gaye-Haake et al., 2005), are between 6 and 8 ‰. This is identical to the signal of sub-thermocline nitrate, which feeds productivity primarily via seasonal deep mixing outside the upwelling areas (Gaye et al., 2013a, b). The δ¹⁵N values > 11 ‰ in the central part of the basin are a result of (i) the offshore advection of ¹⁵N-enriched nitrate from upwelling areas where nitrate is not completely utilized (Naqvi et al., 2003) and (ii) an early diagenetic increase in δ¹⁵N in deep sea sediments (Gaye-Haake et al., 2005; Möbius et al., 2011).

The salient millennial-scale oscillation of the Pleistocene δ¹⁵N records calls for a strong mechanism linking OMZ intensity with Northern Hemisphere climate oscillations. The low δ¹⁵N values between 4 and 6 ‰ from < 1000 m of water depth from the time slice 17–18 ka BP (Fig. 5b) and the LGM (Fig. 6; Table S1) suggest that denitrification was much reduced or absent. The δ¹⁵N values up to 7 ‰ during IS1 and IS2 in the entire basin except in the eastern Arabian Sea and Somali upwelling area (Table S1 and Fig. S9) indicate that denitrification was enhanced but restricted to the northern and northwestern part of the basin during the IS. The δ¹⁵N minima of the YD and H1 are found in all records of high resolution and the average curves (Figs. 6 and S9). The oscillations in the δ¹⁵N records thus follow the primary productivity, which was enhanced during the warm phases due to stronger upwelling and reduced during the LGM, YD and Heinrich events (Leuschner and Sirocko, 2003; Reichart et al., 1997; Schulte et al., 1999a; Suthhof et al., 2001). Exceptions from this productivity pattern were reported from the eastern Arabian Sea, where the LGM had enhanced productivity at some locations (Naik et al., 2017), and from the NE-monsoon-dominated northern Arabian Sea where productivity was enhanced during the YD and Heinrich events due to the strong winter cooling (Reichart et al., 2004).

The OMZ sediments from the northern Arabian Sea were indistinctly laminated during most of the last glacial, indicating suboxic conditions, but were clearly bioturbated, indicating fully oxic conditions only during the YD and Heinrich events (Suthhof et al., 2001). Aragonite preservations and δ¹³C of benthic foraminifera further suggest that enhanced formation of GAAIW ventilated the lower OMZ from about 800 to 1800 m, especially during stadials (Böning and Bard, 2009; Jung et al., 2009; Naidu et al., 2014; Schmiedl and Leuschner, 2005). A similar increase in GAAIW formation was observed in the Atlantic and Pacific and was explained by the strong reduction of NADW formation and a simultaneous strengthening of SAMW and GAAIW formation due to the reduced salinity and density of North Atlantic surface water (Rickaby and Elderfield, 2005; Ronge et al., 2015). A complete breakdown of the Atlantic meridional overturning circulation (AMOC) during the stadials led to the strongest GAAIW production (Rickaby and Elderfield, 2005) and thus to the observed complete OMZ ventilation even in the northern Arabian Sea (Suthhof et al., 2001). The planktonic nitrate source is at 250–300 m depths in upwelling areas and at even shallower sub-thermocline depths outside upwelling areas. The water masses feeding nitrate to the surface are thus the ASHSW and ICW. ASHSW formation was possibly stronger when the NE monsoon was stronger and the climate was more arid so that it better ventilated the OMZ from above. The enhanced formation of SAMW also contributed to better ventilation. At present, SAMW is the major oxygen source to the Arabian Sea OMZ, while PGW, RSW, and IIW are only small contributors of oxygen (Fine et al., 2008; You, 1998). During glacial conditions, the increased SAMW production occurred further north due to the northward shift of the subpolar front similar to the GAAIW (Rickaby and Elderfield, 2005) and better ventilated the upper OMZ (Böning and Bard, 2009). It carried more oxygen due to less mixing with IIW and RSW, and an accelerated circulation (Böning and Bard, 2009) led to a lower residence time in the Arabian Sea. Moreover, it is quite possible that glacial SAMW carried fewer preformed nutrients due to the better nutrient utilization related to increased eolian supply of phosphate and trace metals to the region of SAMW formation (Somes et al., 2017). The lower amount of preformed nutrients further reduced productivity in the Arabian Sea. Better ventilation and reduced upwelling of nutrient-poorer water thus coincided during stadials and explain the complete oxygenation. The observed suboxic conditions discernible from laminated sediments in the OMZ during normal glacial conditions (Suthhof et al., 2001) did not produce enhanced δ¹⁵N signals in the sediments. It is feasible that the oxygen concentrations did not drop below the threshold of denitrification, but it is also possible that conditions in the Arabian Sea were comparable to those in the present Bay of Bengal where nitrate reduction and denitrification occur locally at a low level, but the enriched nitrate is not transported into surface waters due to stratification (Bristow et al., 2017). Enhanced N fixation has been suggested as an alternative reason for the low δ¹⁵N found, especially during stadials, in the Arabian Sea (Altabet et al., 1995; Emeis et al., 1995; Suthhof et al., 2001). It could have been stimulated by the supply of excess phosphate and iron from the more arid continents (Prins, 1999; Sirocko et al., 2000). In this case, N fixation in surface waters provided N with low δ¹⁵N that may have masked the high δ¹⁵N signal from denitrification.

5.3 Nitrogen cycling in the Holocene

During the Holocene the good coherence with the GISP δ¹⁸O record ceases. The global oceanic circulation of the Holocene is stabilized by the permanent salinity and density gradient between NADW and AAIW so that dramatic ocean-wide ventilation changes as in the Pleistocene cannot occur (Keeling and Stephens, 2001). In the Holocene the SAMW production is reduced so that ICW flowing into the Arabian Sea has a stronger contribution of IIW (Naidu and Govil, 2010). The considerable δ¹⁵N fluctuations by up to 1.5 ‰ indicate that existing changes in productivity and circulation can still lead to a pronounced Holocene reorganization of the nitrogen cycle within the basin. The different regional patterns (Fig. 6) can help to constrain the driving mechanisms. The present pattern of the decoupling of the productivity and denitrification maximum evolved in the middle and late Holocene as denitrification intensified in the northern and eastern part of the basin (Fig. 6b, c). The δ¹⁵N minimum between 9 and 5 ka BP is only found in the western part of the basin. It is most prominent in the Oman upwelling area and could be a signal of enhanced early and mid-Holocene OMZ ventilation. Benthic foraminifera indicate that oxygen concentrations were high and denitrification was low during this period despite enhanced productivity (Das et al., 2017). The vigorous upwelling during the Holocene climatic optimum was fed by inflow of ICW from the south, which could have better ventilated the western Arabian Sea and thus compensated for the enhanced respiration (Rixen et al., 2014).

Denitrification has continuously increased during the late Holocene in almost the entire basin but focused in the northern Arabian Sea (Fig. 6b). This trend coincides with Holocene cooling and a strengthening of the NE monsoon. Only in the open western Arabian Sea outside direct upwelling influence do δ¹⁵N values decrease in the late Holocene (Fig. 6d). This may be related to a shift of the Findlater Jet in the offshore direction as modeled by Le Mézo et al. (2017), which may have led to better nutrient availability in the western Arabian Sea. But as there are very few late Holocene data from this region (Table S1) this record has to be interpreted with caution.

Circulation probably changed after the sea level in the Persian Gulf and Red Sea reached its present position at about 6 ka BP, and water masses from these two basins prevented the strong ingression of ICW to the northeastern part of the Arabian Sea (Naidu and Govil, 2010; Pichevin et al., 2007). The ventilation of the OMZ by PGW and RSW today is restricted to the western part where the OMZ is much weaker than in the northeastern part of the basin (Gaye et al., 2013a; Morrison, 1997). The interplay of reduced OMZ ventilation in the north and enhanced NE monsoon productivity is likely the reason for the relocation of the OMZ and denitrification maximum to the NE during the Holocene. Enhanced productivity in the eastern Arabian Sea since the mid-Holocene as reconstructed from sediment cores (Agnihotri et al., 2003; Kessarkar et al., 2013, 2010) could have added to this relocation. It is likely that the inflow of low-density surface water suppressed primary productivity in the eastern Arabian Sea in the early Holocene. After precipitation declined and the sea level difference between the Bay of Bengal and Arabian Sea dropped at about 8 ka BP the inflow of warm, low-saline water with the northeast monsoon current and WICC to the eastern Arabian Sea declined (Mahesh and Banakar, 2014). This coincides with a rise in eastern Arabian Sea δ¹⁵N (Fig. 6c). The OMZ is generally weaker along the west coast of India due to the northward undercurrent, which leads to the oxygenation of subsurface water during the SW monsoon (Resplandy et al., 2012), and its upwelling and convective transport into surface waters along the coast is likely to explain the low δ¹⁵N in the sediments off the west coast of India (Fig. 5a).

Results from the global climate and ocean biogeochemistry model (KCM–PISCES; Sect. 3.4) driven by astronomical forcing over the Holocene suggest that ventilation changes were important drivers of the late Holocene Arabian Sea OMZ intensification (Fig. 8). The model produces a continuous increase in the OMZ volume in the Arabian Sea from 9 ka BP to the present. This is driven mainly by the increasing age (time since contact with the surface) of the water masses in the Arabian Sea OMZ. Arabian Sea export production is fairly constant in the model (Fig. 8) and can thus be ruled out as the driver for deoxygenation. An increase in export production is modeled only in a small area west of the southern Indian coast, indicating that export changes may have played only a local role (not shown). The decelerated circulation allowed more oxygen to be consumed by remineralization; it thus appears to be the main driver of the progressive deoxygenation in the model (Fig. 8) and can explain the increasing water column denitrification in the Arabian Sea in the δ¹⁵N records (Fig. 6).

Total organic carbon mass accumulation rates (TOC MARs; Fig. 6g) reflect productivity, organic matter preservation, and burial efficiency (Cowie et al., 2014; Cowie and Hedges, 1993; Müller and Suess, 1979). Similar to δ¹⁵N, Arabian Sea TOC MAR deviates from the global pattern (Cartapanis et al., 2016). Whereas the global TOC MAR significantly declines during deglaciation and remains low throughout the Holocene, the TOC MAR of the Arabian Sea shows the decline starting at about 20 ka BP but rises during the middle and late Holocene to values similar to those of the glacial (Fig. 6g). This pattern is consistent in the entire Arabian Sea (see data compiled by Cartapanis et al., 2016). As discussed in detail in Cartapanis et al. (2016) the drop in TOC MAR during deglaciation may indicate (i) a reduction of productivity, (ii) a lower transfer and burial efficiency of TOC to the sediments due to the reduced mineral ballast and the temporary storage on the growing continental shelves, and (iii) a reduced oxygen exposure time due to faster burial and reduced bottom water oxygenation of the glacial ocean. As we have discussed above, productivity increased in large parts of the basin due to a strengthening of the SW monsoon during deglacial warming. Productivity thus could not explain the reduced TOC MAR after the LGM. A reduced burial efficiency and increased deep water oxygen content are thus the most likely drivers of the deglacial TOC MAR drop in the basin. The Holocene rise in TOC MAR, inconsistent with the global trend, may likewise be due to better preservation caused by progressive mid-water deoxygenation so that oxygen exposure time again decreased. TOC MAR may have been augmented by enhanced NE monsoon productivity in the northern part of the basin and increasing burial efficiency by rising dust supply from the continents due to aridification after the mid-Holocene (Menzel et al., 2014; Overpeck et al., 1996; Prasad et al., 2014; Sirocko et al., 1993). The results of the KCM model, however, imply that productivity changes are not required as the increasing age of the water mass intensifies the Arabian Sea OMZ during the Holocene.

Figure 8 Simulated volume of the Arabian Sea OMZ (70 µmol L⁻¹ threshold, black) and export production in the entire Arabian Sea (red), western Arabian Sea (green) and eastern Arabian Sea (blue) over the Holocene (a). OMZ volume and age of water masses (time since contact with the surface) averaged over the OMZ (b). Arabian Sea defined as 55–75^∘ E, 8.5–22.5^∘ N; border between west and east defined at 68.5^∘ E.

6 Summary and conclusions

The compilation of SST reconstructions from the Arabian Sea showed up to 4 ^∘C lower glacial SST compared to the Holocene. Glacial ocean surface circulation in the Arabian Sea was generally reduced compared to Holocene circulation. Monsoonal upwelling along the western coasts was very much reduced or absent, as was the inflow of low-salinity water from the Bay of Bengal. Therefore, the general temperature gradient had a stronger insolation-driven N–S trend compared to the circulation-driven NW–SE trend of the Holocene. Upwelling indices calculated from the temperature difference between the northern Arabian Sea and the Somali and Oman upwelling centers reveal depleted or even paused upwelling during the LGM and stadials, while upwelling was enhanced during IS. A shift of the Findlater Jet to a stronger E–W orientation during the glacial could have prevented upwelling off Somalia, while it continued off Oman at a reduced rate. The prevalence of strong upwelling during the mid-Holocene despite a weakening of the SW monsoon could be due to a shift of the Findlater Jet in the offshore direction.

The compilation of δ¹⁵N data shows strong millennial-scale oscillations during the glacial with depleted δ¹⁵N values during the LGM and stadials and enriched δ¹⁵N during IS. These oscillations are caused by changes in the OMZ ventilation driven by millennial-scale fluctuations in the oceanic thermohaline circulation. Complete oxygenation of the OMZ occurred during stadials when SAMW formation was enhanced, and it was more strongly oxygenated due to its reduced residence time. The relative instability in thermohaline circulation led to fast changes in OMZ oxygenation with minima during the IS. In analogy to conditions in the recent Bay of Bengal, δ¹⁵N values of 4–6 ‰ in glacial sediments may not necessarily indicate that denitrification was completely absent. Moderate or occasional denitrification may have taken place in a more oxygenated OMZ, but its δ¹⁵N signal was not recorded in the sediment because the sub-thermocline water mass was isolated from the euphotic zone by stratification. Also, atmospheric N sources could have contributed to the low δ¹⁵N and compensated for the enhanced δ¹⁵N in the OMZ.

Holocene δ¹⁵N fluctuations by up to 1.5 ‰ and different patterns in the Arabian Sea show a strong local reorganization of the nitrogen cycle as global climate and thermohaline circulation provides more stable conditions. Stronger upwelling in the mid-Holocene was accompanied by stronger ventilation of the western part of the basin, which is also ventilated by PGW, especially after sea level reached its present maximum at 6 ka BP. The present denitrification maximum in the northeastern part of the basin was formed during the middle and late Holocene and is induced by a strengthened NE monsoon due to Holocene cooling but also due to reduced ventilation of the northern part of the basin. Results of the KCM–PISCES model simulation show a progressive intensification of the OMZ over the entire model run of 9 kyr. Productivity is constant in the model and the main driver of increasing deoxygenation and denitrification is the prolonged residence time of OMZ waters. OMZ intensification probably explains the increase in TOC MAR throughout the Holocene also deviating from the global trend.