Variability of the Indian Monsoon in the Andaman Sea Across the Miocene‐Pliocene Transition

We reconstructed the variability of the Earth's strongest hydrological system, the Indian monsoon, over the interval 6.24 to 4.91 Ma at International Ocean Discovery Program (IODP) Expedition 353 Site U1448 in the Andaman Sea. We integrated high‐resolution benthic and planktic foraminiferal carbon and oxygen isotopes with Mg/Ca measurements of the mixed layer foraminifer Trilobatus sacculifer to reconstruct the isotopic composition of seawater (δ18Osw) and the gradient between planktic and benthic foraminiferal δ13C. A prominent increase in mixed layer temperatures of ~4°C occurred between 5.55 and 5.28 Ma, accompanied by a change from precession‐ to obliquity‐driven variability in planktic δ18O and δ18Osw. We suggest that an intensified cross‐equatorial transport of heat and moisture, paced by obliquity, led to increased summer monsoon precipitation during warm stages after 5.55 Ma. Transient cold stages were characterized by reduced mixed layer temperatures and summer monsoon failure, thus resembling late Pleistocene stadials. In contrast, an overall cooler background climate state with a strengthened biological pump prevailed prior to 5.55 Ma. These findings highlight the importance of internal feedback processes for the long‐term evolution of the Indian monsoon.


Long-Term Variability of the Indian Monsoon
Monsoonal circulation in the Bay of Bengal and Andaman Sea is primarily induced by seasonal changes in the thermal gradient between the southern Indian Ocean and Asian continental landmasses, in particular the southern Tibetan Plateau (e.g., Webster et al., 1998). Long-term variations in monsoonal intensity are influenced by insolation changes induced by variations in the Earth's precession. During precession minima, when the perihelion occurs during Northern Hemisphere summer, the intensity of summer monsoonal rainfall increases, due to the development of an intense low pressure cell over the Asian continent. Changes in the Earth's obliquity additionally control the seasonality of radiation and the interhemispheric insolation gradient, driving interhemispheric heat and moisture transfer and enhancing summer monsoon rainfall (Bosmans et al., 2015;Mantsis et al., 2014).
The pre-Pleistocene variability of Indian monsoonal circulation and precipitation patterns in the Bay of Bengal, the Earth's strongest monsoonal regime, is still poorly understood (e.g., Clemens et al., 2016a;Kuhnt et al., 2020). In particular, the sensitivity of the East Asian and Indian monsoon systems to internal feedback processes such as changes in atmospheric greenhouse gases and ice volume, which affect zonal and meridional temperature gradients and moisture transport, remains a matter of intense debate (e.g., Clemens et al., 2018;Gebregiorgis et al., 2018). The prevailing view that precessional Northern Hemisphere summer insolation is the main control on the Indian monsoon (Kutzbach, 1981;Ruddiman, 2006) has been challenged by a growing number of studies, suggesting that heat and moisture transfer from the Southern Hemisphere driven by orbital obliquity and ice volume feedbacks also influences the long-term variability of the Indian monsoon (e.g., Bosmans et al., 2015;Clemens & Prell, 2003, 2007Gebregiorgis et al., 2018;Mantsis et al., 2014). This combination of external radiation forcing with internal climate feedbacks results in phase shifts and lagged responses of monsoon proxies that cannot be explained by a straightforward response to Northern Hemisphere radiation (e.g., Clemens & Prell, 2007).
The late Miocene to early Pliocene period is of particular interest, since it allows evaluation of monsoonal sensitivity in relation to changing orbital parameters during an interval when climate boundary conditions showed major shifts. This period is characterized by a major cooling trend between 7.0 and 5.5 Ma (Herbert et al., 2016;Holbourn et al., 2018) associated with aridification of central Asia (Miao et al., 2012) and spreading of C4 grasslands that are better adapted to low pCO 2 and more arid conditions (Cerling et al., 1997). The late Miocene cooling was terminated by a warming spell, starting at 5.5 Ma and culminating at 5.3 Ma, which also marks the final stage of the Messinian Salinity Crisis in the Mediterranean close to the Miocene/Pliocene boundary (e.g., Vidal et al., 2002). Within the early Pliocene and late Miocene, pronounced, transient increases in benthic foraminiferal δ 18 O, corresponding to Thvera (T) and Thvera-Gilbert (TG) cold stages defined by Shackleton et al. (1995), remain highly enigmatic events. The unusually high benthic δ 18 O values during prominent cold isotope stages TG12, TG14, TG20, TG22, and TG34 have been tentatively related to temporary buildup of Northern Hemisphere ice sheets between 5.5 and 6.0 Ma . However, the relationship between these rapid climate oscillations and tropical climate variability is poorly understood, since continuous, high-resolution records from tropical oceans spanning the Miocene-Pliocene transition are especially scarce.
Our new high resolution proxy records, derived from sediment cores recovered during Integrated Ocean Discovery Program (IODP) Expedition 353 at Site U1448 (10°38.03′N, 93°00′E, 1,098 m water depth) in the Andaman Sea ( Figure 1) closely track the variability of the Indian monsoon across the Miocene-Pliocene transition. We present new benthic and planktic foraminiferal stable oxygen (δ 18 O) and carbon (δ 13 C) isotope records together with Mg/Ca-derived mixed layer temperature reconstructions for the interval between 6.25 and 4.91 Ma. We use δ 18 O sw reconstructions from paired stable isotope and Mg/Ca analysis on the mixed layer dwelling foraminifer Trilobatus sacculifer as well as differences between planktic and benthic stable isotopes (Δδ 13 C and Δδ 18 O) to assess regional changes in precipitation and primary productivity. To identify potential forcing mechanisms, we performed spectral analysis on our climate proxy records. Our results provide new insights into the response of the Indian monsoon to orbital forcing and highlight the role of internal feedbacks during an interval of warming climate.

Regional Hydrologic and Oceanographic Setting
Today, summer monsoon precipitation in the Bay of Bengal accounts for up to 75% of the annual rainfall in this region (Dhar & Nandargi, 2003). Massive amounts of water and sediment are transported into the  (Zweng et al., 2013) plotted using Ocean Data View (Schlitzer, 2015).
northeastern Indian Ocean via the Ganges/Brahmaputra and Irrawaddy/Salween river systems during and after the wet season (June to August) (Chapman et al., 2015;Milliman & Farnsworth, 2011). This leads to seasonal salinity differences, which are greatest between the late summer/autumn (August-October) and late winter/spring (February-April) months due to the lagged salinity response to monsoonal freshwater fluxes (Akhil et al., 2016). The resulting buildup of freshwater pools and extreme regional differences in sea surface salinity (Figure 1) promotes stratification of the upper ocean in late summer and autumn, thereby prohibiting upwelling of cold and nutrient rich water into the mixed layer. Weakening of surface stratification in late winter/spring allows vertical transport of nutrients, enhancing primary productivity (Kay et al., 2018). Additionally, the Southwest Monsoon Current contributes to the strong contrast between upper ocean stratification during boreal winter and the eastward transfer of nutrient rich high-salinity upwelling waters from the Arabian Sea during boreal summer (Gordon et al., 2016).
The export of organic matter into the deep ocean in the Bay of Bengal and Andaman Sea is increased by aggregation of organic material to the large amount of fluviatile mineral particles, raising the sinking speed and prohibiting decomposition of organic matter in the water column (e.g., Le Moigne et al., 2013), leading to the development of a weaker oxygen minimum zone compared to the Arabian Sea (Al Azhar et al., 2017). Low dissolved oxygen (<50 μmol/kg) and δ 13 C values below zero in the Andaman Sea intermediate water masses indicate a strong influence of the oxygen minimum zone in the Bay of Bengal, while, in contrast, Andaman Sea surface waters are enriched in oxygen and δ 13 C (Talley, 2013). The deep waters in the Andaman Sea are composed of a mixture of Indian Central and Deep Water from the equatorial and Southern Hemisphere Indian Ocean and local thermocline water including Indonesian Throughflow water generated from North Pacific Intermediate Water modified by intense mixing while flowing through the Indonesian seaways (Gordon, 2005;Tomczak & Godfrey, 2002). Pacific-derived intermediate and deep water masses may have played a more important role in the Andaman Sea during the late Cenozoic when the connection between the Pacific and Indian Oceans was less restricted .

Sampling Strategy
Our study is based on marine sediments from two holes drilled with the half-length advanced piston corer (HLAPC) and the extended core barrel (XCB) drilling system at Site U1448 in the Andaman Sea ( Figure 1). Detailed site location, core recovery, and lithological descriptions can be found in Clemens et al. (2016aClemens et al. ( , 2016b. The shipboard composite depth scale was subsequently revised between 351.65 and 401.69 r-mcd (revised meters composite depth) and extended down to 439.45 r-mcd . Cores were sampled at 15 cm intervals between 353-U1448B-47F-1, 27 cm and 353-U1448B-47F-CC, 29 cm (346.73 to 351.63 mcd), following the original splice and at 10 cm intervals between 353-U1448B-48F-1, 2 cm and 353-U1448A-56X-5, 53 cm (351.76 to 439.37 r-mcd) following the revised splice.

Stable Isotopes
All samples were dried in an oven at 40°C prior to weighing and washing through a sieve with a mesh size of 63 μm. The residue was dried at 40°C on a sheet of filter paper. Dry samples were weighed and separated in four size fractions using sieves with 315, 250, and 150 μm mesh sizes. Stable isotope analysis was performed on three to eight well preserved specimens of the epibenthic foraminifers Cibicidoides wuellerstorfi and/or Cibicidoides mundulus (size fraction >250 μm). In samples where these species were not abundant, we analyzed other Cibicidoides species (C. robertsonianus, C. bradyi, and C. pachyderma in 74 samples) or the endobenthic genus Uvigerina (in 69 samples). We applied the commonly used correction factor of −0.64‰ to Uvigerina δ 18 O measurements to normalize the benthic δ 18 O to C. wuellerstorfi (Shackleton & Opdyke, 1973).
We selected~40 well preserved specimens per sample of the mixed layer foraminifer T. sacculifer (size fraction 250-315 μm) for paired Mg/Ca and stable isotope analysis. Samples were crushed between two glass plates to open all chambers of the tests. Approximately one quarter of the sample material (~10 specimens) was used for stable isotope analysis and three quarters (~30 specimens) for Mg/Ca analysis. In 32 samples with low abundance of T. sacculifer, we used all specimens from the size fraction 250-315 μm for Mg/Ca analysis (minimum 200 μg) and performed stable isotope measurements on tests separately picked from the >315 μm size fraction of the respective samples. After crushing, samples were cleaned with ethanol in an ultrasonic bath and oven dried at 40°C.
Measurements were performed at the Leibniz Laboratory, Kiel University, with a Finnigan MAT 253 mass spectrometer connected to a Kiel IV device for automated CO 2 preparation from carbonate samples. Samples were reacted by individual acid addition (99% H 3 PO 4 at 75°C). A precision better than ±0.09‰ on lab-internal and international standards is reached. Results were calibrated using the NIST (National Institute of Standard and Technology, Gaithersburg, Maryland) carbonate isotope standard, NBS (National Bureau of Standard) 19, and IAEA-603 and are reported in ‰ on the Vienna PeeDee Belemnite (V-PDB) scale.

Mg/Ca Paleothermometry
We measured Mg/Ca on the mixed layer dwelling foraminifer T. sacculifer. The samples were cleaned following the protocol of Martin and Lea (2002). Cleaning steps include methanol treatment for removal of clays, hydrazine treatment for removal of metal oxides, and hydrogen peroxide treatment for removal of organic matter and a weak acid leach (100 μl of 0.01 N HNO 3 ). After dissolution in nitric acid and weight dependent dilution, the samples were measured at the ICP-MS Laboratory, Institute of Geosciences, University of Kiel with a Spectro Ciros CCD SOP ICP-OES and calibrated using the ECRM752-1 standard. Lack of correlation between trace element ratios (Fe/Ca, Al/Ca, and Mn/Ca) and Mg/Ca ratios indicates a successful cleaning procedure. Duplicate measurements of 29 samples show an average standard deviation of 0.148 mmol/mol. To estimate mixed layer temperatures, we used the exponential equation for T. sacculifer (without sac, 350-500 μm) of (Anand et al., 2003) derived from sediment trap samples: A comparison of temperature estimates using the equations of Nürnberg et al. (1996) derived from culture experiments on T. sacculifer and the core top calibration of Dekens et al. (2002) on T. sacculifer without sac (250-350 μm) shows differences of <1°C ( Figure S3 in the supporting information). Several recent studies discussed nonthermal influences, such as pH and salinity, on Mg/Ca variability in foraminiferal calcite (e.g., Gray & Evans, 2019;Holland et al., 2020;Tierney et al., 2019). These authors found a negligible sensitivity of T. sacculifer to changes in pH. We evaluated the influence of salinity on our Mg/Ca based temperature reconstructions following the equation of Gray and Evans (2019) for T. sacculifer ( Figure S4). Assuming a variability of mixed layer salinities between 31 and 33 psu, which is close to the modern variability in the Andaman Sea (Figure 1), the effect on reconstructed temperature changes is between 1°C and 2°C ( Figure S4). Using the equation of Anand et al. (2003) without salinity correction, we aim to directly compare our record to the published Mg/Ca-based temperature reconstruction of ODP Site 1146 , rather than provide precise absolute values.
The uncertainties in mixed layer temperature values originating from Mg/Ca measurements and the Mg/Ca calibration equation are calculated using the error propagation equations of Mohtadi et al. (2014). Following the approach of Steinke et al. (2010), a correction was applied for long term changes in the Mg/Ca composition of seawater between the Pliocene ratio of 4.3 mol/mol (Evans et al., 2016) and the modern ratio of 5.1 mol/mol (Broecker & Peng, 1982). Mg/Ca-derived temperatures are presented on both, a corrected and an uncorrected scale. The covariance between temperature and seawater carbon chemistry (dissolved inorganic carbon [DIC] concentration and carbonate ion ([CO 3 ] 2− ) saturation) also influences temperature estimates on short timescales (Holland et al., 2020). The late Miocene/early Pliocene variability of these parameters, however, is unknown. Furthermore, the calcite saturation state of the bottom water (Ω) has been shown to be after temperature the second most important influence on Mg/Ca of foraminiferal shells (Tierney et al., 2019). Modern values of Ω at 1,000 m water depth close to the core location in the Andaman Sea are approximately 1.4 according to the GLODAPv2 Bottle dataset (Olsen et al., 2019), which suggests that post depositional dissolution alters foraminiferal Mg/Ca only to a minor extent at the shallow location of Site U1448 (1,098 m water depth).

Reconstruction of δ 18 O sw
The δ 18 O value of foraminiferal calcite is influenced by the calcification temperature and the oxygen isotopic composition of the ambient seawater (δ 18 O sw ), which is a function of ice volume, the local 10.1029/2020PA003923

Paleoceanography and Paleoclimatology
precipitation-evaporation budget and the isotopic composition of rain and/or riverine runoff. We used the δ 18 O-temperature relationship of Bemis et al. (1998) and added 0.27‰ to convert to Standard Mean Ocean Water (Hut, 1987): We assumed that the ice volume related component of benthic δ 18 O ranged between~70%, as reconstructed from paired benthic δ 18 O and Mg/Ca measurements for the middle Miocene (Shevenell et al., 2008), and~50% during the late Pleistocene transition (Sosdian & Rosenthal, 2009). Here we applied an ice volume correction of 50%, since we assumed that the temperature component to benthic δ 18 O is increased at the relatively shallow water depth of 1,098 m at Site U1448. A comparison of corrected δ 18 O sw values is shown in Figure S5.

Spectral Analysis
Spectral analysis was performed with the Past3 software (Hammer et al., 2001) on an unevenly spaced time series. We used the REDFIT function, an implementation of the REDFIT procedure of Schulz and Mudelsee (2002) with a Blackman-Harris window of two segments and a frequency oversampling value of 2. Significance lines of 99%, 95%, 89%, and 85% confidence intervals (CIs) are based on parametric approximation (chi-square test). We divided the U1448 planktic δ 18 O and δ 18 O sw records in two distinct intervals (6.20 to 5.55 Ma and 5.55 to 4.90 Ma) for analysis, since a distinct change in cyclicity and signal amplitude occurs at 5.55 Ma. For the evolutive wavelet spectrum, we used the wavelet transform function in Past3 on an evenly spaced time series with a 1 kyr temporal resolution with a Morlet basis function. The p ¼ 0.05 significance level derived by a chi-square test is displayed as a contour line after Torrence and Compo (1998). A cone of influence indicates the area with possible boundary effects. Cross spectral analysis of our planktic and benthic δ 18 O, δ 18 O sw , and mixed layer temperature records against precession and obliquity was performed with AnalySeries 2.0.8 on linear detrended data on a common 1 kyr-spaced scale with the Blackman-Tukey spectral method using a Bartlett window.

Chronology
The age model is derived by correlation of the U1448 benthic foraminiferal stable isotope data to the well dated records of ODP Site 1146  and ODP Site 982 (Drury et al., 2018) (Figures 2, S1, S2a, and S2b and supporting information Text S1). For correlation, we defined eight tie points (Table 1) at distinct maxima in δ 18 O that also coincide with minima in δ 13 C ( Figure 2). These pronounced δ 18 O excursions were identified as isotope stages T2, T4, T8, TG2, TG12, TG14, TG22, and TG34, following the nomenclature of Shackleton et al. (1995), where odd numbers refer to warm stages, even numbers to cold stages and letters to the relevant magnetic polarity chron. The average sedimentation rate between 6.24 and 4.91 Ma is 0.08 m/kyr with a maximum of 0.14 m/kyr and a minimum of 0.03 m/kyr, resulting in a temporal resolution of 1.44 kyr for the stable isotope series ( Figure 2).

Mixed Layer Foraminiferal Stable Isotope Records
Planktic foraminiferal δ 18 O exhibits well-defined cyclic variations between 6.24 and 4.95 Ma (Figure 3c). Highest values are reached during cold stages T4, T8, TG2, TG4, TG20, TG22, and TG34 with a maximum difference of 1.42‰ between cold stage T6 (−2.13‰ at 5.01 Ma) and warm stage T5 (−3.55‰ at 5.04 Ma). The δ 18 O record can be divided into two distinct intervals, based on changes in dominant cyclicity and signal amplitude.
Between 6.24 and 5.55 Ma, mean δ 18 O values fluctuate around −2.69‰ (n ¼ 301, stdev ¼ 0.29) and decrease slightly to −2.90‰ (n ¼ 372, stdev ¼ 0.36) between 5.55 and 4.95 Ma. To assess changes in amplitude variability, we calculated minima and maxima using a 20 kyr moving window on δ 18 O values interpolated to 0.5 kyr intervals using the Stineman function in KaleidaGraph. Amplitude variability is higher between 5.55 and 4.95 Ma in comparison to the older interval between 6.24 and 5.55 Ma. The most extreme fluctuations are evident in the δ 18 O maxima, while the δ 18 O minima variability remains consistent.

Mixed Layer δ 18 O sw Reconstruction
The temperature and ice volume-corrected δ 18 O sw record retains the high-amplitude and cyclic behavior of planktic δ 18 O (Figure 4b). Values are relatively stable between 6.20 and 5.74 Ma with a mean of −1.12‰ (n ¼ 210, stdev ¼ 0.26), increasing from 5.74 to 5.25 Ma to a mean of −0.47‰ (n ¼ 55, stdev ¼ 0.30) between 5.30 and 5.20 Ma, except for a short-term rebound at 5.40 Ma. The highest values above 0‰ occur between 5.30 and 5.20 Ma (0.17‰ at 5.28 Ma and 0.18‰ at 5.24 Ma). From 5.25 to 5.10 Ma, δ 18 O sw exhibits a declining trend and displays a mean of −0.76‰ between 5.10 and 4.95 Ma (n ¼ 108, stdev ¼ 0.29).

. Spectral and Evolutive Wavelet Analysis
A distinct change in cyclicity in the evolutive wavelet power spectrum of T. sacculifer δ 18 O (Figure 5a) occurs at 5.55 Ma, indicating a shift toward higher frequency precession related periodicities and an increase in spectral power after 5.55 Ma. A similar but less pronounced shift toward higher periods is evident in the evolutive power spectrum of δ 18 O sw at 5.55 Ma.

Cross Spectral Analysis
We calculated cross-phase spectra between the precessional (19 and 23 kyr) and obliquity (41 kyr) components and the U1448 mixed layer temperature, benthic and planktic δ 18 O, and δ 18 O sw data for the intervals 6.20 to 5.55 Ma and 5.55 to 4.91 Ma (Table 3 and Figure 6). The phase relationships show fundamental differences in these two time slices. From 5.55 to 4.91 Ma, mixed layer temperature is almost in phase with obliquity maxima (1.0 ± 2.9 kyr), whereas benthic and planktic δ 18 O as well as δ 18 O sw lag obliquity maxima by 3.1 to 4.4 kyr. In the 23 kyr band, mixed layer temperature and benthic δ 18 O lag precession minima by 4.6 to 4.7 kyr, whereas δ 18 O sw lags by 7.6 kyr. By contrast, the phase relationships at the 23 kyr precession band are virtually reversed from 6.20 to 5.55 Ma: Mixed layer temperature and benthic and planktic δ 18 O lead precession minima by −7.1 and −7.9 kyr, respectively. Ice volume corrected δ 18 O sw leads precession minima by −5.1 kyr but only shows a weak coherence of 0.37, which lies below the 80% CI (nonzero coherence is higher than 0.385). Figure S6 shows that the effect of correcting δ 18 O sw for different estimates of the ice volume component (uncorrected, 50% and 70% ice volume component) on the cross spectral analysis of δ 18 O sw with precession and obliquity has a negligible influence on the phase relationship between orbital forcing and the response of δ 18 O sw .
Between 6.20 and 5.55 Ma, the phasing in the obliquity band is also fundamentally different to that in the 5.55 to 4.91 Ma interval: Mixed layer temperature and planktic and benthic δ 18 O lead obliquity maxima by −6.5 (benthic δ 18 O) and −10.8 kyr (mixed layer temperature), with an overall lower coherence. The

Paleoceanography and Paleoclimatology
and river discharge and is, thus, proportional to salinity (Sijinkumar et al., 2016). However, the modern regional relationship between δ 18 O sw and salinity in the northeastern Indian Ocean is highly variable, depending on seasonal changes in rainwater provenance, location, evaporation, and runoff/precipitation (Delaygue et al., 2001;Kumar et al., 2010;Singh et al., 2010). Extrapolation of the modern regional average mixed layer δ 18 O sw -salinity relationship into the geological past is, thus, complex and leads to relatively large errors in mixed layer temperature estimates and the δ 18 O sw -salinity regression slope (Sijinkumar et al., 2016). In particular, changes in the δ 18 O of precipitation and runoff strongly influence δ 18 O sw but not salinity, which is only a function of the amount of precipitation and runoff.
Two main factors control the freshwater δ 18 O in the Andaman area: (1) the δ 18 O composition of freshwater originating from the discharge of the Irrawaddy and Salween Rivers and (2) the δ 18 O composition of local rainfall over the Andaman Islands and adjacent seas, which is strongly dependent on the moisture source of precipitation. The δ 18 O composition of rainfall over the Himalayan mountain chain and Tibetan Plateau, the main source area of river discharge into the Bay of Bengal and Andaman Sea, varies spatially (with generally more depleted values in mountainous areas) and temporally due to summer depletion with prevailing monsoonal summer rainfall (Kumar et al., 2010). Across the eastern Tibetan Plateau, in the vicinity of the Brahmaputra and Irrawaddy drainage areas, local summer precipitation δ 18 O ranges between −13‰ and −18‰ (Bershaw et al., 2012). Summer (July-September, JJAS) rainfall δ 18 O over the southern Tibetan Plateau ranges from −18.26‰ in Lhasa to −16.62‰ in Nyalam (Gao et al., 2013). In contrast, tropical-subtropical summer monsoon rainfall close to the Ganges/Brahmaputra Delta flood plains exhibits considerably higher δ 18 O values (for instance, summer monsoon rainfall in the year 2004 averaged −5.74‰; Dar & Ghosh, 2017). Further south, values of −5.76‰ at Bangalore (Dar & Ghosh, 2017) and −2.16‰ in Sri Lanka (Edirisinghe et al., 2017) are comparable or even higher, as they are less influenced by rainout

10.1029/2020PA003923
Paleoceanography and Paleoclimatology depletion. As a result, δ 18 O sw variations in the Andaman area not only reflect changes in the local precipitation-evaporation budget but also the seasonality of rainfall and the relative contribution of local precipitation versus river discharge of δ 18 O depleted freshwater from mountainous regions.
On longer timescales, the local freshwater δ 18 O strongly depends on the location of the initial moisture source for the local rainwater. Today, modeled backward airmass trajectories indicate a common transport of precipitation into the core area of the Indian monsoon from the Arabian Sea (Dar & Ghosh, 2017). On a local scale, moisture sources for precipitation depend on the prevailing wind directions over the Indian Together, these studies underline the importance of large scale convective systems and local precipitation recycling processes in influencing the isotopic composition of summer monsoon precipitation. Changes in long-term averages of mixed layer δ 18 O sw in the Andaman Sea (Figure 4b) may, thus, be related to large scale changes in monsoonal atmospheric circulation in addition to changes in annual average or seasonal amount of precipitation and runoff.

δ 18 O sw Variability Across the Miocene-Pliocene Transition
High amplitude changes in the ice volume corrected δ 18 O sw record from Site U1448 follow the orbital cyclicity of planktic δ 18 O and are probably related to changes in the amount of precipitation/runoff into the northern Indian Ocean (Sijinkumar et al., 2016). High values associated with increased salinity during cold isotope stages indicate a weakening of Indian summer monsoon rainfall. In addition, an increase in mean amplitude variability (standard deviation increases from 0.25 prior to the warming between 6.20 and 5.74 Ma to 0.29 after the warming between 5.10 and 4.95 Ma) points toward an increased seasonality (monsoonality) of precipitation ( Figure 4b). However, δ 18 O sw also exhibits a distinct increasing trend between 5.74 and 5.25 Ma, coincident with a major increase in mixed layer temperature. Part of this δ 18 O sw increase must be caused by enhanced evaporation in a warming ocean. Furthermore, a change in the provenance of the monsoonal rainfall toward a larger proportion of moisture from sources in the southern equatorial Indian Ocean, where modern surface waters show significantly higher δ 18 O values between 0‰ and 5‰ (Schmidt et al., 1999), may have resulted in a heavier δ 18 O signal of precipitation.

Orbital Forcing of the Indian Monsoon
The covariance between monsoonal precipitation and changes in insolation, driven by variations of the Earth's orbit, is complicated by varying contributions of internal climate processes involving, for instance, Northern and Southern Hemisphere ice sheets, atmospheric and ocean circulation, the topography of monsoonal areas, and the concentrations of greenhouse gases in the atmosphere . The interaction between external insolation forcing and internal feedbacks that operate at different periods can lead in monsoonal areas to complex interference patterns, known as heterodynes (combination of two frequencies by a nonlinear mixer), calculated by adding/subtracting the interacting orbital frequencies (Table 4). Spectral analysis of the planktic δ 18 O and ice volume corrected δ 18 O sw time series from the Andaman Sea reveals several peaks outside the Milankovitch periodicities of 19 and 23 kyr (precession), 41 kyr (obliquity), and 98, 123, and 130 kyr (eccentricity). These additional peaks represent heterodynes resulting from mixing of the usual orbital frequencies ( Figure 5 and Table 4).
The 1/27 kyr heterodyne, which represents a mix of eccentricity and precessional frequencies, often features in Pleistocene monsoonal time series such as in the Indian monsoon record derived from the Bittoo Cave  (Kathayat et al., 2016), the East Asian Pearl River Valley δ 2 H wax record (Thomas et al., 2014), the South China Sea alkenone-derived sea surface temperature and runoff records (Thomas et al., 2014(Thomas et al., , 2016, and the rainfall reconstructions from 10 Be analysis of Chinese loess records (Beck et al., 2018). The occurrence of the 1/27 kyr heterodyne in the planktic δ 18 O and δ 18 O sw records from Site U1448 between 6.20 and 5.55 Ma together with the weaker spectral peaks (<95% CI) at the precessional frequencies 1/23 and 1/19 kyr suggest a predominantly local insolation control on precipitation ( Figure 5). This implies that internal climate processes played a less important role in driving the Indian monsoon prior to 5.55 Ma than external orbital forcing. In contrast, the planktic δ 18 O and δ 18 O sw records exhibit strong spectral power in the heterodyne 1/36 kyr in addition to precession and obliquity between 5.55 and 4.95 Ma. This heterodyne results from the interference of obliquity with precession (Table 4). In the δ 18 O sw record, precession and obliquity are less pronounced (1/23 kyr precession with >89% CI and 1/41 kyr obliquity with <85% CI). The occurrence of the 1/36 kyr heterodyne with a strong obliquity component suggests enhanced response of the Indian monsoon to obliquity forcing, possibly linked to a strengthened temperature gradient between the two hemispheres.
In the East Asian monsoon region, the heterodyne variance in the spectrum of δ 18 O sw during the Pleistocene increases from south to north (Thomas et al., 2016). The prominent occurrence of obliquity influenced heterodynes (1/29 and 1/72 kyr) in this region was explained by an increasing response of multiple environmental parameters to obliquity forcing, such as continental and oceanic moisture sources, terrestrial precipitation recycling and temperature (Thomas et al., 2016). Changes in moisture sources and intensification of the cross-equatorial flow of heat and moisture are also likely candidates for the appearance of the prominent 1/36 kyr heterodyne between 5.55 and 4.91 Ma in our Andaman Sea records.
In the interval 6.20 to 5.55 Ma, the phase relationship of the U1448 planktic and benthic δ 18 O and mixed layer temperature records at the 23 kyr precession band are almost 180°out of phase (Table 3 and Figure 6), suggesting that the ocean-climate system was more responsive to Southern Hemisphere precessional insolation. These changes may be related to variations in ice volume due to the different response time of smaller/larger ice shields during glaciation/deglaciation. In contrast, the phase relationships of monsoonal proxies (mixed layer temperature, planktic δ 18 O, and ice volume corrected δ 18 O sw ) to orbital precession and obliquity after 5.55 Ma are similar to the phase relationships in Pleistocene summer monsoon records from the Arabian Sea (Caley et al., 2011;Clemens & Prell, 2003) and from the tropical eastern Indian Ocean and Andaman Sea (Bolton et al., 2013;Gebregiorgis et al., 2018).
From 5.55 to 4.91 Ma, the phase of the U1448 monsoon proxy records plots on the obliquity band between the Pleistocene Arabian Summer monsoon stacks of Caley et al. (2011) and Clemens and Prell (2003), supporting intensified response to obliquity forcing following warming at 5.55 Ma ( Figure 6). Obliquity was recently suggested to play an important role as driver of late Pleistocene East Asian  and Indian monsoonal variability (Gebregiorgis et al., 2018). These authors proposed that monsoon intensification was induced by the strengthening of the cross-equatorial moisture flow due to the stronger interhemispheric temperature gradient at high obliquity. Our records suggest that the response of the Indian summer monsoon to obliquity forcing was even stronger on a warmer Earth. The increase in the amplitude of the U1448 mixed layer temperature and planktic δ 18 O records after 5.55 Ma together with the enhanced response to obliquity (Figure 4) suggests increasing moisture and heat transfer from the Southern Hemisphere into the Indian monsoonal region, as in the late Pleistocene.

Transient Cooling Events and Summer Monsoon Failure
The interval 6.0 to 5.5 Ma is characterized by transient, intense Northern Hemisphere cooling events with benthic δ 18 O maxima of 2.63‰ and 2.84‰ during the most prominent cold isotope stages TG22 and TG20 ( Table 2), suggesting that climate conditions temporarily reached the threshold for transient Northern Hemisphere ice buildups (Herbert et al., 2016;Holbourn et al., 2018). However, the mixed layer temperature difference between cold and warm isotope stages remains relatively small at Site U1448 (Figure 7c). The δ 18 O sw record shows only low amplitude variability (e.g., maximum difference 0.55‰ between TG23 and TG22), indicating reduced seasonal rainfall between 6.0 and 5.5 Ma (Figure 7c). These results suggest that dry, rather than cold conditions, prevailed in the Andaman Sea and that summer monsoon intensity was reduced.
In contrast, cold isotope stages between 5.50 and 4.91 Ma exhibit larger mixed layer temperature differences of up to 7°C with preceding warm isotope stages (Figure 7b). Variations in planktic δ 18 O and δ 18 O sw are also more pronounced, reaching 1.10‰ between TG5 and TG4 and 1.25‰ between TG3 and TG2, which suggests increased seasonality in temperature and precipitation. The increased temperature gradient can be explained by an enhanced cross-equatorial flow of heat and moisture from southerly sources during the summer monsoon season, also accounting for the presence of heterodynes ( Comparison of cold and warm TG isotope stages during the late Miocene cool period and early Pliocene warm period at Site U1448 with Pleistocene glacial-interglacial cycles (Gebregiorgis et al., 2018) indicates differences in mixed layer temperature of only~2°C in the late Miocene and Pleistocene in contrast tõ 6-7°C in the early Pliocene warm period (Figure 7). Pronounced differences in ice volume corrected δ 18 O sw in the early Pliocene are similar to Pleistocene glacial-interglacial differences, suggesting summer monsoon intensification associated with deglacial warming. Overall, monsoon variability changed from sustained cool, dry conditions prior to 5.50 Ma to much warmer and wetter conditions during the early Pliocene warm stages.

Coupling of Monsoonal Subsystems During Latest Miocene Warming
Temperature and productivity exhibit similar trends in the core region of the Indian monsoon at Site U1448 and in the East Asian Monsoon region at Site 1146 in the South China Sea   (Figure 8). Mg/Ca-derived temperatures at Site 1146 also show an initial warming trend after cold stages TG20 and TG22, interrupted by pronounced, transient cooling during TG12 and TG14 (Figure 8f). Furthermore, declining Δδ 13 C values in the 1146 record, indicating weakening of the biological pump between 5.85 and 5.35 Ma , are also interrupted by a short rebound between 5.65 and 5.50 Ma. This parallel evolution of temperature and productivity proxy records, indicating an intensification of the entire Asian summer monsoon system, suggests a close coupling of the two monsoonal subsystems in their response to global warming. Comparison of spectral characteristics in the planktic δ 18 O records from Sites U1448 and 1146 before and after 5.55 Ma reveals similar changes in cyclicity (Figure 9). A spectral peak of 23 kyr (precession) dominates the 1146 planktic δ 18 O record (<99% CI) between 6.20 and 5.55 Ma. After 5.55 Ma, the South China Sea record shows peaks at heterodyne frequencies of 1/36 (>85% CI) and 1/72 kyr (>95% CI), indicating an increased influence of obliquity. This coupled response highlights the increasing influence of obliquity as a pacemaker for the East Asian and Indian monsoonal subsystems in the latest Miocene.
A recent study showed that Pleistocene Indian monsoon precipitation was phase locked with obliquity and intensified during Southern Hemisphere precessional warming (Gebregiorgis et al., 2018). Precession-band variance accounted for~30% of the total variance of Indian monsoon precipitation. In contrast, precession completely dominates the orbital signal of the East Asian Monsoon, as inferred from Chinese cave δ 18 O records (Gebregiorgis et al., 2018). Yet, the interpretation of speleothem δ 18 O records as representative of East Asian Monsoon intensity remains a matter of intense discussion (Gebregiorgis et al., 2020;Lachniet, 2020;Zhang et al., 2019Zhang et al., , 2020. The different response of individual monsoonal subsystems to orbital forcing in the Pleistocene and the synchronous change in both monsoonal systems at the end of the Miocene highlight the importance of internal processes (ice volume, greenhouse gases, sea surface temperatures, and atmospheric/ocean circulation) in controlling monsoon variability. However, the underlying reason for the end-Miocene change remains an open question. Possible influential factors include changes in latitudinal and interhemispheric temperature gradients associated with the onset of transient Northern Hemisphere glaciations, vegetation shifts on Northern Hemisphere landmasses, changes in ocean circulation and heat budgets linked to modifications in ocean gateway configurations (closure and reopening of the Gibraltar gateway and initial restriction of the Panama and Indonesian gateways) and/or increasing atmospheric greenhouse gas concentrations in the latest Miocene.
A massive reduction in monsoonal precipitation and runoff from North Africa to the Mediterranean and northeast Atlantic also preceded warming at~5.5 Ma (e.g., Colin et al., 2008). The Indian summer monsoon failure between 5.60 and 5.55 Ma in the Andaman Sea record closely corresponds to the interval of near complete desiccation of the Mediterranean during the second phase of the Messinian Salinity Crisis between 5.59 and 5.30 Ma (Krijgsman, Hilgen, et al., 1999;Krijgsman, Langereis, et al., 1999;Vidal et al., 2002). This drying out was linked to a major decrease in river discharge from drainage basins within the North African monsoon belt (Colin et al., 2014). The coeval occurrence of cool and dry conditions between 5.60 and 5.55 Ma in three major monsoonal systems of the Northern Hemisphere regimes indicates a major reduction in tropical convection, possibly associated with southward displacement of the monsoonal rain belt and Intertropical Convergence Zone in the Northern Hemisphere. Monsoonal rainfall rapidly intensified in these three regions with renewed warming after 5.55 Ma. However, the stronger response to obliquity forcing suggests enhanced contribution of cross-equatorial moisture transfer from sources in the southern equatorial Indian Ocean.

Conclusions
At Andaman Sea Site U1448, planktic δ 18 O shows increased amplitude variability and a shift from a dominant precession (19 and 23 kyr) to obliquity (41 kyr) signal, which coincides with the onset of global warming at 5.5 Ma (Herbert et al., 2016;Holbourn et al., 2018). Mixed layer warming between 5.55 and 5.28 Ma follows an initial rise between 5.78 and 5.65 Ma, interrupted during the intense cold isotope stages TG12 and TG14 (5.65 to 5.55 Ma). This warming trend was paralleled by a stepwise decrease in primary productivity, suggesting intensified upper ocean stratification stemming from increased seasonal freshwater input.
Mixed layer temperature and δ 18 O sw reconstructions suggest that late Miocene cold events in the Andaman Sea were characterized by drier rather than colder conditions, thus resembling summer monsoon failures  during late Pleistocene stadial and glacial stages. By contrast, substantial variations of up to 7°C in mixed layer temperatures and large differences in ice volume corrected δ 18 O sw between the early Pliocene cold and warm isotope stages suggest increased transfer of heat and moisture from the Southern Hemisphere into the Bay of Bengal.
Comparison with South China Sea Site 1146 reveals a similar evolution of mixed layer temperature and productivity records, indicating summer monsoon intensification and enhanced response to obliquity forcing associated with global warming in both the Indian and East Asian regions. The late Miocene warming and summer monsoon intensification in the Andaman Sea was also concomitant with the end of the Messinian Salinity Crisis in the Mediterranean, culminating at 5.28 Ma during warm stage TG5, suggesting common drivers of climate change.

Data Availability Statement
Data presented in this study are available from the Data Publisher for Earth and Environmental Science (PANGAEA.de; https://doi.pangaea.de/10.1594/PANGAEA.921091).