Diversity of the Madden-Julian Oscillation

See allHide authors and affiliations

Science Advances  31 Jul 2019:
Vol. 5, no. 7, eaax0220
DOI: 10.1126/sciadv.aax0220


Madden-Julian Oscillation (MJO) is the dominant mode of atmospheric intraseasonal variability and the cornerstone for subseasonal prediction of extreme weather events. Climate modeling and prediction of MJO remain a big challenge, partially due to lack of understanding the MJO diversity. Here, we delineate observed MJO diversity by cluster analysis of propagation patterns of MJO events, which reveals four archetypes: standing, jumping, slow eastward propagation, and fast eastward propagation. Each type exhibits distinctive east-west asymmetric circulation and thermodynamic structures. Tight coupling between the Kelvin wave response and major convection is unique for the propagating events, while the strength and length of Kelvin wave response distinguish slow and fast propagations. The Pacific sea surface temperature anomalies can affect MJO diversity by modifying the Kelvin wave response and its coupling to MJO convection. The results shed light on the mechanisms responsible for MJO diversity and provide potential precursors for foreseeing MJO propagation.


Madden and Julian (1, 2) discovered the most pronounced tropical intraseasonal mode of variability, which is now named as Madden-Julian Oscillation (MJO). The MJO strongly modulates the occurrence of many types of extreme weathers in the global tropics and midlatitudes, including floods, droughts, heavy snow storms, cold surges, heat waves, hurricanes and tropical storms, wildfires, aerosols, tornadoes, hail days, and many more (3, 4). Successful prediction of MJO propagation is vitally important for subseasonal prediction of tropical cyclones, extreme weather events, monsoon, and El Niño–Southern Oscillation (5). The prevalent eastward propagation is the most essential feature of the MJO yet remains a big challenge in numerical modeling and prediction over the past three decades (6).

One of the difficulties in numerical simulation and prediction of MJO arises from the complex behavior in the propagation of the MJO convective envelope that encloses other complex subscale activities. While, as shown by statistical analyses, the convective envelope of an averaged MJO event travels eastward with a mean speed of ~4 to 5 m/s over the warm pool of the Indo-Pacific oceans (6), the propagation, amplitude, and life cycles of the convective envelopes for individual MJO events differ considerably from event to event. This can be seen in Fig. 1. The large-scale intraseasonal convective anomalies [measured by outgoing longwave radiation (OLR) and representing the MJO convective envelope] exhibit various propagation forms along the equatorial belt, ranging from standing oscillation to systematic propagations with different speeds. The MJO diversity can also be seen from the MJO trajectory in the phase-space diagram constructed by the two leading multivariate empirical orthogonal functions that reflect the amplitude and propagation phases of the MJO (7). A better understanding on the causes of the event-to-event differences is critical for improving the simulation and prediction of MJO in global climate models (GCMs) and for a faithful projection of MJO’s future change. Thus far, most of the diagnostic analysis of MJO has been focused on the averaged properties using variety of effective statistical tools (8). The diversified mesoscale structures and the equatorial wave activity imbedded in the MJO envelope have been documented in previous literatures, e.g., (9, 10), but the complex forms and causes of MJO diversity in terms of propagation of MJO convective envelope have not been explored systematically. This work attempts an exploratory study on this issue.

Fig. 1 Example of diversified propagation patterns of the MJO.

Hovmöller diagram of intraseasonal OLR averaged between 10°S and 10°N for the 1983–1984 boreal winter and the 1996–1997 boreal winter is shown. The blue contours indicate −10 W m−2 contours. The thick horizontal dash lines indicate the reference dates for the selected events. The clustered events are labeled. The methods for data filtering, selection of the MJO events, and the MJO classification are described in the Material and Methods section. A zonal three-point running mean was applied to the Hovmöller diagram to remove small-scale noises.

Revealing MJO diversity requires scrutinizing a large number of individual events and treating the MJO as a convectively coupled dynamic system. By tracking intraseasonal anomalies on sequential pentad mean maps of OLR anomalies associated with the MJO convection, Wang and Rui (11) first reveal a variety of distinct propagation tracks that evolve with seasons. Recently, the features of propagating and nonpropagating MJO were compared in some studies (1214). Other studies have noticed that the propagating MJO events exhibit different propagation speed (13, 15). However, a systematic and objective delineation of MJO diversity has not yet emerged in terms of the propagation pattern, intensity, and associated teleconnection and impacts.

In this study, we propose a way to objectively delineate different propagation forms of the observed MJO and to unravel their relationships with MJO circulation and moist thermodynamic structures, as well as to seek for the root causes of the MJO diversity and potential precursors for predicting MJO propagation when the wet phase of the intraseasonal oscillation (ISO) is established in the eastern equatorial Indian Ocean (EIO). The results are expected to advance our understanding of the mechanisms governing complex MJO propagation patterns and to improve MJO simulation and prediction.


Diversity of equatorial MJO propagation and the associated distinct dynamic structures

The K-means cluster analysis is used to objectively classify more than 100 MJO events into four optimal clusters (see Materials and Methods for cluster analysis). The composited propagation patterns for each cluster on Hovmöller diagrams of equatorial OLR anomalies are shown in Fig. 2. The first cluster represents a standing oscillation over the EIO (Fig. 2A). The second cluster shows that an MJO wet phase “jumps” from the EIO to the equatorial western Pacific (EWP) about 5 to 10 days later (Fig. 2B), which can also be considered as a disruption of MJO propagation by the maritime continent (MC) “barrier,” while another convective region emerges over the EWP around 160°E. Both the standing and the jumping clusters are referred to as nonpropagating events, but they are distinguished by the existence of convective anomalies over the EWP. The other two clusters show systematic eastward propagation from the EIO across the MC to EWP, but they differ by their speeds and zonal extents of propagation. The cluster with slower (than average) propagation shows an average speed of 3.5 m/s and the propagation terminates near the dateline (Fig. 2C), whereas the cluster with faster propagation has a mean speed of 5.4 m/s and the propagation extends to ~160°W (Fig. 2D). Note that the canonical MJO obtained by regression of all cases yields a propagation speed (about 4.4 m/s) in between the slow and fast clusters. Thus, the objective cluster analysis unravels four major forms of MJO propagation: standing (24%), jumping (16%), slow propagation (31%), and fast propagation (29%). The number in the parenthesis indicates the percentage of each cluster to the total number of events.

Fig. 2 Four types of MJO propagation patterns along the equator.

Composited longitude-time diagram of intraseasonal OLR anomalies (contour, W m−2) averaged between 10°S and 10°N during boreal winter from November to April for (A) standing, (B) jumping, (C) slow propagation, and (D) fast propagation clusters. The dashed black contour is the −5 W/m2 contour. The shading shows significant signals at 95% confidence level. The black solid lines in (C to D) indicate phase propagation, which is obtained by a least-square fit of the minimum OLR. The phase speed is shown on the top right corner of the panels.

It should be noted that a pure standing ISO cannot be considered as an MJO because the MJO is defined as eastward propagating ISO. By examining the individual events in the standing cluster, we found that some standing events propagate eastward only over the IO and the propagation stops as the convection reaches the MC. These cases may be considered as nonpropagating MJO events discussed in the literature [e.g., (12)]. There are 11 such cases accounting for 50% of the standing cluster cases. Hence, the composited standing cluster shows weak eastward propagation in the IO (Fig. 2A). Given this note, the standing MJO is a loosely defined terminology, acknowledging that some purely standing cases cannot be called MJO events.

To reveal the potential structural differences among the four types of MJO events, we start from diagnosis of the vertical structures of the equatorial MJO circulation anomalies (Fig. 3). Over the EIO between 60°E and 100°E, where the MJO major convection is located, maximum equivalent potential temperature (EPT) anomalies are coupled with a deep updraft, which is a feature common to all clusters. However, there are notable differences between the propagating and nonpropagating clusters. First, the major updraft in the propagating clusters tends to tilt rearward (westward and upward) from the lower to the upper troposphere, but it is nearly upright in the nonpropagating clusters. Second, to the east of the major updraft, the propagating events display notable subsidence, whereas for the nonpropagating clusters, the subsidence is absent (in the standing cluster) or the air is even ascending over the western Pacific (~160°E) (in the jumping cluster). Third, the propagating events feature a deep layer of lower-tropospheric easterlies, whereas in the nonpropagating clusters, the easterlies are either absent (in the jumping cluster) or very weak and shallow (in the standing cluster). In summary, eye-catching–enhanced low-level easterly anomalies are seen to the east of the major convection in the propagating clusters, but they are absent in the nonpropagating clusters.

Fig. 3 Equatorial vertical structures of the four types of MJO.

Composited equatorial (averaged between 10°S and 10°N) vertical structure of circulation (vectors), equivalent potential temperature (EPT, contours in units of K), and precipitation heating (color shading, in units of 1 × 10−2 J kg−1 s−1) for (A) standing, (B) jumping, (C) slow propagation, and (D) fast propagation clusters. Stipples denote where the EPT anomalies are significant at 95% confidence level. The precipitation heating anomalies are only shown for those significant at 95% confidence level. Vectors represent the zonal and vertical velocity (units are m/s for zonal wind and 0.01 Pa s−1 for vertical pressure velocity). Only wind vectors that are significant at 95% confidence level are shown. The vertical lines indicate 100°E. The blue circle indicates the location of maximum upward motion.

There are also notable differences between the slow-propagating and fast-propagating clusters. The subsidence mainly occurs over the MC (120°E to 150°E) for the slow propagating cluster, while it appears over the EWP (150°E to 180°E) for the fast propagating cluster. The deep low-level easterlies blow primarily from 100°E to 160°E in the slow propagating cluster, while they extend from 100°E to the central Pacific (~160°W) in the fast propagating cluster. Thus, the slow and fast propagating clusters display different zonal spans of the leading low-level easterly anomalies, consistent with the theoretical results that the speed of MJO propagation increases with its zonal scale (1618).

The east-west asymmetry in circulation can also be seen from the circulation anomalies at 700 hPa (Fig. 4). The standing cluster features a strong Rossby wave response to the west of the major convection (as evidenced by the Rossby wave cyclonic gyres) and a hardly discernable Kelvin wave response to the east, suggesting that a relatively strong Rossby wave component is a characteristic of standing oscillation. The Kelvin wave is an eastward propagating gravity wave with maximum zonal velocity and geopotential perturbations centered on the equator (19). A heating over the equator could induce a Kelvin wave response to its east and a Rossby wave response to its west (20). In the jumping events, both the Rossby wave westerly and Kelvin wave easterly anomalies are weak. In contrast, in the propagating clusters, the anomalous Kelvin wave easterlies to the east of the major convection are much stronger than those in the nonpropagating events. The fast cluster has also stronger Kelvin wave easterlies, consistent with the results of (16) that a stronger Kelvin wave response produces faster eastward propagation of MJO.

Fig. 4 Horizontal structures of the four types of MJO.

Composited 700-hPa winds (vector, m/s) for (A) standing, (B) jumping, (C) slow propagation, and (D) fast propagation clusters. Only those winds that are statistically significant at the 95% confidence level are plotted. The shading indicates the convective instability index measured by the EPT difference between 850 and 400 hPa (EPT 850 minus EPT 400). The stipples indicate regions where the convective instability indices are significant at 95% confidence level. A nine-point smooth is applied to the convective index field. Blue circles mark the MJO convective center on the equator.

It should be noted that the circulation anomalies to the east of the MJO convection are more complicated than a simple Kelvin wave response, as the circulation is affected by the wave-mean flow interaction (21, 22) and the subsidence over the EWP (12). Since these influences are more significant at the upper troposphere, the circulation anomalies in the lower troposphere can still be largely attributed to the Kelvin wave responses, as the overall configuration of the low pressure and zonal easterly anomalies are largely consistent with the theory.

The thermodynamic structures also show pronounced differences among the four clusters in consistence with the east-west asymmetric circulation patterns. The EPT represents the moist static energy (MSE), a variable that is central to the MJO thermodynamics. The contour in Fig. 3 shows that, in the lower troposphere, positive EPT extends from the maximum center downward and eastward in the propagating clusters, but it tends to be upright in the nonpropagating clusters. The fast propagating events is characterized by a deep layer of MSE extending from the eastern IO to the western Pacific, whereas in the slow propagating events, the deep layer of MSE is confined to the west of the 120°E. The convective instability is closely linked to the deep layer of MSE. The shading in Fig. 4 presents the spatial structure of the lower tropospheric convective instability, which is approximately measured by the differences in the EPT between 850 and 400 hPa (23). A large value in convective instability indicates large potential for the development of congestus convection. Obviously, the propagating clusters are characterized by strong lower tropospheric convective instability located over the MC (105°E to 150°E) for the slow cluster and located further east over the eastern MC-EWP (130°E to 180°E) for the fast cluster. Note that the regions of convective instability in the slow and fast clusters (Fig. 4) coincide with their respective longitudinal locations of enhanced low-level easterlies. In sharp contrast, no significant regions of convectively unstable layer are found in the standing and jumping clusters.

The unstable thermodynamic structure to the east of the major convection favors the development of shallow and congestus clouds. In observations, we use precipitation-induced condensational heating as a surrogate to deduce the corresponding cloud distributions. The MJO precipitation heating is derived from the moisture sink term in the study of Yanai et al. (24). The shadings in Fig. 3 show the condensational heating rate. To the east of the EIO (100°E), the heating in the propagating clusters extends from the major convective center eastward and downward, suggesting the existence of congestus and shallow clouds ahead of the MJO deep convection. The lower tropospheric heating, particularly the congestus cloud heating around 700 hPa over the MC, is considerably stronger in the fast events than in the slow events. On the other hand, the low-level heating to the east of the major convection is absent in the standing cluster and relatively weak in the jumping cluster. Furthermore, the propagating clusters are characterized by a strong midtropospheric cooling over the MC (100°E to 150°E) for the slow propagation and over the EWP (150°E to 180°E) for the fast propagation. In sharp contrast, no midtropospheric cooling is found in the nonpropagating clusters. Rather, there is a warming in the jumping cluster around 160°E to 170°E, signaling a subsequent re-emergence of convective anomaly in the EWP. In summary, both the circulation and thermodynamic structures exhibit distinctive characteristics among the four clusters of MJO propagation, suggesting that the MJO system is a propagation-structure nexus.

Mechanisms of the diversified MJO propagation

A key feature distinguishing propagating and nonpropagating events is that the propagating clusters are characterized by strong Kelvin wave easterly and its tight coupling to the major convection, while the nonpropagating clusters are not (Figs. 3 and 4). We argue that a strong Kelvin wave easterly and its coupling to the major convection could help the MJO eastward propagation. The enhanced low-level equatorial Kelvin wave response can reinforce the boundary layer (BL) moisture convergence (2527) because the equatorial easterly is accompanied by an equatorial low pressure, and in the equatorial region, the BL moisture convergence (BLMC) occurs primarily in a low-pressure area (28). It is the BL convergence that can effectively accumulate moisture in the lower troposphere (14, 23), thus increasing the low-level MSE or EPT (Fig. 3) and destabilizing the lower troposphere (Fig. 4). Note that the horizontal advection of mean moisture can also contribute to this low-level moistening (29), but its amplitude is weaker than the vertical moisture advection associated with the BL convergence near the equator (14) where the MJO’s convection is strongest. In addition to the aforementioned effects, the descending motion to the east of the major convection, which is partly attributed to the Kelvin wave response, can deplete water vapor in the midtroposphere, thus decreasing mid-level MSE and further destabilizing the lower troposphere. As a result, in the propagating clusters, strong convective instability develops to the east of EIO convection (Fig. 4), which promotes the development of congestus convection and lower tropospheric heating (Fig. 3). The congestus convection further preconditions the midtroposphere, promotes development of deep convection in a later time, and leads to the eastward propagation of the MJO.

Observed anomalous surface low pressure and BLMC precede the MJO major convective center (2, 27, 30). To the east of the MJO convective center occurs gradual deepening of the moist BL (31, 32), increased convective instability (23), and a “stepwise” transition from shallow cumulus, congestus clouds to deep convection and anvil stratiform clouds (3335). These structural features are consistent with the enhanced Kelvin wave response and its coupling to the major convection, which can promote eastward propagation.

What are the root causes of the differences of the Kelvin wave response among different clusters? Figure 5 shows the background sea surface temperature (SST) anomalies associated with the four clusters. The most notable feature is that the standing cluster is associated with a significant sea surface cooling over the central and eastern Pacific, reminiscent of a La Niña state, while the fast cluster is associated with a significant warming over the central and eastern Pacific, reminiscent of an El Niño state. The slow cluster is associated with a weak central Pacific cooling, while the jumping cluster is associated with a weak central Pacific warming, but these background SST anomalies are not statistically significant.

Fig. 5 Background SST anomalies associated with the four MJO clusters.

The background SST anomalies are obtained by compositing the 3-month average of monthly SST anomalies (derived from climatological monthly mean) for each MJO cluster, with the central month being the month that contains day 0 of an MJO event.

When the central and eastern Pacific warms (an El Niño condition), the Indo-Pacific warm pool and the associated active convective region expand eastward, which leads to an increased zonal scale of the MJO, as the MJO convective activity is observed over the warm ocean (6) and it tends to amplify on a longer spatial scale (1618). As revealed by the observed wave number–frequency power spectral of MJO signals (36) and theoretical studies (1618), the MJO propagation speed is proportional to its zonal scale. This explains why the fast cluster has a longer zonal scale in circulation and a faster eastward propagation speed than the slow cluster (Fig. 4). In addition, the BL convergence feedback favors enhancing the Kelvin wave response on longer scale (16), leading to a stronger Kelvin wave response on the longer scale.

On the other hand, when the central and eastern Pacific cools (a La Niña condition), the zonal scale of the Indo-Pacific warm pool shrinks, leading to a decrease in the MJO zonal scale. The reduced scale of the convectively coupled system does not favor the development and growth of the Kelvin wave response to the east of the major convection, leading to a weak Kelvin wave response in the standing cluster, which tends to decouple with the major convection (Figs. 3A and 4A).


We have shown that the observed MJO propagation patterns can be objectively delineated into four archetypes of clusters: standing, jumping, slow eastward propagation, and fast eastward propagation. Each type of MJO propagation exhibits distinctive circulation and thermodynamic structures. The Kelvin wave response to the MJO major convection and the associated features of the east-west structural asymmetry provide a set of indicators for predicting the ensuing propagation when MJO major convection center is established in the EIO. A prominent forerunner that distinguishes propagating from nonpropagating events is the strength of the Kelvin wave response and its coupling to the east of the MJO main convection: The Kelvin wave response is strong and tightly coupled with the MJO major convection in the propagating groups while it is very weak and tends to decouple from the major convection in the nonpropagating group. The strength and length scale of the Kelvin wave response also distinguish the slow and fast propagation groups: The fast group has stronger and longer low-level Kelvin wave response. It is further shown that the background SST condition affects the Kelvin wave response: An El Niño state tends to increase the zonal scale of Kelvin wave response, to amplify it, and to enhance its coupling to the convection, while a La Niña state tends to decrease the zonal scale of Kelvin wave response, to suppress it, and to weaken its coupling to the major convection.

The variation of background SST over the Pacific is shown to potentially affect the propagation of the MJO. It is also conceivable that the seasonal variation of background states could significantly affect the propagation patterns of the MJO. For example, the basic SST and convection are more aligned with the equator in November and April, while the basic SST and convection is more off the equator during December to March. A thorough investigation of this issue should be carried on in the future studies.

The knowledge obtained here could yield insights as to why most models do not adequately simulate the MJO and why many of the models exhibit too slow propagation. The observed evidence presented here helps in understanding of the diversified MJO propagation mechanism. It would be of interest to find out how the current GCMs capture the MJO propagation diversity. The observed MJO diversity offers new metrics for validating GCMs and for projecting MJO’s future changes under global warming scenarios.


Data and filtering process

We used daily averaged OLR data from National Centers for Environmental Prediction /National Oceanic and Atmospheric Administration–interpolated OLR dataset (37) for the period of 1979–2013 (35 years) to depict the major deep convection associated with the MJO. For the horizontal and vertical winds, temperature, and specific humidity, the daily averaged European Centre for Medium-Range Weather Forecasts Interim reanalysis data (38) from 1000 to 100 hPa with 50-hPa interval were used. To study the background SST conditions, the SST dataset from the NOAA Extended Reconstructed Sea Surface Temperature (ERSST) version 5 (39) has been adopted. For OLR and atmospheric data, intraseasonal (20- to 70-day) band-pass–filtered anomalies during boreal winter (from November to April) were analyzed. The anomalous fields were obtained by first removing the time mean and climatological annual cycle based on the 1980–2010 reference period and then applying a 20- to 70-day band-pass filter. The OLR and atmospheric data were interpolated to 2.5° × 2.5° grids to match each other. Here, we defined observations to be NOAA OLR, NOAA ERSST, and ERA-Interim combined.

Cluster analysis

A popular and straightforward way of depicting MJO eastward propagation along the equator is to display a lag-longitude correlation map (also known as Hovmöller diagram) of OLR anomalies (8). The OLR propagation along the equatorial belt exhibits various forms, ranging from standing oscillation to systematic propagations with different speeds (Fig. 1). Such a variety of propagation pattern with eye-catching irregularities are difficult to be grouped by subjective analysis without ambiguity. Here, the K-means cluster analysis (40) was used to objectively classify diverse propagation patterns.

The first step was to identify individual MJO events. For this purpose, we focused on the EIO region where the MJO convective variability center is located and the climatological Intertropical Convergence Zone is closest to the equator. During boreal winter from November to April, an MJO event was selected if the box-averaged OLR over the EIO (10°S to 10°N and 75°E to 95°E) was below 1 SD for five successive days. There are a total of 103 selected cases during the 35-year period from 1979 to 2013. The reference date for a selected MJO event was chosen as the day when the EIO (10°S to 10°N and 75°E to 95°E) averaged OLR reaches its minimum.

The second step was to apply cluster analysis to the time-longitude (Hovmöller) propagation patterns of the OLR anomalies averaged along an equatorial belt between 10°S and 10°N for all selected events. The reference date was denoted as day 0. For each qualified MJO event, the cluster analysis domain in the Hovmöller diagram covers a 30-day period from days −10 to 20 and a longitudinal extent from 60°E to 180°E. A zonal three-point running mean was applied to the Hovmöller diagram to remove small-scale noises. To focus on convective anomalies, the cluster analysis focused only on the wet phase of the MJO, i.e., the region where OLR anomalies <−5 W/m2; the OLR value in the region where OLR >−5 W/m2 is set to zero in the cluster analysis.

In cluster analysis, correlation was used to measure the “similarity” between each cluster member and the corresponding cluster centroid. We used the silhouette clustering evaluation criterion to assess the skill of the cluster analysis. The silhouette value for each member ranges from −1 to +1 and indicates how similar a member is to its own cluster compared to other clusters. A high silhouette value indicates that the member is well matched to its own cluster and poorly matched to neighboring clusters (40). The cluster analysis of the 103 MJO events shows that these events can be optimally fitted into four clusters. In our analysis, the members having a silhouette value lower than 0.06 were excluded from the corresponding clusters as they were considered as sufficiently “dissimilar” to the corresponding centroid of the cluster. By removing 13 “outliers,” 90 events remain in the four clusters, which were used to perform composite analysis.

MJO diagnosis

To diagnose the MJO structure, circulation and thermodynamic fields were depicted by compositing 20- to 70-day filtered anomalies of horizontal winds, vertical velocity, EPT, and diabatic heating fields with regard to the phases at the reference location in the EIO (i.e., compositing pentad mean of the selected MJO events with central date being the reference date). This location was chosen because, in the EIO, the mean state tends to be more equatorially symmetric, the corresponding MJO precipitation and circulation are also more symmetric at the equator, and thus it is more suitable for the study of the structure of the MJO with less influence from the mean flows. Because no observed diabatic heating field was available, we estimated it using the budget analysis proposed by Yanai et al. (24). The condensational heating was estimated by the moisture sink defined asQ2=Lv(qt+Vq+ωqp)where Lv is the latent heat at condensation, q is the specific humidity, V refers to the horizontal components of the wind, and ω is the vertical velocity. This quantity was first calculated by daily data, and time filter was applied to the variable to retrieve the intraseasonal component.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We acknowledge the useful discussion with A. Sobel on an early version of the manuscript. Funding: This work is jointly supported by the NSF/Climate Dynamics award no. AGS-1540783, the China National 973 Project (2015CB453200), the National Natural Science Foundation of China (grant no. 41420104002), and Public Science and Technology Research Funds Projects of Ocean (grant no. 201505013). This is the SEOST publication 10705, IPRC publication 1382, and ESMC publication 259. Author contributions: B.W. and G.C. initiated this research and wrote the paper. G.C. made the data analysis and graphics. F.L. contributes to discussion and comments on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article