Long-term trend and variability of ocean heat content in the Indonesian Maritime Continent

The evolution of 64 years of the ocean heat content for the upper 300 meters (OHC300) in the Indonesian Maritime Continent, viz the Banda Sea, the Sulu Sea, the Java Sea, and the Karimata Strait, was investigated using the combination of ocean and atmosphere reanalysis data. The study aims to understand its long-term trend variability and association with El Niño-Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD). The result showed a positive trend in all basins except the Banda Sea. The most significant rise was the Sulu Sea. The IOD fluctuation was more pronounced in the Western part of Indonesia. In contrast, the ENSO variation was more apparent in the Eastern part of Indonesia. The coincided La Nina events and the negative IOD significantly increased OHC—conversely, the coinciding of the El Niño events and the positive IOD caused a decrease in OHC magnitude. The significant factor increasing OHC in the maritime continent is the ENSO fluctuation, followed by the IOD fluctuation.


Introduction
The Earth's climate system has an energy imbalance of nearly 1 W m -2 (Trenberth et al., 2014).The Earth's energy imbalance has been mainly attributed to increased greenhouse gases (GHGs) and aerosols driven forward by anthropogenic activity during the industrial period.The GHGs and aerosols trap more radiation and hence, contribute to global warming (IPCC, 2014).The world's oceans absorb more than 90% of this extra heat, which causes an increase in the ocean heat content (OHC) and sea level rise (Trenberth & Fasullo, 2012).Because of the ocean's prominent role in global heat storage changes, the OHC dynamics significantly impact Earth's energy imbalance on interannual and longer periods.
One of the climate indicators that correctly indicates that the world is warming up is the OHC in the upper 300 m, hereafter OHC300 (Xue et al., 2012).The OHC300 averaged over the regional seas exhibits substantial decadal variability rather than consistent warming.In addition, trends in OHC300 over the past 50 years are spatially quite varied, showing significant effects of oceanic heat redistribution caused by changes in ocean circulation (Doney et al., 2007).The uneven spatial variations in the OHC300 in the tropical oceans will likely impact natural modes of variability and regional climate conditions, particularly in the Indonesian Maritime Continent (IMC).The IMC is a large continent in Southeast Asia that consists of the Indonesian Archipelago, straits, and marginal seas as Figure 1.Observations reveal that seasonal warming variability in the IMC controlled by the heat exchange at the air-sea interface (Ismail et al., 2023, Ismail & Karstensen., 2023).The previous study focused on the OHC derived from Argo data in west Sumatera compared to Indian Ocean Dipole (IOD) in the period 1999-2011 (Radjawane et al., 2015), OHC derived from temperature and salinity of HYCOM model in Lombok area and compared to El Niño-Southern Oscillation (ENSO) and IOD for the period 2011-2015 (Chairuasni et al., 2020).
Considering the role of the maritime continent in the global ocean circulation and limitation of study of OHC, in this study, we investigate the spatial and temporal characteristics of the decadal variability of the OHC300 and estimate the possible relationship and mutual correspondence of the oceanic and atmospheric components such as ENSO and IOD based on the OCEAN5 reanalysis outputs.This work attempts to explore the decadal variability of OHC300 in the IMC using a more extended temporal period from 1958 -2021.

Data
2.1.1.OHC300.We used OHC300 from Copernicus data, the upper Global Ocean Heat Content (0-300m) from Reanalysis & Multi-Observations Reprocessing with the most prolonged period, January 1958 -December 2021 (64 years).The source was https://data.marine.copernicus.eu/product/GLOBAL_OMI_OHC_area_averaged_anomalies_0_300/description.The data was used to see the trend, the variability for 64 years and association with the IOD and the ENSO index.We chose the Sulu Seas, the Banda Seas, Karimata Strait, and the Java Sea in order to represent the Indonesian maritime continent area.

ENSO.
The ENSO data originated from the National Oceanic and Atmospheric Administration (NOAA) (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php).It shows a threshold of 3 months running mean of ERSST.v5SST anomalies in the Niño 3.4 region (5 o N-5 o S, 120 o -170 o W) of +/-0.5 degrees centered 30-year base periods updated every five years.If the threshold is more than 0.5 degrees, it means warm on the west coast of South America.While the threshold is below 0.5 degrees, it means cold on the west coast of South America (Huang et al., 2017) To understand more detail, we included one example of El Niño and La Nina event in the year 1998 in Figure 1 (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensocycle/ensocycle.shtml).During El Niño, the west coast of South America experienced an increased sea surface temperature (about 1-2 degrees above average), which is indicated by red (figure 1).On the contrary, during La Nina, the west coast of South America experienced a decrease in sea surface temperature (about 1-2 degrees below average), indicated by the blue color.
How about the impact on the Indonesia area?During El Niño, Indonesia experiences a decreased temperature (cold) and is dry.On the contrary, during La Nina, Indonesia experiences an increase in 3 SST (warm) and tends to be wet.Regarding this impact, we signed red for warm SST (during La Nina event) and green for cold SST (during the El Niño event) in Figure 3e.

IOD. Intensity of the IOD is represented by anomalous SST gradient between the western equatorial Indian Ocean (50E-70E and 10S-10N) and the south eastern equatorial Indian Ocean (90E-110E and 10S-0N
).This gradient is named as Dipole Mode Index (DMI).When the DMI is positive then, the phenomenon is refereed as the positive IOD and when it is negative, it is refereed as negative IOD (https://psl.noaa.gov/gcos_wgsp/Timeseries/DMI/).To understand the effect of IOD on OHC300, we used IOD in the same period, January 1958 -December 2021, which originated from https://www.cpc.ncep.noaa.gov/products/international/ocean_monitoring/IODMI/DMI_month.html/.The data was calculated using ERSST.V5 (Huang et al., 2017).It was estimated from the subtraction between the western Indian Ocean (WTIO) and the south-eastern Indian Ocean (SETIO).WTIO is SSTA averaged in 50ºE-70ºE, 10ºS-10ºN.SETIO is SSTA averaged in 90ºE-110ºE, 10ºS-0°N.Then, DMI is WTIO subtracted by SETIO.Positive and negative periods are based on a threshold of +/-0.4°C for the DMI and centered 30-year base periods updated every ten years.Currently, the climatology used was from 1991-2020.A positive IOD means warm in the WTIO areas or cold in the SETIO.In other words, the western Indonesia area experiences a decrease in SST and tends to be dry.A negative IOD represents cold in the WTIO areas or warm in the SETIO areas.Western Indonesia area experiences an increase in SST and tends to be wet.Regarding the impact on the Indonesia areas, we signed cold (warm) in the SETIO by green (red) in Figure 3f.

Method
In this study, the OHC was computed using the following equation: ௭ Where ρ is the density of seawater (1026 kg m -3 ), C is the specific heat of sea water (3990 J kg -1 K -1 ), z is the depth limit of the calculation (300 meters), and T is the temperature at each depth ( o K).

Trend
We analyzed a trend using a function trend of the Climate data toolbox (Greene et al., 2019) which calculates the linear trend of a data series by least squares.The sampling rate is 12 months.Because we want to know the trend per decade, we multiply the sampling rate by 10.
In this study, we investigate not only in time series but also in spatial.So, a spatial trend analysis OHC was calculated for all spatial grid points inside the area of interest (Fig. 2) over the entire data period from January 1958 to December 2021.While a long-term trend, we average the OHC of all grids every year for each area, namely Banda Sea, Sulu Sea, Java Sea and Karimata Strait (Fig 3).Note that a seasonal cycle in the data has been removed.The unit is Joule m -2 decade -1 (Jm -2 decade -1 ).

Time series analysis and association with the IOD and the ENSO using correlation.
We calculate the OHC anomalies each year by subtracting the long-term mean of OHC from OHC's mean per year.We applied it to each area, such as Banda, Sulu, Java, and Karimata Strait (Figure 3).The mean of OHC per year is the average of OHC from all grids of each area per year.At the same time, the long-term mean is the average of OHC for all grids over time .
Furthermore, we analyze the OHC anomalies associated with the IOD and the ENSO.Regarding IOD and the ENSO, we calculated the association using the Pearson correlation coefficient (https://nl.mathworks.com/help/matlab/ref/corrcoef.html).The correlation coefficient of two random variables measures their linear dependence.The value is between -1 to 1.A value of -1 means that the two variables have a perfect negative linear correlation, with 0 being no correlation and + 1 indicating a perfect positive correlation (Williams et al., 2020).If each variable has N scalar observations, then the Pearson correlation coefficient is defined as Where μ A and σ A are the mean and standard deviation of A, respectively, and μ B and σ B are the mean and standard deviation of B. In our case, A can be the OHC anomalies, and B can be the IOD or the ENSO index.
Furthermore, we analyze the OHC anomalies associated with the coincided IOD and ENSO (El Niño and La Niña event).

Hovmoller diagram
To understand the evolution of OHC spatially, we used the Hovmoller diagram, such as timelongitude and time-latitude.Time-longitude enables us to see the development of the OHC pattern from one longitude to another over time.In comparison, time-latitude allows us to see the evolution of the OHC pattern from one latitude to another by time.

Trend analysis 3.1.1. Analysis of spatial trend
Figure 2 shows a spatial trend of OHC from 1958-2021 with the removed seasonal cycle.Four boxes represented the Banda Sea (red), the Sulu Sea (yellow), the Java Sea (green), and the Karimata Strait (white), respectively.The multiplier of the contour values is 10 7 Jm -2 decade -1 .The white contour means a spatial trend from 1958-2021 below -1.The dark blue contour contains a value between -1 and 1.The light blue contour represents values above 1 of the trends.The trend in the Banda Sea demonstrates a decreasing value.The Sulu Sea tends up to 1.The Karimata Strait tends to increase, as well as the Java Sea.The Karimata Strait is dominated by trend values ranging from -1 to 1.The Java Sea varies between -1 and more than 1.The Sulu Sea is huge variability about below -1 until above 1.The Banda is dominated by less than -1.The magnitude range of each area is quite similar to the long-term trend in the time series (Figure 3).Karimata Strait (white), Java Sea (green), Sulu Sea (yellow), and Banda Sea (red).White shape contour (< -1), dark blue (between -1 and 1), light blue (> 1).

Analysis of long-term trend from time series
Besides, Figure 3 also shows a trend value for each area displayed in the legend (black text).A positive trend of OHC300 covers all basins except the Banda Sea.The positive trend of the Sulu Sea, the Java Sea, and the Karimata Strait is 1.4693 x 10 7 Jm -2 decade -1 , 1.2877 x 10 7 Jm -2 decade -1 , and 0.21193 x 10 7 Jm -2 decade -1 , respectively.At the same time, the Banda Sea has a negative value (-1.6336 x 10 7 Jm - 2 decade -1 ), meaning that it experiences a decrease.The most significant rise is the Sulu Sea.Generally, the pattern of OHC300 of the Banda Sea and Sulu Sea are similar-likewise, Karimata Strait and Java Sea.Due to the adjacent locations, the Java Sea and the Karimata Strait are identical from 1997 to 2021.

Analysis of time series
Figure 3 shows the time series of the OHC300 from 1958-2021 for four areas: Banda Sea (a) and Sulu Sea (b), Java Sea (c), and Karimata Strait (d).Every time series has varied negative values (green) and positive values (red) with unit Jm -2 .However, a multiplier value is different depending on the calculation result and to plot the graph and analyses easier.A negative OHC and positive cold OHC mean cold and warm OHC, respectively.Figure 3a and Figure 3b show a similarity between warm and cold periods, namely warm years (1962-1980 and 2008-2014) and cold years .Whereas cold years in the Java Sea and Karimata strait (Figure 3c and 3d) took longer .The warm years are 1958-1974 and 1998-2021.

Correlation analysis
We calculated a correlation coefficient from the time series of OHC anomalies in Figure 3 and displayed it in Table 1.The table shows a correlation between the OHC300 anomalies and the ENSO index, the OHC300 anomalies and the IOD.The OHC and the ENSO index in east Indonesia, such as the Banda Sea and Sulu Sea, are more correlated than in west Indonesia, such as the Java Sea and Karimata Strait.On the opposite, the correlation between the OHC and IOD in East Indonesia is less correlated than the correlation in West Indonesia.That was reasonable because the area position, the more distance, the lesser effect.The most significant correlation with the ENSO index is the Sulu Sea, followed by the Banda Sea, Karimata Strait, and Java Sea.In comparison, the most considerable correlation with IOD is the Java Sea, followed by Karimata Strait, Sulu Sea, and Banda Sea, respectively.

Analysis OHC versus ENSO
The Banda is theoretically influenced by ENSO fluctuation because the area is in eastern Indonesia, which is adjacent directly to the Pacific Ocean (Gordon & Susanto, 2001).The study shows that SST in the Banda Sea correlates well with the El Niño and La Niña events.The evidence of interannual SST variation has been observed in some areas of the Indonesian Seas such as in the Banda Sea (Gordon and Susanto, 2001).It is found that lower/higher SST is identified during the El-Nino/La-Nina periods.Interannual SST variability of the Indonesian waters are not only influenced by interannual variability of wind forcing but also interannual variability of thermocline temperatures (Gordon and Susanto, 2001).They suggested that the shallowing/deepening thermocline advected from the Pacific Ocean into the Indonesia Seas is the main factor in inducing lower/higher SST during the El-Nino/La-Nino.The shallowing/deepening thermocline is related to the upwelling/downwelling wave propagation from the Pacific Ocean scattering into the Indonesian Archipelago (Iskandar, 2010).Because the OHC is derived from sea temperature, the association between the OHC and the ENSO index exists.It was also valid for the Sulu Sea.Fluctuation of OHC in the Banda and Sulu Seas tend to be interannual cycle.The OHCs in the Banda and Sulu Seas are colder than normal during the El Niño event.Meanwhile, during the La Niña event, We applied a time-latitude analysis of the OHC anomaly to investigate the impact of a coincidence ENSO with the IOD fluctuation and OHC. Figure 4 shows interannual variability along latitude for the period 2002-2021.We only display the last 20 years to analyse quickly, meaning that spatial variability in a short period can be seen clearly.Previously, we hypothesized that the possible warmest OHC happened when the negative ENSO (La Niña) coincided with the negative IOD event.From Figure 3, the coincided event between the La Niña and the Negative IOD is 2011, 2013, 2017, and 2021.Compared to Figure 4, the warmest is the year 2011, during the La Niña event coincided with the Negative IOD, although it did not cover the whole months of 2011.It occurred in all areas, with eastern Indonesia being more influenced than western Indonesia.Although the magnitude of negative IOD is quite large, if the magnitude of the ENSO negative is small, the impact only happened in the Java Sea and Karimata Strait.In addition, during the negative IOD (2011, 2013, 2017, and 2021), as in Figure 4, an increasing OHC lies on the north Karimata Strait.We suggested it is the impact of La Niña.
In contrast, we argued that this coldest OHC occurred due to the influence of the period of positive ENSO (El Niño) and positive IOD events.Figure 3 shows the coincided event in 2004, 2007, 2015, and 2019.The coldest was in 2015 during the El Niño coincided with the Positive IOD.Same as above, the Karimata Strait does not experience much because the magnitude of the Positive IOD is relatively tiny.The other coincided event is 2019; the magnitude of positive IOD is quite significant, and the positive ENSO value is relatively small, but the OHC in Karimata Strait is still warm.The situation's impact is only in Banda, Sulu, and Java.The positive IOD value is relatively high.However, it did not have a strong influence.We conclude that IOD significantly affects the Java Sea and Karimata Strait.The ENSO influenced the Banda Sea and Karimata Strait.However, when they coincide, the ENSO event is more controlling than the IOD.

Conclusion
Based on analysis above, we conclude that the trend of both spatial and time series shows a similar pattern which is a positive trend in all basins except the Banda Sea.The most significant rise was the Sulu Sea.The IOD fluctuation more influenced west Indonesia.In contrast, the ENSO variation more influenced east Indonesia.The coincided La Niña event and Negative IOD significantly increased OHC-conversely, the coinciding of the El Niño event and the Positive IOD caused a decrease in OHC magnitude.The significant factor increasing OHC in the maritime continent is the ENSO fluctuation, followed by the IOD fluctuation.

Figure 3 . 1 .
Figure 3.Comparison of time series of OHC300 for each area against the ENSO and IOD and the trend.

Figure 4 .
Figure 4. Anomaly of OHC by a longitude averaged.OHC anomaly based on period from 1958 to 2021 (64 years).

Table 1 .
Correlation between OHC300 per area against the ENSO and the IOD