Anthropogenic CO₂, air–sea CO₂ fluxes, and acidification in the Southern Ocean: results from a time-series analysis at station OISO-KERFIX (51° S–68° E)

The temporal variation of the carbonate system, air–sea CO₂ fluxes, and pH is analyzed in the southern Indian Ocean, south of the polar front, based on in situ data obtained from 1985 to 2021 at a fixed station (50°40^′ S–68°25^′ E) and results from a neural network model that reconstructs the fugacity of CO₂ (fCO₂) and fluxes at monthly scale. Anthropogenic CO₂ (C_ant) is estimated in the water column and is detected down to the bottom (1600 m) in 1985, resulting in an aragonite saturation horizon at 600 m that migrated up to 400 m in 2021 due to the accumulation of C_ant. At the subsurface, the trend of C_ant is estimated at $+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹ with a detectable increase in the trend in recent years. At the surface during austral winter the oceanic fCO₂ increased at a rate close to or slightly lower than in the atmosphere. To the contrary, in summer, we observed contrasting fCO₂ and dissolved inorganic carbon (C_T) trends depending on the decade and emphasizing the role of biological drivers on air–sea CO₂ fluxes and pH inter-annual variability. The regional air–sea CO₂ fluxes evolved from an annual source to the atmosphere of 0.8 molC m⁻² yr⁻¹ in 1985 to a sink of −0.5 molC m⁻² yr⁻¹ in 2020. Over 1985–2020, the annual pH trend in surface waters of $-\mathrm{0.0165}\pm \mathrm{0.0040}$ per decade was mainly controlled by the accumulation of anthropogenic CO₂, but the summer pH trends were modulated by natural processes that reduced the acidification rate in the last decade. Using historical data from November 1962, we estimated the long-term trend for fCO₂, C_T, and pH, confirming that the progressive acidification was driven by the atmospheric CO₂ increase. In 59 years this led to a diminution of 11 % for both aragonite and calcite saturation state. As atmospheric CO₂ is expected to increase in the future, the pH and carbonate saturation state will decrease at a faster rate than observed in recent years. A projection of future C_T concentrations for a high emission scenario (SSP5-8.5) indicates that the surface pH in 2100 would decrease to 7.32 in winter. This is up to −0.86 lower than pre-industrial pH and −0.71 lower than pH observed in 2020. The aragonite undersaturation in surface waters would be reached as soon as 2050 (scenario SSP5-8.5) and 20 years later for a stabilization scenario (SSP2-4.5) with potential impacts on phytoplankton species and higher trophic levels in the rich ecosystems of the Kerguelen Islands area.

1 Introduction

The ocean plays an important role in mitigating climate change by taking up a large part of the excess of heat (Cheng et al., 2020a; Fox-Kemper et al., 2021) and of CO₂ released by human activities (Sabine et al., 2004; Gruber et al., 2019a; Canadell et al., 2021). Since 1750, the global ocean has captured 185±35 PgC (petagram of carbon) from a total of 700±75 PgC of anthropogenic carbon emissions from fossils fuels and land-use changes (Friedlingstein et al., 2022). The oceanic sink for anthropogenic CO₂ increased progressively from 1.1±0.4 PgC yr⁻¹ in the 1960s to 2.3±0.4 PgC yr⁻¹ in the 2000s. Over the decade 2012–2021, the partitioning of the anthropogenic CO₂ uptake was roughly equal between the ocean (2.9±0.4 PgC yr⁻¹) and the land (3.1±0.6 PgC yr⁻¹) (Friedlingstein et al., 2022). This partitioning has been confirmed for the decade 2013–2022 (Friedlingstein et al., 2023).

Ocean observations indicate that the Southern Ocean (SO) south of 45° S has been accumulating each year about 0.5 PgC yr⁻¹ since the 1990s (e.g., Takahashi et al., 2009a, b; Lenton et al., 2013; Rödenbeck et al., 2013; Long et al., 2021; Fay et al., 2024; Gray, 2024). Results based on BGC-Argo floats (Southern Ocean Carbon and Climate Observations and Modeling project, SOCCOM) suggest that the CO₂ sink in the SO might be much lower (0.16 PgC yr⁻¹ south of 44° S for the period 2015–2017; Gray et al., 2018; Bushinsky et al., 2019), but there is an ongoing debate on the size of the carbon sink in this region depending the periods and methods (Long et al., 2021; Sutton et al., 2021; Hauck et al., 2023b; Gray, 2024). It is also well established that the CO₂ sink in the SO undergoes substantial decadal variability first documented for the 1990s (Le Quéré et al., 2007; Metzl, 2009; Lenton et al., 2013) and subsequently identified for the period 1982–2018 (Landschützer et al., 2015; Keppler and Landschützer, 2019; Mackay et al., 2022; Hauck et al., 2023a, b). However as for the mean state, there are also uncertainties on both the magnitude and phasing of decadal variability in the SO carbon sink mainly due to insufficient sampling (Gloege et al., 2021; Hauck et al., 2023a, b). A recent extension of the period to 1957–2020 suggests that the inter-annual to decadal variability of the SO CO₂ sink was most pronounced after the 1980s (Rödenbeck et al., 2022; Bennington et al., 2022). Whatever the variability of the SO CO₂ sink since the 1960s, the ocean continuously absorbs atmospheric CO₂ ,and the distribution of anthropogenic CO₂ (C_ant) in the SO is now relatively well documented (e.g., Pardo et al., 2014; Gruber et al., 2019a) thanks to the Global Ocean Data Analysis Project (GLODAP) data synthesis effort for the global ocean (Olsen et al., 2016, 2019, 2020). The SO takes up about 40 % of the total anthropogenic carbon that enters the ocean (Khatiwala et al., 2013; Gruber et al., 2019a).

The anthropogenic CO₂ uptake in the ocean results in lowering carbonate ion concentrations and pH, a chemical process termed “ocean acidification” (OA) (Caldeira and Wickett 2003; Doney et al., 2009). This decreases the saturation state with respect to carbonate minerals (aragonite, ΩAr, and calcite, ΩCa), a process most pronounced in the cold waters at high latitudes where the saturation state is naturally low (Orr et al., 2005; Takahashi et al., 2014; Jiang et al., 2015). The first estimate of C_ant distribution in the global ocean (for a nominal year 1994; Sabine et al., 2004) shows that the accumulation of C_ant led to an upward migration of the ΩAr and ΩCa saturation horizon in all ocean basins (Feely et al., 2004). This change is particularly pronounced south of the polar front (PF) in the SO due to both C_ant uptake and the enhanced upwelling of carbon-rich deep waters (e.g., Hauck et al., 2010; Pardo et al., 2017). It has been suggested, through numerical studies, that depending on future CO₂ emission levels, surface waters in the SO could reach undersaturation state for aragonite by 2030–2050 (Orr et al., 2005; Gangstø et al., 2008; McNeil and Matear, 2008; Negrete-Garcia et al., 2019). Such a change would have multiple and detrimental impacts on marine ecosystems (Fabry et al., 2008; Doney et al., 2012; Bopp et al., 2013), in particular calcifying marine organisms, especially aragonite producers such as pteropods (Hunt et al., 2008; Gardner et al., 2023), as well as calcite-producing planktonic foraminifera (Moy et al., 2009), coccolithophorids (Beaufort et al., 2011), and non-calcifying species such as the abundant SO diatoms (e.g., Benoiston et al., 2017; Petrou et al., 2019; Weir et al., 2020; Duncan et al., 2022) and krill (Kawaguchi et al., 2013).

Hindcast simulations with global ocean biogeochemical models (GOBMs), as well as projections with Earth system models (ESMs), have been used to evaluate the ocean carbon cycle over the past decades and future changes in C_ant storage, ocean acidification, or impacts of global changes on marine ecosystems. However, current model-based estimates of the contemporary SO CO₂ sink are subject to relatively large uncertainties (e.g., Long et al., 2013; Gooya et al., 2023; Hauck et al., 2020, 2023a, b; Mayot et al., 2023; DeVries et al., 2023). Differences between GOBMs can reach up to 0.7 PgC yr⁻¹ in the SO (Hauck et al., 2020), which is roughly equivalent to the mean climatological flux of 0.5 PgC yr⁻¹ (McNeil et al., 2007; Takahashi et al., 2009a, b; Lenton et al., 2013). At the high latitudes of the SO (>50° S) for the 2010s, ESMs from the Coupled Model Intercomparison Project Phase 6 (CMIP6) simulated either a large sink or a modest source of CO₂ (McKinley et al., 2023). This is mainly due to incorrect or missing physical and/or biological processes in the models (e.g., Pilcher et al., 2015; Kessler and Tjiputra, 2016; Mongwe et al., 2018; Lerner et al., 2021), leading to biases in the seasonality of temperature, dissolved inorganic carbon C_T, partial pressure of CO₂ (pCO₂), air–sea CO₂ fluxes, pH, or Ω (e.g., McNeil and Sasse, 2016; Rodgers et al., 2023; Rustogi et al., 2023; Joos et al., 2023). Such model imperfections should be resolved to gain reliability in future projections of CO₂ uptake, OA, productivity, and the responses of the marine ecosystems (Frölicher et al., 2015; Hauck et al., 2015; Sasse et al., 2015; Kessler and Tjiputra, 2016; McNeil and Sasse 2016; Kwiatkowski and Orr, 2018; Negrete-Garcia et al., 2019; Burger et al., 2020; Krumhardt et al., 2022; Jiang et al., 2023; Mongwe et al., 2024). In this context, long-term biogeochemical observations are particularly valuable to quantify and understand recent past and current changes, and ultimately evaluate model simulations, as often concluded in modeling studies (e.g., Kessler and Tjiputra, 2016; Gooya et al., 2023; Wright et al., 2023; Hauck et al., 2023a; Mayot et al., 2023; Rodgers et al., 2023).

Table 1 Trends of oceanic fCO₂ (µatm yr⁻¹) and pH (per decade) in the Southern Ocean south of the polar front based on observations. IO: Indian Ocean sector. PO: Pacific Ocean sector. AO: Atlantic Ocean sector. SO SPSS: Southern Ocean subpolar seasonally stratified biome (around 50–60° S). PZ: polar zone. NR: not reported. Standard deviations when available are given in brackets.

Although the SO south of the polar front remains much less observed than other oceanic regions, several observations-based studies have estimated the decrease in pH in surface waters in response to the increase in oceanic CO₂ fugacity, fCO₂ (Midorikawa et al., 2012; Takahashi et al., 2014; Lauvset et al., 2015; Munro et al., 2015; Xue et al., 2018; Iida et al., 2021; Leseurre et al., 2022; Brandon et al., 2022). Results showed a large range in the pH trends from −0.008 −0.035 per decade depending on the period and the region of interest. Most of these analyses were based on summer observations (Table 1), and some studies highlighted contrasting pH trends on a 5–10-year timescale probably linked to large-scale climate variability such as the Southern Annular Mode (SAM) (e.g., Xue et al., 2018). Given such variability, it is important to continue monitoring fCO₂ and pH trends and, if possible, at different seasons as future change in CO₂ uptake and potential tipping points of the carbonate saturation state also depend on seasonality (Sasse et al., 2015). The above observational studies were dedicated to pH changes in surface waters. In contrast to northern high latitudes (e.g., Olafsson et al., 2009, 2010; Franco et al., 2021; Skjelvan et al., 2022), few studies in the SO evaluated decadal changes of carbonate system properties and acidification in the water column based on time-series stations.These changes in the SO water column were investigated from data collected during cruises generally 3 to 15 years apart (e.g., Hauck et al., 2010; Van Heuven et al., 2011; Pardo et al., 2017; Tanhua et al., 2017; Carter et al., 2019).

Figure 1 Left: annual air–sea CO₂ flux (molC m⁻² yr⁻¹) in the southern Indian Ocean for year 2020 from the FFNN model (negative flux for ocean sink, positive flux for ocean source). The black box identifies the location of the study southwest of Kerguelen Islands. Right: track of cruises with underway fCO₂ data southwest of Kerguelen Islands. The station at 50°40^′ S–68°25^′ E occupied in 1985, 1992–1993, and 1998–2021 is indicated by a yellow circle. The yellow square is the region selected to calculate the mean values from the underway surface observations and from the FFNN model. Figures produced with Ocean Data View (ODV; Schlitzer, 2018).

The present study complements in time, seasons, and in the water column the surface fCO₂ and pH trends investigated by Leseurre et al. (2022) in different regions of the southern Indian Ocean for the period 1998–2019 during austral summer. South of the PF around 50° S, Leseurre et al. (2022) showed that in summer the surface fCO₂ increase and pH decrease over 20 years were mainly driven by the accumulation of anthropogenic CO₂ by about $+\mathrm{0.6}\pm \mathrm{0.2}$ µmol kg⁻¹ yr⁻¹ and by a small warming of $+\mathrm{0.03}\pm \mathrm{0.02}$ °C yr⁻¹. In addition Leseurre et al. (2022) showed that in the recent decade, 2007–2019, the fCO₂ trend was low ($+\mathrm{0.3}\pm \mathrm{0.2}$ µatm yr⁻¹) compared to the previous decade ($+\mathrm{5.3}\pm \mathrm{0.4}$ µatm yr⁻¹ over 1998–2007), highlighting the sensitivity of the fCO₂ and pH trends to the selected time period (especially during summer). In particular, they observed relatively stable pH values over 2010–2019 (i.e., no decrease in pH) with no clear explanation on the origin of the slow-down of the fCO₂ and pH trends in surface waters south of the PF in recent years. To complement the analysis by Leseurre et al. (2022) based on summer observations over the period 1998–2019, this study focuses on one location regularly visited south of the polar front (around 50° S–68° E southwest of Kerguelen Island; Fig. 1). The analysis period is first extended back to 1985 and forward to 2021 to investigate the recent status of fCO₂ and pH. We also evaluate the trends during late winter using sparse data in October/November. The combination of in situ observations and monthly estimates from a neural network model over the period 1985–2020 (Chau et al., 2022) enables us to assess potential changes in seasonality of the surface ocean carbonate system (including fCO₂, C_T, pH, Ω) as suggested in recent decades or in future scenarios (Hauck and Völker, 2015; Gallego et al., 2018; Landschützer et al., 2018; Kwiatkowski and Orr, 2018; Kwiatkowski et al., 2020; Lerner et al., 2021; Fassbender et al., 2022; Yun et al., 2022; Rodgers et al., 2023; Joos et al., 2023). The changes observed in surface waters will be related to changes in C_ant concentrations estimated in the water column and will be complemented by an analysis of OA at depth between 1985 and 2021. Finally we will explore the long-term change of surface fCO₂ and pH since the 1960s and potential future changes of the carbonate system at this time-series site.

2 Data selection, methods, and quality control

2.1 Study area and data selection

This study focuses on a high-nutrient low-chlorophyll area (HNLC; Minas and Minas, 1992) in the Indian sector of the Southern Ocean (SO) in the Permanent Open-Ocean Zone (POOZ) south of the polar front (PF) and southwest of Kerguelen Islands (around 50° S–68° E, Fig. 1). The Kerguelen Plateau is an extended topographic feature that controls part of the Antarctic Circumpolar Current (ACC) and generates eddies (Daniault and Ménard, 1985) and the northward deflection of the PF just east of the island (Pauthenet et al., 2018). The plateau is also a region of relatively high chlorophyll-a (Chl-a) concentration (Moore and Abbott, 2000; Mongin et al., 2008) and strong CO₂ uptake during austral spring–summer that contrasts with the weaker sink over the POOZ/HNLC (Metzl et al., 2006; Jouandet et al., 2008, 2011; Lo Monaco et al., 2014; Leseurre et al., 2022). The POOZ/HNLC region west (upstream) of the Kerguelen Plateau is characterized by rather stable water mass properties (temperature, salinity, oxygen, or nutrients) over time and low eddy activity compared to the plateau (Daniault and Ménard, 1985; Chapman et al., 2015; Dove et al., 2022). In this region, located in the deep Enderby Basin, the flow is not constrained by topography and there is no local upwelling that would import C_T-rich waters to the surface layers as observed on the eastern side of the Kerguelen Plateau (Brady et al., 2021).

The Indian sector of the SO is also recognized to host the strongest winds in the SO leading to year-round high gas transfer coefficients (Wanninkhof and Trinanes, 2017). As a result, and in contrast to the Atlantic sector of the SO, the Indian region south of 45° S was a periodic annual CO₂ source, especially in the 1960s to the 1980s (Rödenbeck et al., 2022; Bennington et al., 2022; Prend et al., 2022; Gray, 2024). In the POOZ–HNLC region, high winter wind speed (monthly average up to 16 m s⁻¹) and associated heat loss drive deep mixing. Deep winter mixing entrains subsurface properties to the surface layer and increases surface C_T concentrations, leading to wintertime outgassing of CO₂ (Metzl et al., 2006). This combination of characteristics makes the region an ideal test bed for 1-D modeling studies investigating the temporal dynamics and drivers of biogeochemical processes including nutrients, iron, phytoplankton, and carbon (Pondaven et al., 1998, 2000; Louanchi et al., 1999, 2001; Jabaud-Jan et al., 2004; Metzl et al., 2006; Mongin et al., 2006, 2007; Kane et al., 2011; Pasquer et al., 2015; Demuynck et al., 2020).

Here we used surface and water column observations around location 50° 40^′ S–68° 25^′ E (Fig. 1, Table S1 in the Supplement), historically called KERFIX station (KERguelen FIXed station) sampled from 1990 to 1995 in the framework of the WOCE and JGOFS programs (Jeandel et al., 1998). The station was first occupied in March 1985 during the INDIGO-1 cruise (Indian Ocean Geochemistry; Poisson, 1985; Poisson et al., 1988), and since 1998 it has been regularly visited during the OISO cruises (Océan Indien Service d'Observations; Metzl and Lo Monaco, 1998, https://doi.org/10.18142/228). The regular occupation from 1985 to 2021 makes it the longest time-series station in the Southern Ocean POOZ/HNLC area for investigating the inter-annual to decadal trends of carbonate properties in surface waters and across the water column (0–1600 m). Despite the occasional large anomalies in surface waters properties (e.g., lower temperature in December 1998, lower salinity in February 2013), we consider all observations selected for this study both in surface waters and the water column to be representative of the water masses in this POOZ/HNLC region upstream of the Kerguelen Plateau.

Data for the period 1985–2011 were extracted from the GLODAP data product, version V2.2021 (Lauvset et al., 2021a, b; Table S1a). Observations collected during OISO cruises from 2012 to 2021 will be included in GLODAP-V3. For the surface water properties, all available underway fCO₂ data were selected (Fig. 1). This includes one cruise in November 1962 (Keeling and Waterman, 1968) and 41 cruises from 1991 to 2021 (Table S1b). All surface temperature, salinity, and fCO₂ data were extracted from the SOCAT (Surface Ocean CO₂ Atlas) data-product version v2022 (Bakker et al., 2016, 2022) and have an accuracy of fCO₂ between 2 and 5 µatm.

2.2 Methods

The methods for surface underway fCO₂ and biogeochemical properties (oxygen, C_T, total alkalinity A_T, nutrients) in the water column for the INDIGO-1, KERFIX, and OISO cruises were described in previous studies (e.g., Poisson et al., 1993; Louanchi et al., 2001; Metzl et al., 2006; Metzl, 2009; Mahieu et al., 2020; Leseurre et al., 2022). Here we briefly recall the methods for underway fCO₂ and water column observations.

2.2.1 Surface fCO₂ data

For fCO₂ measurements in 1991–2021, sea-surface water was continuously equilibrated with a “thin film” type equilibrator thermostated with surface seawater (Poisson et al., 1993). The xCO₂ in the dried gas was measured with a non-dispersive infrared analyzer (NDIR, Siemens ULTRAMAT 5F or 6F). Standard gases for calibration (around 270, 350, and 490 ppm) were measured every 6 h. To correct xCO₂ dry measurements to fCO₂ in situ data, we used polynomials from Weiss and Price (1980) for vapor pressure and from Copin-Montégut (1988, 1989) for temperature. Note that when incorporated in the SOCAT database, the original fCO₂ data are recomputed (Pfeil et al., 2013) using temperature correction from Takahashi et al. (1993). Given the small difference between equilibrium temperature and sea surface temperature ($+\mathrm{0.56}\pm \mathrm{0.30}$ °C on average for the cruises in 1998–2021), the fCO₂ data from SOCAT used in this analysis (Bakker et al., 2022) are almost identical (within 1 µatm) to the original fCO₂ values from our cruises (http://www.ncei.noaa.gov/access/ocean-carbon-data-system/oceans/VOS_Program/OISO.html, last access: 15 January 2024).

2.2.2 Water column data

Over the period 1990–1995, water samples were collected during the KERFIX program on the ship La Curieuse at standard depths using 8 L Niskin bottles mounted on a stainless-steel cable and equipped with reversing SIS pressure and temperature probes. Methods and accuracy for the geochemical measurements used in this analysis (A_T, C_T, oxygen, nutrients) are detailed by Jeandel et al. (1998) and by Louanchi et al. (2001). From 1998 onwards, the station was occupied within the framework of the OISO long-term monitoring program on board the R/V Marion Dufresne. We used conductivity–temperature–depth (CTD) sensors mounted on 24 rosette samplers equipped with 12 L Niskin bottles. Temperature and salinity measurements have an accuracy of 0.002 °C and 0.005 respectively (Mahieu et al., 2020). Samples for A_T and C_T were filled in 500 mL glass bottles and poisoned with 300 µL of saturated mercuric chloride solution to halt biological activity. Discrete C_T and A_T samples were analyzed onboard by potentiometric titration derived from the method developed by Edmond (1970) using a closed cell. Based on replicate samples from the surface or depth, the repeatability for A_T and C_T varies from 1 to 3.5 µmol kg⁻¹ depending on the cruise. The accuracy of ±3 µmol kg⁻¹ was ensured by daily analyses of certified reference materials (CRMs) provided by Andrew Dickson's laboratory (Scripps Institute of Oceanography).

Dissolved oxygen (O₂) concentration was determined by a sensor fixed on the rosette, and values were adjusted based on discrete measurements (Winkler method, Carpenter, 1965) using a potentiometric titration system. Accuracy for O₂ is ±2 µmol kg⁻¹ (Mahieu et al., 2020). Although long-term deoxygenation in the Southern Ocean has been suggested (Ito et al., 2017; Schmidtko et al., 2017; Oschlies et al., 2018), no significant trend in O₂ was identified over 1985–2021 at this station around 50° S in both the surface and the subsurface layer (at the depth of the temperature minimum representing winter water, a layer used for C_ant calculations as described later). However, in the station data a small O₂ decrease was detected around 800 m in the O₂ minimum layer over 36 years ($-\mathrm{0.22}\pm \mathrm{0.07}$ µmol kg⁻¹ yr⁻¹). As this has no impact on the interpretation for pH and Ω trends for this analysis, the observed change of O₂ at depth will not be discussed further. Here the O₂ data are mainly used for the calculation of anthropogenic CO₂ concentrations, and the observed O₂ change at depth is too small to have an impact on temporal variations of C_ant concentrations given the uncertainty of the calculation.

Nitrate (NO₃) and silicate (DSi) were analyzed on board or at LOCEAN/Paris by colorimetry following the methods described by Tréguer and Le Corre (1975) for 1998–2008 or from Coverly et al. (2009) for 2009–2021. The uncertainty of NO₃ and DSi measurements is ±0.1 µmol kg⁻¹. Based on replicate measurements on deep samples, we estimate an error of about 0.3 % for both nutrients. Phosphate (PO₄) samples were analyzed from a few cruises following the method of Murphy and Riley (1962) revised by Strickland and Parsons (1972) with an uncertainty of ±0.02 µmol kg⁻¹. When nutrient data are not available for a cruise, we used climatological values based on the seasonal nutrient cycles inferred from data from 1990 to 2021. This method has a very small impact on the carbonate system calculations and the trend analysis as we did not detect any significant trends in nutrients in surface or at depth since 1985 (not shown) as opposed to what has been observed at higher latitudes of the SO (Iida et al., 2013; Hoppema et al., 2015). However, we will see in Sect. 3.1 that the inter-annual variability of nutrients (especially DSi in the HNLC region) might inform on potential changes in biological processes.

Samples were collected in the top layers (0–150 m) for chlorophyll a (Chl a). For that, 1–2 L of seawater was filtered onto 0.7 µm glass microfiber filters (GF/F, Whatman), and filters were stored at −80 °C onboard. Back at the LOCEAN/Paris laboratory, samples were extracted in 90 % acetone (Strickland and Parsons, 1972), and the fluorescence of Chl a was measured on a Turner Type 450 fluorometer for the period 1998–2007 and since 2009 at 670 nm on a Hitachi F-4500 spectrofluorometer (Neveux and Lantoine, 1993).

2.2.3 Data quality control and data consistency

When exploring the trends of ocean properties based on different cruises more than 35 years apart, it is important to first verify the consistency of the data and to correct for any bias or drift. The INDIGO data from 1985 (i.e., prior to CRMs available for A_T and C_T) were first controlled prior to their incorporation into the original GLODAP product (Sabine et al., 1999; Key et al., 2004), and corrections for A_T and C_T were revisited within the framework of the CARINA project (CARbon IN the Atlantic; Lo Monaco et al., 2010) and the GLODAPv2 synthesis (Olsen et al., 2016). A secondary quality control was performed on the data from the OISO cruises collected between 1998 and 2011 within the CARINA and GLODAP-v2 initiatives (Lo Monaco et al., 2010; Olsen et al., 2016). Significant offsets were identified for A_T and C_T in samples from the KERFIX cruises (1990–1993) compared to INDIGO and OISO data, and it was proposed to correct the original values by −35 µmol kg⁻¹ for C_T and −49 µmol kg⁻¹ for A_T (Metzl et al., 2006). These corrections were applied in GLODAP version v2.2019 (Olsen et al., 2019) and resulted in coherent A_T and C_T concentrations for KERFIX in the deep layers compared to other cruises (Supplement Table S2, Fig. S1). The same data quality control protocol as for GLODAP-v2 was applied to data from OISO cruises for the period 2012–2021 (Mahieu et al., 2020). Given the accuracy of the data no systematic bias (except in 2014) was found for the properties measured in 2012–2021. The time series of A_T and C_T at depths below 1450 m for all cruises in 1985–2021 show some variability but no trend over 36 years as expected in the bottom waters in this region (Fig. S1). However, we identified a small bias for C_T in 2014 (cruise OISO-23) where C_T concentrations in the deep water appeared slightly lower (2228–2234 µmol kg⁻¹ in 2014 compared to the mean value of 2240.7±3.7 µmol kg⁻¹; Table S2, Fig. S1). When compared to fCO₂ in surface waters, we also suspect the C_T data in the mixed layer in 2014 to be too low by about 10 µmol kg⁻¹ (Figs. S2, S3). Therefore we applied a WOCE/GLODAP flag 3 for C_T data of this cruise and will not use the station data in 2014 for the C_ant calculations and the trend analysis described in this study.

2.2.4 CMEMS-LSCE-FFNN model

As most of the cruises took place during austral summer and data are not available each year, we completed the observations with the results from an ensemble of feed-forward neural network models (CMEMS-LSCE-FFNN or FFNN for simplicity here; Chau et al., 2022). The FFNN model allows mapping at global-scale monthly surface fCO₂ from the SOCAT gridded datasets and ancillary variables. The reconstructed fCO₂ is then used to derive monthly surface C_T and pH fields as well as air–sea CO₂ fluxes. This data product is used to investigate the trends for different seasons and to derive estimates of annual air–sea CO₂ fluxes to interpret the change in CO₂ uptake, if any. For a full description of the model, access to the data, and a statistical evaluation of fCO₂ reconstructions, please refer to Chau et al. (2022). Within this study, we compared the FFNN fCO₂ with observations from 35 cruises for the years between 1991 and 2020 (Table S3, Fig. S2a). Except for a few periods (January 1993 and January 2002), model–data differences are generally within ±10 µatm with a mean difference of 2.1±7 µatm for the 35 co-located periods. Note that, as opposed to sea surface fCO₂, no temporal trend was identified for the differences between the observed and reconstructed fCO₂ (Fig. S2b); i.e., the trends of sea surface fCO₂ derived from the observations and from the FFNN model should be the same. Aside from the fCO₂ reconstructions, surface ocean alkalinity (A_T) fields are also provided by using the multivariate linear regression model LIAR (Carter et al., 2016, 2018) based on sea surface temperature, salinity, and nutrient concentration.

2.2.5 Calculations of carbonate properties

Based on the data available for each cruise (fCO₂, or A_T and C_T) or from the FFNN model (fCO₂ and A_T), other carbonate system properties (pH, [H⁺], [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], and Ω) were calculated using the CO2sys program (version CO2sys_v2.5; Orr et al., 2018) developed by Lewis and Wallace (1998) and adapted by Pierrot et al. (2006) with K1 and K2 dissociation constants from Lueker et al. (2000) as recommended (Dickson et al., 2007; Orr et al., 2015; Wanninkhof et al., 2015). The total boron concentration was calculated according to Uppström (1974) and KSO₄ from Dickson (1990). To calculate the properties with the underway surface fCO₂ dataset, we used the A_T–S relationship based on A_T and C_T data from the OISO cruises over the period 1998–2019 in the southern Indian sector as described by Leseurre et al. (2022):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{A}_{\mathrm{T}}& =\mathrm{64.341}\times S+\mathrm{106.764}(\mathrm{rmse}=\mathrm{7.5}{\mathrm{\mu }\mathrm{mol}\mathrm{kg}}^{-\mathrm{1}},\\ & n=\mathrm{4775}).\end{array}\end{array}$$

The use of other A_T–S relationships (e.g., Millero et al., 1998; Jabaud-Jan et al., 2004; Lee et al., 2006; Carter et al., 2018) would change slightly the A_T concentrations but neither the A_T trend nor the interpretation of the C_T, pH, or Ω trends. However, as salinity is an important predictor in the calculation of A_T, C_T, or pH from fCO₂ data, we have assessed the original underway salinity data and found biases for a few cruises in 1992, 1993, and 1995 (Table S1b). For these cruises or when salinity was not measured, we used the salinity from the World Ocean Atlas (WOA; Antonov et al., 2006) in the SOCAT datasets (Pfeil et al., 2013, identified “WOA” in Table S1b). Monthly fCO₂ and A_T data extracted from the CMEMS-LSCE-FFNN datasets at the station location (50.5° S–68.5° E) over 1985–2020 were used to calculate the carbonate properties in the same way as from observations.

2.2.6 Comparisons of different datasets and the FFNN model

To validate the properties calculated using the fCO₂ data for 1991–2021 or from the FFNN model over 1985–2020, we compared the calculated values (A_T, C_T, pH, [H⁺], [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], Ω) with those calculated from A_T and C_T data measured in the mixed layer at the OISO-KERFIX station occupied in 1985 and between 1993 and 2021. For this comparison, we averaged the continuous underway fCO₂ data selected in a box around the station location (50–51.5° S, 67.5–69° E; yellow box in Fig. 1). Results of the comparisons between various datasets are detailed in the Supplement (Tables S3 and S4). During the period 1993–2021, there are 22 station occupations with co-located underway fCO₂ data for different seasons (but mainly in summer). Since we found a close agreement between measured fCO₂ and the FFNN model (Table S3, Fig. S2), mismatches in all calculated carbonate system properties between the underway fCO₂ dataset and the FFNN model are small, falling within the range of the errors associated with the calculations (Orr et al., 2018). For example, for 35 co-located periods, the mean differences in calculated C_T of 1.5±5 µmol kg⁻¹ or pH of $-\mathrm{0.002}\pm \mathrm{0.008}$ are in the range of the theoretical error of about 5 µmol kg⁻¹ and 0.007 respectively when taking into account measurements errors on salinity, temperature, nutrients, fCO₂, and A_T (Orr et al., 2018). On the other hand, compared to the station data in the mixed layer (Table S4), the calculated A_T using Eq. (1) is slightly higher by about 5 µmol kg⁻¹. This explains the relatively high differences for C_T (mean difference around 8 µmol kg⁻¹) and for pH (mean difference around 0.008) calculated with fCO₂ and the A_T–S relationship. The differences of calculated values with observations are, on average, in the range of uncertainties of the carbonate system calculations using A_T–C_T pairs (error for fCO₂ around 13 µatm and for pH around 0.0144). Importantly, there is no temporal trend for the differences between calculated and observed properties (Fig. S3b). We are thus confident using the selected fCO₂ data for the trend analysis presented in this study. The independent comparison with A_T and C_T measurements in the mixed layer also indicates that the FFNN model results for A_T and C_T are close to the observations (Tables S4, S5, Fig. S4) as well as for calculated pH, [H⁺], [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], Ω_Ca, and Ω_Ar. This somehow validates the use of the FFNN data for the trend analysis over the period 1985–2020 and for different seasons, although the FFNN model was not constrained by in situ fCO₂ before 1991, few data in austral winter since 1991, and no Chl-a satellite data available before 1998. Nevertheless, the model shows a good agreement with observations collected in March 1985 (Table S5, Fig. S4).

Table 2 Trends of oceanic fCO₂ (µatm yr⁻¹), pH total scale (TS per decade), and C_T (µmol kg⁻¹ yr⁻¹) at the OISO-KERFIX location (50°40^′ S–68°25^′ E) in the southern Indian Ocean for different periods based on observations (Obs.) and the FFNN model (FFNN). Standard deviations are given in brackets.

Figure 2 (a) Time series of sea surface fCO₂ observations (µatm) southwest of Kerguelen Islands in 1985–2021 (inset map shows the location of observations selected around station OISO-KERFIX at 50°40^′ S–68°25^′ E). The color dots correspond to five periods of the year (January–February, March–April, July–August, October, and December), and triangles show the average for each month. The monthly sea surface fCO₂ from the FFNN model is presented for the period 1985–2020 (grey line), and the atmospheric fCO₂ is represented by the red line. In March 1985 there was no underway fCO₂ observation, and the triangle corresponds to fCO₂ calculated with A_T and C_T measured in the mixed layer. (b) Trends of atmospheric and oceanic fCO₂ (µatm yr⁻¹) in summer and winter over four different periods based on observations (January) and the FFNN model (January and August). Inset map produced with ODV (Schlitzer, 2018).

3 Results and discussion

3.1 Variability and trend of sea surface fCO₂ and air–sea CO₂ fluxes: 1985–2021

The fCO₂ observations around 50° S–68° E and their mean values for each cruise are shown in Fig. 2a. fCO₂ measurements are available for different seasons since 1991, though most of them stem from austral summer (January–February). During austral summer, the ocean fCO₂ was generally lower than in the atmosphere (i.e., the ocean was a CO₂ sink), whereas from July to October it was near equilibrium. The same seasonal change is obtained from the FFNN model for the period 1991–2020 (Fig. 2a). The model also indicates that between 1985 and the mid-1990s the fCO₂ during austral winter (May–September) was always higher than the atmospheric fCO₂, leading to an annual CO₂ source during this period (Fig. 3). In 1985 the oceanic fCO₂ from the FFNN model was higher than in the atmosphere from March to October (Fig. S4), resulting in an annual CO₂ source of +0.8 molC m⁻² yr⁻¹. The model estimates a decrease of the annual CO₂ source until the end of the 1990s followed by an increase of the source over the following decade (Fig. 3). Around the year 2010, the annual CO₂ flux was around +0.5 molC m⁻² yr⁻¹ and then decreased over the last decade to change into an annual CO₂ sink that increased to reach −0.5 molC m⁻² yr⁻¹ in 2020. For this reason and given the data available since 1991, we evaluated the summer and winter trends in fCO₂, C_T, and pH from the FFNN model over the three periods 1991–2001, 2001–2010, and 2010–2020 and compared the summer trends with those deduced from observations (Table 2). The analysis of trends and their associated drivers for different seasons and periods will allow us to explore links with the variability of primary production and/or the Southern Annual Mode (SAM). Shifts from a negative to a positive SAM index (Fig. 3) may have strengthened the upwelling of deep waters and could therefore impact ocean properties throughout the water column including C_T, nutrients, primary production, or pH (e.g., Lovenduski and Gruber, 2005; Lenton et al., 2009; Hauck et al., 2013; Hoppema et al., 2015; Pardo et al., 2017).

Figure 3 Time series of the SAM index in the Southern Ocean (in grey), wind speed (in black, m s⁻¹), and air–sea CO₂ flux (molC m⁻² yr⁻¹) from the FFNN model (in red) at location 50.5° S–68.5° E. A positive (negative) flux represents a CO₂ source (sink). Wind speed and SAM are presented for a 24-month running mean based on monthly values. Note the positive SAM (>0.5) in 1998–2002 and 2010–2019. SAM data from Marshall (2003), http://www.nerc-bas.ac.uk/icd/gjma/sam.html (last access: 14 August 2021). Wind speed data from ERA5 (Hersbach et al., 2020).

From the first underway measurements obtained at the OISO-KERFIX site in February 2021 to the last measurements used in this study in February 1991, the average oceanic fCO₂ increased by +50.5 µatm (from 344.4±1.5 to 394.9±1.5 µatm; Fig. 2a). During the same period, the atmospheric CO₂ increased by 57 µatm in this region (recorded at the Crozet Islands; Dlugokencky and Tans, 2022). This first comparison of two cruises 30 years apart indicates that the oceanic fCO₂ increase was close to that of the atmosphere. During the same period, we observed small variations in A_T (average $A_{\mathrm{T}}=\mathrm{2276.5}\pm \mathrm{4.5}$ µmol kg⁻¹) and a clear increase in C_T (Figs. 4a and S5). This suggests that most of the change observed in oceanic fCO₂ and C_T over the last 30 years is due to the uptake of anthropogenic CO₂. However, the evolution of air–sea CO₂ fluxes (Fig. 3) suggests that other mechanisms were at play over shorter periods, and changes in the air–sea fCO₂ disequilibrium (Fig. 2a) suggest that different drivers may be involved in summer and in winter.

Figure 4 (a) Time series of surface C_T (µmol kg⁻¹) around station OISO-KERFIX at 50°40^′ S–68°25^′ E calculated from fCO₂ data (Fig. 2) using the A_T–S relation (see Sect. 2.2.5). The color dots correspond to five periods of the year (January–February, March–April, July–August, October, and December), and triangles show the average for each month. The monthly sea surface C_T from the FFNN model is presented for the period 1985–2020 (grey line). The annual C_T trend of $+\mathrm{0.58}\pm \mathrm{0.05}$ µmol kg⁻¹ yr⁻¹ (dashed line) is derived from the FFNN monthly data. In March 1985 the triangle corresponds to the observed C_T in the mixed layer. (b) Trends of sea surface C_T (µmol kg⁻¹ yr⁻¹) in summer and winter over four different periods based on observations (for January) and the FFNN model (for January and August). The trend for C_ant (µmol kg⁻¹ yr⁻¹) is also shown (red bars) based on estimates in the winter water.

Summer data are characterized by a strong inter-annual variability in both fCO₂ and C_T (Figs. 2a and 4a) with the ocean being a CO₂ source in January 2002, but a strong sink in January 1993, 1998, 2014, 2016, and 2019. In January 1998, when the surface ocean experienced a warm anomaly (Jabaud-Jan et al., 2004), the low fCO₂ of 337 µatm and the low C_T of 2110 µmol kg⁻¹ (Figs. 4a and S5) co-occurred with intense primary production (Fig. 5), probably supported by diatoms as suggested by very low DSi concentrations (<2 µmol kg⁻¹ down to 100 m, Fig. S6). In January 2014 and 2016, mixed-layer DSi concentrations were also remarkably small (<5 µmol kg⁻¹ down to 75 m; Fig. S6). In 2014 low DSi coincided with Chl-a levels that started to increase in mid-November 2013 and stayed at a high level until February 2014 (surface Chl a >0.3 mg m⁻³; Figs. 5 and S7). The intense primary production contributed to the low fCO₂ of 365 µatm reached by mid-January 2014, a value as low as 10 years earlier (Fig. 2a). To the contrary, in 2002 relatively low Chl a (mean Chl a <0.2 mg m⁻³; Fig. 5) was associated with higher levels of fCO₂ (373 µatm), C_T (2128 µmol kg⁻¹; Figs. 4a, S5a), and DSi (Fig. S6). This was also associated with higher salinity indicative of entrainment that might be related to storm events that would have occurred few days before the measurements leading to brief positive fCO₂ anomaly as recently observed from glider data in the subpolar South Atlantic (Nicholson et al., 2022). As opposed to the other periods the ocean was a source of CO₂ in summer 2002 (this particular year was not well reconstructed by the FFNN model; Figs. 2a and S2b). The important inter-annual variability observed in summer indicates that in this region historically referred to as HNLC (Minas and Minas, 1992), primary production could significantly impact fCO₂ levels in summer (Jabaud-Jan et al., 2004; Pasquer et al., 2015; Gregor et al., 2018), a result that needs to be taken into account when evaluating drivers of inter-annual variability (Rustogi et al., 2023) and the decadal trends of fCO₂ or pH.

Figure 5 Time series (1998–2021) of sea surface Chl a (mg m⁻³) in summer (December–February) from weekly satellite data (SeaWiFS and MODIS, open triangles) and associated mean (black triangles). The trends in 1998–2010 and 2010–2021 of respectively $-\mathrm{0.0099}\pm \mathrm{0.0041}$ and $+\mathrm{0.0078}\pm \mathrm{0.0032}$ mg m⁻³ yr⁻¹ (dashed grey) indicate a decrease or increase of the primary production. The full Chl-a record is shown in Fig. S7.

The Chl-a time series derived from MODIS suggests higher concentrations in recent years compared to 2002–2013, with Chl-a peaks identified in 2014, 2016, 2018, 2019, and 2021 (Figs. 5 and S7) when the oceanic fCO₂ in summer was well below the atmospheric level (Fig. 2a).

The primary production lowers C_T concentrations and fCO₂, i.e., opposite to the C_T increase from anthropogenic CO₂ uptake. These counteracting processes might explain the relatively stable fCO₂ previously observed in the Indian POOZ in summer 2007–2019 with an annual fCO₂ rate of increase of only $+\mathrm{0.3}\pm \mathrm{0.2}$ µatm yr⁻¹ (Leseurre et al., 2022). This low rate is confirmed here with the recent data obtained in 2020–2021 (Figs. 2b and S8). For the period 2010–2021, the oceanic fCO₂ trend in summer derived from observations and the FFNN model is lower than +1 µatm yr⁻¹ (Table 2), i.e., much lower than the atmospheric fCO₂ rate of +2.4 µatm yr⁻¹ and the oceanic fCO₂ trend of $+\mathrm{2.21}\pm \mathrm{0.17}$ µatm yr⁻¹ estimated in winter by the FFNN model (Table 2, Fig. 2b). This rate is also lower compared to the change observed in October (+2.9 µatm yr⁻¹) albeit being only based on two cruises in October 2011 and 2016 (Fig. 2a). As the low fCO₂ trend in recent years is detected for summer only this is likely linked to an increase in primary production, as suggested by Chl-a records (Fig. 5). From 1998 to 2010 the summer Chl-a concentrations decreased at a rate of $-\mathrm{0.099}\pm \mathrm{0.041}$ mg m⁻³ per decade, whereas from 2020 to 2021 Chl a increased by $+\mathrm{0.078}\pm \mathrm{0.032}$ mg m⁻³ per decade (Fig. 5). These trends are coherent with previous studies, e.g., the reduced net primary productivity reported in the Indian Antarctic zone over 1997–2007 (e.g., Arrigo et al., 2008; Takao et al., 2012) and the shift of the Chl-a trend in 2010 also reported at large scale in the HNLC region of the Southern Ocean (Basterretxea et al., 2023). As a consequence, after 2010 the difference between oceanic and atmospheric fCO₂ (ΔfCO₂= fCO${}_{\mathrm{2}}^{\mathrm{oce}}$–fCO${}_{\mathrm{2}}^{\mathrm{atm}}$) decreased in summer (−1.4 µatm yr⁻¹), and as it remains relatively steady during winter, the annual CO₂ flux progressively varied from a source of +0.45 molC m⁻² yr⁻¹ in 2010 to a sink of −0.63 molC m⁻² yr⁻¹ in 2020 (Fig. 3). In addition, as the wind speed was stable during this period (12.0±0.9 m s⁻¹ on average in 2010–2020; Fig. 3), the variation of the air–sea CO₂ flux was mainly controlled by ΔfCO₂ (e.g., Gu et al., 2023), and the decadal variation of primary production imprinted a significant change on the fCO₂ trend and air–sea CO₂ flux in this HNLC region. In the region investigated here, increasing Chl-a levels co-occurred with shifts of the SAM index to a positive state (Fig. 3), a link previously suggested south of the polar front in the SO but for a short period over 1997–2004 (Lovenduski and Gruber, 2005). Modeling studies also suggest that summertime biological activity could play an important role for the variability of the CO₂ sink in the SO in response to the SAM (Hauck et al., 2013).

Another process to take into account for interpreting fCO₂ trends is the change in temperature in surface waters. Previous analysis suggested a progressive warming in the region investigated here (Auger et al., 2021 for summer 1993–2017). Over 1998–2019 Leseurre et al. (2022) estimated a warming of Indian POOZ surface waters of $+\mathrm{0.03}\pm \mathrm{0.02}$ °C yr⁻¹. Extending the time series for the period 1991–2021 (Fig. S9a), we note that the surface temperature presents sub-decadal variability and that the ocean cooled after 2018 with a trend of $-\mathrm{0.47}\pm \mathrm{0.16}$ °C yr⁻¹ over 2018–2021 based on the monthly sea surface temperature (SST, Fig. S9b). The trend derived from our in situ observations in summer over this period was $-\mathrm{0.25}\pm \mathrm{0.09}$ °C yr⁻¹.

In 2019, the lower temperature and relatively high Chl a led to low fCO₂ (380 µatm, Fig. 2a) and low C_T (2128 µmol kg⁻¹) compared to 2018 (fCO₂=386 µatm; C_T=2137 µmol kg⁻¹, Fig. 4a). The decrease in observed fCO₂ from summer 2018 to 2019, also reconstructed by the FFNN model (Fig. 2a), is contrary to the expected fCO₂ and C_T increase due to anthropogenic uptake. In 2020, although the temperature was also lower than in 2019, the oceanic fCO₂ was higher (392 µatm) probably due to lower primary production as suggested by higher DSi (Fig. S6), as well as from C_T (2135 µmol kg⁻¹, Fig. 4a) and Chl-a records (Fig. 5). In January 2021 the temperature was close to that in January 2020, and both fCO₂ and C_T were slightly higher (395 µatm, 2139 µmol kg⁻¹). A_T concentrations were stable between 2018 and 2021 (2278.9±1.8 µmol kg⁻¹, Fig. S5), indicating no effect of A_T on the observed fCO₂ change in this region as opposed to the areas north of the polar front in the Indian Ocean where A_T variations are often linked to coccolithophore blooms (Balch et al., 2016; Smith et al., 2017).

The inter-annual and pluri-annual variability observed over 1991–2021 highlights the competitive processes that drive C_T, fCO₂, or pH temporal variations. In order to separate natural and anthropogenic contributions, the anthropogenic CO₂ signal is estimated in the following section.

3.2 Anthropogenic CO₂

3.2.1 Anthropogenic CO₂ in the water column

To calculate anthropogenic CO₂ concentrations (C_ant), we used the TrOCA method developed by Touratier et al. (2007) and previously applied in the southern Indian Ocean (Mahieu et al., 2020; Leseurre et al., 2022). Such an indirect method is not suitable for evaluating C_ant concentrations in surface waters due to biological activity and gas exchange, and we restrict the C_ant calculations below the productive layer around 150 m. In the region south of the polar front, a well-defined subsurface temperature minimum is observed each year characterizing the winter water (WW) at a depth range of 150–250 m (Fig. 6a).

Figure 6 Hovmöller section (depth–time) of (a) potential temperature (°C) and (b) anthropogenic CO₂ (C_ant, µmol kg⁻¹) over 1985–2021 at station OISO-KERFIX (50°40^′ S–68°25^′ E). The section for temperature is presented in the layer 0–500 m and for summer to highlight the temperature minimum around 200 m (winter water, WW). The section for C_ant is limited below 200 m. Section produced with ODV (Schlitzer, 2018).

The C_T and C_ant concentrations increased over time in the water column, a signal that is most pronounced in the top layers (200–400 m, Fig. 6b). In the deep layer, the presence of the Indo-Pacific Deep Water (IPDW) around 600–800 m is identified by a maximum of C_T (C_T>2250 µmol kg⁻¹) and a minimum of O₂ (O₂ close to or <180 µmol kg⁻¹, Fig. S10) (Talley, 2013; Chen et al., 2022). In the IPDW layer restricted to the neutral density (ND) range 27.75–27.85 kg m⁻³, there is no significant change in C_T over time (Fig. S10). In that layer the C_ant concentrations in 1985 (17.3 µmol kg⁻¹) were almost identical to those evaluated in 2021 (21.2 µmol kg⁻¹), considering the uncertainty in the C_ant calculations (±6.5 µmol kg⁻¹; Touratier et al., 2007). Below 800 m, the C_ant concentrations were small but not null (Fig. 6b). The average C_ant concentration below 800 m for all years and seasons was 8.0±5.3 µmol kg⁻¹ (n=123) with a very small change detected over time (C${}_{\mathrm{ant}}=\mathrm{7.7}\pm \mathrm{1.3}$ µmol kg⁻¹ in 1985 and C${}_{\mathrm{ant}}=\mathrm{10.4}\pm \mathrm{0.6}$ µmol kg⁻¹ in 2021). As discussed above (Sect. 2.2.3) the C_T and A_T concentrations in the bottom layer (>1450 m) were stable over 1985–2021 (Table S2, Fig. S1).

Figure 7 (a) Time series of anthropogenic CO₂ (C_ant µmol kg⁻¹) estimated in the winter water layer (WW around 200 m; see Fig. 6) from 1985 to 2021 at station OISO-KERFIX (50°40^′ S–68° 25^′ E). Black dots are the individual data in the WW and the grey squares the average for the 1990s, 2000s, and 2010s (anomalies in 1998, 2005, and 2020 discarded). The C_ant trend of $+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹ is represented (dashed line). (b) Same data for C_ant versus atmospheric fCO₂ (the slope is $+\mathrm{0.263}\pm \mathrm{0.042}$ µmol kg⁻¹ µatm⁻¹).

3.2.2 Anthropogenic CO₂ trend in the winter water

To separate the natural and anthropogenic signals in surface waters for the driver analysis, we assume that C_ant in the WW is representative of C_ant in the mixed layer (ML). This is confirmed with few stations occupied during winter showing that C_ant concentrations in the WW in summer are almost equal to C_ant in the ML during the preceding winter (Fig. S11). The evolution of C_ant in the WW from 1985 to 2021 is presented in Fig. 7a for all seasons. In 1985 the C_ant concentration in the WW was 47.1 µmol kg⁻¹, and C_ant reached a maximum of 71.7 µmol kg⁻¹ in 2021. The data selected at 200 m present some inter-annual variability such as the relatively low C_ant in 1998, 2005, and 2020 probably related to natural variability. In 1998 and in 2020 the O₂ concentrations were slightly lower in the WW (<300 µmol kg⁻¹), explaining the lower C_ant concentration (44.8 µmol kg⁻¹ in 1998 and 53.8 µmol kg⁻¹ in 2020), but no anomaly was observed for C_T. This suggests that the biological contribution may have been overestimated (lower O₂ is interpreted by the TrOCA method as more organic matter remineralization which should be associated with higher C_T). This could be instead related to a change in mixing or circulation. In 2005 anomalies of C_T, A_T, and O₂ concur to explain the lower C_ant (43.9 µmol kg⁻¹).

From 1985 to 2021, we estimated a C_ant trend in WW of $+\mathrm{0.49}\pm \mathrm{0.09}$ µmol kg⁻¹ yr⁻¹. When the C_ant anomalies in 1998, 2005, and 2020 were discarded, this C_ant trend was $+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹ (Fig. 7a). As expected, the C_ant concentrations in the ocean are positively related to atmospheric CO₂ (slope $+\mathrm{0.26}\pm \mathrm{0.04}$ µmol kg⁻¹ µatm⁻¹; Fig. 7b). Interestingly, the slope observed south of the PF in the Indian Ocean is close to that observed in the Antarctic Intermediate Water (AAIW) in the South Atlantic ($+\mathrm{0.23}\pm \mathrm{0.05}$ µmol kg⁻¹ µatm⁻¹; Fontela et al., 2021). Gruber et al. (2019a, b) evaluated C_ant changes between 1994 and 2007 in the global ocean. In the southern Indian sector, they estimated a mean C_ant accumulation at the surface of $+\mathrm{6.0}\pm \mathrm{1.1}$ µmol kg⁻¹ in the band 50–55° S south of the PF. At our station location (50–52° S, 68° E) in the layer 0–250 m, the C_ant accumulated from 1994 to 2007 was $+\mathrm{5.7}\pm \mathrm{1.5}$ µmol kg⁻¹. In 13 years, this corresponds to a trend of $+\mathrm{0.44}\pm \mathrm{0.11}$ µmol kg⁻¹ yr⁻¹. Gruber et al. (2019a, b) did not use the data presented here, allowing for an independent comparison to the present study. Estimates of C_ant accumulation by Gruber et al. (2019a, b) are in agreement with ours for the period 1994–2007 ($+\mathrm{0.49}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹) but lower than reported here in recent years ($+\mathrm{0.61}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹ over 2008–2021). Indeed our estimates over 3 decades indicate an increase in the uptake of anthropogenic CO₂ with time (Fig. 4b).

3.2.3 Anthropogenic and surface C_T seasonal trends

The C_ant trend in the WW over 1985–2021 ($+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹) is slightly lower than the annual surface C_T trend derived from the FFNN model for 1985–2020 ($+\mathrm{0.58}\pm \mathrm{0.05}$ µmol kg⁻¹ yr⁻¹; Fig. 4a, Table 2) suggesting that anthropogenic CO₂ uptake explains 86 % of the C_T increase in surface waters. Over 1991–2020 the surface C_T trend appears slightly higher in January than in August (Fig. 4b, Table 2). This suggests that in addition to the increase of C_T due to anthropogenic CO₂, other processes such as the variability of the biological activity, vertical mixing, or upwelling contributed to the observed trend. Indeed, as for fCO₂ (Fig. 2b), the C_T growth rate also depends on seasons and decades (Fig. 4b). Over 1991–2001 the C_T trend from the observations ($+\mathrm{0.05}\pm \mathrm{0.64}$ µmol kg⁻¹ yr⁻¹; Table 2) is highly uncertain due to few data and the large variability (Fig. 4a, b). The FFNN model showed that the C_T trend in summer was faster than the trend in C_ant (Fig. 4b), suggesting that natural processes would have increased C_T. This could be explained by an increase in vertical mixing due to the increase in wind speed (Fig. 3). On the contrary, the winter C_T trend was lower than the C_ant trend estimated in subsurface waters (Fig. 4b).

Over 2001–2010 the C_T trends were much faster than over the previous decade, and they were the same for both seasons (Fig. 4b, Table 2). For this decade the summer C_T trends from the observations and the FFNN model are coherent. They were also twice the C_ant rate in the WW, which could be explained by enhanced upwelling of C_T-rich deep waters during this period after the SAM reached a high positive index (Fig. 3; Lenton and Matear, 2007; Le Quéré et al., 2007; Hauck et al., 2013). However, over this period we did not detect any clear change at depth for ocean properties (except for C_T and C_ant) that would support this assumption (enhanced upwelling). The rapid C_T (and fCO₂) trend for this decade is probably due to processes occurring at the surface (e.g., biological activity, as discussed later) rather than changes in the water column (vertical mixing or upwelling). Over the last decade C_T trends were lower than over the previous one (Fig. 4b). For summer, this is identified from both observations and the FFNN model. In winter the C_T trend (from FFNN) is close to C_ant, indicative of the anthropogenic CO₂ accumulation. The low C_T trend at the surface in summer, about half the C_ant trend, is likely due to the increase of primary production after 2010 as described above (Fig. 5). Thus, it appears that the impact of biological activity and its variability in summer could counteract that of anthropogenic CO₂ and explain the low temporal change of the carbonate system at the surface in recent years.

Figure 8 (a) Time series of the seasonal amplitude (August minus January) for surface C_T (black; µmol kg⁻¹) and fCO₂ (grey; µatm) from the FFNN model at station OISO-KERFIX (50°40^′ S–68°25^′ E). Also shown are the mean surface Chl a values (green; mg m⁻³) in summer from 1998 to 2021. (b) Seasonal amplitude of fCO₂ and C_T versus summer Chl a over 1998–2020. The dashed lines indicate that the seasonal amplitude (August–January) increases when Chl a is higher.

Given the differences of the fCO₂ and C_T trends in summer and winter (Figs. 2b and 4b, Table 2), we explored the temporal variations of the seasonality. For each year we estimated the differences between August and January (Fig. 8a). The seasonal amplitude for C_T was on average 26.1±3.4 µmol kg⁻¹ and for fCO₂ 15.1±5.6 µatm. Some large inter-annual variations appear related to the variability of Chl a in summer (Fig. 8a). Interestingly, the fCO₂ seasonal amplitude reached a minimum around 2008–2010 and then increased over 2010–2020. This signal also appears correlated with the evolution of surface Chl a in summer (Fig. 8). This supports the conclusion that low phytoplanktonic biomass between 2008 and 2010 reduced the seasonal amplitude of fCO₂.

The inter-annual variability of the seasonality is clearly identified when comparing C_T with C_T calculated due only to C_ant accumulation (Fig. S12). This supports the conclusion that in addition to the C_ant accumulation, the variations of phytoplanktonic biomass imprinted inter-annual variability on C_T and fCO₂ in summer. This holds for the seasonal amplitude as the results for winter follows the C_ant trend (Figs. 4b, S12a). The same is true for pH for which reduced seasonal amplitude was found when the production was low (not shown). However, over 36 years (1985–2020) we did not identify a long-term trend of the seasonal amplitude for C_T or for fCO₂ as suggested by other studies (Landschützer et al., 2018; Rodgers et al., 2023; Shadwick et al., 2023). Our results highlight a variability over 5–10 years (Fig. 8a) and suggest a potential change in seasonality and annual CO₂ sink if primary production changes in the future (e.g., Bopp et al., 2013; Leung et al., 2015; Fu et al., 2016; Kwiatkowski et al., 2020; Krumhardt et al., 2022; Seifert et al., 2023).

Figure 9 (a) Time series of surface pH (total scale, TS) around station OISO-KERFIX(50°40^′ S–68°25^′ E) calculated from fCO₂ data (Fig. 2) using the A_T–S relation (see Sect. 2.2.5). The color dots correspond to five periods of the year (January–February, March–April, July–August, October, and December), and triangles show the average for each month. The monthly sea surface pH from the FFNN model is presented for the period 1985–2020 (grey line). The annual pH trend in 1985–2020 of $-\mathrm{0.0165}\pm \mathrm{0.0004}$ per decade (dashed line) is derived from the FFNN monthly data (the same figure for [H⁺] concentrations is presented in Fig. S13). (b) Trends of pH (TS per decade) in summer and winter over four different periods based on observations (January) and the FFNN model (January and August).

3.3 Anthropogenic CO₂ drives acidification in surface waters and in the water column

3.3.1 Surface pH trend

To explore the temporal change of pH in surface waters, we used the fCO₂ observations and the monthly results from the FFNN model. For both datasets pH was calculated from fCO₂ and A_T reconstructed as described in Sect. 2.2.5. Figure 9a presents the time series of pH at the surface (the same time series for [H⁺] concentrations is shown in Fig. S13). For the full period, 1985–2020, the annual pH trend derived from the FFNN model is $-\mathrm{0.0165}\pm \mathrm{0.0004}$ per decade (Table 2), exactly the same as derived at large scale in the Southern Ocean (south of 44° S) for the period 1993–2018 (Iida et al., 2021, Table 1), but when restricted to this period, 1993–2018, the trend from the FFNN model appears slightly faster ($-\mathrm{0.0182}\pm \mathrm{0.0006}$ per decade). This is less than the pH trend derived from pCO₂ data in the SO subpolar seasonally stratified biome around 40–50° S (SO-SPSS) for 1981–2011 ($-\mathrm{0.020}\pm \mathrm{0.002}$ per decade; Table 1; Lauvset et al., 2015) and close to the pH trend based on OceanSODA-ETH reconstructed fields in the SO-SPSS for the period 1982–2021 ($-\mathrm{0.0189}\pm \mathrm{0.0010}$ per decade; Ma et al., 2023). However, as for fCO₂ and C_T, we estimated different pH trends in summer and winter, depending on the periods (Fig. 9b, Table 2).

The winter pH decrease estimated over the last 2 decades was twice as fast as estimated during the previous one, mirroring the winter fCO₂ trends (Table 2). In summer, the pH trend presents a large variability at decadal scale as it was 3 times faster over 2001–2010 than during the previous and following decades (Fig. 9b, Table 2). Although the trends based on the observations are less robust because the cruises were not conducted each year, the reduced pH trend in summer after 2010 is confirmed from in situ data (Fig. 9b, Table 2).

Figure 10 (a) Time series of pH (red dots) and pre-industrial pH (pH-PI; black dots) estimated in the winter water layer (WW around 200 m; see Fig. 6) over 1985–2021 at station OISO-KERFIX (50°40^′ S–68°25^′ E). pH-PI for each sample was calculated after subtracting C_ant to C_T. The pH trend from the present day is $-\mathrm{0.0161}\pm \mathrm{0.0033}$ per decade (dashed red line). No trend is observed for pH-PI (dashed black). The mean pH-PI in the WW is 8.173±0.020 (n=45). (b) Profiles of pH and pH-PI evaluated from March 1985 (black symbols) and February 2021 data (red symbols). The profiles for pH-PI are shown below 150 m only as C_ant estimates are not available in the surface layer. Note that the pH-PI profiles are the same either using either the 1985 or 2021 data.

Our results show that the pH trend varied significantly from decade to decade and that part of the variations could be explained by the evolution of phytoplanktonic biomass, but overall the decrease of pH since 1985 was mainly driven by the accumulation of anthropogenic CO₂. This is revealed in the winter water when comparing pH and pre-industrial pH (Fig. 10a). Here, the pre-industrial pH (pH-PI) was calculated after subtracting C_ant values from the observed C_T concentrations for each sample in the WW layer. Interestingly, the pH trend in the WW of $-\mathrm{0.0161}\pm \mathrm{0.0033}$ per decade (here deduced from the station A_T and C_T data over 1985–2021) is very close to the annual trend at the surface deduced from the FFNN model over 1985–2020 ($-\mathrm{0.0165}\pm \mathrm{0.0004}$ per decade). This trend is slightly faster than the pH trends of $-\mathrm{0.0134}\pm \mathrm{0.001}$ per decade recently estimated in subsurface waters (100–210 m) of the Southern Ocean south of the PF and derived for years 1994–2017 from historical data and BGC-Argo floats (Mazloff et al., 2023). For the same period, 1994–2017, at the OISO-KERFIX station we estimated a pH trend in the WW of $-\mathrm{0.0168}\pm \mathrm{0.0043}$ and of $-\mathrm{0.0186}\pm \mathrm{0.0006}$ per decade in surface waters from the FFNN model.

As for other properties (A_T, O₂, temperature, salinity, and nutrients), the pre-industrial pH (pH-PI) does not change over time in the WW (mean pH-PI $=\mathrm{8.173}\pm \mathrm{0.020}$, n=45; Fig. 10a). The pH-PI in the WW is in the range of the pre-industrial surface pH value in the Southern Ocean (8.2 for year 1750 and 8.18 for year 1850) derived from Earth system models (Jiang et al., 2023, their Table S9). In the WW at our location the modern pH (1985–2021) was on average $-\mathrm{0.147}\pm \mathrm{0.021}$ lower than pre-industrial pH. In 1985 pH in the WW was −0.119 lower than pH-PI, and in 2021 it was −0.184 lower than pH-PI (Fig. 10a). The progressive decrease of pH was clearly linked to C_ant concentrations in the WW layer and the pH decrease identified below that layer in the water column (Fig. 10b).

Figure 11 Anthropogenic pH (pH_ant) and anthropogenic ΩAr (ΩAr_ant) versus anthropogenic CO₂ concentrations (C_ant, µmol kg⁻¹) at station OISO-KERFIX (50°40^′ S–68° 25^′ E). The data are selected in the layer 150–1600 m for the periods 1985 (blue), 1993–2020 (black), and 2021 (red). The green arrow identifies the data in the layer 150–500 m (for C_ant>30 µmol kg⁻¹). Below 500 m (brown arrow) no change of C_ant was observed from 1985 to 2021 and thus for pH_ant and ΩAr_ant.

3.3.2 Temporal change in the water column

From 1985 to 2021, signals of decreasing pH and increasing C_T in surface waters are propagated in the water column down to about 500 m. As mentioned above the data in 1985 (first occupation of the station) reveal significant C_ant levels across the water column (Fig. 6b). Therefore the pH down to the bottom was already lower in 1985 than at pre-industrial times (Fig. 10b). However, the largest C_ant increases were found in the top layers, and changes in pH from 1985 to 2021 were small below 500 m (Figs. 10b, S14). While observations for all years fall on a common linear relationship between C_ant and pH_ant for depths greater than 500 m, the change in pH for a given level of C_ant increases with time for layers shallower than 500 m (Fig. 11).

The increase in C_ant concentrations over time (Fig. 6b) also leads to a decrease of carbonate ion concentrations [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] and of ΩAr and ΩCa (Figs. S14, S15). These decreases are well identified since the pre-industrial era in the whole water column, but in the last 36 years observations do not show any appreciable changes below 500 m (Fig. 11). The aragonite saturation horizon (ΩAr=1) was found around 600 m in 1985 and around 400 m in recent years (2015–2021; Figs. S14, S15). Moreover, during the period covered by observations (1985–2021), we did not detect abrupt change of the aragonite saturation horizon from one year to the next (nor between winter and summer; Fig. S16). This contrasts with previous regional studies in the SO and most notably with results from the layers close to the deep minimum of carbonate ion concentrations (Hauri et al., 2015; Negrete-Garcia et al., 2019). At our station the [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] minimum lies around 500–600 m (Figs. S14, S15) and, along with the superimposed C_ant accumulation, explains the upward shift of the aragonite and calcite saturation horizon between the pre-industrial and modern periods (Fig. S15). At pre-industrial time undersaturation with regard to aragonite (ΩAr <1) was found at the bottom only (1600 m), whereas between 1985 and 2021 it was found in the water column below 600 or 400 m (Fig. S15). The subsurface pre-industrial ΩAr value was around 1.9–2 (Fig. S15) and in the range of ΩAr value in the Southern Ocean at pre-industrial time from ESMs (Jiang et al., 2023, their Fig. 4).

The aragonite undersaturation already occurred in 1985 at 600–700 m, a layer corresponding to the [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] minimum (Fig. S15), and a small increase of C_T just above this layer (via C_ant accumulation) would rapidly shift the aragonite saturation horizon above 600 m. This might have already occurred and could explain that ΩAr value was 1.02 at 350 m in 2021 (Fig. S15). These results suggest that for pelagic calcifiers living in subsurface waters (150 m or deeper) such as pteropods and foraminifera (e.g., Hunt et al., 2008; Meilland et al., 2018), the impact of acidification might occur sooner than at the surface.

For the interpretation of the trend analysis based on observations, only data below 150 m could be used as C_ant was not evaluated in the surface layer. At 200 m, based on A_T and C_T data, we estimated a decrease in pH from 1985 to 2021 by −0.059 (Fig. 10b), corresponding to an increase by +1.1 nmol kg⁻¹ in [H+] (Fig. S13) and a decrease by −0.16 in ΩAr (Fig. S15). Over 36 years, this represents about 30 % of the total change since the pre-industrial era for pH (−0.184), [H⁺] (+3.5 nmol kg⁻¹), and ΩAr (−0.6). This is mainly linked to the C_ant change that also represents over 36 years 30 % of the total accumulation (+24.6 µmol kg⁻¹ from 1985 to 2021 for a total concentration of +71.7 µmol kg⁻¹ at 200 m in 2021; Fig. 7). We conclude that the accumulation of anthropogenic CO₂ drives the change of the carbonate system in subsurface waters and probably also in surface waters.

Figure 12 Profiles (0–400 m, left axis) of observed and calculated properties (C_T, pH, ΩAr, [CO${}_{\mathrm{3}}^{\mathrm{2}-}$]) at station OISO-KERFIX (50°40^′ S–68°25^′ E) in March 1985, January 2020, and January 2021 along with surface time series in 1985–2020 (right axis) of the same properties in January (yellow line) and August (blue line) from the FFNN model. The FFNN values in January 2020 are coherent with January 2020 or January 2021 observations in the mixed layer and in January 1985 are close to the observations in March 1985. Note that the differences of properties between 2020–2021 and 1985 have a similar magnitude as the seasonal amplitude (illustrated by the FFNN values for January and August).

In order to quantify the propagation of surface trends to depth, the temporal variations of carbonate system properties at the surface for both summer and winter derived from the FFNN model are compared to the changes observed across the water column (Fig. 12). The comparison shows that the seasonal amplitude of surface waters properties was of a similar magnitude to the observed changes in the mixed layer between 1985 and 2021. For example, the C_T and ΩAr seasonal amplitude, respectively around 20 µmol kg⁻¹ and 0.2, correspond to the C_T increase and ΩAr decrease from 1985 to 2021. The comparisons also highlight that in summer the FFNN results were close to observations in the mixed layer (e.g., C_T was 2120 µmol kg⁻¹ in 1985 and 2140 µmol kg⁻¹ in 2021). In winter, at the surface, C_T was higher and pH, [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], and ΩAr were lower (from the FFNN model; blue line in Fig. 12). The winter surface values in 1985 and 2020/2021 are in good agreement with observations at depth in the winter water (150–200 m). As an example, in 1985 surface C_T in winter was 2145.5 µmol kg⁻¹, which corresponds to the concentration measured at 150 m during summer (purple line in Fig. 12). In 2020, the winter C_T at the surface (2168.3 µmol kg⁻¹) is equal to C_T concentrations observed at 150–180 m in January 2020 or in 2021. For ΩAr, the surface value derived from the FFNN model in winter 1985 (1.6) was equal to the ΩAr observed at 125 m in March 1985. In 2020, the surface winter estimate of ΩAr (1.42) was equal to ΩAr observed at 100–150 m in January 2020 or 2021. The same correspondences between winter surface and WW data were identified for pH and [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] (Fig. 12). This supports the use of winter and summer surface data from the FFNN model to investigate the seasonal ΩAr trends and their projection in the future.

The surface water ΩAr (ΩCa) trend from the FFNN model in summer of −0.0059 yr⁻¹ (−0.0094 yr⁻¹) was stronger than the winter of −0.0050 yr⁻¹ (−0.0079 yr⁻¹) and also higher than the trend derived from observations in the WW (−0.0043 yr⁻¹ for ΩAr and −0.0069 yr⁻¹ for ΩCa). Our results indicate that the change of carbonate properties in the years 1985–2021 were mainly driven by C_ant accumulation in surface waters and across the water column. However, a potential increase in primary productivity after 2010 mitigated the effects of increasing C_ant accumulation in response to increasing atmospheric CO₂ leading to relatively stable summer C_T and fCO₂ and to a stronger CO₂ sink (Fig. 3). Consequently, when restricted to the period 2010–2020, the trend of ΩAr in surface waters in summer was much smaller, $-\mathrm{0.024}\pm \mathrm{0.027}$ per decade than during the preceding period. This was much smaller than derived from all the data over 1985–2021 (−0.048 per decade) or estimated from reconstructed fields in the SO-SPSS over 1982–2021 (−0.0616 per decade; Ma et al., 2023). It underscores the uncertainty in extrapolating time series to the future depending on the selection of data and periods.

3.4 Long-term change in surface water: from the 1960s to the future

The data described above allowed evaluating the temporal variations of the properties of the carbonate system and C_ant over 1985–2021 along with a comparison to the pre-industrial state in the water column. The results over 36 years informed on the recent changes, inter-annual variations, and trends, but the time series appears somehow short to extrapolate the trends over time. What was the change of the carbonate system in surface waters before 1985 and what will be its future evolution?

3.4.1 Back to the 1960s: observed trends since 1962

To explore the long-term change, we start by comparing our recent data with observations from the LUSIAD cruise conducted in 1962–1963 (Keeling and Waterman, 1968). Some data from this cruise were obtained in mid-November 1962 south of the polar front, in the region southwest off Kerguelen Islands. Because of the seasonality, we compared the November 1962 data with our observations obtained in October–November in 1995, 2011, and 2016 and with the FFNN model results for November (Fig. 13). The C_T concentration, pH, ΩAr, and ΩCa for 1962 were calculated using fCO₂ data and A_T (from the A_T–S relationship in Eq. 1) with salinity from the World Ocean Atlas (Antonov et al., 2006).

Figure 13 Observed (black dots) sea surface temperature (°C), fCO₂ (µatm), C_T (µmol kg⁻¹), pH (TS), ΩAr, and ΩCa around station OISO-KERFIX at 50°40^′ S–68°25^′ E for October–November. Also shown are the results for the FFNN model for November in 1985–2020 (purple). The C_T concentrations, pH, ΩAr, and ΩCa were calculated from fCO₂ data using the A_T–S relation (Eq. 1). The red line is the atmospheric fCO₂, and dashed red lines in each plot are the evolution of properties since 1960 corrected for C_ant where fCO₂, pH, ΩAr, and ΩCa were recalculated using C_T+C_ant, A_T constant at 2290 µmol kg⁻¹, and SST at 2 °C. Grey triangles identify the mean values for C_T and pH.

First, we note that the measured SST in November 1962 (1.7 °C) was slightly lower compared to recent years (on average by about −0.6 °C), but SST as low as 1.8 °C for this season was also found in other years (e.g., November 1995, 2014). The change in SST is unlikely to explain the long-term increase in fCO₂ or decrease in pH since 1962 (Fig. 13). In 1962, the oceanic fCO₂ was 324 µatm, which is slightly higher than in the atmosphere (ΔfCO${}_{\mathrm{2}}=+\mathrm{8}$ µatm, a small source), whereas in November 1985–2020 the ocean was a small CO₂ sink on average (ΔfCO${}_{\mathrm{2}}=-\mathrm{3.3}\pm \mathrm{4.5}$ µatm). The C_T concentration in 1962 (2135 µmol kg⁻¹) was much lower than observed since 1995, and the pH (8.115) was much higher than in the last 3 decades (Fig. 13). Compared to 1962, pH in 2016 was −0.078 lower, i.e., representing 70 % of the pH decrease of −0.11 in the global ocean since the beginning of the industrial era (Jiang et al., 2019). In November 1962, surface C_T was lower by −15.1 µmol kg⁻¹ compared to the data in October 1995, i.e., a trend of +0.46 µmol kg⁻¹ yr⁻¹ over 33 years close to the C_ant trend observed in the WW over 1985–2021 as described above ($+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹). Having the C_T value in 1962, we can project the C_T in time by adding the C_ant concentration based on the relationship observed between C_ant and atmospheric CO₂ (Fig. 7b) assuming that the anthropogenic CO₂ uptake since the 1960s is representative of the C_T change (i.e., the change of C_T due to natural variability was small). This projection is shown for all properties (dashed red lines in Fig. 13) and confirms that the progressive C_ant accumulation explained most of the C_T and fCO₂ increase in surface waters since 1962. We note that the C_T derived from the FFNN model suggests slightly lower C_T compared to the C_ant projection especially before 2006. The difference of projected C_T and the FFNN model (on average $-\mathrm{2.2}\pm \mathrm{2.7}$ µmol kg⁻¹) is within the uncertainty of C_T calculations (error is ±5 µmol kg⁻¹ when using the A_T/fCO₂ pairs), and the trend of the difference over 1985–2020 (−0.15 µmol kg⁻¹ yr⁻¹) is too small to be related with confidence to changes associated with natural processes. On the other hand, the oceanic fCO₂ recalculated with the projected C_ant trend suggested that for this season (November) the ocean moved from a CO₂ source in 1962–1985 (ΔfCO₂>0) to a sink in 1986–2021 (ΔfCO₂<0) in line with results from the FFNN model. The recalculated fCO₂ with C_ant (dashed red line in Fig. 13) was close to that observed in 1995 or from the FFNN model in 1985–2014 (mean difference over 1985–2014 is $-\mathrm{1.2}\pm \mathrm{5.2}$ µatm). After 2016, the recalculated fCO₂ suggests a stronger sink, and the difference with observations in 2011 and 2016 or the FFNN model is slightly higher (mean difference over 2016–2020 is $-\mathrm{8.8}\pm \mathrm{1.5}$ µatm). Although the differences are in the range of the error in fCO₂ calculation using A_T–C_T pairs (±13 µatm), this might indicate that after 2016 a process could contribute to increase fCO₂ faster than the effect of C_ant only. This difference could be due to the warming that occurred after 2016 when SST was higher than 2 and up to 3 °C in November 2017 (Figs. 13 and S9). The same could be applied for pH that was slightly lower than the pH recalculated from C_ant trend after 2015 (the mean difference between recalculated pH and FFNN-pH over 1985–2020 is only 0.002±0.006). Therefore, we conclude that for November the pH decrease since 1962 was mainly driven by the accumulation of anthropogenic CO₂. Aragonite and calcite saturation states also show a clear decrease since 1962 (Fig. 13), a diminution of 11 % over 59 years for both ΩAr and ΩCa. Based on these results over almost 60 years that confirm the conclusions from the observations in 1985–2021, we now evaluate the long-term change of the carbonate system in surface waters in the future.

3.4.2 Projecting the observed trends in the future

The trends of the properties based on observations in 1962–2021 and the FFNN model in 1985–2020 indicate relatively linear trends linked to C_ant uptake albeit with some decadal variability in summer (Fig. 4). A simple linear extrapolation of the trends in the future suggests that aragonite undersaturation in surface waters would be reached in year 2110 for the winter season and 2120 for summer (Fig. S17), whereas the subsurface trend suggests undersaturation in 2090. In year 2100, surface pH and [H⁺] would be around 7.9 and 12 nmol kg⁻¹ (Fig. S17). However, ESM CMIP6 models suggest that under a high emission scenario (SSP5-8.5), pH in 2100 in the Southern Ocean near 50° S would be around 7.65 and [H⁺] around 22 nmol kg⁻¹ (Jiang et al., 2023, their Fig. 4). This shows that the simple linear extrapolation based on recent observed trends (Fig. S17) underestimated the future change of the carbonate system for a high emission scenario as previously shown in the southeastern Indian Ocean based on summer trends derived from observations in 1969–2003 (Midorikawa et al., 2012, their Fig. 4).

Table 3 Results of the simulated properties for year 2020, 2050, and 2100 for two emission scenarios (SSP5-8.5 and SSP2-4.5). For 2020 the results based on observations in January (Obs.) and the FFNN model in January and August also listed. Sensitivity tests: “SSP85 W-T” is for winter with constant temperature, and “SSP85 W-A-T” is for winter with constant A_T and temperature.

Figure 14 Evolution of atmospheric CO₂ (ppm) and sea surface C_T (µmol kg⁻¹) between 1960 and 2110 evaluated for two scenarios (SSP2-4.5 dashed black and SSP5-8.5 dashed red), for summer (red line for SSP5-8.5) and winter (blue line for SSP5-8.5 and green line for SSP2-4.5). Grey lines are the high and low C_T for winter SSP5-8.5 based on the error in the C_ant/fCO₂ relationship (Fig. 7b). Also shown are the results for the FFNN model in 1985–2020 for summer (red diamonds) and winter (blue diamonds) and C_T in 1962 (black diamonds). The C_T values for different seasons and scenarios were used to calculate the carbonate properties in the future (Fig. 15).

To better investigate the changes for the next decades, we assumed that the C_ant trend for the modern period (Fig. 7) that experienced a “business as usual” scenario after the 1960s is representative of the future changes in the surface ocean carbonate system. For this analysis, we used two emissions scenarios (shared socioeconomic pathways, SSPs; Meinshausen et al., 2020) with atmospheric xCO₂ reaching 1135 ppm in 2100 (a “high” emission scenario SSP5-8.5) or xCO₂ reaching 603 ppm in 2100 after a stabilization around 2080 (scenario SSP2-4.5). This enables us to simulate future C_T concentrations for summer or winter (Fig. 14) and to calculate other carbonate properties using C_T and A_T (Fig. 15, Table 3) in response to approximated future changes in physical and geochemical properties excluding impacts of changes in atmospheric and oceanic circulation. As the calculated properties are sensitive to A_T values, we used a fixed A_T of 2280 µmol kg⁻¹ or applied a correction based on the long-term change of sea surface salinity observed in the last 6 decades (1960–2017), i.e., a freshening in the Southern Ocean of around −0.01 to −0.02 per decade (Durack and Wijffels, 2010; Cheng et al., 2020b). The decrease in salinity in the southern Indian Ocean ($-\mathrm{0.02}\pm \mathrm{0.01}$ per decade) was recently analyzed by Akhoudas et al. (2023), who showed that in the years 1993–2021 the freshening was mainly due to an increase in precipitation linked to the acceleration of the atmospheric hydrological cycle. From our data in the mixed layer over 1985–2021, we estimated a trend in salinity of $-\mathrm{0.021}\pm \mathrm{0.004}$ per decade. For the A_T sensitivity test, we thus selected a salinity trend of −0.01 per decade in 1962–1985 and −0.02 per decade after 1985 and applied these trends to simulate A_T over 1960–2100 using the A_T–salinity relationship (Eq. 1). This leads to a salinity of 33.650 and A_T of 2272 µmol kg⁻¹ in the year 2100, about 8 µmol kg⁻¹ lower than observed in 2021 (2280 µmol kg⁻¹). Compared to the C_T change from 2021 to 2100 (+50 and +193 µmol kg⁻¹ for the “low” and “high” emissions scenario, Fig. 14), the impact of the A_T decrease has a minor effect on the future change for pH, [CO${}_{\mathrm{3}}^{\mathrm{2}-}$], and Ω (Table 3). For example, in winter for the SSP5-8.5 scenario, when the A_T decrease is taken into account, pH in 2100 is 7.316 and Ω_Ar is 0.33 against 7.372 and 0.34 when A_T is constant (Table 3). In both cases, the surface aragonite undersaturation (Ω_Ar=1) in winter occurred in 2055, whereas in summer it is identified in 2070. The effect of lower A_T in the future appeared also small compared to the seasonal differences of pH and Ω in 2100.

As noted above, the Southern Ocean experienced a warming in recent decades (e.g., Auger et al., 2021), and it is projected that warming will continue in the future (IPCC, 2022). Therefore, to test the sensitivity of calculated properties to warming we applied a warming of +0.0125 °C yr⁻¹ in 1985–2020 and +0.025 °C yr⁻¹ after 2020 (Azarian et al., 2023). As for A_T, these results are compared for winter using constant SST (Table 3). The effect of the long-term warming does mainly impact the projection of [H⁺] and pH (Table 3).

Figure 15 Evolution of sea surface fCO₂ (µatm), pH (TS), [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] (µmol kg⁻¹), [H⁺] (nmol kg⁻¹), ΩAr, and ΩCa between 1960 and 2110 evaluated for the SSP5-8.5 scenario for winter (blue line) and summer (red line) taking into account both A_T and SST future trends. For winter the results are also presented using the SSP2-4.5 scenario (grey lines). Also shown are the results for the FFNN model in 1985–2020 for summer (red diamonds) and winter (blue diamonds). Atmospheric fCO₂ is also shown for SSP5-8.5 (dashed red) and SSP2-4.5 (dashed grey). Values in 2020, 2050, and 2100 for different sensitivity tests are listed in Table 3.

These sensitivity tests for temperature and A_T showed that as for the observed period 1962–2021 (Fig. 13), the projection in the future depends mainly on the anthropogenic CO₂ accumulation. Here, the C_T concentrations were calculated using the C_ant versus atmospheric CO₂ relationship (Fig. 7b). We thus tested the results for winter based on the error associated with this relationship (Fig. S18). This leads to either higher or lower C_T compared to the original calculation (Fig. 14). For the SSP5-8.5 scenario, the winter C_T concentrations in 2100 range between 2328 and 2378 µmol kg⁻¹, higher than simulated in the ESM CMIP6 models around 50° S (2300 µmol kg⁻¹; Jiang et al., 2023). As in the ESMs, the projected C_T concentration in 2100 at our location for the SSP2-4.5 scenario is much lower, 2217 µmol kg⁻¹ (Fig. 14). The future change of the carbonate system is not significantly different using low or high C_ant accumulation (Fig. S18), but this test gives a range of years to reach aragonite and calcite undersaturation. In winter (SSP5-8.5 scenario), aragonite (calcite) would reach undersaturation between year 2050 and 2060 (between year 2070 and 2080). Note that for summer we derived undersaturation for Ω_Ar in year 2065 and for Ω_Ca in year 2085. For the SSP2-4.5 scenario, where C_T is 143 µmol kg⁻¹ lower in 2100 compared to SSP5-8.5, aragonite undersaturation would not be reached before 2070 (Fig. 15).

4 Summary and concluding remarks

The times series of high-quality observations collected between 1985 and 2021 and the results from the FFNN model at one location, south of the polar front in the southern Indian Ocean (50° S–68° E), presented in this analysis offered new results on the inter-annual variability, decadal to long-term trends of the carbonate system in surface waters, air–sea CO₂ fluxes, and associated drivers. The evaluation of anthropogenic CO₂ concentrations in the water column indicated that the trends of the carbonate species are mainly driven by the anthropogenic CO₂ uptake leading to a progressive acidification in surface waters and at depth.

In 1985, the C_ant concentrations were approaching 50 µmol kg⁻¹ at 200 m, and C_ant was detected in the water column down to the bottom (1600 m). This explains why aragonite undersaturation was observed at 600–700 m in 1985, where [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] concentration was at its minimum, whereas for the pre-industrial era the whole water column was super-saturated (this study Fig. S15; Lauvset et al., 2020, their Fig. S15). Because of the anthropogenic CO₂ accumulation, 36 years later, we observed an upward migration of the aragonite saturation horizon that was found around 400 m in 2021 (a shoaling rate of around −6 m yr⁻¹).

At the subsurface, in the winter water layer, the C_ant trend is estimated at $+\mathrm{0.53}\pm \mathrm{0.01}$ µmol kg⁻¹ yr⁻¹ over 1985–2021 with a detectable increase of the trend in recent years. The C_ant concentrations in the ocean are closely related to the atmospheric CO₂ concentrations, and the slope we observed south of the PF in the Indian sector of $+\mathrm{0.263}\pm \mathrm{0.042}$ µmol kg⁻¹ µatm⁻¹ is close to that observed in the AAIW in the South Atlantic (Fontela et al., 2021). This suggests that local observations in the southern Indian POOZ captured the link between C_ant and atmospheric CO₂ at a larger scale.

In surface waters, over 1991–2020 the oceanic fCO₂ increased at a rate close to or slightly lower than in the atmosphere (Fig. 2b), and the C_T trend followed the C_ant accumulation (Figs. 4b, S12a). However in the last decade both observations and the FFNN model showed low fCO₂ trends in summer (less than 1 µatm yr⁻¹). The change in summer trend appears related to primary production as revealed by a decrease of Chl a in 1998–2010 followed by an increase after 2010. Biological activity counteracts the C_T increase due to C_ant, resulting in rather stable C_T and fCO₂ in summer during the last decade. As a result, the region moved from an annual source of +0.8 molC m⁻² yr⁻¹ in 1985 to a sink of −0.5 molC m⁻² yr⁻¹ in 2020. Adding historical data from November 1962 that indicate an annual source of 2.2 molC m⁻² yr⁻¹ and extrapolating to the entire southern Indian POOZ (50–58° S, 20–120° E; 6.5 M km²) suggest that this region changed from a CO₂ source of 0.17 PgC yr⁻¹ in 1962, reduced to 0.06 PgC yr⁻¹ in 1985, to a CO₂ sink of −0.04 PgC yr⁻¹ in 2020. This can be compared with reconstructed fluxes from a data-based model that produced a CO₂ source around 0.10 PgC yr⁻¹ in 1960–1990 and a sink around −0.05 PgC yr⁻¹ in 2020 in the southern Indian sector (Rödenbeck et al., 2022, their Fig. 6). Based on the FFNN reconstructions, the increase of the ocean CO₂ sink was particularly pronounced after 2011 (Fig. 3) when phytoplankton biomass was increasing in this HNLC region and occurred when the SAM index was in a positive state. Although observations in the water column do not suggest enhanced upwelling, we cannot eliminate the possibility that the SAM influenced changes in primary production.

For October/November, the estimated increase in C_T concentration in surface waters over 54 years (+21 µmol kg⁻¹) was almost equal to the increase of C_ant (+22.3 µmol kg⁻¹). As a result, surface ocean pH dropped from 8.11 in 1962 to 8.044 in 2020. Over a multi-decadal timescale (30 years or more), acidification in the southern Indian POOZ was mainly controlled by the uptake of anthropogenic CO₂. However, our data also indicate a modulation of the summer pH trend by natural processes. After 2010, a very small pH trend was estimated in summer (−0.0098 ±0.0042 per decade) when the region experienced an increase in primary productivity. On the opposite, in winter, the pH trends continuously increased with time. At the subsurface (winter water layer), the trend of pH based on A_T and C_T data over 1985–2021 ($-\mathrm{0.0161}\pm \mathrm{0.0033}$ per decade) is also almost equal to the annual surface trend from the FFNN model. A simple extrapolation of the trends in the WW indicated that undersaturation (Ω<1) would be reached at year 2090 for aragonite and year 2180 for calcite. However, as atmospheric CO₂ is expected to increase and ocean C_T will increase in the future, pH and Ω will decrease at a faster rate than observed in the last decades. A projection of future C_T concentrations based on two emission scenarios, excluding changes in ocean circulation, indicated that the winter surface pH in 2100 would decrease to 7.32 for a high emission scenario (SSP5-8.5) or to 7.782 for a low emission scenario (SSP2-4.5). This is up to −0.86 lower than pre-industrial pH and −0.71 lower than pH observed in 2020. For the winter season the aragonite undersaturation in surface would be reached around 2050 for a high emission scenario and 2070 for a low emission scenario.

The time series presented here for the Southern Ocean, along with other historical time series of A_T and C_T in the water column (BATS, HOT, ESTOC, KNOT, Iceland, or Irminger seas; Bates et al., 2014; Lange et al., 2024) or the recent BGC-Argo floats in the Southern Ocean (Mazloff et al., 2023), offer useful data for the evaluation of biogeochemical and Earth system models, especially for the physical and biological drivers of the carbonate system not well represented in current models at seasonal to decadal scales in the Southern Ocean (e.g., Hauck et al., 2023a; Rodgers et al., 2023; Joos et al., 2023). Observing the decadal changes of the carbonate system in the water column is also an important step to extend the evaluation of biogeochemical and ESMs below the surface (Jiang et al., 2023). It is important to maintain such time series for monitoring the future evolution of the ocean CO₂ sink, of the acidification, and its impact on phytoplankton species and higher trophic levels. This is especially the case in marine protected areas such as the French sub-Antarctic islands including the Kerguelen Archipelago, which was listed as a UNESCO World Heritage Site in 2019.