the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# The rate of information transfer as a measure of ocean–atmosphere interactions

### Stéphane Vannitsem

### Alessio Bellucci

Exchanges of mass, momentum and energy between the ocean and atmosphere are of large importance in regulating the climate system. Here, we apply for the first time a relatively novel approach, the rate of information transfer, to quantify interactions between the ocean surface and the lower atmosphere over the period 1988–2017 at a monthly timescale. More specifically, we investigate dynamical dependencies between sea surface temperature (SST), SST tendency and turbulent heat flux in satellite observations. We find a strong two-way influence between SST and/or SST tendency and turbulent heat flux in many regions of the world, with the largest values in the eastern tropical Pacific and Atlantic oceans, as well as in western boundary currents. The total number of regions with a significant influence by turbulent heat flux on SST and on SST tendency is reduced when considering the three variables (this case should be privileged, as it provides additional sources of information), while it remains large for the information transfer from SST and SST tendency to turbulent heat flux, suggesting an overall stronger ocean influence compared to the atmosphere. We also find a relatively strong influence by turbulent heat flux taken 1 month before on SST. Additionally, an increase in the magnitude of the rate of information transfer and in the number of regions with significant influence is observed when looking at interannual and decadal timescales compared to monthly timescales.

The climate on Earth is strongly affected by exchanges of mass, momentum and energy between the ocean and atmosphere. The ocean absorbs a large amount of solar energy and releases part of this energy to the atmosphere. In turn, the atmosphere modifies the ocean state through changes in wind, humidity and temperature. The classical view is that the slowly changing upper ocean is modulated by the high-frequency atmospheric variability (Hasselmann, 1976; Frankignoul and Hasselmann, 1977). While this paradigm has been successful in explaining the variability in sea surface temperature (SST) and surface heat flux over large parts of the ocean, it has been challenged over ocean regions characterized by intense mesoscale activity, such as western boundary currents and the Antarctic Circumpolar Current (Chelton et al., 2004; Brachet et al., 2012; Kirtman et al., 2012; Bishop et al., 2017; Roberts et al., 2017; Small et al., 2020; Bellucci et al., 2021).

Chelton et al. (2004) used 25 km resolution satellite radar scatterometer measurements over 1999–2003 and revealed the existence of persistent small-scale features in wind stress. According to Chelton et al. (2004), much of that mesoscale variability is attributable to SST modification. In particular, they found that surface wind speed is locally higher over warm water and lower over cool water, i.e., a positive correlation that is opposite to the one found at a large scale. Bishop et al. (2017) found that monthly scale lead–lag correlations between SST, SST tendency and turbulent heat flux (THF, positive upwards) allow one to discriminate between atmospheric-driven variability and ocean-led variability using results from a stochastic energy balance model (Wu et al., 2006) designed to represent ocean–atmosphere interactions, as well as satellite observations over 1985–2013. In their analysis, when the SST variability is dominated by atmospheric weather, SST tendency is negatively correlated with THF anomalies; on the other hand, when the SST fluctuations are driven by intrinsic ocean processes (ocean weather), the correlation between SST and THF is positive. Following this approach, they showed that eddy-rich regions associated with pronounced SST gradients, such as western boundary currents, are characterized by ocean-driven SST variability, while less dynamically active open-ocean regions, characterized by weaker SST gradients, exhibit an atmosphere-driven SST variability. These two regimes are reproduced by both eddy-parameterized (∼1^{∘} spatial resolution) and eddy-permitting (∼0.25^{∘}) coupled global climate models (Bellucci et al., 2021), with an increased ocean resolution leading to a substantially improved representation of SST and THF cross-covariance patterns.

According to a commonly accepted interpretation, the positive SST–THF zero-lag correlation (identifying an ocean-driven regime) indicates a damping role of THF on the SST anomalies generated by ocean dynamics, while the negative SST-tendency–THF correlation (identifying an atmosphere-driven regime) is attributed to the ocean surface cooling determined by the release of heat from the ocean to the atmosphere. However, the presence of such correlations does not firmly demonstrate causal influences between these variables, as correlation does not mean causation. Thus, the use of a dedicated causal method is of crucial importance to corroborate these findings. Several causal inference frameworks have been developed in the past years to identify such causal links (Granger, 1969; Liang and Kleeman, 2005; Sugihara et al., 2012; Krakovská et al., 2018; Paluš et al., 2018; Runge et al., 2019).

The Liang–Kleeman information flow method allows one to identify the direction and magnitude of the cause–effect relationships between variables (Liang and Kleeman, 2005). It is based on the rate of information transfer in dynamical systems and is rigorously derived from the propagation of information entropy between variables (Liang, 2016). It was initially developed for two-variable systems (Liang, 2014) and has recently been extended to multivariate systems (Liang, 2021). Compared to other causal inference frameworks, the rate of information transfer is a relatively simple index to compute one-way and/or two-way dependencies between variables. This novel method has been successfully applied to several climate studies, e.g., causal influences between greenhouse gases and global mean surface temperature (Stips et al., 2016; Jiang et al., 2019; Hagan et al., 2022), dynamical dependencies between a set of observables and the Antarctic surface mass balance (Vannitsem et al., 2019), soil moisture and air temperature interactions in China (Hagan et al., 2019), prediction of El Niño Modoki (Liang et al., 2021), causal links between climate indices in the North Pacific and Atlantic regions (Vannitsem and Liang, 2022), and identification of potential drivers of Arctic sea-ice changes (Docquier et al., 2022).

In our study, we analyze upper-ocean–lower-atmosphere interactions using the rate of information transfer developed by Liang (2021). More specifically, we check the two-way influences between SST and/or SST tendency and THF at the air–sea interface in satellite observations. Thus, our study allows one to go one step further than previous studies (Bishop et al., 2017; Bellucci et al., 2021), which have mainly focused on lead–lag correlation analyses, by identifying causal links between these variables. Section 2 presents the data and methods used in this analysis. Section 3 provides the main results of our study and places them in the overall context. Our conclusions are presented in Sect. 4.

## 2.1 Data

We use version 3 of the Japanese Ocean Flux Data Sets with the Use of Remote-Sensing Observations (J-OFURO3; Tomita et al., 2019; Tomita, 2020). This dataset uses the data of multiple satellites to estimate surface fluxes between the ocean and atmosphere over sea-ice-free regions with a resolution of 0.25^{∘}. It makes use of passive microwave radiometers and scatterometers available from 1988 to 2017 (see Tomita et al., 2019 for further details). From this dataset, we extract monthly mean latent and sensible heat fluxes, as well as SST. The latter is computed as an ensemble median obtained from various global SST products.

We also use the SeaFlux Data Products to estimate the observational uncertainty. It consists of estimates of ocean surface latent and sensible heat fluxes, among other variables (Roberts et al., 2020). It relies on the use of the Special Sensor Microwave Imager (SSM/I) and the Special Sensor Microwave Imager Sounder (SSMIS) over the period 1988–2018 (we use 1988–2017 to be consistent with J-OFURO3). The SST is also available for this dataset and is computed using the Reynolds Optimally Interpolated version 2.0 (Reynolds et al., 2007). As for J-OFURO3, SeaFlux data are available on a 0.25^{∘} grid, and we extract monthly mean latent and sensible heat fluxes and SST.

In the main text, we only show results from J-OFURO3, as results obtained with SeaFlux are largely consistent. The latter are presented in Appendix B.

## 2.2 Methods

Our analysis involves three variables, namely SST, SST tendency and THF, following the approach of Bishop et al. (2017) and Bellucci et al. (2021). The choice of these three specific variables is based on the stochastic energy balance model of Wu et al. (2006). As explained in Sect. 1, lead–lag covariances between these three variables are used as a way to diagnose ocean-driven and atmosphere-driven regimes. The goal of our analysis is to go beyond the correlation–covariance relationships identified by Bishop et al. (2017) and Bellucci et al. (2021) and to check the causal links between SST, SST tendency and THF. THF is defined as the sum of latent heat flux and sensible heat flux and is expressed in W m^{−2} (positive upwards), as in Bishop et al. (2017) and Bellucci et al. (2021). SST tendency is computed via a central difference approximation of SST (expressed in ^{∘}C) using a time step of 1 month (following Bishop et al., 2017) and is expressed in ^{∘}C month^{−1}.

We compute the rate of information transfer between SST, SST tendency and THF at each grid point of the globe using monthly data from 1988 to 2017. As we would like to extract statistical information from time series independent of any specific trend or cycle, we remove the trend and seasonality of all three variables using a linear regression and additive decomposition. The absolute rate of information transfer from variable *X*_{j} to variable *X*_{i} is computed assuming linearity following Liang (2021) as follows:

where **C** is the covariance matrix, *d* is the number of variables, Δ_{jk} are the cofactors of **C**, *C*_{k,di} is the sample covariance between *X*_{k} and the Euler forward difference approximation of d${X}_{i}/$d*t* (d*t* is the time step and equals 1 month in our study), *C*_{ij} is the sample covariance between *X*_{i} and *X*_{j}, and *C*_{ii} is the sample variance of *X*_{i}.

To assess the relative importance of the different cause–effect relationships, we compute the relative rate of information transfer from variable *X*_{j} to variable *X*_{i} following Liang (2021) as follows:

where *Z*_{i} is the normalizer, computed as follows:

where the first term on the right-hand side represents the information flowing from all the *X*_{k} to *X*_{i} (including the influence of *X*_{i} on itself), and the last term is the effect of noise and unobserved processes, computed following Liang (2021).

When *τ*_{j→i} is statistically different from 0 (either positive or negative), *X*_{j} has an influence on *X*_{i}, while if ${\mathit{\tau}}_{j\to i}=\mathrm{0}$, there is no influence. A value of |τ|=100 *%* indicates that *X*_{j} has the maximum influence on *X*_{i}. A positive (negative) value of *τ*_{j→i} means that the variability in *X*_{j} makes the variability in *X*_{i} more uncertain (certain; Liang, 2014); i.e., it increases (decreases) the variability in *X*_{i} (Appendix A; Fig. A1). The statistical significance of *τ*_{j→i} is computed via bootstrap resampling with the replacement of all the terms included in Eqs. (1)–(3) using 500 realizations. These bootstrap realizations are combined together using the false discovery rate (FDR) from Wilks (2016) with a significance level of 5 *%* to account for the multiplicity of tests.

Bishop et al. (2017) identified a strong positive zero-lag covariance between SST and THF at a monthly timescale (over 1985–2013) in western boundary currents, the Agulhas Return Current and eastern tropical Pacific using the OAFlux dataset (1^{∘} resolution) for THF and the NOAA OISST dataset (0.25^{∘} resolution) for SST (see their Fig. 3b). They also found a strong negative zero-lag covariance between SST tendency and THF over many regions of the globe, with the largest values at mid-latitudes (see their Fig. 3e). Using the J-OFURO3 dataset (0.25^{∘} resolution), we find similar results in terms of zero-lag covariance (Fig. B1). Additionally, when mapping the Pearson correlation coefficient instead of the covariance, we find a strong positive correlation between SST and THF in many regions of the world, with the largest values in the eastern tropical Pacific and Atlantic regions and in western boundary currents (Fig. 1a). A strong negative correlation between SST tendency and THF is also identified in most parts of the world, with the exception of a relatively narrow band along the Equator and in western boundary currents (Fig. 1b).

Bishop et al. (2017) and Bellucci et al. (2021) showed that regions of high SST gradient and THF (such as the Gulf Stream) are characterized by an ocean-driven regime. In these regions, the SST–THF zero-lag covariance is positive, suggesting that ocean processes drive SST anomalies, and the SST-tendency–THF lead–lag covariance is anti-symmetric (positive covariance at lag −1 and negative covariance at lag +1). On the contrary, an anti-symmetric SST–THF lead–lag covariance and a negative SST-tendency–THF zero-lag covariance are typical of an atmosphere-driven regime (such as in the North Atlantic subtropical gyre). In the latter case, the release of heat flux from the ocean to the atmosphere acts to cool the upper ocean. While this approach is interesting in terms of identifying whether the SST variability is driven by ocean or atmosphere processes, it does not precisely indicate whether the SST causally influences THF or the other way round. In our study, we quantify the causal relationships between SST and/or SST tendency and THF using the rate of information transfer from Liang (2021).

## 3.1 The two-dimensional (2D) case

If we only take into account SST and THF (i.e., two-dimensional system, hereafter referred to as 2D) in the computation of the rate of information transfer, we find that many regions are characterized by relatively strong two-way influences between SST and THF (Fig. 2). The spatial distribution of the rate of information transfer is relatively similar to the one of the correlation coefficient between SST and THF (Fig. 1a), with the largest values (either positive or negative) in the eastern tropical Pacific and Atlantic regions, western boundary currents and many parts of the Southern Hemisphere. Interestingly, the rate of information transfer is mainly positive for the influence of SST on THF (Fig. 2a), while negative values dominate for the reverse influence (Fig. 2b). This suggests that SST variability generally increases THF variability, while THF variability mainly constrains SST variability. Also, the regions where the influence of SST on THF is strongest are also characterized by a strong influence of THF on SST (in absolute values).

Regarding the SST-tendency–THF relationship, using only these two variables in the computation of the rate of information transfer also provides relatively strong two-way influences (either positive or negative) in many regions of the world (Fig. 3). The spatial distribution is more contrasted than with the SST–THF relationship, as both positive and negative values are now present in both directions, especially for the influence of SST tendency on THF (Fig. 3a). The information transfer from SST tendency to THF is characterized by positive southwest–northeast bands in the North Atlantic and North Pacific, positive northwest–southeast bands in the South Atlantic and South Pacific, and negative values between these regions (Fig. 3a). The information transfer from THF to SST tendency shows relatively symmetrical behavior to that from SST tendency to THF, with positive (negative) values where the reverse information transfer is negative (positive; Fig. 3b). This indicates that, in regions where the variability in SST tendency increases (decreases) the variability in THF, the variability in THF decreases (increases) the variability in SST tendency.

In summary, the 2D analysis shows additional information compared to previous lead–lag correlation studies (Bishop et al., 2017; Bellucci et al., 2021). In particular, we find that the ocean surface influences the lower atmosphere not only in strong boundary currents but also in many other regions of the world (Figs. 2a and 3a). In turn, the lower atmosphere (via surface heat fluxes) influences the ocean surface not only at mid-latitudes but also in tropical regions and western boundary currents (Figs. 2b and 3b). Our results are in line with Bach et al. (2019), who also find significant two-way influences between the upper ocean and lower atmosphere in many regions of the world using the Granger causality but with some methodological differences (daily timescale and other atmospheric fields). This shows that the lead–lag covariance analysis, while interesting in terms of identifying a particular ocean-driven or atmospheric-led regime, is not sufficient to accurately quantify causal links between the upper ocean and lower atmosphere. Importantly, our analysis (like the ones from Bishop et al., 2017 and Bellucci et al., 2021) applies to the monthly timescale, and we will show results beyond this specific timescale in Sect. 3.4.

However, this 2D analysis excludes one of the three key variables considered in the stochastic energy balance model of Wu et al. (2006), either SST tendency in Fig. 2 or SST in Fig. 3. Thus, it may provide a false impression of a two-way influence emerging due to the absence of a hidden variable. Also, the symmetry between the two directions for both the SST–THF (Fig. 2) and SST-tendency–THF (Fig. 3) relationships is questionable due to inaccurate results that were found using a theoretical example (see Sect. 3.5 for further discussion of this aspect). That is why we repeated the computation of the rate of information transfer including all three variables in the next section (Sect. 3.2).

## 3.2 The three-dimensional (3D) case

We computed the rate of information transfer based on the three variables analyzed in Bishop et al. (2017) and Bellucci et al. (2021), namely SST, SST tendency and THF (hereafter referred to as 3D). The 3D case provides additional sources of information compared to the 2D case and should thus be preferred in terms of result interpretation.

In the 3D case, the influence of SST on THF (Fig. 4a) is very similar to that in the 2D case (Fig. 2a). This is logical, since SST has no significant influence on SST tendency (Fig. B2), so all the information from SST goes to THF, demonstrating the robustness of the approach. However, a much-reduced number of regions shows a significant rate of information transfer from THF to SST (Fig. 4b) compared to that in the 2D case (Fig. 2b). This reduction is due to the fact that we now take SST tendency into account in the computation of the rate of information transfer; part of the information transfer from THF also goes into SST tendency, as we will see below. Despite this reduction in the number of regions with a significant transfer of information, the eastern tropical Pacific and Atlantic regions still show a strong negative rate of information transfer, suggesting that THF variability constrains SST variability in these regions (Fig. 4b). Interestingly, some regions (such as the Agulhas Return Current) show a positive rate of information transfer from THF to SST in the 3D case (Fig. 4b), while it is negative or close to 0 in the 2D case (Fig. 2b).

In the 3D case, the rate of information transfer from SST tendency to THF (Fig. 5a) is very similar to that in the 2D case (Fig. 3a) for the same reason as for the SST–THF influence. However, a much-reduced number of regions shows a significant transfer of information from THF to SST tendency (Fig. 5b) compared to that in the 2D case (Fig. 3b). Similarly as for the influence of THF on SST, this is due to the inclusion of a third variable; specifically, part of the information transfer from THF also goes into SST. Nevertheless, we still find regions of significant influence, e.g., negative values in the North Atlantic and northeastern Pacific and positive values in tropical regions.

Thus, computing the rate of information transfer based on the three variables of interest (SST, SST tendency and THF) somehow partitions the total influence of the lower atmosphere (THF) into a contribution to SST (Fig. 4b) and another contribution to SST tendency (Fig. 5b). Additionally, the total number of regions with a significant influence by THF on SST and by THF on SST tendency (both combined) clearly decreases compared to the 2D case. As the influences of SST on THF and of SST tendency on THF remain strong in the 3D case (almost unchanged compared to the 2D case), this is in favor of a stronger ocean influence compared to the atmosphere influence at monthly timescale, especially in extra-tropical regions.

In western boundary currents, we find large values of the rate of information transfer from SST to THF (Fig. 4a), suggesting a strong ocean influence, in agreement with Bishop et al. (2017). However, in extra-tropical regions far away from western boundary currents, such as in the central North Atlantic, we also find a strong transfer of information from SST (Fig. 4a) and from SST tendency (Fig. 5a) to THF, which is generally stronger than the reverse influence (from THF to SST and to SST tendency). This somewhat brings into question previous findings that suggest an atmospheric-driven SST variability in these regions (Bishop et al., 2017; Bellucci et al., 2021), and this shows that lead–lag covariance analyses should be supplemented by causality studies. The use of the SeaFlux observational dataset provides results that are in broad agreement with J-OFURO3 (Figs. B3–B4), which confirms the robustness of our findings.

The extension to additional variables can of course be performed to refine the analysis further, making sure that the additional variables are not nearly parallel (otherwise, the singularity of the covariance matrix could numerically deteriorate the results). We can do this by including other fields, higher-order tendencies or lagged fields. A first step toward an analysis of such a type is provided below. In our study, we prefer to keep the three fields used in the stochastic energy balance model from Wu et al. (2006) and analyzed by Bishop et al. (2017) and Bellucci et al. (2021), i.e., SST, SST tendency and THF. We will therefore focus on the use of a lagged field (Sect. 3.3), as well as on the analysis of interannual to decadal variability (Sect. 3.4).

## 3.3 Lagged transfer of information

Due to the inertia of the ocean mixed layer, the SST does not necessarily respond directly to changes in THF (Deser et al., 2003; Shi et al., 2022). To take this effect into account, we added a fourth variable to our analysis, namely THF leading SST by 1 month, hereafter referred to as THF(-1). Thus, the following four variables are considered here: SST, SST tendency, THF and THF(-1). The rate of information transfer has been applied to lagged variables in a previous study to predict El Niño Modoki based on solar activity (Liang et al., 2021). Note that the rate of information is not entirely free of time lag, as it involves a tendency term (d${X}_{i}/$d*t* in Eq. 1), which is estimated in a discrete way.

We find that there is a significant positive rate of information transfer from THF(-1) to SST in eastern tropical Pacific and Atlantic regions, western boundary currents (Gulf Stream and Kuroshio Extension) and the Agulhas Return Current, as well as negative values in other parts of the world (Fig. 6a). Thus, the lagged analysis shows that THF taken 1 month before strongly controls SST variability, especially in northern extra-tropical regions where this influence is almost absent in the original 3D case (Fig. 4b). There is also additional information provided by this fourth variable in the transfer of information from THF(-1) to SST tendency (Fig. 6b), but this is mostly restricted to the eastern tropical Pacific and Atlantic (negative values). The influence of SST and/or SST tendency on THF(-1) is found to be very close to 0 and to be not significant almost everywhere (Fig. B5), which confirms the robustness of the method, as causality cannot go back in time.

Interestingly, most regions showing a positive rate of information transfer from THF(-1) to SST (Fig. 6a) also have a positive lead–lag covariance between THF(-1) and SST (Fig. 3c of Bishop et al., 2017). Bishop et al. (2017) show that these regions have a symmetric lead–lag structure between SST and THF, which is characteristic of an ocean-driven regime. However, these similarities between the rate of information transfer and the lead–lag covariance disappear when we look at the relationships between THF(-1) and SST tendency. According to Bishop et al. (2017), the strongest values of lead–lag covariance between THF(-1) and SST tendency appear in western boundary currents (negative values; Fig. 3f of Bishop et al., 2017), while the rate of information transfer is significant only in the eastern tropical Pacific and Atlantic (negative values; Fig. 6b).

## 3.4 Interannual to decadal variability

All previous results are based on monthly mean outputs. Thus, our results are valid at a monthly timescale. In order to figure out what happens at interannual and decadal timescales, we take the 12- and 120-month (respectively) running-mean SST and THF (moving average with the same number of samples), and we re-compute the rate of information transfer in the 3D case (SST, SST tendency and THF). This approach is similar to the one used in Vannitsem and Liang (2022), who found differences in the rate of information transfer between climate indices depending on the timescale used.

At the interannual timescale (12-month running mean), the influence of SST on THF encompasses approximately the same regions as at a monthly timescale but with a general increase in the magnitude of the rate of information transfer (Fig. 7a). We also find regions that have a negative rate of information transfer, such as in the western North Atlantic Ocean, whereas such regions are not present with monthly means. The reverse influence of THF on SST provides a more contrasted pattern compared to that in the original 3D case, with a reduced number of regions with negative values along the Equator but also the additional presence of regions with positive values (Fig. 7b).

At a decadal timescale (120-month running mean), almost the whole globe is covered by significant information transfer between SST and THF in the two directions (Fig. 8). Also, the magnitude of the rate of information transfer clearly increases at this timescale compared to at interannual and decadal timescales. These results suggest that ocean-atmosphere interactions become more pronounced at larger timescales.

## 3.5 Limitations of the method

The rate of information transfer used in our study has major advantages compared to other causal methods, including its derivation from the first principles of information entropy (Liang, 2016) and its relative simplicity (Eq. 1). However, several limitations exist, which need to be kept in mind.

The first limitation arises from the linearity assumption (Liang, 2014). While the method has been tested and validated with highly nonlinear synthetic examples (e.g., Liang, 2014, 2018, 2021), it provides an approximated solution, and a generalization to the fully nonlinear case is still to be developed to get a more accurate solution (Liang, 2021). However, recent studies applying the rate of information transfer to climate data show that the method is successful at representing some key causal influences, such as the interactions between El Niño and the Indian Ocean Dipole (Liang, 2014), CO_{2} and global mean temperature (Stips et al., 2016), and Arctic sea ice and its drivers (Docquier et al., 2022).

The second limitation is linked to the problem of hidden variables. If we omit an important variable in our system, the rates of information transfer might be biased. In our specific study, some differences appear between the 2D case and the 3D case (Sect. 3.2). As we do not include one of the key variables of the stochastic energy balance model designed by Wu et al. (2006) in the 2D case, we end up with a biased result. In particular, the rate of information transfer from THF to SST is overestimated in northern and southern extratropical regions (compare Fig. 2b to Fig. 4b). This overestimation also appears in the rate of information transfer from THF to SST tendency (compare Fig. 3b to Fig. 5b). Also, a relatively symmetric behavior appears in the 2D case between both causal directions for the SST–THF relationship (Fig. 2) and for the SST-tendency–THF relationship (Fig. 3), although a number of regions do not show such a symmetry. This somehow goes in hand with the inaccurate symmetry obtained in computing the rate of information transfer in the unidirectionally coupled Rössler systems (Paluš et al., 2018) when using only two variables (Paluš, 2022; Fig. C1f). Thus, the results obtained in the 2D case in our study must be taken with great care, and 3D results should be privileged. Further research should be carried out to better understand the impact of hidden variables in computing the rate of information transfer.

Finally, the information flow method depends on the sampling frequency when applied to unidirectionally coupled Rössler systems and possibly to other systems and real-case studies. In the case of the Rössler systems, the method becomes inaccurate when the sampling frequency is too coarse (Appendix C; Fig. C1a–e). In our real-case study, it is difficult to accurately quantify this effect due to the shortness of the time period (360 months). Further research is needed to better understand the effect of the sampling frequency on the computation of the rate of information transfer.

In summary, we find that the rate of information transfer provides a more detailed quantification of dependencies between SST, SST tendency and turbulent heat flux (THF) than previous classical correlation–covariance studies. We do not argue that causal methods should replace covariance analyses, but they should rather be used as a complement in order to get a better understanding of physical interactions between variables. We show that the ocean surface (SST and SST tendency) strongly drives changes in the lower atmosphere (THF) and that the lower atmosphere also has an important influence on the ocean surface in many regions of the world. This result is somewhat different from what has been found in covariance analyses, in which ocean-driven regimes exist in western boundary currents and where atmospheric-led regimes dominate in the open ocean (Bishop et al., 2017; Bellucci et al., 2021). It is, however, supported by another recent analysis using the Granger causality, which shows that many regions of the world present a significant two-way influence between the lower atmosphere and the upper ocean, with a stronger ocean influence in tropical regions and a stronger atmospheric impact in the extra-tropics (Bach et al., 2019). As the latter study presents methodological differences – i.e., it focuses on the daily timescale and uses different atmospheric variables – we need to be cautious in comparing it to our analysis. In any case, our interpretation of these results is that ocean–atmosphere interactions are more complex than presented by classical covariance analyses.

Furthermore, we find that the influence of THF is partitioned between SST and SST tendency if we consider the three variables together (3D case, which should be privileged due to the inclusion of additional sources of information) so that the single impact of either THF on SST or of THF on SST tendency is decreased in the 3D case compared to that in the 2D case. Also, the number of regions with a significant rate of information transfer from THF to SST and from THF to SST tendency (combined) is smaller than the one with a significant rate of information transfer from SST to THF and from SST tendency to THF (combined) in the 3D case. This suggests an overall stronger upper-ocean influence compared to the lower atmosphere. However, when adding THF taken 1 month before to take the lagged effect into account, we find that this variable has a relatively strong influence on SST in a large portion of the globe. Finally, a larger timescale (going from monthly to interannual and decadal) provides larger values in terms of the rates of information transfer between the three variables.

In our study, we only considered the ocean surface, but several studies have shown that variations in the ocean heat content are controlled by both air–sea fluxes and ocean heat transport convergence, with a more important role for the latter with a deeper integration of ocean heat and with higher-resolution climate models (Roberts et al., 2017; Small et al., 2020). Also, only observations were considered here. Thus, extending our analysis to the ocean heat budget terms in climate models would provide further insights into the causal influences between the ocean and atmosphere. This is of large importance, as ocean–atmosphere interactions constitute an important regulator of our climate. Finally, from a theoretical perspective, additional investigations of the role of hidden and lagged variables should be performed, as well as a comparison of the rate of information transfer to other causal methods.

As explained in the main text (Sect. 2.2), a positive value of the relative rate of information transfer *τ*_{j→i} means that the variability in *X*_{j} increases the variability in *X*_{i}, while a negative value means that the variability in *X*_{j} decreases the variability in *X*_{i}. More broadly, an increase in the rate of information transfer from *X*_{j} to *X*_{i} leads to an increase in the variability of *X*_{i}. We demonstrate this by computing the rate of information transfer *τ*_{j→i} from variable *X*_{j} to variable *X*_{i} (based on Eq. 3 in the main text) using a three-dimensional stochastic linear system of equations as follows:

where *X*_{1}, *X*_{2} and *X*_{3} are the three variables; *a*_{kl} are the different coefficients; *t* is time and varies between 0 and 100 with 100 000 time steps (Δ*t*=0.001); and *W*_{1}, *W*_{2} and *W*_{3} represent normal random noises (standard Wiener process). We set ${a}_{\mathrm{11}}={a}_{\mathrm{22}}={a}_{\mathrm{33}}=-\mathrm{1}$, and we vary the six other coefficients one by one with five different values between −1 and 1 (−1, −0.5, 0, 0.5 and 1). When varying one of the six coefficients, we set the other five coefficients to a fixed value (${a}_{\mathrm{12}}={a}_{\mathrm{13}}=\mathrm{0.5}$ and ${a}_{\mathrm{21}}={a}_{\mathrm{23}}={a}_{\mathrm{31}}={a}_{\mathrm{32}}=\mathrm{0}$).

We solve the linear system (A1) using the Euler–Maruyama method, and 40 different values of the random noise are taken in order to take into account the uncertainty related to the rate of information transfer. The variance in each variable is compared to the rate of information transfer from any other variable to this variable to test our hypothesis. Results show that when the rate of information transfer from *X*_{j} to *X*_{i} increases, the variance in *X*_{i} increases (Fig. A1).

In order to test the dependence of the information flow method on the sampling frequency, we compute the absolute rate of information transfer (Eq. 1) on the two unidirectionally coupled Rössler systems using the same parameters as Paluš and Vejmelka (2007) and Paluš et al. (2018) as follows:

where *ω*_{1}=1.015, *ω*_{2}=0.985, ${a}_{\mathrm{1}}={a}_{\mathrm{2}}=\mathrm{0.15}$, ${b}_{\mathrm{1}}={b}_{\mathrm{2}}=\mathrm{0.2}$ and ${c}_{\mathrm{1}}={c}_{\mathrm{2}}=\mathrm{10}$. We use the same time step for integration as Paluš et al. (2018); i.e., d*t*=0.0785. We use the same initialization method as Paluš et al. (2018), solve the equations using the fourth-order Runge–Kutta scheme, and run the model for 7.5 × 10^{6} time steps. We discard the first 1.5 × 10^{6} time steps from the analysis. We use a coupling strength *ϵ* that varies between 0 and 0.25 with an increment of 0.025. We test different numbers of samples per pseudo-period and retain here the five of them (20, 60, 80, 200 and 6000) that best illustrate our results.

As illustrated in Fig. C1, results depend on the number of samples per pseudo-period. When using 20 samples per pseudo-period, ${T}_{y\mathrm{1}\to x\mathrm{1}}>{T}_{x\mathrm{1}\to y\mathrm{1}}$ when the coupling strength *ϵ*≥0.2, which is physically not correct (Fig. C1a). For 60 samples per pseudo-period, both influences are comparable for larger coupling strengths (Fig. C1b). From 80 samples per pseudo-period, ${T}_{x\mathrm{1}\to y\mathrm{1}}>{T}_{y\mathrm{1}\to x\mathrm{1}}$ (Fig. C1c–e), which is physically correct, and one needs to reach ∼ 6000 samples per pseudo-period to have ${T}_{y\mathrm{1}\to x\mathrm{1}}=\mathrm{0}$ (Fig. C1e). Importantly, all six variables need to be considered in the computation of the rate of information transfer (multivariate formula; Liang, 2021) to reach this result. If we only consider *X*_{1} and *Y*_{1} and compute the rate of information transfer using the bivariate formula (Liang, 2014), we obtain an inaccurate result with a symmetry between the two directions whatever the number of samples per pseudo-period (Fig. C1f).

J-OFURO3 observational data (Tomita et al., 2019) are available on https://doi.org/10.20783/DIAS.612 (Tomita, 2020). SeaFlux observational data (Roberts et al., 2020) are accessible from NASA (https://doi.org/10.5067/SEAFLUX/DATA101). The Python scripts to compute the rate of information transfer and to produce the figures of this article are available on Zenodo: https://doi.org/10.5281/zenodo.7547073 (Docquier, 2023).

DD wrote the paper with contributions from all the co-authors. DD, SV and AB designed the science plan. DD analyzed the satellite data, ran the Rössler systems and produced the figures. All the authors participated in the interpretation of the results and provided useful comments to help improve the analysis.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We would like to thank Milan Paluš and an anonymous reviewer for providing insightful comments, which helped to improve our paper. We also thank the Editor Rui A. P. Perdigão for his encouragements in preparing an improved version of the paper. We thank X. San Liang for the useful discussions related to the rate of information transfer. We thank Claude Frankignoul for the very helpful feedback provided following one of our presentations. David Docquier, Stéphane Vannitsem and Alessio Bellucci are supported by ROADMAP (Role of ocean dynamics and Ocean–Atmosphere interactions in Driving cliMAte variations and future Projections of impact-relevant extreme events; https://jpi-climate.eu/project/roadmap/, last access: 5 May 2023), a coordinated JPI-Climate/JPI-Oceans project.

David Docquier and Stéphane Vannitsem have received funding from the Belgian Federal Science Policy Office (contract B2/20E/P1/ROADMAP) and Alessio Bellucci has been supported by MUR – Italian Ministry of University and Research (D.M. 593/2016).

This paper was edited by Rui A. P. Perdigão and reviewed by Milan Palus and one anonymous referee.

Bach, E., Motesharrei, S., Kalnay, E., and Ruiz-Barradas, A.: Local atmosphere-ocean predictability: Dynamical origins, lead times, and seasonality, J. Clim., 32, 7507–7519, https://doi.org/10.1175/JCLI-D-18-0817.1, 2019. a, b

Bellucci, A., Athanasiadis, P. J., Scoccimarro, E., Ruggieri, P., Gualdi, S., Fedele, G., Haarsma, R. J., Garcia‑Serrano, J., Castrillo, M., Putrahasan, D., Sanchez‑Gomez, E., Moine, M., Roberts, C. D., Roberts, M. J., Seddon, J., and Vidale, P. L.: Air‑Sea interaction over the Gulf Stream in an ensemble of HighResMIP present climate simulations, Clim. Dynam., 56, 2093–2111, https://doi.org/10.1007/s00382-020-05573-z, 2021. a, b, c, d, e, f, g, h, i, j, k, l, m

Bishop, S. P., Small, R. J., Bryan, F. O., and Tomas, R. A.: Scale dependence of midlatitude air-sea interaction, J. Clim., 30, 8207–8221, https://doi.org/10.1175/JCLI-D-17-0159.1, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t

Brachet, S., Codron, F., Feliks, Y., Ghil, M., Le Treut, H., and Simonnet, E.: Atmospheric circulations induced by a midlatitude SST front: A GCM study, J. Clim., 25, 1847–1853, https://doi.org/10.1175/JCLI-D-11-00329.1, 2012. a

Chelton, D. B., Schlax, M. G., Freilich, M. H., and Milliff, R. F.: Satellite measurements reveal persistent small-scale features in ocean winds, Science, 303, 978–983, https://doi.org/10.1126/science.1091901, 2004. a, b, c

Deser, C., Alexander, M. A., and Timlin, M. S.: Understanding the persistence of sea surface temperature anomalies in midlatitudes, J. Clim., 16, 57–72, https://doi.org/10.1175/1520-0442(2003)016<0057:UTPOSS>2.0.CO;2, 2003. a

Docquier, D.: Liang Index to quantify ocean-atmosphere interactions (v2), Zenodo [code], https://doi.org/10.5281/zenodo.7547073, 2023. a

Docquier, D., Vannitsem, S., Ragone, F., Wyser, K., and Liang, X. S.: Causal links between Arctic sea ice and its potential drivers based on the rate of information transfer, Geophys. Res. Lett., 49, e2021GL095892, https://doi.org/10.1029/2021GL095892, 2022. a, b

Frankignoul, C. and Hasselmann, K.: Stochastic climate models, Part II. Application to sea-surface temperature anomalies and thermocline variability, Tellus, 29, 289–305, https://doi.org/10.3402/tellusa.v29i4.11362, 1977. a

Granger, C. W.: Investigating causal relations by econometric models and cross-spectral methods, Econometrica, 37, 424–438, https://doi.org/10.2307/1912791, 1969. a

Hagan, D. F. T., Wang, G., Liang, X. S., and Dolman, H. A. J.: A time-varying causality formalism based on the Liang–Kleeman information flow for analyzing directed interactions in nonstationary climate systems, J. Clim., 32, 7521–7537, https://doi.org/10.1175/JCLI-D-18-0881.1, 2019. a

Hagan, D. F. T., Dolman, H. A. J., Wang, G., Sian, K. T. C. L. K., and Yang, K.: Contrasting ecosystem constraints on seasonal terrestrial CO2 and mean surface air temperature causality projections by the end of the 21st century, Environ. Res. Lett., 17, 124019, https://doi.org/10.1088/1748-9326/aca551, 2022. a

Hasselmann, K.: Stochastic climate models Part I. Theory, Tellus, 28, 473–485, https://doi.org/10.1111/j.2153-3490.1976.tb00696.x, 1976. a

Jiang, S., Hu, H., Lei, L., and Bai, H.: Multi-source forcing effects analysis using Liang–Kleeman information flow method and the community atmosphere model (CAM4.0), Clim. Dynam., 53, 6035–6053, https://doi.org/10.1007/s00382-019-04914-x, 2019. a

Kirtman, B. P., Bitz, C., Bryan, F., Collins, W., Dennis, J., Hearn, N., Kinter, J. L., Loft, R., Rousset, C., Siqueira, L., Stan, C., Tomas, R., and Vertenstein, M.: Impact of ocean model resolution on CCSM climate simulations, Clim. Dynam., 39, 1303–1328, https://doi.org/10.1007/s00382-012-1500-3, 2012. a

Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., and Paluš, M.: Comparison of six methods for the detection of causality in a bivariate time series, Phys. Rev. E, 97, 042207, https://doi.org/10.1103/PhysRevE.97.042207, 2018. a

Liang, X. S.: Unraveling the cause-effect relation between time series, Physical Review E, 90, 052150, https://doi.org/10.1103/PhysRevE.90.052150, 2014. a, b, c, d, e, f

Liang, X. S.: Information flow and causality as rigorous notions ab initio, Physical Review E, 94, 052201, https://doi.org/10.1103/PhysRevE.94.052201, 2016. a, b

Liang, X. S.: Causation and information flow with respect to relative entropy, Chaos, 28, 1–8, https://doi.org/10.1063/1.5010253, 2018. a

Liang, X. S.: Normalized multivariate time series causality analysis and causal graph reconstruction, Entropy, 23, 679, https://doi.org/10.3390/e23060679, 2021. a, b, c, d, e, f, g, h, i

Liang, X. S. and Kleeman, R.: Information transfer between dynamical system components, Phys. Rev. Lett., 95, 244101, https://doi.org/10.1103/PhysRevLett.95.244101, 2005. a, b

Liang, X. S., Xu, F., Rong, Y., Zhang, R., Tang, X., and Zhang, F.: El Niño Modoki can be mostly predicted more than 10 years ahead of time, Scientific Reports, 11, 17860, https://doi.org/10.1038/s41598-021-97111-y, 2021. a, b

Paluš, M.: RC2: 'Comment on egusphere-2022-942', EGUsphere, https://doi.org/10.5194/egusphere-2022-942-RC2, 2022. a

Paluš, M. and Vejmelka, M.: Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections, Phys. Rev. E, 75, 075307, https://doi.org/10.1063/1.5019944, 2007. a

Paluš, M., Krakovská, A., Jakubík, J., and Chvosteková, M.: Causality, dynamical systems and the arrow of time, Chaos, 28, 056211, https://doi.org/10.1103/PhysRevE.75.056211, 2018. a, b, c, d, e

Reynolds, R. W., Smith, T. M., Liu, C., Chelton, D. B., Casey, K. S., and Schlax, M. G.: Daily High-Resolution-Blended Analyses for Sea Surface Temperature, J. Clim., 20, 5473–5496, https://doi.org/10.1175/2007JCLI1824.1, 2007. a

Roberts, C. D., Palmer, M. D., Allan, R. P., Desbruyeres, D. G., Hyder, P., Liu, C., and Smith, D.: Surface flux and ocean heat transport convergence contributions to seasonal and interannual variations of ocean heat content, J. Geophys. Res.-Oceans, 122, 726–744, https://doi.org/10.1002/2016JC012278, 2017. a, b

Roberts, J. B., Clayson, C. A., and Robertson, F. R.: SeaFlux Data Products V1, NASA Global Hydrometeorology Resource Center DAAC, Huntsville, Alabama, USA [data set] https://doi.org/10.5067/SEAFLUX/DATA101, 2020. a, b

Runge, J., Bathiany, S., Bollt, E., Camps-Valls, G., Coumou, D., Deyle, E., Glymour, C., Kretschmer, M., Mahecha, M. D., Munoz-Mari, J., van Nes, E. H., Peters, J., Quax, R., Reichstein, M., Scheffer, M., Scholkopf, B., Spirtes, P., Sugihara, G., Sun, J., Zhang, K., and Zscheischler, J.: Inferring causation from time series in Earth system sciences, Nat. Commun., 10, 2553, https://doi.org/10.1038/s41467-019-10105-3, 2019. a

Shi, H., Jin, F.-F., Jacox, M. G., Amaya, D. J., Black, B. A., Rykaczewski, R. R., Bograd, S. J., Garcia-Reyes, M., and Sydeman, W. J.: Global decline in ocean memory over the 21st century, Sci. Adv., 8, https://doi.org/10.1126/sciadv.abm3468, 2022. a

Small, R. J., Bryan, F. O., Bishop, S. P., Larson, S., and Tomas, R. A.: What drives upper-ocean temperature variability in coupled climate models and observations, J. Clim., 33, 577–596, https://doi.org/10.1175/JCLI-D-19-0295.1, 2020. a, b

Stips, A., Macias, D., Coughlan, C., Garcia-Gorriz, E., and Liang, X. S.: On the causal structure between CO2 and global temperature, Scientific Reports, 6, 21691, https://doi.org/10.1038/srep21691, 2016. a, b

Sugihara, G., May, R., Ye, H., Hsieh, C.-H., Deyle, E., Fogarty, M., and Munch, S.: Detecting causality in complex ecosystems, Science, 338, 496–500, https://doi.org/10.1126/science.1227079, 2012. a

Tomita, H.: J-OFURO3, Data Integration and Analysis System (DIAS) [data set], https://doi.org/10.20783/DIAS.612, 2020. a, b

Tomita, H., Hihara, T., Kako, S., Kubota, M., and Kutsuwada, K.: An introduction to J-OFURO3, a third-generation Japanese ocean flux data set using remote-sensing observations, J. Oceanogr., 75, 171–194, https://doi.org/10.1007/s10872-018-0493-x, 2019. a, b, c

Vannitsem, S. and Liang, X. S.: Dynamical dependencies at monthly and interannual time scales in the climate system: Study of the North Pacific and Atlantic regions, Tellus A, 74, 141–158, https://doi.org/10.16993/tellusa.44, 2022. a, b

Vannitsem, S., Dalaiden, Q., and Goosse, H.: Testing for dynamical dependence: Application to the surface mass balance over Antarctica, Geophys. Res. Lett., 46, 12125–12135, https://doi.org/10.1029/2019GL084329, 2019. a

Wilks, D. S.: `The stippling shows statistically significant grid points': How research results are routinely overstated and overinterpreted, and what to do about it, B. Am. Meteorol. Soc., 97, 2263–2273, https://doi.org/10.1175/BAMS-D-15-00267.1, 2016. a

Wu, R., Kirtman, B. P., and Pegion, K.: Local air-sea relationship in observations and model simulations, J. Clim., 19, 4914–4932, https://doi.org/10.1175/JCLI3904.1, 2006. a, b, c, d, e

- Abstract
- Introduction
- Data and methods
- Results and discussion
- Conclusions
- Appendix A: Impact of the rate of information transfer on the variability
- Appendix B: Additional figures
- Appendix C: Dependence on the sampling frequency
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Data and methods
- Results and discussion
- Conclusions
- Appendix A: Impact of the rate of information transfer on the variability
- Appendix B: Additional figures
- Appendix C: Dependence on the sampling frequency
- Code and data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References