Ocean feedback to pulses of the Madden–Julian Oscillation in the equatorial Indian Ocean

Dynamical understanding of the Madden–Julian Oscillation (MJO) has been elusive, and predictive capabilities therefore limited. New measurements of the ocean's response to the intense surface winds and cooling by two successive MJO pulses, separated by several weeks, show persistent ocean currents and subsurface mixing after pulse passage, thereby reducing ocean heat energy available for later pulses by an amount significantly greater than via atmospheric surface cooling alone. This suggests that thermal mixing in the upper ocean from a particular pulse might affect the amplitude of the following pulse. Here we test this hypothesis by comparing 18 pulse pairs, each separated by <55 days, measured over a 33-year period. We find a significant tendency for weak (strong) pulses, associated with low (high) cooling rates, to be followed by stronger (weaker) pulses. We therefore propose that the ocean introduces a memory effect into the MJO, whereby each event is governed in part by the previous event.

T he Madden-Julian Oscillation (MJO) is a tropical intraseasonal oscillation that forms in the Indian Ocean and whose dominant surface signature in the equatorial Indian and western Pacific Oceans is a pulse of intense winds and heavy precipitation coupled with deep atmospheric convection lasting 2-3 days at a given location, covering an area of 50,000 km 2 and recurring at intervals of 30-90 days 1 . Its signal propagates around the globe in the equatorial band, where it influences weather phenomena as disparate as Atlantic hurricanes and the Pineapple Express along the west coast of North America 2 . Here, we use the term 'pulse' to represent the westerly wind burst at the sea surface, beginning at the time when the surface wind stress exceeds a threshold value and ending when the stress subsequently falls below this value.
The effect of sea surface temperature (SST) on some forms of weather phenomena is well established. For example, hurricanes derive most of their energy from heat stored in the upper ocean 3 . The significance of diurnal warm layers on the intensity of atmospheric convection has been shown through simulations 4 and observations 5 . It has been suggested that, in March 2015, anomalously warm SST in the western Pacific provided the heat energy for MJO amplification to record amplitudes 6 . However, neither the role of air-sea fluxes 7,8 nor that of internal ocean processes 9,10 in MJO development is clear at present.
By virtue of its high heat capacity (4 Â that of air) and density (10 3 Â air), the ocean responds slowly to energetic, rapidly moving and evolving atmospheric disturbances like the MJO. These same factors also mean that the oceanic response continues after the atmospheric disturbance has passed. They raise the issue that there may be an effect of this slowly evolving response on subsequent atmospheric disturbances.
Here, we link insight gained from recent short-term field measurements with a statistical analysis of the existing longer-term record to show an unanticipated dependence of MJO intensity on SST. Dynamics of the Madden-Julian Oscillation (DYNAMO) 11,12 was a large-scale air-sea interaction experiment in the central equatorial Indian Ocean in boreal autumn of 2011 that included local and coincident measurements of surface fluxes and subsurface currents, stratification and mixing. Early results from shipboard measurements made during this experiment showed the persistence of the subsurface disturbance following passage of the MJO in the atmosphere 13 . Here we include longer-term moored records, also part of DYNAMO, that demonstrate the full subsurface response to two consecutive MJO pulses, a strong pulse following a weak pulse. Since it was at least possible that the greater upper ocean heat content (HC) remaining after the weak pulse might have contributed to the stronger following pulse, this anecdotal result prompted us to pose a hypothesis that stronger (weaker) pulses always follow weak (strong) pulses. To address this hypothesis we exploited a set of reanalysis data products that begins in 1980 and includes sufficient information to identify MJO pulses, and to quantify their intensity as well as pre-and post-pulse SST. These data revealed 84 MJO pulses between 1980 and 2013. Of these, 18 pulse pairs separated by o55 days were identified. These pulse pairs showed a significant correlation between SST cooling produced by the initial pulse and the intensity of the following pulse. The change in SST also correlated negatively with the intensity of each MJO pulse. Taken together, these suggest negative feedback between successive MJO pulses, thus supporting our stated hypothesis.

Results
Details of upper ocean response to two successive MJO pulses. Subsurface measurements during and following pulses of the MJO in the Indian Ocean suggest that the ocean's response may influence subsequent pulses. These measurements have led to distinctions in both form and intensity between observed pulses [12][13][14] . Two pulses in October and November 2011 were particularly well-sampled and helped in illustrating differences between relatively weak and relatively strong pulses. We next assess the generality of these patterns using historical data.
A large-scale quantification of the intensity of each pulse is provided by the real-time multivariate MJO (RMM) index 15 . RMM provides an assessment based on outgoing longwave radiation (OLR) and winds at 200 and 850 hPa. RMM has both amplitude and phase, where phase roughly references geographic regions around Earth at the equator. The central equatorial Indian Ocean, the region of detailed measurements discussed here, roughly corresponds to phase 3 (ref. 16), and we refer to the value of RMM computed there as 3 RMM. A value of 3 RMM41 is considered to be a significant event. The examples in Fig. 1 are characterized by 3 RMM ¼ 1.3 (significant but relatively weak pulse) and 3 RMM ¼ 2.3 (strong pulse), where these values represent three-day averages centred at the peak daily value.
The two pulses were sensed at three locations along the equator where detailed subsurface measurements were made (0°, 78°E; 0°, 80°E; 0°, 90°E) from oceanographic moorings instrumented with turbulence sensors 17 . An averaged oceanic response was derived by shifting the time series using the nominal MJO pulse propagation speed (5 m s À 1 ) to reference time at each location to the arrival of the pulse. Time series of surface wind stress (t), SST, zonal velocity (u), heat content anomaly (dHC), net heat flux at the sea surface (J 0 q ) and subsurface turbulent heat flux (J t q ) were averaged daily.
Significant distinctions are apparent between the relatively weak early MJO pulse and the later, stronger pulse ( Fig. 1). During the stronger pulse surface forcing was greater and lasted longer, represented by large t (Fig. 1a) and large negative J 0 q (Fig. 1e). Before pulse 2, SST was 0.5 K warmer than it was before pulse 1 (Fig. 1b), after which it cooled by 1.0 K, compared with a net cooling of o0.5 K for pulse 1. The high SST preceding the MJO pulse has been referred to as the positive intraseasonal SST anomaly 18,19 . Upper ocean HC (integrated over 0-40 m) before pulse 1 was smaller than it was before pulse 2 by 44 MJ m À 2 . Following pulse 1, cooling led to a minimum value dHC ¼ À 55 MJ m À 2 4 days after the wind burst when J 0 q reversed sign to heating, and the upper ocean was heated at a rate of 55 W m À 2 , recovering the full heat lost by the pulse after 16 days (Fig. 1c). In contrast, seven days after pulse 2, the ocean had lost 117 MJ m À 2 of heat energy. The initial recovery coincided with the sign reversal of J 0 q , but was weaker (31 W m À 2 ) until day 13, followed by rapid heating (158 W m À 2 ). After 16 days, the HC was still smaller by 60 MJ m À 2 than its initial value and indeed, smaller than the pulse 1 minimum. Note that the cooling of SST during the disturbed (cloudy) state before the WWB is not reflected in dHC. It is primarily due to a deeper and weaker diurnal warm layer.
The greater wind stress of pulse 2 accelerated an eastward jet in the upper 100 m of the ocean by 0.8 m s À 1 over 8 days (Fig. 1d), causing a doubling of eastward mass transport 13 , a deepening of the jet and an enhancement of shear at its base. The combination of strong surface forcing and enhanced shear is presumably responsible for the increased and persistent subsurface turbulence, expressed here as J t q (Fig. 1e), the average turbulent heat flux 20 over the upper 79 m (see 'Methods' section). The subsurface cooling inferred from J t q through down-gradient, Fickian-like diffusion (deeper waters are cooler) was large for at least an additional 7 days past the cessation of the wind burst. The slope change in dHC after pulse 2 coincides with diminished J t q . Larger and extended cooling of the upper ocean during the passage of pulse 2 is consistent with larger and extended surface forcing. But after the passage of the wind burst, dHC and SST recovery rates were considerably less than following pulse 1, despite comparable net surface heating. Figure 1e suggests that cooling from below via subsurface mixing is the likely cause of this. The net result is that the stronger MJO pulse reduces upper ocean HC more than the weaker pulse does, and this cooling continues after the passage of the atmospheric disturbance. We next assess the generality of these patterns using historical data.
Properties of MJO pulses 1980-2013. The differences in oceanic response to the varying intensity of the pulses is important if upper ocean HC is a significant contributor to MJO pulse intensity. This would be the case if the greater heat available following a weak pulse feeds back into a stronger following pulse. We test this idea using a longer record (1980-2013) of atmospheric and surface variables derived from reanalysis products of daily, gridded estimates of OLR (ref. 21), t, SST and J 0 q (ref. 22) in the tropics. Data were averaged over 3 days in time, and spatially over 6°of latitude and 6°of longitude centred at 0°, 80°E, roughly the centre of RMM phase 3. These data provide a broader spatial perspective than is obtained from the set of point measurements shown in Fig. 1. From these data, MJO pulse pairs were selected under the following criteria: we required that 3 RMM41, that pulses propagated eastward in OLR at 3-9 m s À 1 , that zonal winds at 850 hPa (National Centers for Environmental Prediction (NCEP) reanalysis) were westerly and that the pulse pair separation was o55 days. The resulting pulse durations ranged from 3 days (which is our averaging period, hence minimum length) to 14 days with mean and median of 6 days.
In our analysis, we use SST as a proxy for HC. While it is SST that directly communicates heat from ocean to atmosphere, a thin layer cannot do so for very long. A more representative measure of the ocean's capacity to effect significant influence on the atmospheric boundary layer above is obtained by integrating over the upper ocean to determine its HC. To assess the correspondence of SST to HC, subsurface temperature and SST data from the Research Moored Array for African-Asian-Australian Monsoon Analysis and Prediction (RAMA) mooring at 0°, 80°E between late 2008 and through 2012 was averaged daily and HC computed over the upper 40 m. The comparison of HC and SST is shown in Fig. 2. The strong correlation between the two provides a rationale for using SST as a proxy for HC.
From the complete records starting in 1980, our selection criteria delivered 18 MJO pulse pairs. The beginning of each pulse was identified as the time at which westerly t increased above 0.025 N m À 2 , accompanied by a sign change in J 0 q from surface heating to surface cooling. The end of each pulse coincides with t decreasing below the same threshold value. The data series generated from reanalysis products were referenced to the beginning of the principal wind burst, so that resulting series are in terms of time since the beginning of the wind burst. They were then simply averaged into composite weak ( 3 RMMr1.8) and strong ( 3 RMM41.8) time series; their characteristics are depicted in Fig. 3. The stronger (blue) MJO pulses are characterized by larger negative OLR anomalies indicating cooler cloud tops, hence deeper convection (Fig. 3a), greater cooling of the sea surface (Fig. 3b) and larger and longer-lasting surface wind stress (Fig. 3c). SST cooling (Fig. 3d) begins before the principal wind burst (as in Fig. 1) presumably during the disturbed state preceding the active MJO. As in the DYNAMO pulses represented in Fig. 1, there is no reason to suspect that upper ocean HC decreases since J 0 q 40 during this period. This early decrease in SST is likely due to suppression of the near-surface diurnal warm layer. Net SST cooling is, on average, 0.25 K for weak pulses, the minimum occurring 5 days after the ARTICLE wind burst begins and 0.5 K for strong pulses, with the minimum occurring 10 days after the beginning of the wind burst (t ¼ 0). Consequently, after 30 days, SST has fully recovered to its pre-pulse value after the weak pulses but is still 0.3 K smaller 30 days after the strong pulses.
A particular feature to note in Fig. 3c is the twin peaked structure in wind stress during strong MJO pulses. The peaks are separated by 5 days, as in the stronger, well-sampled DYNAMO pulse in Fig. 1. The signal in Fig. 1 has been associated with a pair of atmospheric Kelvin waves tracked in satellite precipitation records and embedded within the greater envelope of the propagating MJO pulse 13 . Is it possible that stronger MJO pulses are always associated with convectively coupled Kelvin waves in the central equatorial Indian Ocean?
Stronger MJO pulses cause greater SST cooling. For each individual pulse, the net sea surface cooling was determined as the SST difference between post-wind burst minima (for to15 days) and the value at t ¼ 0. In Fig. 4, paired MJO pulses are colour coded, the first in the pair (pulse n ) denoted by a circle and the second (pulse n þ 1 ) by a triangle. The stronger pulses (quantified by larger 3 RMM, an index independent of SST) are associated with stronger SST cooling (dSST in Fig. 4). If SST is a reliable proxy of HC, then Fig. 4 indicates that, on average, pulses following strong initial pulses pass over an upper ocean with reduced HC compared with those following weak pulses.
Negative feedback between pulse pairs. Apparently, the atmosphere responds to these changes in ocean HC, at least over relatively short time scales (o55 days). When we compare the intensities of following pulses ( 3 RMM n þ 1 ) to SST changes caused by the corresponding initial pulses (dSST n ), there is significant correlation (Fig. 5). This correlation progressively disappears when the time between pulses exceeds 55 days. The sense of the correlation in Fig. 5 is such that weak (strong) SST cooling is followed by a stronger (weaker) pulse. This is consistent with a scenario in which short-term changes in upper ocean HC feed back (negatively) to the intensity of MJO pulses, acting as a short-term governor on the system. Stronger MJO pulses generally follow larger pre-pulse SST. The correlation shown in Fig. 5 is also consistent with a broader potential association of high SST and high 3 RMM (Fig. 6). The occurrence and intensity of tropical cyclones (TCs, including hurricanes and typhoons) is dynamically associated with warm SST through its effect on the moisture content of the atmospheric boundary layer. However, showing this statistically (important for projecting occurrence frequencies and intensities in a changing climate) has not been straightforward. Tropical SST has been shown to be correlated with the maximum cubed surface wind speed within hurricanes, a direct measure of the power in individual TCs (ref. 23). A strong relationship also exists between the frequency of TCs with SST (ref. 24), however indirect the SST/TC number relationship. These studies benefited from relatively long records (480 years) and direct measurements of cyclone winds and occurrences.
By comparison, our record of 3 RMM is short and measurements are less direct. However, a relatively unsophisticated analysis using 3 RMM of the 84 MJO pulses identified in our record compared with SST averaged locally at 0°, 80°E is suggestive. We sorted the 84 pulses into three ranges of preexisting SST and then counted the number of pulses in six ranges of 3 RMM (Fig. 6). This shows a progressively greater propensity for large values of 3 RMM to occur with progressively larger SST.
It is possible that 3 RMM is not the best metric for this analysis. 3 RMM is not a direct measure of energy while SST, as a proxy for HC, is. A more insightful analysis might potentially use SST over a broader equatorial expanse of the Indian Ocean to include the origins of the MJO pulses.

Discussion
For the two MJO pulses we have been able to observe in detail the total heat extracted from the upper ocean is lost in roughly equal parts to the atmosphere via J 0 q and to the deeper stratified ocean via J t q (Fig. 1). Through the first MJO, 19 MJ m À 2 was lost to surface fluxes over the first 3 days of the wind burst and 24 MJ m À 2 lost to subsurface fluxes over 5 days. In the second MJO pulse, 69 MJ m À 2 of heat energy was lost to surface fluxes over 2 two-day periods associated with each individual wind burst and 85 MJ m À 2 to subsurface fluxes over 12 days. The extended cooling from below is due to turbulence generated via instabilities by the persistent shear at the base of the Yoshida-Wyrtki Jet that is excited by the wind bursts. Roughly 25% greater cooling is attributed to subsurface fluxes compared with surface fluxes.
Contrast this to two extreme cases. In the open ocean in moderate winds and away from strong currents, nighttime convection is principally due to cooling of the ocean from above by the atmosphere. A little more than 10% of the cooling of the mixed layer is effected by mixing of cooler waters across the mixed layer base through penetrative convection 25 . At the other extreme, TCs extract their energy from the ocean, which contributes several thousands of W m À 2 to cooling the upper ocean during cyclone passage 26 . In this case subsurface fluxes are 5-10 Â greater than surface fluxes. These two extremes can be characterized by J 0 q =J t q 4 41 (nighttime convection) and J 0 q =J t q o o1 (TC). The MJO pulses observed here represent J 0 q =J t q ffi 1.

ARTICLE
In the case of the TCs, the redistribution of heat between near-surface and deeper ocean layers is so great that it has been suggested (via modelling) that the resultant deep warm anomalies, or deep mixed layers persist long enough to act as positive feedback to subsequent TCs 27 . That is, the sense of the correlation is opposite to that shown in Fig. 5. The deeper mixed layers induce smaller surface cooling as subsequent TCs pass, thus promoting stronger cyclone intensity. That the opposite correlation is found in Fig. 5 for MJO feedback must be due to the far less extreme values of J t q during and following the passage of MJO pulses.
The correlation in Fig. 5 suggests a significant influence of MJO pulses on those immediately following, indicating that the MJO leaves a memory of its passage behind as it propagates around the globe. The 55-day time scale is consistent with the relatively slow response time of the ocean, suggesting that the memory effect reflects heat storage in the upper ocean.
If the change in upper ocean HC effected by a pulse influences the following pulse, and the ocean contributes as much to reducing ocean HC via subsurface fluxes as does the atmosphere via surface fluxes, then there ought to be a general improvement in forecast skill by proper inclusion of subsurface mixing. The influence of the ocean on the atmosphere is via surface fluxes through the air-sea interface, which depend in part on SST. We have shown that subsurface fluxes contribute significantly to SST. Therefore, surface fluxes and internal ocean processes must both contribute to MJO intensity.

Methods
Identification of the pulse pairs in Figs 4 and 5. The Wheeler-Hendon computation of 3 RMM uses OLR and winds at 200 and 850 hPa to make a fairly good predictor of MJO-like activity in the equatorial Indian Ocean. However, it does not guarantee eastward propagation of the OLR anomaly nor does it guarantee surface westerlies. These are further factors we have flagged for in our definition of MJO pulses. We select significant 3 RMM values occurring between September and May and also require that: eastward propagation of the OLR anomaly with a phase speed of 3-9 m s À 1 , winds at 200 hPa and 850 hPa are out of phase and winds are westerly at 850 hPa.

Confidence limits in Figs 4 and 5.
For each MJO pulse, dSST was computed as the difference in SST at its minimum (within 15 days following pulse onset) and the value at pulse onset. Confidence limits were estimated as the sum of the confidence limits of the minimum and maximum values of SST over the averaging domain (3 Â 3 degree box centred at 0°N, 82.5°E). In turn, these are 1.96s/On, where s is the s.d., n the number of data points (36 one-degree resolution data points per day) and the factor 1.96 represents 95% confidence for a normal distribution.
The same approach was used to estimate confidence limits for 3 RMM. We determined s for 3 RMM values within the averaging time window of 3 days.
Definition of J t q . The vertical (more strictly, diapycnal) flux of heat by turbulence is represented by Fickian diffusion, enhanced by a turbulence diffusion coefficient, 28). The temperature-variance dissipation rate (w T ) is computed by scaling temperature-gradient spectrum measured by fast thermistors on small autonomous instruments (wpods 20 ) attached to oceanographic moorings and packaged with inertial navigation units to quantify the component of flow speed past the sensor that is due to wave-induced motions of the surface float 29 . The other component of the flow speed is the ocean current speed, measured separately by velocity sensors on the mooring. The vertical temperature gradient, T z , is defined locally by the vertical motion of the sensor through the water or, on a larger scale, by additional temperature sensors on the mooring. The depthdependent vertical heat flux caused by turbulence is J t q ¼ À rC p K T T z , where r is the density of seawater and C p its heat capacity. For this experiment, wpods were distributed differently in depth for each of the three moorings and we simply averaged J t q over all depths (in Table 1) at each mooring to provide relative estimates of subsurface heat flux between the weak and strong pulses in Fig. 1 Data availability. The data that support the findings of this study are available from the corresponding author upon request.  30 . The mooring at 0°, 78°E was deployed specifically for the DYNAMO experiment by the Applied Physics Lab at University of Washington. All three moorings were equipped with meteorological sensors from which the wind stress and net surface heat flux were derived from bulk formulae 31 . This table shows depths of velocity, temperature and wpod measurements on each mooring. HC; heat content.