To investigate year-to-year variability in IOSTMW, we used the latest version of the “In Situ Analysis System” (ISAS; Gaillard et al., 2016; Kolodziejczyk et al., 2021). ISAS was created to construct gridded temperature and salinity fields with the aim of fully utilizing the temporal and spatial resolution provided by the Argo network via the Optimal Interpolation method. In this study, we specifically used the latest version (ISAS20), which analyzes temperature and salinity data from the Argo and Deep-Argo programs for the period from 2002 to 2020. Previous studies have demonstrated the utility of ISAS for detecting signals from mode waters in other ocean basins (Kolodziejczyk et al., 2019; Bernardo and Sato, 2020). In this study, the ISAS20 dataset was used to identify layers with weak temperature gradients in the formation and distribution region of IOSTMW.

Additionally, ISAS provides the percentage of a priori variance (PCTVAR), a metric that reflects the impact of data sampling and coverage, the covariance scale (which determines the region influenced by the data), and measurement errors for the resulting gridded fields. Thus, PCTVAR is a useful parameter by which the coverage of Argo profiles may be assessed. At any given location, a PCTVAR value of 100% indicates the absence of data at or near that location; in such cases, the temperature and salinity values are relaxed toward the climatology (Kolodziejczyk et al., 2021). Figure 1A shows Argo profiles in the western end of the Indian Ocean subtropical gyre (15–65°E, 30–45°S; see inset map in Figure 1A ) used in ISAS-20 from 2002 to 2020. Very few profiles are present at the beginning of 2002, but since then, the number of profiles has steadily grown. Since 2005, approximately 200 profiles have been recorded each month and the number has rapidly increased since 2008. Accordingly, the PCTVAR parameter significantly decreased from almost 100% in 2002 to approximately 70% in 2005 ( Figure 1B ), and the trend continues to decline to almost 50% toward 2011, indicating improved data coverage in the Indian Ocean subtropical gyre region ( Figures 1B, C ). Based on the results shown in Figure 1 , 2005 was chosen as the starting year of the analysis herein.

To compute the total heat flux, we used the Objectively Analyzed Air–Sea Fluxes (OAFlux; Yu et al., 2008) monthly latent and sensible heat dataset combined with surface radiative flux data from the Clouds and the Earth’s Radiant Energy System (CERES) project (Kato et al., 2018). We utilized the Energy Balanced and Filled (EBAF)-Surface product as an appropriate means by which the net radiative heat flux emitted from the ocean may be determined. We examined data obtained during the period from 2005 to 2020. All previously mentioned datasets were transformed onto a 0.5° grid using the Python package xESMF (Zhuang, 2018). For the OAFlux and CERES datasets, bilinear interpolation was applied, whereas for the ISAS20 gridded fields, a conservative interpolation method was employed. We utilized the Gibbs-SeaWater (GSW) Oceanographic Toolbox (McDougall and Barker, 2011) to compute density. All computations and data processing were performed using the xarray package (Hoyer and Hamman, 2017; Hoyer et al., 2023).

To understand the primary reason for the heat flux anomaly, we analyzed each component of the net heat flux. Net heat flux (Qnet) was calculated as follows:

where SW is shortwave radiation, LW is longwave radiation, LH is latent heat, and SH is sensible heat. A positive value of the net heat flux and its components indicates that heat is released from the ocean to the atmosphere (positive upward). Furthermore, we decomposed the latent heat (LH) anomaly to examine its controlling factors, such as surface wind speed and humidity, and LH anomalies were decomposed based on the following equation (following Tanimoto et al., 2003 and Takahashi et al., 2021):

where ρ_a is atmospheric density, L is latent heat of vaporization, C_e is the bulk coefficient, U_a is wind speed at 10 m above the sea surface, q_s is specific humidity at the sea surface, and q_a is specific humidity 2 m above the sea surface. The overbar indicates the climatological mean for each month, whereas the prime denotes the anomaly from the mean. The first term on the right-hand side (RHS) of the equation represents the contribution of differences in specific humidity at the sea surface and humidity at a height of 2 m (q′ = q s ' − q a ' ), while the second term represents the contribution of anomalous wind speed at 10 m height (( U a ' ). The third term describes the combined effect of changes in wind speed anomaly and humidity difference anomaly ( U a ' q′), and the fourth term is the monthly climatology of the third term. The last two terms have very small values and can therefore be disregarded.

A mode water is distinguished by a layer with low stratification, and there are two commonly used approaches to assess its vertical stratification: vertical temperature gradient (dT/dz) and potential vorticity (PV). Potential vorticity is defined as follows:

where f is the Coriolis parameter, ρ is in situ density, and σ is potential density. Relative vorticity was deemed negligible since it is much smaller than planetary vorticity. We calculated the vertical potential density gradient (dσ/dz) by taking the potential density difference between adjacent grid points above and below (with a depth difference of 10 m). We used the same approach to compute the vertical temperature gradient (dT/dz), where T represents potential temperature.

To examine the upper layer stratification in summer, we used the Brunt-Vaisala frequency squared (N ²) as a proxy of stratification. N ² is calculated as follows:

where N is the Brunt-Vaisala frequency, g is the acceleration due to gravity, ρ is in situ density, and σ is potential density.

We defined the MLD as the depth at which the density increases by 0.125 kg m⁻³ and the temperature changes by 0.5°C relative to the surface, with temperature and density values at 10 m considered as the surface values. The MLD was determined as the shallowest depth calculated between the two definitions at each grid point. Using only temperature as a criterion to determine MLD can lead to overestimation in areas of strong salinity stratification. Similarly, relying solely on density criteria can result in overestimation in the region of compensated layers, where the density is vertically uniform, while the temperature and salinity change rapidly with depth (de Boyer Montégut et al., 2004; Oka et al., 2007). To avoid such potential overestimation, we utilized a combination of temperature and density criteria, taking advantage of the temperature and salinity data provided by the Argo float system.

The MLD distribution in September is shown in Figure 2A . September was selected because it represents the late winter period when the MLD is deepest. The dense mixed layer temperature contours around 40°S indicate the presence of the Agulhas Return Current. As shown in Figure 2A , thick mixed layers are observed on the equatorial flank of the Agulhas Return Current. The red dashed box in Figure 2A represents the formation region of IOSTMW. We intentionally excluded the coast of South Africa and the Agulhas Current region because we found that winter heat loss and variability in summer stratification in these areas were not clearly associated with the formation of IOSTMW formation. The characteristics of the winter mixed layer play a crucial role in the study of IOSTMW. During the spring, the seasonal thermocline forms and caps these deep winter mixed layers, resulting in the formation of IOSTMW. The properties of the winter mixed layer are important because they are closely tied to those of IOSTMW. This newly formed IOSTMW is then transported away from its formation zones, leaving behind a distinctive signature of winter mixing across vast regions of the Indian Ocean subtropical gyre.

Figures 2B, C shows the properties of the deep winter mixed layer north of the Agulhas Return Current. According to Tsubouchi et al. (2010), IOSTMW is estimated to have a temperature range of 16°C–17°C and a density range of 25.9–26.1 σ_θ , while Feucher et al. (2019) found that the average temperature of IOSTMW lies between 16.5°C–17.8°C with an average density of 25.9–26.0 σ_θ . Although the results of these studies may vary slightly, they are nonetheless within the range of mixed layer properties presented in Figure 2A , suggesting that the IOSTMW is linked to the late winter mixed layer. In this study, we excluded the deep mixed layers located east of 55°E with a temperature below 15°C and a density greater than 26.4 σ_θ. These properties of this water mass bear greater similarity to the lighter variants of Subantarctic Mode Water (SAMW), as described in prior studies (Fine, 1993; Wong, 2005; Koch-Larrouy et al., 2010).

Figure 5 displays scatter plots of IOSTMW thickness in September plotted against surface heat flux during winter (June–August) ( Figure 5A ) and summer (January–March) N ² values from the surface to 150 m ( Figure 5B ). Figure 5A shows that the thickness of IOSTMW is lowest in 2018, 2012, and 2015; this can be attributed to weak heat loss during those years. However, summer stratification varied significantly over these years. 2012 exhibited a high stratification, which is expected to inhibit the formation of a thick IOSTMW layer. In contrast, 2015 exhibited a relatively low stratification, whereas 2018 had the lowest summer stratification among all the data evaluated in Figure 5B . 2015 and 2018 share a notable common characteristic: they have the weakest winter heat loss among all the years examined. Hence, this finding suggests that even if summer stratification is weak, a thick IOSTMW cannot form in late winter if the winter heat loss is extremely weak. On the other hand, the year 2011 featured strong heat loss, but its summer stratification was well-stratified, which could have hindered the formation of a thick mode water layer. The years 2009, 2010, and 2019 are examples of thick years for IOSTMW, which provide insight into how summer stratification and winter heat loss combine to influence the formation of IOSTMW. In 2019 and 2009, moderate heat loss with weak summer stratification resulted in the formation of thick mode waters. The year 2010 exhibited both intense heat loss and weak summer stratification, naturally leading to the formation of thick mode water during late winter. These examples illustrate how winter heat loss and summer stratification work together to form IOSTMW each year. However, the cases of years 2015 and 2018 highlight that winter heat loss is the primary factor determining the IOSTMW thickness. Although a poorly stratified upper layer is favorable for the formation of a thick IOSTMW, extremely weak heat loss makes it difficult for a thick IOSTMW layer to form.

In the previous section, we showed that winter heat loss is the primary factor responsible for determining the annual thickness of IOSTMW. In this section, we explore in more detail how winter heat loss contributes to the formation of IOSTMW. To begin, we break down the net heat flux into its components, as illustrated in Figure 6A . It is noticeable that in some years, sensible heat (represented by the dashed light blue line) and latent heat (dashed dark-blue line) work together to produce an anomalous heat flux, such as observed in 2008, 2010, 2011, and 2018. However, we found that, in 2012 and 2015, latent heat alone was responsible for anomalous weak winter heat loss. We therefore further investigated the variability in latent heat.

Figure 6B presents the time series of the first three LH components on the RHS of Equation 2, highlighting that the three thinnest years of IOSTMW are caused by different factors. Specifically, the thin IOSTMW in 2012 is the result of an anomalously small humidity difference at the sea surface and at a height of 2 m, whereas the thin IOSTMWs in 2015 and 2018 are due to anomalously weak wind speed over the formation region. Conversely, the thick IOSTMW in 2010 was the consequence of a large difference in humidity leading to a strong winter heat loss. However, the wind was anomalously weak in 2010, resulting in decreased total latent heat. These findings demonstrate that thin IOSTMW years are associated with anomalously weak latent heat flux.

Furthermore, to elucidate the formation of the thick and thin years, we selected years wherein the mean thickness value exceeds 1 standard deviation (σ = 10.35 m; standard deviation calculated from thickness anomaly timeseries). The thick years include 2009, 2010, and 2019, while thin years include 2012, 2015, and 2018. Figure 7 displays the average winter net heat flux and summer stratification for both the thick and thin years. Figure 7B reveals a negative anomaly in the net heat flux, which indicates a weak heat loss during the thin years. In contrast, the formation region of IOSTMW, particularly east of 40°E, experiences intense heat loss during thick years ( Figure 7A ). However, this strong heat loss distribution is not as prominent as the weak heat loss distribution observed during thin years because the net heat loss during 2009 and 2019 is not particularly strong, as seen in Figures 4 , 6 . During these years, the formation of thick IOSTMW layers is aided by the weakly stratified upper layer during the summer season, as shown in Figure 7C . In thin years, the summer average N ² shows a positive anomaly value inside the formation region ( Figure 7D ). However, it is important to note that this is mainly due to the averaging process. The negative anomaly value during the thick years ( Figure 7C ) is due to the relatively weak summer stratification for 2009, 2010, and 2019 ( Figure 5 ). On the other hand, during the thin years, there appears to be an inconsistency between the weak stratification in 2018 and the very well-stratified summer of 2012 ( Figure 5 ). The presence of thin IOSTMW layers during thin years is attributed to the extremely weak winter heat loss, as indicated in section 4.1.

To better understand the impact of anomalous latent heat on the heat flux anomaly, Figure 8 displays the spatial distribution of U a ' and ( q s ' − q a ' ) calculated using Equation 2. During thin years, the wind speed distribution shows a weak anomaly ( Figure 8B ) coupled with an anomalously small difference in humidity between the sea surface and the air ( Figure 8D ). This combination provides the ideal conditions for creating an anomalously thin IOSTMW layer. Figure 6 demonstrates that the latent heat anomaly contributions from the wind speed ( U a ' ) and humidity difference ( q s ' − q a ' ) have negative values for 2012, 2015, and 2018, i.e., weakened heat loss from the ocean to the atmosphere. On the other hand, the opposing effects of the wind and humidity contributions during thick years can be observed ( Figure 8A, C ). This can be explained by examining the latent heat decomposition in Figure 6B . In particular, the wind contribution in 2019 has a positive effect on latent heat, enhancing heat loss to the atmosphere, but this is counteracted by the negative wind contribution in 2009 and 2010 (green dashed line in Figure 6B ). Therefore, the distribution of the wind speed anomaly is weakly negative ( Figure 8A ). In Figure 8C , a positive anomaly value in the humidity difference suggests an intense heat loss to the atmosphere. The ( q s ' − q a ' ) anomalies in 2009 and 2019 are small and close to zero, as shown by the red dashed line in Figure 6B . Therefore, the remaining positive contribution observed in 2010 dominates when averaging the three thick years, showing positive anomalies around the formation region.

The lower panels of Figure 8 display a similar pattern to that shown in the upper panels of Figure 7 , indicating the significant role played by ( q s ' − q a ' ) in the total winter heat flux. During thin years, the negative contribution from ( q s ' − q a ' ) is further augmented by the negative contribution from U a ' , resulting in a weak net heat loss anomaly. Figure 6B shows that the latent heat anomaly is the primary cause of the weak heat loss, with the exception of 2018, when sensible heat also exhibited a negative contribution. These findings demonstrate how ( q s ' − q a ' ) and U a ' combine to engender a weak heat loss anomaly, which hinders the formation of thick IOSTMW. Conversely, during thick years, ( q s ' − q a ' ) displays positive anomalous values ( Figure 8C ) but is slightly counteracted by the contribution of U a ' ( Figure 8A ). However, the net heat flux values are similar or greater at some locations within the formation region ( Figure 7A ), implying that other contributions from net heat flux components are also present. During the thick years, contributions from U a ' and ( q s ' − q a ' ) act in opposition (e.g., 2010 and 2019 in Figure 6B ); as such, the total latent heat anomaly is not strongly positive.