Assessing the potential of Individual Foraminifera Analyses (δ¹⁸O) to reconstruct variability, seasonality and extremes in the tropical Indian Ocean

The Indian Ocean plays a critical role in the global climate system by regulating heat and moisture transport, influencing major fluxes that drive atmospheric circulation. Individual foraminiferal analyses (IFA) of the stable oxygen isotope composition (δ¹⁸O) of their test provides a powerful approach to reconstruct past climate variability beyond changes in mean conditions, with the potential to capture seasonal to interannual climatic extremes associated with phenomena such as the Indian monsoon and the Indian Ocean Dipole (IOD). However, field-based calibrations of this proxy remain sparse in the Indian Ocean, especially in the western part of the basin, limiting our understanding of how IFA records regional hydrographic signals. We combine a forward-modelling approach with IFA in Indian Ocean core-top sediment samples to evaluate how δ¹⁸O variability responds to changes in temperature and salinity linked to seasonal and interannual variability. Using ocean reanalysis data, we simulate surface and thermocline conditions at several sites representative of different oceanographic provinces across the Indian Ocean and generate synthetic core-top IFA datasets through random sampling. We evaluate the sensitivity of results by artificially amplifying or reducing seasonal and interannual climate variability. Forward-modelled data are then compared with IFA data from core-tops at four key sites across the Indian Ocean. Our results show that surface-ocean individual foraminiferal δ¹⁸O changes in the western and central Indian Ocean are mainly controlled by temperature seasonality related to monsoon dynamics, whereas in the eastern Indian Ocean, seasonality is less pronounced and interannual variability is more dominant. On the other hand, thermocline δ¹⁸O variability for all sites is primarily associated with interannual variability and temperature changes, highlighting the potential of thermocline species to record interannual temperature variability. Forward-modelled IFA match the SD of real datasets from Late Holocene core-tops at most sites for surface and thermocline depths. This framework provides a basis for interpreting IFA in the Indian Ocean in terms of underlying climate processes and offers perspectives for reconstructing past seasonal and interannual climate variability in this important region.

1 Introduction

The evolution of climate variability at seasonal and interannual timescales plays a key role in shaping the occurrence and intensity of extreme weather events. In order to better understand the sensitivity of seasonal to interannual climate variability and its consequences on the hydrological cycle (Swain et al., 2025), it is essential to reconstruct variability at these timescales under climatic background conditions different from those observed during the instrumental period.

Seasonal to interannual sea surface temperature (SST) variability in the tropical oceans disrupts the global hydrological cycle, with strong repercussions for tropical countries (Rasmusson and Carpenter, 1982; Gadgil, 2003). Ocean-atmosphere interactions give rise to phenomena such as El Niño–Southern Oscillation (ENSO) in the equatorial Pacific Ocean and the Indian Ocean Dipole (IOD) in the tropical Indian Ocean. ENSO is the primary mode of interannual variability on Earth and leads to strong modifications in equatorial upwelling and easterlies. Whilst the IOD presents relatively weaker variability than ENSO, it manifests on the eastern and western sides of the Indian Ocean basin (Saji et al., 1999). Its consequences for regional hydroclimate are profound, driving anomalous rainfall, drought, and flood events across East Africa, the Indian subcontinent, and Southeast Asia (Ummenhofer et al., 2009; Anderson et al., 2019), making it one of the most socially impactful climate modes in the world, affecting the livelihoods of over a billion people across the basin. IOD variability is mainly influenced by the equatorial SST gradient and the direction of prevailing winds. On the seasonal timescale, the northern Indian Ocean is under the influence of seasonally-reversing monsoon winds, with the southwest monsoon bringing high precipitation to southern Asia and the Bay of Bengal during boreal summer and a northeast monsoon during boreal winter, bringing moisture to southern India, Sri Lanka, and parts of Southeast Asia, resulting in drier conditions over much of the northern Indian Ocean and South Asia. The monsoon impacts the hydrological cycle and in situ populations living around the basin, causing extreme events.

Understanding the variability of the Indian Ocean climate and extreme events remains challenging, especially in the context of climate change. Future projections from modelling studies suggest that the Indian Ocean could harbor enhanced variability with more intense and frequent extreme events arising due to an intensified zonal SST gradient (Abram et al., 2008; Cai et al., 2013, 2014; Roxy et al., 2014; Zheng et al., 2013). Indeed, recent observations suggest that these changes are already significantly altering the hydrological cycle across the Indian Ocean basin (Ummenhofer et al., 2018). Yet, the future evolution of rainfall dynamics at regional scales remains model-dependent, and their impact on hydrological cycle extremes may depend on how well the future evolution of Indo-Pacific SST variability is reproduced by climate models (Park et al., 2022). Documenting how large climate perturbations may have affected the magnitude and frequency of climate extremes in south to south-east Asia in the past could help evaluate the ability of climate models to reproduce the ocean-atmosphere interactions that underpin climate extremes. To reconstruct past climate variability, individual planktonic foraminiferal analyses (IFA) are a well-established technique that can provide such diagnostics under contrasting climate backgrounds (Fehrenbacher et al., 2024; Leduc et al., 2009). Traditional δ¹⁸O paleoclimate analysis methods focus on reconstructing average climate conditions, while IFA offers the advantage of recording the full climate spectrum within a time slice rather than mean climate (Thirumalai et al., 2013). Due to low (generally <5 cm kyr⁻¹) accumulation rates and bioturbation in open-ocean marine sediments, multi-test “bulk” foraminiferal analyses cannot typically record seasonal or interannual climate signatures. Yet they contain interesting snapshots on total variability that, once correctly calibrated with respect to regional ocean characteristics, can help decipher the mode of climate variability that influenced the recorded changes in IFA variability. IFA has been applied in many studies to reconstruct past modes of ENSO or IOD variability (Leduc et al., 2009; Ford et al., 2015; Thirumalai et al., 2019, 2024) and to capture monsoon dynamics (Ganssen et al., 2011; Thirumalai and Clemens, 2020; Thirumalai et al., 2025).

The application and accurate interpretation of IFA in the Indian Ocean requires careful consideration to account for region-specific climate, ocean dynamics, and foraminiferal ecology using modern datasets (Thirumalai et al., 2013; Bienzobas Montávez et al., 2024). We aim to understand how modern climatic and oceanographic variability in the Indian Ocean is recorded by forward-modelling IFA-δ¹⁸O and testing the sensitivity of this proxy using modern reanalysis datasets (ORAS5). We do this at four distinct sites representative of distinct oceanographic provinces, in northwestern, western, central and eastern Indian Ocean (Fig. 1). Forward modelled data are then compared with new and published core-top (late Holocene) IFA measurements performed on the planktonic species Globigerinoides ruber and Neogloboquadrina dutertrei, to constrain surface and thermocline variability, respectively. We will discuss the applicability of IFA to reconstruct variability and extremes in the Indian Ocean by addressing the following specific objectives in three steps. (1) Describe seasonal and interannual variability patterns in the surface and thermocline based on reanalysis in our key areas; (2) Illustrate how an increase in seasonality and interannual variability would manifest itself in modelled IFA by modifying these patterns based on reanalysis; and (3) Compare IFA measured in real core-tops with forward-modelled IFA from the reanalysis.

2 Study sites and climatic setting of the tropical Indian Ocean

We focus on four sites distributed across the western and central tropical Indian Ocean to investigate the controls on planktonic carbonate δ¹⁸O (δ¹⁸O_c) variability using both IFA forward-modelling and measurements: International Ocean Discovery Program (IODP) Site U1467 (73.17° E, 4.51° N, 487.4 m water depth), Site MD96-2060 (40.46° E, 8.30° S, 2026 m), Site NIOP 905 (51.95° E, 10.76° N, 1575 m) and Site SO189-39KL (99.90° E, 0.78° S, 517 m) (Fig. 1). These sites span distinct oceanographic settings influenced by different modes of variability and regional water-mass dynamics.

IODP Site U1467 located in the central tropical Indian Ocean (CTIO) is governed by the seasonal reversal of monsoonal winds under the strong influence of the south Asian monsoon system (Wang et al., 2005). These wind reversals strongly modulate surface ocean currents, SST, salinity, and upwelling intensity in the western Arabian Sea. During boreal summer, strong southwesterly winds drive the eastward-flowing Summer Monsoon Current. In contrast, during boreal winter, the Winter Monsoon Current flows westward, transporting low-salinity waters from the Bay of Bengal to the Arabian Sea (Fousiya et al., 2016; Behara and Vinayachandran, 2016). Thermocline depth is highly sensitive to IOD phases especially in the western Indian Ocean, shoaling during negative events and deepening during positive years (Sankar et al., 2021).

Site NIOP 905, in the western Arabian Sea (WAS) is located in the most intense upwelling system in the Indian Ocean. The upwelling system is primarily driven by the Somali Jet during the southwest (summer) monsoon, which enhances alongshore winds off Oman and induces coastal and offshore upwelling via Ekman transport and the Somali coastal current (Schott and McCreary, 2001; Hood et al., 2017). During the northeast monsoon, upwelling is reduced promoting higher SST (Izumo et al., 2008), leading to high seasonality in SST (9–11 °C). As compared to other sites, this site is most affected by upwelling, with intense primary productivity occurring during most of the year culminating during summer (Wiggert et al., 2005) and a shallower thermocline depth.

Site MD96-2060, in western tropical Indian Ocean (WTIO), offshore Tanzania, lies within the East African Coastal Current (EACC). The EACC exhibits clear seasonal modulation by monsoonal winds: it strengthens during the southwest monsoon (April–September) and weakens during the northeast monsoon months (November–February) (Mayorga-Adame et al., 2017; Vinayachandran et al., 2021). Thermocline temperatures exhibit a seasonal range of ∼1–2 °C.

Site SO189-39KL, located in the eastern tropical Indian Ocean (ETIO), is situated within the Sumatra coastal upwelling system, which exhibits strong seasonal intensification during the southwest monsoon (Susanto et al., 2001). Seasonal variability in SST and Sea surface salinity (SSS) is nonetheless recorded and is controlled by the migration of the Intertropical Convergence Zone (ITCZ), not by the monsoon or by ENSO directly (Aldrian and Susanto, 2003). Within the thermocline, temperature fluctuations are mainly influenced by IOD events (Yu et al., 2005; Qiu et al., 2012).

Together, these sites provide a comprehensive view of the spatial heterogeneity of Indian Ocean variability, encompassing monsoon-dominated northern regions, the western boundary current system, and the IOD and upwelling-influenced eastern equatorial zone. Their contrasting oceanographic settings allow us to test the sensitivity of δ¹⁸O_c variability to seasonal and interannual forcing and to evaluate the performance of IFA forward-modelling across diverse climatic regimes.

Figure 1 Modern seasonality in SST and SSS, and estimated δ¹⁸O_c in the tropical Indian Ocean at surface and thermocline depths. (A)–(C) Annual standard deviation of sea surface temperature (SST), sea surface salinity (SSS), and δ¹⁸O_c at the surface (0.5 m depth). (D)–(F) Annual standard deviation of temperature, salinity and δ¹⁸O_c for thermocline depth (97 m). δ¹⁸O_c was calculated using species-specific paleotemperature equations (Farmer et al., 2007; Eqs. 1–2) and δ¹⁸O_w–S relationships (Eqs. 3–4). Calculations were performed based on the ORAS5 reanalysis dataset from 1958 to 2018.

3 Methods

3.1 Forward-modelling of planktonic IFA δ¹⁸O_c

3.1.1 Calculating surface and thermocline δ¹⁸O_c time series from ocean reanalysis

To simulate δ¹⁸O_c variability at the selected sites, we used the ORAS5 ocean reanalysis dataset, which provides monthly fields of temperature and salinity at 1° resolution for the period 1958–2018. For each site, we extracted time series at two depths representative of the surface mixed-layer and thermocline conditions in the Indian Ocean (Fig. 2A): 0.5 m for surface conditions, and 97 m for the thermocline for all sites. The choice of 97 m for thermocline depth was made based on available depths in the ORAS5 analysis dataset, and which of these most closely matched the depth of the highest SD extracted from temperature profiles at our four sites (∼100 m depth; Fig. 2B). Depth profiles of temperature and salinity show similar thermocline depths for all four sites, despite being influenced by contrasting oceanographical regimes (Fig. 2A). We also compared ORAS5 to three other products, HadISST1.1 (Rayner et al., 2003), GODAS (Behringer and Xue, 2004) reanalyses, and the EN4 (Good et al., 2013) instrumental dataset for surface and thermocline temperature and salinity at the four proxy sites from 1980 to 2018 (Fig. S1 in the Supplement). ORAS5 reproduces the seasonal cycle of surface temperature with near-perfect to high-confidence agreement at all sites, and captures interannual SST variability well (r²=0.11–0.68 across sites, p<0.001). Thermocline-depth temperature likewise shows a significant correlation with all tested products at every site (p<0.001). Interannual salinity correlations are generally low, but consistent with known reanalysis limitations in the tropical Indian Ocean (Copernicus Climate Change Service, 2021) and subsurface salinity contributes only marginally to modelled δ¹⁸O_c variance. Overall, this validation confirms that ORAS5 faithfully reproduces large-scale temperature and salinity conditions at the four proxy sites, supporting its use in our forward-modelling framework.

While reanalysis products such as ORAS5 provide temperature and salinity at multiple discrete depth levels, we favour a layer-integrated extraction (0–50 m for the surface mixed layer, 70–120 m for the thermocline) over a single fixed depth where possible, as it better represents the vertical range over which G. ruber and N. dutertrei are known to calcify. We verified that this choice has a negligible effect on the simulated δ¹⁸O_c variability at three of our four sites (ΔSD < 0.065 ‰; Fig. S2). At SO189-39KL, however, layer-integration measurably dampens the near-surface interannual signal that is central to detecting extreme events at this upwelling-influenced site; we therefore retain the fixed 0.5 m depth specifically at this site to preserve this signal, and discuss both approaches where relevant.

The mean thermocline depth, defined as the depth of maximum vertical temperature gradient (|d$T/$ dz|) derived from ORAS5 was 83.2 ± 30.7 m at MD96-2060, 103.3 ± 37.9 m at NIOP 905, 85.1 ± 17.6 m at U1467, and 90.2 ± 27.2 m at SO189-39KL. The depth of 97 m falls within one standard deviation of the mean thermocline depth at all four sites, justifying its use as the calcification depth for N. dutertrei.

Simulated δ¹⁸O_c values were computed using species-specific paleotemperature equations: the G. ruber equation of Farmer et al. (2007) (1) for surface simulations, and the N. dutertrei equation of Farmer et al. (2007) (2) for thermocline simulations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }G.\mathit{ruber}:T=\mathrm{15.4}-\mathrm{4.78}\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{w}}\right)\text{(2)}& {\displaystyle }& {\displaystyle }N.\mathit{dutertrei}:T=\mathrm{14.6}-\mathrm{5.09}\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{w}}\right)\end{array}$$

The corresponding δ¹⁸O sea water (δ¹⁸O_w) values were estimated from salinity using depth-specific δ¹⁸O_w–S relationships. For the surface, we chose to apply the same equation for all sites to produce comparable changes under an altered climate. We applied the equation of Thirumalai and Clemens (2020), which has been adapted from Singh et al. (2010), which is similar to the δ¹⁸O_w–S relationships commonly found in tropical surface waters (Schmidt et al., 2007) (Eq. 3):

$$\begin{array}{}\text{(3)}& \mathrm{Surface}:{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{w}}=\mathrm{0.26}S-\mathrm{8}\end{array}$$

This relationship was adapted from the coastal Arabian Sea regression of Singh et al. (2010) as applied by Thirumalai and Clemens (2020) (R²=0.62, p<0.01, n=152; slope: 0.26 ± 0.02 ‰ psu⁻¹; intercept: $-\mathrm{8.0}\pm \mathrm{0.5}$ ‰). It is applied uniformly across all sites in the variability sensitivity analysis (step 1 of the forward model) to ensure that differences in simulated δ¹⁸O_c reflect contrasting local oceanographic conditions alone.

For the subsurface (75–125 m), we derived a new relationship based on the global dataset of LeGrande and Schmidt (2006), (R²=0.931, p<0.001, n=2073) restricted to this depth range (Fig. S3; Eq. 4):

$$\begin{array}{}\text{(4)}& \mathrm{Thermocline}:{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{w}}=\left(\mathrm{0.450}\pm \mathrm{0.003}\right)S-\left(\mathrm{15.513}\pm \mathrm{0.104}\right)\end{array}$$

Figure 2 Modern vertical oceanographic profiles and temperature variability at the four study sites. Top panels: Mean vertical profiles of temperature (red) and salinity (blue) for each site. The shaded area indicates the depth range of the thermocline (defined on the increase of temperature gradient). The dashed line indicates the thermocline depth and the dotted line the extraction depth of 97 m. Bottom panel: Temperature standard deviation (°C) versus depth (0–300 m) derived from ORAS5 reanalysis with the location of the sites represented by coloured stars.

3.1.2 Simulating synthetic planktonic foraminifera IFA core-tops

The temperature and salinity timeseries extracted from ORAS5 are used to identify how increasing variability will be imprinted in the surface and subsurface Indian Ocean when the corresponding δ¹⁸O_c values are calculated. We therefore opted to use a constant δ¹⁸O_w–S relationship to better illustrate how modifying seasonal and interannual temperature and salinity values would affect an idealized, continuous synthetic δ¹⁸O_c timeseries. However, for the simulation of synthetic IFA data, to allow for direct comparison with measured data from the same locations, we use region-specific δ¹⁸O_w–S relationships. In the second step of the forward model, where synthetic IFA distributions are compared directly with measured core-top IFA data, site-specific regional equations are applied as reported in the original publications when available, to minimise systematic biases. This differs from the first step (Sect. 3.1.1), in which a uniform δ¹⁸O_w–S relationship is applied across all sites to ensure comparability of the sensitivity analysis (Fig. S4). The use of a constant relationship is further justified by our sensitivity analysis (Fig. S5), which demonstrates that the choice of δ¹⁸O_w–S equation has a negligible effect on the simulated δ¹⁸O_c SD (ΔSD < 0.023 ‰) across all sites and equations tested.

Based on the modelled surface and thermocline δ¹⁸O_c timeseries using the ORAS5 monthly reanalysis data for any coordinate and depth, we carried out random IFA sampling (60 individuals) across the 1958–2018 period (732 months) following the procedure of Thirumalai et al. (2013). Each individual foraminifera is assumed to record the average monthly environmental conditions of its calcification period. This assumption is justified by the 2–4-week lifespan of both G. ruber and N. dutertrei (Spero, 1998; Schiebel and Hemleben, 2017), consistent with lunar-paced reproduction documented for both species (Bijma et al., 1990; Jonkers et al., 2015), which is shorter than or comparable to the monthly resolution of the ORAS5 reanalysis, such that each individual is considered to have calcified within a single monthly interval. By sampling randomly throughout the 60-year timeseries, we simulate a continuous and constant shell flux, effectively integrating the climatic signal while avoiding any assumption of a seasonally phase-locked life cycle, meaning we do not prescribe a fixed monthly calcification window tied to the reanalysis calendar. To assess the influence of altered climate variability on δ¹⁸O_c, we artificially increased (or decreased) the amplitude of seasonal variability in surface waters by ± 75 %. This range of amplitude was selected as a sensitivity benchmark to simulate extreme monsoonal or IOD/ENSO regimes, consistent with the values used in a previous study (Thirumalai and Clemens, 2020). For the thermocline, we modified the amplitude and frequency of interannual variability by injecting random events similar to potential IOD/ENSO events that create anomalies in which temperature and salinity are modified by ± 75 %. Although an increase in SD is expected by construction, this test quantifies in δ¹⁸O with propagated uncertainty, how much of an imposed change in variability is recovered by a 60-specimen IFA sample at each site, and provides the baseline for the model-data comparison of Step 2.

This choice is further supported by analysis of peak ENSO and IOD events in the ORAS5 reanalysis (1958–2018), which shows that thermocline temperature anomalies during the strongest events range from 2.78 to 7.11 °C across the four sites, approximately 50 %–75 % of the mean seasonal amplitude at each site (Figs. 3, S6). The ± 75 % perturbation therefore likely represents a physically meaningful upper bound for altered variability states. We tested three configurations at the surface (0–50 m) across all four sites: seasonal-only, interannual-only, and combined perturbations. At MD96-2060 and NIOP 905, interannual changes produce nearly identical or lower IFA SD than seasonal (0.57 ‰ vs. 0.58 ‰ and 0.57 ‰ vs. 0.69 ‰ respectively), confirming that seasonal variability dominates surface IFA dispersion at these western sites. At U1467 and SO189-39KL, interannual variability dominates over seasonality (0.54 ‰ vs. 0.43 ‰ and 0.48 ‰ vs. 0.30 ‰), reflecting a stronger ENSO/IOD surface signature at these central and eastern sites. A combined seasonal and interannual perturbation increases IFA SD by 0.17 ‰–0.25 ‰ across sites (Fig. S7).

For each altered and original dataset, we generate a single representative draw of 60 values randomly sampled with replacement to simulate an IFA-style core-top sampling with a statistically robust SD (Leduc et al., 2009). To assess whether this SD value reliably represents the underlying population rather than an artefact of that particular draw, we additionally applied a Monte-Carlo bootstrap procedure, repeating the same 60-individual sampling 1000 times for each dataset and scenario and recording the resulting SD at each repetition (Fig. S8). We consider that the resulting SD estimate is statistically robust in the sense that the uncertainty is small relative to the effect sizes interpreted, the SD standard error across the 1000 bootstrap repetitions is typically ± 0.02 ‰ to 0.05 ‰, much smaller than the SD values themselves (0.19 ‰ to 0.70 ‰). To propagate uncertainty in the forward-modelling framework itself, the paleotemperature equation coefficients and δ¹⁸O_w–S equation parameters (slope and intercept) were resampled from a normal distribution defined by their reported standard errors at each Monte Carlo iteration, rather than held fixed at their mean values. For each site and depth, the 1000 SD estimates obtained under the original and altered scenarios differ significantly (paired Wilcoxon signed-rank test, p<0.001), and the difference between their means exceeds the spread of the bootstrapped SD, indicating that IFA can detect such changes in variability.

We finally evaluate the effect of these scenarios using Quantile–Quantile (Q–Q) plots, which jointly compare the distributions of δ¹⁸O_c values of the original and altered datasets. Q–Q plots are useful for visualizing both systematic and non-linear deviations caused by changes in temperature and salinity. An alignment along the 1:1 reference is expected if two IFA distributions are statistically the same. Any deviations from this line will reflect changes in the background (without changes in the slope), magnitude (alignment along a different slope) or with other recognizable structures of changing total variability, such as skewness (arcuate-shaped Q-Q plot). This method has already proven useful for visualizing how interannual variability was altered in previous IFA studies (Ford et al., 2015; Thirumalai et al., 2019).

3.2 IFA Measurements and Statistics

We present new measurements of IFA δ¹⁸O in core-top samples (Late Holocene, 0–4 ka) from Sites MD96-2060, IODP U1467, and NIOP 905. We selected two species with distinct ecological niches (Chaabane et al., 2026): Globigerinoides ruber is a surface-dwelling, symbiont-bearing species that calcifies throughout the year within the upper mixed layer (Fairbanks et al., 1982; Stainbank et al., 2019) and Neogloboquadrina dutertrei is a thermocline-dwelling species commonly found in the deep chlorophyll maximum (DCM), associated with the pycnocline (Chaabane et al., 2026; Kemle-von Mücke and Oberhänsli, 1999). N. dutertrei is particularly useful in reconstructions of interannual variability, due to its sensitivity to thermocline depth changes linked to phenomena such as ENSO and IOD (Leduc et al., 2009; Thirumalai et al., 2013; Ford et al., 2015). New δ¹⁸O data are measured on G. ruber and N. dutertrei for site MD96-2060, N. dutertrei for site U1467 and N. dutertrei for site NIOP 905 (Table 1). We also use published G. ruber data from Site NIOP 905 (Ganssen et al., 2011) and IODP Site U1467 (Stainbank et al., 2021).

Planktonic foraminifera were picked from the 250–315 µm size range. New δ¹⁸O_c analyses were conducted on individual specimens (15–20 µg, up to 40 µg for N. dutertrei) using a Thermo-Finnigan Delta V-Plus IRMS coupled with a Carbo-Kiel IV carbonate preparation device at CEREGE. The international standard NBS-19 was used for calibration (30–50 µg), and δ¹⁸O values are reported in per mil (‰) relative to the Vienna Pee Dee Belemnite (VPDB) standard. The mean external reproducibility (1σ) was better than 0.08 ‰ for δ¹⁸O. Measurements yielding a mass-44 signal below 1500 mV, or an internal δ¹⁸O standard deviation above 0.05 ‰ over the six dual-inlet cycles, were rejected.

To assess the similarity between the IFA distributions observed in core-tops and those generated by forward-modelling, we applied several statistical tests. These included assessments of normality (e.g., Shapiro-Wilk), and comparisons of the means, SD, and distribution (skewness, kurtosis). Classical significance tests were also applied to assess whether differences between datasets (modelled vs. measured) were statistically significant. Outliers were identified using the interquartile range (IQR) method, as values outside Q1 − 1.5 × IQR and Q3 + 1.5 × IQR, where IQR = Q3 − Q1. This method has been successfully applied in Ganssen et al. (2011) on core NIOP 905, and only 5 individual analyses among a total of 139 measurements for the newly generated IFA data were identified as outliers (Table S1 in the Supplement). An independent criterion based on the median absolute deviation (modified z-score > 3.5; Iglewicz and Hoaglin, 1993) identifies the same five values. Converted into temperature, the most extreme one, at U1467, would imply thermocline temperatures ∼12 °C colder than the site mean, well beyond any interannual anomaly at this site. We therefore consider these values unlikely to reflect real environmental extremes and excluded them from the statistical analyses. Core-tops dating and age models applied are reported in Table 1 and described in Ivanochko et al. (2005) for NIOP 905, Mohtadi et al. (2010, 2011, 2014) for SO189-39KL and Stainbank et al. (2020) for Site U1467, and Dauvier et al. (2026) for MD96-2060 (Table 1).

Another critical consideration in interpreting IFA variability is the potential for dissolution of the calcite tests analysed. First, the bathymetric setting of the cores ranging from 487 m at Site U1467 to 2026 m at Site MD96-2060 places all samples well above the modern Indian Ocean lysocline, typically situated between 3500 and 3800 m (Feely et al., 2004; Peterson and Prell, 1985). Second, visual inspection of G. ruber and N. dutertrei tests reveals good preservation, consistent with SEM observations of pristine calcite and clearly apparent spines at other well-preserved sites in the region (Wang et al., 2013). This is further corroborated by an independent dissolution-proxy dataset for the tropical Indian Ocean (Adebayo et al., 2025), which reports a Fragment Variance Index of 0.92 and a fragmentation percentage of only 3 % at MD96-2060 among the best-preserved core-tops of their 62 sites compilation. Consequently, we conclude that the observed IFA distributions represent primary climate signals rather than a signature of differential preservation across the varying water depths of the sites that would lead to difference in homogenisation and altered variability.

Table 1 Summary of core-top sediment sites and associated Individual Foraminifera Analysis (IFA) datasets: location, water depth, age models, and data provenance (new vs. published) for G. ruber and N. dutertrei.

4 Results

4.1 Forward modelled synthetic time-series and IFA core-top δ¹⁸O_c under modern conditions

4.1.1 Time series extracted from reanalysis dataset

At the seasonal timescale, SST exhibits the highest variability at sites NIOP 905 and MD96-2060, with SST amplitude ranging from 23 °C in boreal summer to 29 °C in spring and 25 °C in summer to 29 °C in spring, respectively. In contrast, the other two sites show a weaker SST seasonality, with Site U1467 showing seasonal variability of 2–3 °C and SO189-39KL less than 1 °C (Fig. 3B). Seasonal SSS variability is relatively low for Sites MD96-2060 and NIOP 905 (<1 psu), and higher for Sites SO189-39KL (1.5 psu) and U1467 (2 psu) (Fig. 3B). This variability translates into a modelled δ¹⁸O_c variability that increases in line with SST seasonality, exhibiting maximum amplitude of 1.5 ‰ at Site NIOP 905. Only the eastern Indian Ocean site shows seasonal surface δ¹⁸O_c variability of less than 0.5 ‰.

SST and δ¹⁸O_c variability are higher at the seasonal scale for the surface than the thermocline for all sites. Overall, results highlight that the drivers of surface δ¹⁸O_c are mainly related to temperature variability in the western basin whereas salinity becomes more significant for the Maldives and in the east. Interannual variability at the surface is relatively low (<1.5 °C in SST and <1 psu in SSS), resulting in δ¹⁸O_c fluctuations that do not exceed 0.5 ‰. Maximum interannual variability is observed at eastern Indian Ocean Site SO189-39KL (Fig. 3B). At most of the sites, surface water variability is mainly driven by seasonal variability and theoretically the most sizeable interannual signal is found within the thermocline where we find the highest temperature SD (Fig. 2). The interannual standard deviation of monthly δ¹⁸O_c at the surface confirms these patterns (Fig. 3B). Site NIOP 905 displays the largest interannual SD in δ¹⁸O_c (up to 0.21 ‰), with site U1467 (up to 0.19 ‰) and site SO189-39KL (up to 0.17 ‰). Site MD96-2060 shows intermediate SD values (0.06 ‰–0.12 ‰), primarily driven by interannual SST variability. Overall, interannual SD values at the surface remain well below seasonal amplitudes at all sites, confirming that seasonal variability is the dominant driver of surface δ¹⁸O_c distributions.

Interannual thermocline (97 m) temperatures at interannual timescale (thick curves on Fig. 3C) show the strongest amplitudes at sites MD96-2060 (4.8 °C) and SO189-39KL (6.3 °C), whereas Sites U1467 and NIOP 905 display lower amplitudes (3.6 and 2.7 °C). Salinity at thermocline depth varies between 0.2 and 0.5 psu at all sites, with a slightly higher amplitude (0.52 psu) at Site SO189-39KL. The simulated thermocline δ¹⁸O_c exhibits a maximum amplitude of 1 ‰ at MD96-2060 and 1.3 ‰ at SO189-39KL, followed by 0.64 ‰ at U1467 and 0.48 ‰ at NIOP 905. This indicates that sites SO189-39KL (east) and MD96-2060 (west) are the most sensitive to interannual thermocline variability. At the seasonal timescale (thin curves, Fig. 3D), the maximum amplitudes are found in sites NIOP 905 and SO189-39KL. Over the instrumental period (1958–2018), the superposition of interannual variability on the seasonal cycle increases the temperatures range to 11 °C at site NIOP 905. This strong total variability results in a calculated δ¹⁸O_c amplitude of up to 2.78 ‰. The interannual standard deviation of δ¹⁸O_c at thermocline depth further illustrates the contrast between sites (Fig. 3D). Site SO189-39KL displays the largest SD (∼ 0.55 ‰), consistent with the strong interannual modulation of thermocline temperature by the Indonesian Throughflow and ENSO forcing. Site MD96-2060 also shows elevated thermocline δ¹⁸O_c SD (∼ 0.35 ‰). In contrast, Sites NIOP 905 and U1467 display lower interannual SD values (∼ 0.13 ‰–0.25 ‰ and ∼0.19 ‰–0.27 ‰ respectively).

Figure 3 Modern surface and subsurface (thermocline) variability at the four studied sites (U1467, NIOP 905, MD96-2060 and SO189-39KL), derived from ORAS5 reanalysis datasets (1958–2018). (A) and (C): Monthly time series of δ¹⁸O_c, salinity, and temperature at both the surface (top) and thermocline depth (bottom). (B) and (D) The right panels show the corresponding monthly mean values (averages for each month over the 1958–2018 period) for each variable.

At all sites, modern climatology data show that seasonal variations dominate in the surface ocean, with relatively lower-amplitude interannual timescale variability (Fig. 3A). At the depth of the thermocline, all sites show increased sensitivity to interannual variability relative to at the surface (Fig. 3C). Site SO189-39KL shows the strongest interannual variability in temperature (and therefore in modelled δ¹⁸O_c), followed by site MD96-2060. At sites NIOP 905 and U1467, interannual thermocline variability is much smaller than the dominant seasonal cycle, with maximum interannual SD during boreal summer reaching 1.6 °C at NIOP 905 (ORAS5, 1980–2018). The seasonal amplitude (∼ 11 °C at NIOP 905) is larger than the interannual signal, and seasonal variations in SST and SSS at thermocline depths dominate the total variability (Fig. 3D). These variations are interpreted in the discussion below by linking our data with modes of variability and measured δ¹⁸O-IFA.

4.1.2 Artificial modern amplification of seasonal and interannual variability and its effect on forward-modelled IFA

The temporal evolution of the calculated δ¹⁸O_c from the reanalysis datasets (1958–2018) illustrates the theoretical proxy response to an increase in high-frequency variability by 75 % (Fig. 4). For each site, we compare δ¹⁸O_c timeseries of ORAS5 (blue curve), and the altered timeseries scenario (red curve), where seasonal or interannual variability was artificially amplified or reduced to test IFA sensitivity.

At the surface (Fig. 4A–D), we only change the seasonal variability by amplifying the annual oscillation of δ¹⁸O_c, which is highly regular and characterised by a relatively high amplitude increase of 75 % (sinusoidal pattern). The “Altered” scenario (Fig. 4A–D) clearly shows an amplification of seasonal extremes. Conversely, at the thermocline depth (Fig. 4E–H), the seasonal signal remains unchanged and interannual anomalies appear, increasing/decreasing the total amplitude of variability. The “Original” and “Altered” curves overlap in a more complex manner, highlighting interannual anomalies.

Figure 4 Comparison of monthly reanalysis dataset for an “original” (blue curve) and “altered” (red curve) climate scenarios for the surface and thermocline. δ¹⁸O_c calculated from monthly temperature and salinity from 1958 to 2018. (A)–(D) Surface reanalysis dataset for Sites MD96-2060, NIOP 905, U1467 and SO189-39KL respectively. (E)–(H) Thermocline reanalysis dataset for the same sites.

We turn then to the results of the IFA picking under the different scenarios to simulate synthetic IFA core-tops. We simulate an increase of seasonality for surface, that we will associate with G. ruber in our IFA analysis and an increase of the interannual variability for the thermocline that we associate with N. dutertrei for our IFA. Forward-modelled IFA δ¹⁸O show an increase in the SD at each site for both depths in an altered climate (Fig. 5). A higher SD indicates a greater dispersion of values and increased variability in the simulated δ¹⁸O-IFA values. At the surface, higher SD is observed at Sites NIOP 905 and MD96-2060, with 0.38 ‰ and 0.35 ‰ respectively for the original scenario, versus 0.69 ‰ and 0.61 ‰ respectively for the altered scenario. Site U1467 shows an increase in SD from 0.28 ‰ to 0.49 ‰, and Site SO189-39KL from 0.19 ‰ to 0.29 ‰. At thermocline depth, SD at Site SO189-39KL increases from 0.63 ‰ for the original dataset to 0.76 ‰ for the altered dataset. Site NIOP 905 shows an increase from 0.37 ‰ to 0.49 ‰, Site MD96-2060 from 0.34 ‰ to 0.46 ‰ and Site U1467 from 0.28 ‰ to 0.42 ‰.

In summary, forward-modelled IFA δ¹⁸O values respond coherently with an increase in SD under altered scenario, when seasonality at the surface or interannual variability at the thermocline is increased. A dedicated sensitivity test further confirms that adding interannual perturbations at the surface produces effects of comparable or greater magnitude than seasonal perturbations at the central and eastern sites U1467 and SO189-39KL (see Sect. 3.1.2), while seasonal perturbations remain dominant at the western sites MD96-2060 and NIOP 905 (Fig. S7). Also, extreme IFA values become more frequent when overall variability increases, indicating that stronger interannual modes generate more anomalous conditions. We observe a significant difference between the western and eastern Indian Ocean for all scenarios with a higher surface variability for the western site and reduced variability for the eastern site, and the opposite tendency for thermocline variability with an increase for the east and a reduction for the west. Sites SO189-39KL and NIOP 905 seem to be the most sensitive to δ¹⁸O changes at thermocline depth.

Figure 5 Modelled IFA δ¹⁸O_c (synthetic core top) values under original and altered climate variability scenarios at four sites in the Indian Ocean. Left panel: Surface modelled IFA variability associated with surface-dwelling G. ruber; right panel: Thermocline-depth modelled IFA variability associated with N. dutertrei. Two distributions are shown: the original modelled IFA δ¹⁸O values (labelled “Original”) and values obtained after artificially enhancing climate variability (“Altered”). The red and blue arrows highlight the relationship between δ¹⁸O values and temperature (warmer = more negative values, colder = more positive values). SD for each distribution is reported for each dataset.

4.2 Core-top IFA analysis

New data include an IFA dataset for the surface at Site MD96-2060 and three IFA datasets for the thermocline at Sites MD96-2060, U1467 and NIOP 905. δ¹⁸O measurements on G. ruber IFA from Site MD96-2060 (n=46) yield a mean value of −1.75 ‰ and a SD of 0.33 ‰ (Fig. 6). At the same site, new N. dutertrei IFA δ¹⁸O values (n=54) display a mean of −1.34 ‰ and a SD of 0.30 ‰ (Fig. 6A), with two outliers identified and removed prior to statistical analysis. Similarly, N. dutertrei measurements from Site U1467 (n=44) have a mean value of −1.68 ‰ and a SD of 0.30 ‰, with one outlier excluded from the dataset and G. ruber data from this site give a mean of −2.75 ‰ and a SD of 0.27 ‰.

At Site NIOP 905, we use the G. ruber dataset published in Ganssen et al. (2011) along with new N. dutertrei data measured on the same core-top (Gerald Ganssen, unpublished data). The G. ruber IFA isotopic values have a mean δ¹⁸O of −1.54 ‰ and a SD of 0.56 ‰, excluding 2 outliers. The N. dutertrei data (n=38) are characterised by a mean δ¹⁸O value of −0.01 ‰ and a SD of 0.53 ‰ (Fig. 6B). Finally, the SO189-39KL core reports a mean G. ruber IFA δ¹⁸O of −2.94 ‰, and a SD of 0.21 ‰ (Thirumalai et al., 2019). All dataset statistics are available in Table S1.

5 Discussion

5.1 Interpretation of IFA- δ¹⁸O_c datasets total variability signal

At the surface, the G. ruber IFA SD at sites MD96-2060 (0.33 ‰) and U1467 (0.27 ‰) do not differ significantly (F-test, p=0.12; Fig. 6). SO189-39KL (0.21 ‰) differs significantly from MD96-2060 (p=0.002) but only marginally from U1467 (p=0.055). NIOP 905 (0.56 ‰) differs significantly from all three other sites (p<0.001 in each case). SO189-39KL is located in the warmest water and shows the lowest variability, whereas NIOP 905 is the coldest and has the highest variability. The opposite pattern emerges when considering salinity, where the highest salinity is found in WAS. Thus, the most and least variable IFA datasets appear to reflect the regional extremes of temperature and salinity observed in the ORAS5 reanalysis. WAS has a SD significantly higher than all the other sites (Fig. 6), indicating that this area is somehow successfully capturing the expected total variability (Fig. 1C). This likely reflects the influence of a highly dynamic seasonal system influenced by monsoonal winds (Schott and McCreary, 2001) and its associated Ekman pumping during the summer upwelling season, as originally suggested in Ganssen et al. (2011). The winter monsoon counterpart of this dynamical upwelling system may also affect the CTIO, in the Maldives archipelago (Bassinot et al., 2011). However, Site U1467 is supposed to be impacted in the same way as the WTIO (Fig. 1), which is less influenced by seasonal wind-induced upwelling and/or mixing. This suggests that, even though they show a similar SD the processes explaining their variability are not the same and reflect regional influences. The permanent upwelling in ETIO that produces the lowest IFA variability (site SO189-39KL, SD of 0.21 ‰) may indicate that the modern thermocline depth is so deep (Fig. 6D) that, despite potential changes in upwelling activity, the surface IFA SD is very low at that site. The ETIO off Sumatra is less influenced by monsoonal dynamics, as suggested by the stable and warm SST (Fig. 3). Here, seasonal and interannual movements of the ITCZ in the Indonesian archipelago induce significant changes in the SSS (Fig. 3), however these are not sufficient to significantly increase the IFA SD in this location (Figs. 4, 6).

Figure 6 Vertical δ¹⁸O_c profiles and IFA measured on coretops across the Indian Ocean. The vertical δ¹⁸O_c profiles (dashed lines) represent synthetic calcite values derived from modern instrumental reanalysis data, calculated using specific δ¹⁸O_c equations (1) and (2) from Farmer et al. (2007) and regional δ¹⁸O_sw-salinity relationships, using the relationship presented in equations (3) and (4) for the surface and subsurface water depth, respectively. The shaded envelope around each profile shows ± 1 standard deviation of the modelled δ¹⁸O_c, reflecting its interannual and seasonal variability at each depth (1958–2018). These hydrographic profiles are compared with Individual Foraminifera Analysis (Boxplots) results obtained from core-top sediments, representing the surface (0.5 m, G. ruber) and the thermocline (97 m, N. dutertrei). Boxes span the first to third quartiles with the median as a black line; whiskers extend to 1.5 × IQR; red dots are outliers. The gray shaded area represents the thermocline determined by a change in the temperature gradient. Yellow stars spot newly generated IFA datasets for this study.

At thermocline depth, IFA on N. dutertrei reveal a similar amplitude in the WTIO and CTIO at Sites MD96-2060 and U1467 (SD of 0.30 ‰, Fig. 6A, C). As in the surface, WAS records the highest subsurface SD (0.53 ‰). The pattern of variability at this site (Fig. 4A, C) confirms that IFA results are mainly driven by seasonality even in subsurface. This suggests a deeply integrated response to seasonal (monsoonal) forcing potentially suggesting a well-mixed water column from the surface to the thermocline depth. In the WAS and CTIO, interannual thermocline variability is smaller relative to the dominant seasonal cycle. At NIOP 905, thermocline temperature anomalies reach up to 5.4 °C during strong IOD and ENSO events with maximum interannual SD during boreal summer (JJA: 1.6 °C), consistent with the known modulation of Somali Current intensity and coastal upwelling by IOD forcing (Schott and McCreary, 2001; Beal and Donohue, 2013).

We note that the most negative mean and the lowest IFA variability results come from the warm, low-salinity in ETIO (Figs. 1, 3), despite the equatorial upwelling that prevails year-round (Fig. 6). In contrast, the highly dynamic upwelling system offshore Somalia coincides with the core-top where the highest SD in IFA is reported (Fig. 6). Overall, these deviations from the mean δ¹⁸O_c profile as well as the differences in the SD at different sites imply that multiple dynamical and ecological processes are at play such as monsoon-driven upwelling intensity (WAS/NIOP 905), thermocline depth (ETIO/SO189-39KL) and species-specific depth habitat and seasonal production timing, which need to be accounted for in the interpretation of planktonic foraminifera IFA.

5.2 Assessing proxy sensitivity through Forward-modelled IFA “Original” and “Altered” variability scenarios

To evaluate the sensitivity and spatial variability of IFA-derived δ¹⁸O_c signals, we modelled four synthetic scenarios in which only salinity or temperature variability was either amplified or reduced. The original scenario aims to represent a synthetic core-top based on a random sampling of the reanalysis datasets presented in Fig. 4. The altered scenario is a sensitivity test that simulates a synthetic sample representing a time window when increased and decreased seasonality or interannual variability is imposed to the reanalysis dataset as shown in Fig. 5. This sensitivity test translates these relationships into proxy δ¹⁸O_c with propagated uncertainty which provide the baseline needed to interpret the Step 2 model-data comparison in Fig. 8.

At the surface, the modelled results reveal a distinct spatial pattern of sensitivity to seasonal and interannual variability (Fig. 1). CTIO at site U1467, shows a coherent difference in the slopes between increased and reduced scenarios (Fig. 7A). Both changes in the SSS and SST annual cycle influence the IFA-δ¹⁸O, especially in the extreme high quantile of the distribution. This site, located in the path of the seasonally reversing West India Coastal Current (WICC) that transports fresher surface water from the Bay of Bengal to the Arabian Sea during winter and saline water from the Arabian Sea eastwards during summer (Shankar et al., 2002), is sensitive to seasonal changes in SSS. The amplitude of SSS seasonality is set by the monsoonal system (Talley et al., 2011). Moreover, it also drives seasonal changes in temperature linked to upwelling dynamics during the southwest monsoon when Ekman pumping that advects waters upwelled from sites situated in the WAS region is stronger, reducing SST in CTIO.

The Quantile-Quantile plot of Site NIOP 905, located offshore Somalia, south of the Arabian Sea, shows a dominant temperature influence on IFA-δ¹⁸O (Fig. 7B). The seasonal variation in SSS is less than 1 psu, whereas SST that reaches a maximum amplitude of ∼ 5.6 °C (Fig. 3B). This major temperature variability is associated with the strong seasonal upwelling during summer monsoon that leads to the shoaling of cold water-masses in this area (Vinayachandran et al., 2021). A second SST drop during the winter monsoon is also observed and is associated with wind-induced mixing (Fig. 3B), interspaced by periods of warm-water advection toward the NIOP 905 site. The IFA-δ¹⁸O surface variability in WTIO at Site MD96-2060 is also mainly temperature-driven (Fig. 7C). The site is under the influence of a seasonally reversing surface current, the East African Coastal Current (Schott and McCreary, 2001; Schott et al., 2009), which is situated near the northern limit of the South Equatorial Current. Its boundary oscillates between 6 and 9° S in response to the South Indian Ocean trade winds. In ETIO, site SO189-39KL, IFA-δ¹⁸O shows both invariant temperature and salinity at seasonal timescales (Fig. 1). Not only is the total IFA SD of the synthetic core-top IFA halved compared to other sites (Fig. 7D), but also an artificial increase or decrease in its variability at the seasonal scale fails to clearly modify the slope of the Q-Q plots distributions.

Decomposition of the IFA-δ¹⁸O signals at thermocline depth shows that they are mainly driven by temperature changes at all sites (Fig. 7E–H), except for CTIO, where salinity has a greater influence (Fig. 3, 7E). Although the seasonal temperature cycle remains significant in the thermocline, the interannual variability manifests much more in the thermocline variability as compared to that of surface waters (Fig. 3). Based on the climatological data, most of the IFA-δ¹⁸O variations at interannual scale are driven by temperature variability (Fig. 3). Given that salinity fluctuations remain minimal across all sites (less than 0.5 psu), we conclude that the influence of salinity on the interannual IFA-δ¹⁸O signal is secondary, allowing thermocline-dwelling species to be used primarily as a proxy for temperature-driven variability at these timescales. WTIO appears to be the most promising site for recording interannual variability with IFA at the thermocline depth. Many studies show that this area is very sensitive to interannual events such as IOD, influencing upper water column temperature and sea surface salinity (Yu et al., 2005; Qiu et al., 2012). Also, strong positive IOD events lead to significant seasonal upwelling, which reduces SST. This is apparent in the Q-Q plots as a clear deviation in the slopes of the IFA distributions for site MD96-2060 (Fig. 7). Interestingly, the ETIO site SO189-39KL is also very sensitive to interannual variability within the subsurface (Fig. 7). However, the surface seasonality is so muted in that core compared to other sites that any changes in the altered IFA in the surface waters will be almost exclusively associated with changes in interannual variability (Thirumalai et al., 2019).

Figure 7 Quantile–Quantile (Q-Q) plots comparing original and altered modelled IFA δ¹⁸O distributions at four sites (U1467, NIOP 905, MD96-2060 and SO189-39KL) for both surface (top row) and thermocline (bottom row) depths. Each panel shows the relationship between the original modelled IFA δ¹⁸O values (x-axis) and the altered distributions (y-axis), generated by applying synthetic perturbations to temperature and salinity. Coloured points represent different types of variability scenarios: temperature increase (red), temperature decrease (blue), salinity increase (orange), and salinity decrease (green). The 1:1 line (dashed) indicates perfect agreement between original and altered distributions.

Overall, our forward-modelling results through Q-Q plots demonstrate that IFA-δ¹⁸O distributions coherently capture changes in seasonal climate variability at the surface and interannual climate variability at thermocline depth across the Indian Ocean (Fig. 7). Sensitivity tests additionally show that interannual surface perturbations produce site-dependent responses, with western sites (MD96-2060, NIOP 905) dominated by seasonality and eastern/central Indian Ocean sites (U1467, SO189-39KL) showing a stronger interannual surface signature (Fig. S7). Specifically, surface-dwelling species such as G. ruber are excellent candidates to record monsoon-driven seasonality, while thermocline-dwelling species will generally primarily reflect temperature-driven interannual anomalies associated with IOD/ENSO dynamics. One exception to this is in regions where large freshwater inputs or mixing processes occur, imprinting an important additional salinity influence (e.g., site U1467). This robust response at both depths validates the use of IFA to identify the prevalence of extreme climatic events in the sedimentary record, even where salinity exerts a secondary influence.

While these sensitivity experiments clarify how IFA signals theoretically respond to changing climate variability, the reliability of our interpretations depends on the model's ability to replicate observed sedimentary data. In the following section, we benchmark our forward-modelled results against the measured core-top distributions to assess the reproducibility of both isotopic variance and mean signals.

5.3 Benchmarking the forward model: Reproducibility of measured IFA-δ¹⁸O variance and mean

Previous studies have generated artificial IFA data in order to assess the ability of foraminifera to record Indian monsoon variability (Thirumalai and Clemens, 2020), ENSO dynamics in the eastern equatorial Pacific (Thirumalai et al., 2013), or for global IFA analysis with the Quantifa algorithm (Glaubke et al., 2021). Our study, which focuses on a streamlined script using modern reanalysis datasets (ORAS5), provides a useful tool to simulate random foraminiferal picking that can effectively reproduce the δ¹⁸O distributions observed in measured Late Holocene core-top samples. We now compare our new IFA dataset with Forward-modelled IFA using this new script as a tool to reproduce Late Holocene variability conditions for all sites.

Comparing forward-modelled (synthetic) and measured (real) core-top distributions (Fig. 8), using a two-sample Kolmogorov-Smirnov test (α=0.05) indicates significant differences at 5 of 7 site/species combinations (all except G. ruber at NIOP 905 and U1467; Table S1), often reflecting the mean offset discussed below rather than SD alone. Shapiro-Wilk tests confirm normality for all distributions except the synthetic G. ruber at MD96-2060 (p=0.007). Site NIOP 905 remains the clearest exception in dispersion specifically, with real SD higher than synthetic (0.565 ‰ vs. 0.377 ‰ surface; 0.529 ‰ vs. 0.431 ‰ thermocline). We also observe a systematic offset in absolute values across all sites: synthetic δ¹⁸O distributions are consistently shifted towards more negative (warmer/fresher) values compared to the real core-top measurements. This offset, which can reach ∼ 0.5 ‰ when comparing mean values, does not appear to affect the model's ability to simulate the range of climate variability and likely stems from factors such as foraminiferal vital effects, depth habitat uncertainties, or the older ages of the real Holocene core-top samples (which range from 0.1 to 3.8 kyr) relative to the synthetic modern ones.

5.3.1 Influence of anthropogenic warming and mean calcification depth on absolute offsets

These differences in mean δ¹⁸O values could be partly linked to the fact that modelled IFA δ¹⁸O values reflect conditions over the period 1958 to 2018 linked to anthropogenic climate change influenced conditions (Chaabane et al., 2024). The Indian Ocean is the region that experienced the fastest warming (0.12 °C per decade) over recent decades (Roxy et al., 2014). This warmer mean state can at least partly explain why synthetic core-top means are systematically lighter than the core-top measurements (Fig. 8).

Nevertheless, forward-modelling using ORAS5 reanalysis data produces similar SD values to core-tops where IFA-δ¹⁸O was measured, except for the very dynamic NIOP 905 site.

Beyond the physical structure of the water column and historical warming trends, the biological complexities of the proxy, specifically the seasonal timing of shell production and the ecological preferences of each species, should be considered to explain the remaining differences in signal recording. Deviation from the synthetic core-top means could also represent a vital effect on the δ¹⁸O, which is traceable when planktonic foraminifera calcite δ¹⁸O at equilibrium are compared (Fairbanks et al., 1982).

However, we note that the amplitude and frequency of extreme events such as ENSO, IOD and monsoon have already increased in the modern era due to anthropogenic global warming (Cai et al., 2013, 2014), and have been shown to intensify upwelling dynamics in regions such as the Arabian Sea and the Sumatra coast (Praveen et al., 2016). Our model-data comparison does not detect such an effect in the SD. It cannot, however, be excluded at a level below the precision of a 60-specimen estimate (± 0.02 ‰–0.05 ‰).

This offset could also partly reflect the choice of paleotemperature equation, as calibrations for a given species tend to differ more in their intercept than in their slope. Replacing Farmer et al. (2007) by Mulitza et al. (2003) for G. ruber shifts the simulated mean by −0.43 ‰ to −0.49 ‰, with only a minor effect on SD (+2 % to +7 %). For N. dutertrei, the calibration of Bouvier-Soumagnac and Duplessy (1985), known to yield unrealistically deep calcification depths (Wejnert et al., 2013), shifts the mean by −0.24 ‰ to −0.29 ‰ and reduces SD by 16 %–23 % (Fig. S9). In both cases, the synthetic values become more negative and the offset with the measured core-tops increases, so we retain the equations from Farmer et al. (2007).

Figure 8 Comparison between measured and modelled IFA δ¹⁸O_c distributions across multiple sites in the Indian Ocean for the surface (G. ruber-like) IFA (A) and the thermocline (N. dutertrei-like) IFA (B). For each site, panel shows the individual foraminiferal δ¹⁸O distributions as boxplots (real data: in green; forward-modelled data: in orange). The Boxplots span the first to third quartiles; whiskers extend to 1.5 × IQR. Red dots from real core-tops distribution are identified as outliers, defined with the 1.5 IQR method explained in Sect. 3.2. The black line in each box represents the median value of each dataset.

5.3.2 Impact of seasonality, shell flux, and ecology of planktonic foraminifera

One hypothesis about these offsets could be linked to the complexity of simulating foraminiferal ecology (Scott et al., 2025), with uncertainties related to ecological niches and living depths of planktonic foraminifera that could explain the difference between synthetic and real core-tops for Site NIOP 905. G. ruber is found in the Indian Ocean all year round and should represent mean annual conditions. On the other hand, simulation of foraminiferal growth shows that species such as G. ruber have a preferential maximum growth rate depending on the season and the region where they are found (Lombard et al., 2011). For the Somali margin, G. ruber shows vertical fluxes all year round, while peaking during the intermonsoon and NE monsoon periods (Curry et al., 1992) contrasting with the other sites where maximum growth rates occur during boreal summer (August and September). Our own reanalysis of FORCIS sediment trap data (Chaabane et al., 2023) at NIOP 905 also shows peak G. ruber flux in February (Fig. S10), during the same periods. To further assess this seasonal production bias, we additionally tested a mechanistic ecosystem-model estimate of monthly G. ruber production (PLAFOM2.0; Kretschmer et al., 2018), extracted at all four sites following Saenger et al. (2024). At site NIOP 905, this estimate converges with the February peak identified from FORCIS, although the two approaches rely on entirely different data. This agreement does not extend to the other three sites. Despite this disagreement, reweighting the Monte Carlo sampling by PLAFOM2.0-derived production changes the simulated G. ruber IFA SD by only 0.002 ‰–0.028 ‰, comparable to the FORCIS-based weighting, and does not alter the site ranking. PLAFOM2.0 does not simulate N. dutertrei, so this comparison is restricted to the surface species (Fig. S11).

Moreover, Stainbank et al. (2019) reported deeper apparent calcification depths in the Maldives (∼ 77 m for G. ruber, ∼ 97 m for N. dutertrei), which they attribute for G. ruber to the last chambers being formed deeper rather than to its living depth. Drawing the calcification depth of each simulated individual around these values, at all sites (Fig. S12), increases simulated IFA SD by 8 %–159 % in seven of eight cases, but reduces it by 36 % for G. ruber at MD96-2060 (p<0.001 in all cases). The G. ruber site ranking is largely reversed, with SO189-39KL showing the highest SD and MD96-2060 the lowest. However, this scheme also shifts the simulated G. ruber means by +0.43 ‰ to +0.87 ‰, well beyond the offset observed between synthetic and measured core-tops (+0.08 ‰ to +0.36 ‰; Fig. 8), and shifts N. dutertrei means in the direction that increases this offset. We therefore retain the layer-integrated approach and consider depth sampling, for instance using the vertical habitats simulated by PLAFOM2.0, a priority for future work.

In addition, the calcification depth of G. ruber, our target species, likely extends down to 50 m depth in the water column (Chaabane et al., 2026) possibly leading to a bias in our modelled IFA-δ¹⁸O. However, as this species remains in the mixed layer depth, there should be limited changes between both estimates. Site NIOP 905 shows minor changes in the mean SD values at different depth for the surface mixed layer (Fig. S13). For N. dutertrei, the FORCIS database shows an average occurrence depth in the Indian Ocean of ∼ 40 m (Chaabane et al., 2026), consistent with a broad vertical range spanning the base of the mixed layer to the sub-thermocline; this estimate is likely biased towards shallower depths given that non-stratified net tows integrate over large depth intervals. Yet, we urge caution while assigning calcification depth to N. dutertrei because the DCM depth is known to vary by tens of meters on local and seasonal scales, as recently shown in the northern Mozambique channel, a region with similar surface hydrological patterns as WTIO (Berthod et al., 2021).

Seasonal flux was characterised using FORCIS sediment-trap data (Chaabane et al., 2023) within a ± 10° box around each core site (Fig. S10). NIOP 905 benefits from a near trap (AS8, < 1 km); the other three sites rely on traps located 946 km (MD96-2060), 1312 km (U1467) and 1223 km (SO189-39KL) away. The MD96-2060 record comprises only 39 sediment-cup samples spanning 2.3 years. A sensitivity test weighting monthly IFA sampling by FORCIS flux weights changes simulated IFA SD by −0.037 ‰ to +0.039 ‰ across sites and species (Fig. S11). The largest change occurs at U1467 for N. dutertrei (−0.037 ‰), and NIOP 905 N. dutertrei is the only case where flux weighting increases rather than reduces IFA SD (+0.039 ‰), consistent with this site's near-trap (< 1 km) resolving a locally-distinct flux seasonality relative to the more distant traps used at the other sites, without changing the overall site ranking by variability.

N. dutertrei, our target thermocline species, is thought to grow mostly in September and October at site MD96-2060 and U1467, in February for the eastern Indian Ocean and in August for NIOP 905 (Lombard et al., 2011). The peak in August corresponds to the summer monsoon and the maximum intensity of the Somalian upwelling with strong advection of high-nutrient waters (Schott and McCreary, 2001), while the second-order mixing events during the winter monsoon trigger a second decrease in SST (Fig. 4). Interestingly, the thermocline δ¹⁸O recorded in warmest and coldest months is of the same magnitude as the ones recorded during extreme interannual events (Fig. 4), so it is difficult to attribute a shift in N. dutertrei IFA SD to a change in seasonal vs. interannual shift in the thermocline variability regime. N. dutertrei is commonly found above the DCM in tropical zones (Schmuker and Schiebel, 2002), so this species is most often reported to adopt an opportunistic behaviour (Kroon and Ganssen, 1989; Schiebel et al., 2004). Multiple studies show that the chlorophyll maximum off the Somalia coast is found within the first 40 m from the surface (Lakshmi et al., 2020) during the peak of the summer upwelling when the highest proportion of N. dutertrei may calcify and be directly influenced by temperature changes associated with cold water advection, analogously to what has been recorded in the G. bulloides fluxes in sediment traps (Curry et al., 1992). Compared to G. ruber, N. dutertrei could dwell over a large range of water depths (Chaabane et al., 2026; Fairbanks et al., 1982). However, the simulated thermocline SD is very close to the core-top IFA distribution (Fig. 8), except for NIOP 905. This highlights the fact that studying the variability of N. dutertrei could help to reconstruct the potential impact of the summer monsoon and the upwelling intensity at both seasonal and interannual scales. The extreme SD value associated with the real core-top of NIOP 905 (Fig. 8B) may represent extreme cold or warm thermocline temperatures (up to 1 ‰ warmer or colder than the mean), suggesting that the SD of N. dutertrei is meaningful for tracking extreme events. Similar seasonal dynamics are observed in the eastern Indian Ocean south of Java, where sediment traps also report that N. dutertrei and G. bulloides are more abundant in boreal summer (August–September) when intense upwelling occurs (Mohtadi et al., 2009), as described for the Somalian upwelling.

Overall, the IFA mean values of both species at all sites deviate from almost zero to values up to ± 1 ‰ when compared to the calculated δ¹⁸O_c profiles (Fig. 6). This likely indicates that the calcification depth or the seasonal variability of foraminiferal fluxes differs from those expected under mean-annual surface and thermocline temperature and salinity conditions, which could produce the observed deviations.

5.3.3 The role of the Salinity to δ¹⁸O_w relationship

Additional uncertainties could be related to the choice of the salinity-δ¹⁸O_w equation for each site to forward-model the synthetic core-tops. For the subsurface water, we used an equation derived from the available dataset from LeGrande and Schmidt (2006) restricted to the tropical Indian Ocean (Fig. S3) that provides an average relationship (Eq. 4). However, for the surface, multiple equations are used in studies focused on regional areas and some of which differ from the global Indian Ocean equation proposed by LeGrande and Schmidt (2006) (δ¹⁸O_w=0.16S–5.31). For each site we used regional and derived equations adapted for each different area following Stainbank et al. (2021) for site U1467 (δ¹⁸O_w=0.28S–9.24), the LeGrande and Schmidt (2006) equation for SO189-39KL, as in Thirumalai et al. (2019) for the eastern Indian Ocean (δ¹⁸O_w=0.16S–5.31), and (Eq. 3) for the two other western Indian Ocean sites, following Thirumalai and Clemens (2020).

For NIOP 905 the few reported data of δ¹⁸O_w and surface salinity presented in Ganssen et al. (2011) represent a relation of δ¹⁸O_w=0.12S–4.03, defining a very shallow slope. One possibility is that such seawater δ¹⁸O and salinity measurements are representative of a limited area over a very short time window and do not reflect the mean oceanographic isotopic and salinity conditions in a regional, long-term context. Moreover, the nearest equation to this site established for the Arabian Sea by Delaygue et al., 2001 (i.e. δ¹⁸O_w=0.26S–8.9) indicates that the slope may be at least twice as large as that calculated using the in-situ data published in Ganssen et al. (2011). This discrepancy is particularly relevant in a region strongly affected by upwelling and mixing that may add a deeper component into the equation that may tend to make the slope even steeper than in the surface. We note that the mean IFA-δ¹⁸O_c generated is very sensitive to the δ¹⁸O_w-S relationship in the forward model, and could be shifted to higher or lower δ¹⁸O values with differences as large as 0.5 ‰ between two equations.

5.3.4 Impact of temporal integration, sedimentation rates, and core-top age

Another source of uncertainty resides in using a framework of cores that have different sedimentation rates reflecting their distinct environments. WAS and ETIO sites are located on continental margins, while CTIO sits on a carbonate platform (Maldives), and WTIO site is under pelagic sedimentation conditions. Our synthetic IFA distributions are generated from a 60-year instrumental window (1958–2018), but our core-top samples represent a time-averaged signal determined by the local sedimentation rate and the amount of time averaged within the 1 cm wide samples after it has been homogenized by bioturbation. Based on regional sedimentation rates, the temporal window integrated into each 1 cm sample varies significantly across the basin. For sites SO189-39KL, NIOP 905, and U1467, sediment rates are very high and range from 20 to 40 cm kyr⁻¹, which represent approximately 25 to 50 years for this thickness (Mohtadi et al., 2014; Ivanochko et al., 2005; Stainbank et al., 2020). Site MD96-2060 has a much lower sedimentation rate (6–8 cm kyr⁻¹), so that a 1 cm sample may integrate 125 to 166 years without considering potential bioturbation effects. A scenario in which the bioturbated layer would reach 10 cm thickness would suggest that the MD96-2060 core-top integrates oceanic conditions over the last millennium. This may lead to an overestimation of the total variability simulated in our core-top because foraminifera potentially lived under conditions when the climate background was different. However, a higher SD in the real vs synthetic core-top is not seen in our data. A more direct test of this hypothesis would require an ensemble of mid-Holocene or last-millennium climate simulations, which is beyond the scope of the present study. The absence of an elevated real SD at this site nonetheless argues against a strong effect of this kind. Our core-top age at site MD96-2060 is ∼ 3.8 kyr, which is much older than the core-tops at other sites (Table 1), indicating that the IFA here recorded environmental conditions during the late Holocene, not the last few centuries. Observing the results from real and synthetic core-tops (Fig. 8) shows that even with a different sedimentation rate and core-top age, site MD96-2060 exhibits the same variability with the same SD as most of the other sites except for NIOP 905. The range of seasonal and interannual climate hence likely remained relatively stable between the modern and Late Holocene periods. The dating of the MD96-2060 core-top (0–1 cm) is likely affected by bioturbation and/or incomplete core-top recovery with the CALYPSO piston corer, making its mean planktonic foraminifera age older than at the other sites. We are aware of a nearby piston core with analogous hydrological conditions (MD96-2062), but onboard preliminary analysis (core description, magnetic susceptibility, sedimentary reflectance) indicates a sedimentary gap at ∼ 4 m depth, precluding its use as an alternative. As the Holocene δ¹⁸O values of G. ruber remain largely stable at more dynamically active sites situated North (Ganssen et al., 2011) and South (Wang et al., 2013) of MD96-2060, we have confidence in the reliability of core-top IFA values at this site for testing the applicability of IFA under contrasting past climate backgrounds.

Overall, the comparison of SD between forward-modelled and core-top IFA reveals good agreement, though a systematic offset in mean values is observed across all sites. This offset likely reflects a potential combination of anthropogenic warming trends embedded in the ORAS5 reanalysis and uncertainties in the regional δ¹⁸O_w-S relationship, both of which would shift synthetic distributions toward lighter values. At site NIOP 905, the higher SD in real core-tops point to additional factors, related to the strong upwelling dynamics of the Somali margin. However, the most probable primary source of discrepancy in both mean values and SD across all sites lies in the seasonal variability of foraminiferal fluxes and ecological factors, as discussed in Sect. 5.3.2. Temporal integration of the IFA signal and sedimentation rate may introduce additional uncertainties but are expected to have a comparatively minor influence on the observed offsets.

6 Summary & Conclusion

We evaluated the sensitivity and spatial variability of individual foraminifera δ¹⁸O_c signals across the Indian Ocean by combining forward-modelled synthetic core-tops with real core-top measurements. Our approach shows that forward-modelling synthetic IFA with modern reanalysis data can reproduce the SD of δ¹⁸O_c observed in Late Holocene core-tops, providing a robust framework to simulate natural variability. The simulated SD are in good agreement with those derived from real IFA data for both surface and thermocline-dwelling species.

Forward-modelled IFA enables statistical analyses and spatial comparisons of modern variability across the Indian Ocean, as well as testing the effects of altered climate variability scenarios. When statistically modified, the seasonality or interannual δ¹⁸O_c variability show coherent response to higher amplitude showing that this method provides a tool to record this type of high frequency variability. Amplification or reduction of seasonal and interannual variability highlights the ability of IFA δ¹⁸O_c to record extreme events.

Our results indicate that at the surface, δ¹⁸O_c variability from western and northern Indian Ocean (MD96-2060, NIOP 905, and U1467) shows strong sensitivity to seasonal variability associated with monsoon dynamics and upwelling processes, whereas the equatorial eastern site (SO189-39KL) is less affected by seasonality and more influenced by interannual events associated with phenomena such as IOD and ENSO influencing the temperature variability. At thermocline depth, δ¹⁸O_c variability is largely driven by temperature changes, with salinity exerting a secondary influence except in regions affected by strong freshwater inputs or mixing such as site U1467. Overall, the good agreement between forward-modelled and measured IFA SD demonstrates that the modern climate variability captured by IFA depends primarily on the depth habitat of the target species and the regional oceanographic setting. In the western Indian Ocean, surface-dwelling G. ruber predominantly records monsoon-driven seasonal variability, while thermocline-dwelling N. dutertrei appears more sensitive to interannual forcing, likely reflecting IOD-driven temperature anomalies at thermocline depth. Site NIOP 905 stands as a notable exception, where the shallow thermocline and intense seasonal upwelling dynamics impart a strong seasonal imprint on the N. dutertrei and lead to a mixed signal of the total variability.

Several sources of uncertainty remain, especially at site NIOP 905 that highlights a significant difference in SD between real and synthetic core-top. These differences may arise from ecological factors, including the actual depth range and seasonal growth patterns of foraminiferal species and their opportunistic behaviour, as well as from the choice of salinity–δ¹⁸O_w equations used for forward-modelling. Also, biases in reanalysis products, especially near coastal regions, may influence the mean δ¹⁸O_c signal, though their effect on variability metrics appears limited and unclear.

The method introduces a template for interpreting IFA-δ¹⁸O_c and its associated spatial fingerprints under different climate backgrounds, to be applied to downcore measurements. Our results allow for a clear distinction between regions dominated by seasonal forcings such as in monsoon systems and upwelling regions, as well as those governed by interannual climate modes. This framework establishes a first robust starting point for future studies aiming to reconstruct past seasonal and interannual climate variability using individual foraminifera analysis δ¹⁸O for distinct periods such as the Last Glacial Maximum or other warm intervals and to test a wider range of perturbations, for instance derived from large ensembles of historical and future climate simulations.