Network science approaches for analyzing regional interactions of sea surface temperature

2.1. Sea Surface Temperature Dataset and Preprocessing

We utilize the optimum interpolation Sea Surface Temperature (SST) dataset version 5 from National Oceanic and Atmospheric Administration (NOAA, USA) (NOAA). This dataset is a high-resolution (2° × 2°) global analysis of SST, produced by blending satellite observations, ship, and buoy data. The SST dataset provides daily SST values, which are optimally interpolated to fill data gaps, and this preprocess makes the dataset suitable for climate monitoring and long-term trend analysis. Version 5 of the dataset includes updated algorithms for better bias correction and improved accuracy in regions with sparse in-situ observations. This dataset is widely used in oceanographic and climate studies due to its consistency and high temporal resolution, making it an ideal choice for analyzing regional SST interactions. In this study, data from 1979 to 2019 were utilized for the analysis, and to analyze the temperature dataset, we calculated the average temperature for each cell over the target period, then subtracted this average from all temperature data points within that cell.

2.2. Heat Equation Approach

In statistical physics, the temporal evolution of a spatially distributed temperature field $$ T\left(\overrightarrow{r},t\right)$$ is often described by the heat equation, where $$ \overrightarrow{r}$$ denotes the position vector, and t represents time. The mathematical form of the heat equation is as follows:

Here, ∇² is the Laplacian, the second-order spatial derivatives of temperature, and D represents thermal diffusivity, indicating the sensitivity of a point to temperature differences in its vicinity. The Laplacian yields negative values when the temperature is higher than the average temperature of the surroundings and positive values when the temperature is lower.

In discrete spaces, Eq (1) can be rewritten as

where index i denotes spatial positions, and τ is time interval of the observed temperature. The M matrix represents the diffusivity of each cell, and the elements of the Laplacian matrix L, calculated as K − A, where K is a diagonal matrix containing the number of neighboring locations for each cell, and A is the adjacency matrix, thus analogously representing the Laplacian operator in continuous space.

The heat equation is typically considered in uniform spaces where interactions among neighboring nodes are equal, similar to an unweighted network. However, in the real world, non-uniform interactions between regions are expected in phenomena such as ocean currents and atmospheric teleconnections. This leads us to consider a weighted and directed network scenario, where the off-diagonal elements of M_iL_ij can take any real value. Here, the real-valued elements represent the weighting of how much a position of a node j influences node i, based on temperature. The previously discussed M_iL_ij is extended to a generalized Laplacian matrix, W_ij, which represents the relative influence of j on i. These weights need not be symmetric. In this case, the diagonal elements W_ii become −$$ -\sum _j{W}_{ij}$$. For example, in a network with uniform temperature T_i = T_j =T, there would be no temporal variation in temperature at any node, leading to $$ \sum _j W_{ij}T=T\sum _j W_{ij}=0$$.

To determine W_ij which best explains the given spatio-temporal SST data T under the aforementioned conditions, the equation suggests that a suitable W_ij should satisfy the following criteria:

and Eq. 3 can be applied to an objective function of $$ E\left(W_{ij}\right)={\left(T_i\left(t+\tau \right)-T_i\left(t\right)-\sum _j W_{ij}{T}_j\left(t\right)\right)}^2$$. By applying the Adam optimization to minimize the error, one can begin with the Laplacian matrix of the fully connected network as the initial W_ij. Since the update rule for W_ij must also account for changes in W_ii due to constraints that preserve the sum of W_ij in W_ii, the equation can be rewritten as follows:

where $$ \frac{\partial W_{ii}}{\partial W_{ij}}=-1$$. With a learning rate β, one can update W_ij as

Here, s denotes the updating step, and δs represents the interval between the updating time steps. Please note that, during the simulation, any change that would result in a positive W_ij will be used. It means off-diagonal components of W_ij were kept positive by setting any negative values to zero.

2.3. Construction of Backbone Network

To extract the backbone of the SST network, we applied the method proposed by Serrano et al. (Serrano et al., 2009). This technique identifies the most significant connections in a weighted network by statistically filtering out less relevant links. Specifically, for each node, the algorithm evaluates the weights of its edges and compares them to a null model of random edge assignments. This method allows us to retain only those edges that exhibit statistically significant interactions, effectively reducing the complexity of the network while preserving its essential structure. In our analysis, we set the significance level to β = 0.1, meaning that only links satisfying the significance level were retained.

Starting from the weighted matrix W, derived from applying the heat equation approach to the SST data, we extracted the SST backbone network. This filtered network highlights the significant interactions between sea regions, offering a clearer understanding of the underlying patterns in SST variations. After having the backbone network with β = 0.1, the total number of nodes in original backbone network is 9,958, and with 1,555,158 directed edges.

2.4. Louvain Community Detection

The Louvain community detection algorithm is a widely used method for identifying communities in large networks by optimizing modularity. Modularity is a measure that compares the density of edges within communities to the density of edges between communities (Blondel et al., 2008; Dugué et al., 2015). The modularity Q for a directed weighted network is mathematically expressed as:

where A_ij is the adjacency matrix of the network, $$ k_i^{\mathrm{i}\mathrm{n}}$$ in and$$ k_j^{\mathrm{o}\mathrm{u}\mathrm{t}}$$ out are the incoming degrees of nodes i and outgoing degrees of node j, m is the total number of edges. In this equation, c_i and c_j represent the communities to which nodes i and j belong, and δ(c_i, cj) is 1 if nodes i and j are in the same community and 0 otherwise. The Louvain algorithm works in two phases. In the first phase, each node is assigned to its own community, and nodes are iteratively moved to neighboring communities to locally maximize the modularity. In the second phase, nodes belonging to the same community are aggregated into a single super-node, and the algorithm is reapplied to this reduced network. The process continues until no further improvements in modularity can be made, allowing the detection of hierarchical community structures at multiple scales. In the current study, we employed the Louvain algorithm from Python’s NetworkX package, using the default resolution parameter of 1, and a convergence threshold of 10^-7.

2.5. Building the Community Network

Each detected community with the Louvain algorithm represents a group of nodes, that have stronger interactions with one another compared to nodes outside the community. In this study, a group of nodes in the backbone network represents the positions of sea surface areas divided into a grid.

To simplify the interaction structures between sea regions, we constructed a new network from the backbone with community indexes, where each node in this new network represents an entire community identified in the backbone network. The connectivity between these new community nodes was established by aggregating the original links between nodes in the backbone that belonged to each pair of communities. In this aggregation, we only considered strong link weights greater than 10^-3. This process resulted in a community-level network with 12 nodes, where the edges represent the total weight of the connections between communities. This community-level network provides insights into the interactions among different sea region communities, enabling us to observe higher-order patterns and relationships between key regions within the backbone structure.

2.6. Network Centralities

3.1. Regions with strong inward and outward weights from heat equation approach

Fig. 1 demonstrates significant spatial variation in the adjusted inward and outward flow weights derived from the heat equation approach. Regions A (the North Pacific) and B (the Indian Ocean) in panel (a) are characterized by high inward flow, suggesting a possible strong heat concentration of influences within these areas. Conversely, regions A, C (ENSO), D (NAWH), and E (the Southern Atlantic Ocean) in panel (b) exhibit high outward flow, indicating that these areas might act as major sources of heat transfer, distributing their heat flow across the network. The black arrows in both panels illustrate the sum of directional flows for each point, providing insights into the overall movement dynamics. Specifically, the directional vector at each location i, denoted as $$ \overrightarrow{V_i}$$, is determined by its relationship with surrounding nodes positioned in the eight directions around i—namely, below, left, right, and the four diagonal directions. For each pair (i,j), the location vectors for i, and j are represented as $$ \overrightarrow{r_i}$$, and $$ \overrightarrow{r_j}$$, respectively. From the location vectors, the unit vector representing the direction from i to j is denoted as $$ {\overrightarrow{v}}_{ij}=\left(\overrightarrow{r_j}-\overrightarrow{r_i}\right)$$. By using W_ij as the weight associated with the direction from j to i, we can define the weighted sum at location i as follows:

Fig. 1.

Adjusted inward and outward weights derived from the heat equation approach. The weights Wij from the heat equation method have been adjusted to range from 0 (blue) to 2 (red). Panel (a) illustrates the adjusted inward weights, while panel (b) shows the adjusted outward weights. The black arrows represent the sum of flow directions at each point’s surrounding positions. High inward flow regions are indicated as regions A and B in panel (a), while high outward flow regions are labeled as regions A, C, D, and E in panel (b).

These directional relationships highlight key regions that play crucial roles in the inward and outward heat transfer, which may have implications for understanding underlying processes such as oceanic circulation or climate phenomena.

For example, regions B, D, and E, located in the Indian Ocean, the North, and South Atlantic Ocean, respectively, are known for their interrelated heat transfer interactions (Hu and Fedorov, 2020; Yang et al., 2022). These regions are highlighted in Fig. 1a and Fig. 1b, showing strong inward flow for region B and a strong outward flow for region E. Additionally, the North Atlantic Warming Hole (NAWH) is marked as region D in Fig. 1b, showing strong outward flows. Given that the causes of the NAWH are still under debate, this result offers a new perspective on the oceanic heat transfer interactions related to this phenomenon and its connection to other sea regions.

Moreover, the primary region of El Niño-Southern Oscillation (ENSO) in the Pacific Ocean is highlighted as region C, indicating a strong outward flow in Fig. 1b. This result demonstrates that the weights derived from the heat equation can capture directional heat transfer relationships, which have been limited in conventional analyses.

3.2. Structural and Spatial Characteristics of Backbone Network

Since the heat equation approach provides directional weights of heat transfer flows, one can analyze directional relationships by transforming them into a heat transfer weight network. While the structural analysis of directional and weighted network can reveal crucial insights into the connectivity patterns and spatial relationships within the oceanic heat transfer network, its dense connectivity can obscure the identification of primary structures. Therefore, extracting the backbone network is essential to focus on the most significant connections and better understand the key components of the system. To this end, we applied a backboning algorithm (see Methods, (Serrano et al., 2009)) to extract statistically significant weighted links.

The in-degree and out-degree distributions of the original network show a characteristic pattern that can be well-described by a log-normal distribution (Fig. 2a). This type of distribution may indicate highly heterogeneous connectivity, with a few nodes having substantially more links than others, depending on the mean and standard deviation. To assess the heterogeneity of in-degree and out-degree distributions, we performed power-law fitting on their tails (Fig. 2b). As observed, the in-degree distribution exhibits an exponent close to 2, demonstrating a heavy-tail tendency. In contrast, the out-degree tail shows an exponent close to 9, which is relatively large and less characteristic of a heavy-tail property. A distribution like in-degree distribution is typical in natural systems, where certain regions play a dominant role in either collecting or distributing influence, further emphasizing the complex dynamics of heat transfer within the ocean.

Fig. 2.

Structural characteristics of the backbone network derived from directional weights. (a) The in-degree and out-degree distributions of the original network derived from the directional and weighted interactions shown in Fig. 1. Circle markers indicate the empirically observed frequency of the number of outgoing links, while triangle markers represent the number of incoming links. The grey and red solid lines show the log-normal fits for the outgoing and incoming links, respectively. (b) Degree distributions of the original correlation networks (on log-log axes), with tail-part fitting applied to the distributions. (c)The distribution of the number of in-degree and out-degree for the backbone network (α = 0.1, see Methods), with markers indicating the relative frequency of each type of degrees, are shown. The fitted tail distributions are shown with dashed lines in red (for in-degree) and grey (for out-degree). (d) Geodesic distance distribution between nodes, calculated based on the center of each cell within the backbone network, highlighting the spatial relationships between nodes. The distances range from 10 km to 20,000 km. The vertical dotted line indicates the starting point of the second peak in long-range distances between nodes (≥3,500 km).

On the other hand, the in-degree and out-degree distributions of the backbone network exhibit power-law behavior, with exponent of 5.74 and 4.37, respectively, as shown in Fig. 2c. Due to the large exponent values, it is challenging to conclude that the distribution exhibits a heavy-tail property. However, this distribution still demonstrates that, after extracting statistically significant links, a few nodes retain a large number of connections (more than links), while most nodes have fewer than 200 links.

In terms of spatial characteristics of the backbone network, Fig. 2d shows the geodesic distance distribution between nodes, revealing that connections between nodes span a wide range of distances, from localized interactions (10 km) to broader, ocean-basin scale connections (up to 20,000 km). This diversity in connection lengths suggests a multiscale structure of heat transfer in sea regions, where both local and global interactions are essential for maintaining the heat transfer backbone network's overall connectivity.

3.3. Communities in Backbone Network

To identify key regions and their connectivity patterns in the backbone network, we performed community detection using the Louvain algorithm (see Methods). In particular, we used links spanning more than 3,500 km to focus on the long-range interactions between sea regions, as 3,500 km represents a starting point (a mid point between 3,000 km and 4,000 km) of the second peak in the distribution of long-range interactions (Fig. 2d). The community detection results in Fig. 3a reveal a clear structural organization within the backbone network, highlighting the distinct communities formed as a result of long-range directional heat transfer interactions.This structure is maintained when the distance criterion is 3,000 km. Panel (a) shows the complete set of communities within the backbone network, with the previously marked regions (A, B, C, D, and E) overlaid on the community structure. The detected communities, each represented by a different color, exhibit unique spatial patterns, indicating the presence of highly interconnected regions across multiple scales of distance. This differentiation allows us to better understand how various parts of the ocean interact and identify the dominant areas of influence in oceanic heat transfer.

Fig. 3.

Community detection in the long-range backbone network and regional categorization. (a) Community structure of the entire long-range backbone network derived from the directional weights, with marks of previous regions labeled as A, B, C, D, and E. The detected communities are represented by different colors based on their community ID (from 0 to 11). (b-f) panels show detailed community distribution for each marked region: (b) Region A, (c) Region B, (d) Region C, (e) Region D, and (f) Region E. Each panel highlights the communities within the corresponding region, revealing the spatial distribution and boundaries of the identified communities.

In the detailed views provided in panels (b) to (f), each region is analyzed individually to highlight its community structure. Region A, as shown in panel (b), contains nine communities, including most of the communities except for those in parts of the Indian and Pacific area (community 9) and the North Pole (community 10), each with distinct spatial boundaries. On the other hand, regions B and C, as shown in Fig. 3c and 3d, include a relatively smaller number of communities compared to region A. Notably, region B, which is primarily located in the Indian Ocean, includes communities 6, 9, and 11, which are part of the ENSO system. Additionally, this region contains community 2, and 6, which is largely related to the North Atlantic Ocean. Region C also highlights the ENSO system, including its relationships with the Indian and Atlantic Oceans.

Region D, which mostly represents the NAWH, is characterized by more diffuse community boundaries involving communities 2, 6, and 7, possibly indicating a broader influence from all sea areas, including the near the South Pole area, and the Indian Ocean. In contrast, region E resembles region A, as it contains a large portion of the communities, including the North Pole area, suggesting that the Southern Atlantic Ocean has a highly mixed influence from other sea communities. Overall, the spatial categorization of these communities reveals the multiscale nature of heat transfer in the ocean, highlighting previously reported crucial connections among the Indian, Atlantic, and Pacific Oceans.

3.4. Centralities of Communities with Total Weights and Averaged Weights

To understand the directional and weighted flow between communities, we reconstructed the community network by aggregating the directional weighted links between them (see Methods). In this network, each community is represented as a node, and all connections within a community have been merged into a single directed weighted link to another community. When merging these links, we used both the average weights and the total weights as they represent different aspects of the relationships between communities (see Methods). On this networkwe applied centrality measures to both nodes and edges, using network science techniques to identify the functional role of each community in sea surface heat transfer.

3.4.1. Betweenness Centrality

Fig. 4 highlights the identification of key communities within the network of oceanic communities using betweenness centrality as a measure of importance (see Methods). Specifically, in panel (a), betweenness centrality was calculated using average weights, which is normalized by the number of links, allowing us to determine which communities play a central bridging role per unit link between different parts of the network. As shown in Fig. 4b, community 8, 11, and 9 play crucial roles as bridges connecting various other communities, with their representative region being the ENSO system. This result indicates the possibility that a significant portion of the ENSO system may act as channels for heat transfer between other communities.

Fig. 4.

Identification of key communities using betweenness centrality in the network of sea region communities. (a) A network between communities shows betweenness centrality measured for the communities (nodes), using average weights normalized by the number of links between them. Node colors represent the strength of each node (each community)'s betweenness centrality, with the top 10% of average weights (13 edges) highlighted in darker grey. (b) The top 3 communities with the highest betweenness centrality from (a) are shown. (c) Betweenness centrality on the same network as in (a) is calculated using total weights instead of average weights, with nodes colored according to their betweenness centrality. The top 10% of edges by total weight (13 edges) are also highlighted in darker grey. (d) The map highlights the top 3 regions with the highest betweenness centrality based on (c). The node size is determined based on the number of links, both incoming and outgoing.

In panel (c), betweenness centrality was recalculated using total weights rather than average weights, providing insight into the overall influence of each oceanic community by considering the strength of all connections. Fig. 4d maps the top three regions (communities 3, 4, and 11) with the highest betweenness centrality, which are primarily located in the South Polar areas and the periphery parts of ENSO system. This result emphasizes their crucial function in global oceanic connectivity, particularly when considering their overall influence. These findings highlight the specific roles of oceanic communities, namely the ENSO system and polar regions, in heat distribution throughout the ocean. Considering the structure of community networks in Figs. 4a, and 4c, the ENSO system-related communities (4, 9, and 11) have relatively weak but diverse inward connections from most other communities in both cases, with particularly strong connections between community 4 and communities 2, 5, 6, and 7 (closely related to the South Pole, Indian Ocean, and parts of the Pacific Ocean) in terms of total weights. Regarding the average weights, the ENSO-related communities receive inward flows from communities 2, 5, 6, and 8 (which are more diffuse), and direct outflows to other communities (0, 6, and 7), which are located near the periphery of the ENSO system, as well as the Indian and Atlantic Oceans. This provides insight into the role of the ENSO system, suggesting that its crucial impact on oceanic and atmospheric dynamics may arise from its role as a bridge connecting various oceanic regions.

3.4.2. Outbound Eigenvector Centrality

To understand the importance of sea communities in outgoing heat transfer, we analyzed the outbound eigenvector centrality. Fig. 5 highlights the key communities using outbound eigenvector centrality in the oceanic community network. Fig. 5a shows the outbound eigenvector centrality of each community with their interactions, calculated using average weights. This approach allows us to determine which communities are most influential in transmitting heat flows to other parts of the network. Considering the measurement’s characteristic, one can see that communities 7, 2, and 8 have high influence to the system in terms of their outgoing heat transfer. These communities are related to a diffuse community spread across the globe, primarily in the Indian and Atlantic Oceans, with partial coverage in the Pacific Ocean, indicating the importance of these oceanic systems in heat transfer (Fig. 5b).

Fig. 5.

Identification of key communities using outbound eigenvector centrality in the network of communities. (a, c) Networks illustrate outbound eigenvector centrality using average weights and total weights, respectively, with node colors representing the centrality values. The top 10% of average weights are highlighted in darker grey. (b, d) The top 3 communities based on the outbound eigenvector centrality from (a), and (c) are shown. The node size is determined based on the number of links, both incoming and outgoing.

In Fig. 5c, outbound eigenvector centrality was recalculated using total weights, providing insight into each oceanic community's comprehensive influence when considering all outgoing connections. The nodes are colored according to their recalculated centrality values, revealing how different regions contribute to the overall heat flow in the ocean. In the network with total weights shown in Fig. 5c, the communities with high outbound eigenvector centrality are communities 6, 4, and 2 (Fig. 5d). These communities consist of two types of regional characteristics: two communities (communities 6 and 4) are close to the ENSO system, while another community (community 2) represents a more dispersed community around the globe, with stronger connections to the North Atlantic Ocean. This suggests that the North Atlantic Ocean and the outer part of the ENSO system may play a primary role in distributing heat flows throughout the oceanic system.

Examining the flow directions between these communities through Figs. 5c, it appears that community 6 directs its flow towards community 2 and 7, which are located near the ENSO system and also spread across the entire system. The fact that the community 11, representing the main parts of the ENSO system, does not rank high in outbound eigenvector centrality again suggests that its functional role may lie in bridging other oceanic regions, with the main outgoing flows likely originating from other oceanic regions near the ENSO system.

3.4.3. Inbound Eigenvector Centrality

In contrast to the previous analysis, inbound eigenvector centrality can highlight important sea regions in terms of their incoming flows, and Fig. 6 depicts the key communities by applying inbound eigenvector centrality within the oceanic community network. Panels (a) and (c) illustrate the inbound eigenvector centrality of each community, calculated using average weights and total weights, respectively. In panel (a), the average weights were used to determine which communities receive the most influence from their connected nodes, with nodes colored based on their centrality values. Communities with the highest inbound influence, such as communities 2, 6, and 7, are primarily associated with parts of the North Pole, the periphery of ENSO region, and the Atlantic Ocean, suggesting that these areas can play a major role in receiving and concentrating heat flows from other regions, considering their weight per unit link.

Fig. 6.

Identification of key communities using inbound eigenvector centrality in the network of communities. (a, c) Networks illustrate inbound eigenvector centrality using average weights and total weights, respectively, with node colors representing the centrality values. The top 10% of average weights are highlighted in darker grey. (b, d) The top 3 communities based on inbound eigenvector centrality from (a), and (c) are shown. The node size is determined based on the number of links, both incoming and outgoing.

In panel (c), inbound eigenvector centrality was recalculated using total weights to provide a broader understanding of the influence each community receives when considering the overall strength of all incoming connections. Communities with high inbound centrality in this analysis are more comprehensively linked, revealing which parts of the oceanic system act as major recipients of heat flows. The analysis shows that communities 7, 2, and 6 exhibit high inbound eigenvector centrality, with most of these communities being widely distributed across the Atlantic and the near ENSO regions. One can notice that these key communities are the same as those identified by inbound eigenvector centrality with average weights, as shown in Figs. a and b, with only the order of the top communities being different. In particular, the reciprocal relationship between community 6 and 7 in Fig. 6d can be crucial, as this mutual heat transfer involves the Indian Ocean, the Atlantic Ocean, and the peripheral regions of the ENSO system as major contributors to heat accumulation, driven by connections from various oceanic sources. Especially, the composition of communities 2 and 7 requires further study, as they include three oceanic regions within the same community, which may exhibit synchronous heat dynamics. Given that previous studies have reported the crucial role of interactions between the Indian and Atlantic Oceans in ocean warming, along with a partial role of the Pacific Ocean (Yang et al., 2022), these alignments of major inbound eigenvector centrality for key oceanic regions suggest that the global oceanic heat transfer network analyzed in this study could provide valuable insights into the functional roles of oceanic regions within the context of their interactions.

3.5. Functional Roles of Key Sea Regions

Fig. 7a focuses on the community network encompassing major oceanic regions—ENSO (community 11), Indian (community 8) and Atlantic Ocean (community 7), and the periphery of the ENSO system (community 6 and 9). The bidirectional links between these communities represent their mutual interactions (total weights), with asymmetric weights on the relationships highlighted through distinct color. Red arrows indicate links where one direction has significantly stronger influence—at least 1.5 times greater than the other—while the strongest mutual links are highlighted in light blue. The results underscore the complex interplay among these regions, particularly emphasizing the dominant influences from the main ENSO system (community 11) to the periphery of the ENSO system (community 9), from the Indian Ocean (community 8) to the Atlantic Ocean (community 7), and from the Atlantic Ocean (community 7) to the periphery regions of the ENSO system (community 6). In addition to these connections, interactions between communities 6 (peripheral regions of the ENSO system) and 11 (the main parts of the ENSO system) maintain strong mutual influences. These major connections demonstrate the strong correlations among the ENSO, Atlantic, and Indian Oceans, suggesting that these ocean regions act as critical hubs for balancing heat exchanges, with significant asymmetric influences involving the ENSO system and Indian Ocean.

Fig. 7.

Functional roles of major sea regions in community networks. (a) Five communities representing five ocean regions are extracted from sea community network with total weights: 6 (ENSO), 7 (Indian and Atlantic Oceans), 8 (Indian Ocean), 9 (Indian and ENSO), and 11 (ENSO). Each pair of nodes has bidirectional links, and links with asymmetric weights are highlighted with different colors. A link with a stronger weight that is at least 1.5 times greater than the weaker link’s weight is shown in red (weaker is in orange). The strongest bidirectional links are colored in light blue. (b) Edges with the top 10% edge betweenness values (ranges from 0 to 0.44) are displayed. A hierarchical layout from Cytoscape is used, with node colors and sizes are based on outbound eigenvector centrality. The node size is determined based on the number of links, both incoming and outgoing.

Fig. 7b provides a deeper examination of the network by identifying edges with the top 10% edge betweenness values, representing the most significant pathways within the oceanic community network. The nodes are colored based on their outbound eigenvector centrality, and the layout is constructed using the hierarchical algorithm in Cytoscape, illustrating the heat flow among communities from top to bottom in the heat transfer network. The high edge betweenness values indicate the primary shortest paths that facilitate major heat exchanges between different ocean regions. When we extract the edges with high edge betweenness, communities 4 (mostly in the Pacific Ocean), 7 (distributed in the Indian, Pacific, and the North Atlantic Ocean), and 8 (mainly in the Indian Ocean) are particularly prominent, suggesting their roles as essential conduits for thermal energy transfer. The hierarchical layout used in this visualization further clarifies the network structure, emphasizing how regions like the peripheral areas of the ENSO system (e.g., community 9, 1, and 6), which link multiple major basins, act as central nodes in the oceanic heat transfer network. These findings suggest that the outer parts of the ENSO system and diffusive oceanic communities (e.g., community 4, 7, and 9) not only accumulate heat but also play a crucial role in redistributing it globally, thereby contributing significantly to the regulation of global oceanic temperatures.