## Abstract

Recent studies have suggested that the Madden–Julian oscillation is a result of an instability driven mainly by cloud–radiation feedbacks, similar in character to self-aggregation of convection in nonrotating, cloud-permitting simulations of radiative–convective equilibrium (RCE). Here we bolster that inference by simulating radiative–convective equilibrium states on a rotating sphere with constant sea surface temperature, using the cloud-permitting System for Atmospheric Modeling (SAM) with 20-km grid spacing and extending to walls at 46° latitude in each hemisphere. Mechanism-denial experiments reveal that cloud–radiation interaction is the quintessential driving mechanism of the simulated MJO-like disturbances, but wind-induced surface heat exchange (WISHE) feedbacks are the primary driver of its eastward propagation. WISHE may also explain the faster Kelvin-like modes in the simulations. These conclusions are supported by a linear stability analysis of RCE states on an equatorial beta plane.

## 1. Introduction

Since its discovery in the early 1970s (Madden and Julian 1971), the MJO has challenged our basic understanding of tropical atmospheric dynamics and thermodynamics [see the review by Zhang (2005)]. Early on, it was noticed that the structure and propagation of the MJO in some ways resemble those of a classical equatorially trapped Kelvin wave (Matsuno 1966), but the observed wave travels at a small fraction of the predicted Kelvin wave speed, given its first baroclinic-mode-like vertical structure. This observation confirmed Madden and Julian’s supposition that coupling to deep convection is fundamental to the oscillation. Initial formulations used the wave–conditional instability of the second kind (CISK) formalism (e.g., Lindzen 1974) to account for latent heat release, but these also produced Kelvin waves that were too fast and suffered the ultraviolet catastrophe of largest growth at the smallest scales. This defect was to some extent remedied by tying latent heat release specifically to frictionally induced convergence (Moskowitz and Bretherton 2000), but the CISK formalism depends on ambient conditional instability, which is small in the tropics (Betts 1982; Xu and Emanuel 1989).

Thus, while convective coupling might help explain the discrepancy between observed and dry Kelvin-mode propagation speeds, it could not explain the energy source of the waves. This led Neelin et al. (1987) and Emanuel (1987) to propose that intraseasonal variability was driven by a feedback between surface winds and surface enthalpy fluxes, a feedback known variously as evaporation-–wind feedback and wind-induced surface heat exchange (WISHE). Simple linear models incorporating this process produce eastward-propagating modes that have largest growth rates at small scales unless coupling to the stratosphere is included (Yano and Emanuel 1991) or small lags in convective response are accounted for (Emanuel 1993). These more advanced theories produced a promising explanation for the higher-frequency equatorially trapped modes beautifully revealed by the wavenumber–frequency decomposition of satellite-observed outgoing longwave radiation (OLR) by Wheeler and Kiladis (1999), but they did not capture the low, nearly constant frequency of the observed MJO.

The strong modulation of OLR by the MJO suggests that the interaction between radiation and clouds and water vapor may be a key process driving the MJO, as first suggested by Hu and Randall (1994). By the late 1990s, a somewhat separate research track developed to understand the properties of radiative–convective equilibrium states using three-dimensional, doubly periodic, cloud-permitting models run under homogeneous boundary conditions. In the first of these, Tompkins and Craig (1998) showed that the initially random convection spontaneously organized into banded structures, driven by cloud–radiation interactions and enhanced surface fluxes associated with convection. The essential role of WISHE and cloud–radiation interaction in this self-aggregation of convection was reinforced by the work of Bretherton et al. (2005), who likewise used a cloud-permitting, doubly periodic model with constant sea surface temperature, though the aggregation took the form of clusters rather than bands. Experiments without one or the other (or both) of WISHE and cloud–radiation interactions failed to self-aggregate. Detailed analyses of the physics of self-aggregation by, for example, Wing and Emanuel (2014), left little doubt that interactions between both longwave and to a lesser extent shortwave radiation with clouds and water vapor are essential to the instability that leads to self-aggregation, that longwave–cloud interactions are the main energy source for mature aggregates, and that wind-dependent surface fluxes have a strong influence on the process as well.

The observation of self-aggregation in nonrotating models leads naturally to the idea that the MJO is a result of self-aggregation in the equatorial waveguide of a rotating planet. Using a simplified numerical model with idealized boundary conditions, Raymond (2001) showed that cloud–radiation interactions could destabilize the radiative–convective equilibrium (RCE) state and lead to MJO-like disturbances that propagate eastward at speeds close to the observed propagation speed of about 5 m s^{−1}. The physics of such modes were explored in simple linear model frameworks by, for example, Fuchs and Raymond (2005), Bony and Emanuel (2005), Zurovac-Jevtić et al. (2006), and Fuchs and Raymond (2017). All of these studies showed that cloud–radiation interactions and WISHE are important elements of the MJO-like modes that emerged. In particular, Fuchs and Raymond (2005) concluded that the main driver of the MJO is cloud–radiation interaction, while WISHE is responsible for its eastward propagation. On the other hand, their updated analysis (Fuchs and Raymond 2017) identified WISHE as the main mechanism for driving the MJO, though cloud–radiation interactions further destabilized the mode.

An important step forward was taken recently by Arnold and Randall (2015), who analyzed simulations on an aquaplanet with uniform sea surface temperature using the superparameterized (SP) Community Atmosphere Model [CAM (SP-CAM)]. This model uses a superparameterization of moist convection wherein convection is explicitly simulated in two-dimensional domains embedded within each three-dimensional grid box. Convection aggregates in this model whether or not the planet rotates. In the rotating case, the aggregation takes the form of an MJO with properties very similar to those of the observed MJO. As with the earlier studies, the aggregation disappears in the nonrotating case when the longwave radiation is horizontally homogenized, but there is still a weak MJO-like mode on the rotating planet. As found by Fuchs and Raymond (2005), the rotating MJO mode is driven eastward by WISHE.

We here extend the work of Arnold and Randall (2015) to a model in which the convection is explicitly simulated everywhere rather than being superparameterized. To make this feasible, we have to use a comparatively low horizontal resolution and limit our domain in the north–south direction. These compromises allow us to conduct mechanism-denial experiments to cast some light on the physical mechanisms.

## 2. The model and case setup

Detailed information about the System for Atmospheric Modeling, version 6.10.6 (SAM), is given by Khairoutdinov and Randall (2003). SAM is a nonhydrostatic cloud-resolving model (CRM) that solves an anelastic approximation to the momentum equations, which filters out sound waves. There is a choice of several types of single-moment and double-moment microphysics. In this study, we use SAM’s original computationally efficient single-moment microphysics, which is based on two prognostic water variables, the mixing ratio of total nonprecipitating (sum of water vapor, cloud water, and ice) and precipitating (sum of rain, graupel, and snow) water. In addition to two prognostic water variables, the model uses the liquid water/ice static energy as a prognostic variable. The cloud and precipitating-water categories are simply diagnosed from the prognostic variables assuming prescribed partitioning functions depending on temperature only. The momentum equations are integrated in time using the third-order Adams–Bashforth scheme and in space using the second-order centered differences in flux form and with enforced conservation of the kinetic energy. All scalars are advected using the second-order positive-definite and monotonic MPDATA algorithm. The subgrid-scale fluxes are evaluated using a Smagorinsky-type prognostic closure. The interactive radiation code is taken from NCAR’s CAM3 climate model.

The computational setup of this study follows the so-called near-global modeling approach (Bretherton and Khairoutdinov 2015; Narenpitak et al. 2017). The grid is a staggered Arakawa C type, with uniform horizontal and variable vertical grid spacing. The computational domain extends 10 200 km in the north–south (N–S) direction, equivalent to latitudes extending to about 46°N/S, where solid-wall lateral boundaries are placed. (We will henceforth refer to north–south distances as “latitude,” which should be understood as distance from the equator, and east–west (E–W) distances as “longitude,” equivalent to degrees along the equator.) In the zonal direction, the domain is periodic and 40 000-km long, which is equal to Earth’s circumference at the equator. Because of the high computational cost of this study, almost all simulations have been performed using rather coarse “cloud permitting” 20-km grid spacing, except for one relatively short simulation, which used a “cloud resolving” 4-km grid spacing. The vertical grid has 38 levels with the top at 34 km. The grid spacing gradually increases from about 70 m near the surface to much coarser spacing of up to 1000 m near the tropopause and up to 1500 m above 22 km. The lowest grid level is at 37 m. There is a Newtonian damping layer in the upper third of the domain to minimize the reflection of gravity waves from the domain top. The time step is 20 s (10 s for 4-km run). The insolation at the model top is constant (no diurnal cycle) with a solar constant of 651 W m^{−2} and a fixed zenith angle of 50.5° (following Tompkins and Craig 1998), while the surface enthalpy and momentum fluxes are interactive. Unlike the aforementioned near-global studies, which used a meridionally varying, but zonally uniform sea surface temperature (SST), the SST in this study is set to 300 K everywhere. The Coriolis parameter depends realistically on the latitude (not a beta plane) and is zero at the N–S center point of the domain. Arguably, the framework of this study is the simplest possible one to study equatorially trapped disturbances.

The initial temperature and humidity profiles are horizontally uniform and are obtained from a small-domain nonrotating RCE simulation with the same SST, resolution, and physics settings as the large-domain run. Initially, there is no wind. To initiate the convection, small-amplitude (about 0.5 K) random noise is added everywhere to the temperature field at a few levels near the surface. The simulations are run for 260–280 days.

## 3. Results

### a. Control case

Figure 1 shows the evolution of the domain-mean precipitation and precipitable water (column-integrated water vapor). For precipitable water, the evolutions of the domain minimum and maximum daily values are also shown. It takes about 80–100 days to reach a well-developed, statistically steady state. An interesting feature of the simulation is the rapid development, over an initial period of about 20 days, of very dry regions, which is consistent with both small-domain CRMs (Wing and Emanuel 2014) and superparameterized GCMs’ (Arnold and Randall 2015) RCE results. The dry patches are illustrated by snapshots of precipitable water at both the beginning and the end of the simulation (Fig. 2). Tropical cyclones develop outside of the narrow equatorial belt and then migrate toward higher latitudes to their eventual demise near the domain’s northern and southern walls.

Despite the constant and spatially uniform SST and insolation, a mean meridional circulation still develops. It is characterized by a pair of weak Hadley cell–like circulations in the vicinity of the equator and Ferrel cell–like circulations further away from the equator, as illustrated by latitude–pressure plots of zonal-mean zonal and mass streamfunction in Fig. 3. The tropical tropospheric easterlies extend to about 20°N/S near the surface and farther away from the equator in the upper troposphere, surrounded by a pair of “extratropical” westerlies with jets near the surface. Such a large-scale circulation pattern is consistent with the GCM simulations of rotating RCE over constant SST and constant insolation by Shi and Bretherton (2014). The zonal-mean precipitation (Fig. 3) has a double-ITCZ structure with a minimum of about 3.5 mm day^{−1} at the equator and surrounding maxima of about 4.5 mm day^{−1} at about 7°N/S and minima near 20°N/S. There are also two maxima of about 6.5 mm day^{−1} right near the walls at 40°N/S, which are probably related to the tropical cyclones that develop in these regions. The meridional distribution of zonal-mean precipitable water (Fig. 3) qualitatively follows the pattern of precipitation, with the equator being the driest place, where mean precipitable water is about 31 kg m^{−2}, surrounded by maxima of about 34 kg m^{−2} at 5°–10°N/S.

Variability in the narrow equatorial belt (10°S–10°N) is illustrated by the Hovmöller diagram of precipitable water anomalies (Fig. 4). One can see an eastward-propagating, planetary-scale, that is, zonal wavenumber-1, MJO-like disturbance (Fig. 4a), which is even better revealed when all the modes with scales shorter than zonal wavenumber-1 mode are filtered away (Fig. 4b). The period of the MJO-like disturbance is about 50–55 days, corresponding to the speed of about 8–9 m s^{−1}. There are also westward-propagating small-scale disturbances that could be the equatorially trapped Rossby waves and also disturbances advected by the equatorial easterlies.

The observed MJO propagates somewhat more slowly over the Indo-Pacific warm pool (about 5 m s^{−1}; e.g., Zhang 2005) but then accelerates over cooler waters, having an overall period similar to those in our simulations. We leave an analysis of this disparity in phase speed to future work but speculate that the coarse resolution of the model may lead to a gross moist stability that is too large and to other problems.

Wavenumber–frequency symmetric and antisymmetric signal-to-background variance spectra of OLR in the narrow equatorial belt, following the methodology of Wheeler and Kiladis (1999), are shown in Fig. 5. The symmetric spectrum (Fig. 5a) is dominated by a zonal wavenumber-1 signal with a period longer than 40 days, which is consistent with an MJO-like disturbance. There is also a relatively strong Kelvin wave signal, which seems to be spectrally separate from the MJO-like signal, as often seen in observed OLR spectra. There is also some indication of the presence of equatorial Rossby waves and high-frequency inertia–gravity waves although their signals are rather weak. The antisymmetric spectrum (Fig. 5b) is dominated by mixed Rossby–gravity waves with periods shorter than 6 days. Overall, the spectra look quite reasonable, resembling the observed spectra of equatorially trapped disturbances even though the zonally symmetric circulation and distribution of moisture and precipitation are quite different from observed.

Figure 6 shows the time evolution of fractional growth rates of normalized spatial variance of the vertically integrated moist static energy (iMSE) for the MJO-like disturbance averaged over the 10°N–10°S equatorial belt computed using the methodology described by Wing and Emanuel (2014). The growth rate is dominated by the positive covariance between the net longwave radiative heating and the iMSE and mostly opposed by the negative covariance with advective convergence of iMSE. Note that the convergence of iMSE itself has been calculated as a residual of the moist static energy (MSE) budget. The contribution of the shortwave radiative heating appears negligible. The contribution of the covariance of the surface enthalpy flux to the iMSE variance growth appears to be rather insignificant with the exception of the first 10 days, when the surface fluxes tend to be nearly as important as the contribution of the longwave heating to the initial growth of the MJO-like disturbance. The dominance of the longwave heating in the iMSE variance maintenance is broadly consistent with the self-aggregation study by Wing and Emanuel (2014) and the results of superparameterized GCM simulations by Arnold and Randall (2015). The role of the advection appears to predominantly damp the iMSE variance, as appears to be true in nature (Kiranmayi and Maloney 2011) and in other numerical simulations (e.g., Arnold and Randall 2015).

Composites of the MJO-like disturbance are created by tracking the minimum of the surface-pressure anomaly in the 10°N–10°S equatorial belt of the eastward-propagating zonal wavenumber-1 signal with a period of 40–80 days, starting from day 80. The surface pressure anomaly is chosen as the tracking anomaly because pressure is one of the least noisy fields. The time series of the longitude of the pressure minimum is then used to reposition the other fields of interest so that they are all centered for each time sample at a chosen central longitude. In addition, to reduce the noise associated with shorter spatial and temporal scales, the composite fields are filtered to retain only zonal wavenumbers from 1 to 4 and periods from 30 to 80 days before averaging in time.

The resultant anomaly composites of the MJO-like disturbance for several fields are shown in Fig. 7. The precipitation and iMSE anomalies show a characteristic “swallowtail” pattern often seen in the composites of the observed MJO (e.g., Zhang and Ling 2012), when the anomalies extend eastward from the center along the equator and split into two off-equatorial branches west of the center. The circulation around the core of the MJO-like disturbance shows a strong quadrupole pattern of vorticity anomalies, with a pair of low-level cyclonic gyres to the west and a pair of anticyclonic gyres to the east of the core. The vorticity patterns reverse sign in the upper troposphere. Such a quadrupole gyre structure is a well-known characteristic of the observed MJO.

It is instructive to look at the phase relationships and amplitudes of the various contributors to the budget of vertically integrated moist static energy (here denoted as ). In a coordinate system moving with the MJO, in which we assume the fields are steady, this budget may be written

where the angle brackets denote vertical integrals through the domain, in pressure coordinates; is the surface turbulent enthalpy flux; and and are the total upward radiative fluxes (longwave and shortwave) at the surface and top of the troposphere, respectively. Because of high-frequency variations in velocities and moist static energy, computing the terms on the left side of (1) can be problematic, and here we calculate them as a budget residual [i.e., the sum of the terms on the right side of (1)]. This does not permit the separation of the advective terms into vertical and horizontal components, though. We next transform (1) into an Earth-fixed coordinate frame by adding in the effect of the coordinate eastward translation at speed *c*, which will be compensated by a local time derivative in the Earth-fixed frame:

Since the term involving *c* is linear, it is easy to calculate from model output. The quantity , where is the unit vector in the *x* direction, is just the ground-relative horizontal velocity. Remember that the nonlinear terms in (2) are calculated as a budget residual in the MJO-relative coordinate frame and that **V** is the MJO-relative horizontal velocity. We will henceforth refer to the sum of first three terms on the right of (2) as the “advective terms.”

The composite net radiative heating and the surface enthalpy flux as well as the advective contribution to the vertically integrated moist static energy budget are shown in Fig. 8. As expected, the radiative heating is collocated with the column-integrated iMSE anomaly, contributing to the maintenance of the MJO-like disturbance. On the other hand, the positive surface flux anomaly is shifted one-quarter-wave zonal distance to the east of the core, which optimally facilitates eastward propagation of the MJO-like disturbance. The positive surface flux anomaly is collocated with the enhanced easterly surface wind anomaly. At the same time, the negative surface flux anomaly west of the core correlates well with the enhanced westerly surface wind. These anomalies, when added to the background equatorial easterlies, can explain the local surface enthalpy flux anomalies as driven predominantly by WISHE. In this simulation, the advective terms mostly act to damp the moist static energy perturbations, consistent with the moist static energy budget calculated from reanalysis data (Kiranmayi and Maloney 2011). On the other hand, the combination of the observed radiative and surface fluxes in Kiranmayi and Maloney (2011) tends to counter the eastward propagation of the MJO, contrary to what is seen in our simulations. Note that the dominant role of the surface fluxes in the eastward propagation of the MJO-like disturbance in our simulation, while consistent with the theory by Emanuel (1987), goes against the observed dominant role of transport of water vapor in the eastward propagation of the MJO (Maloney 2009; Pritchard and Bretherton 2014; Sobel et al. 2014; Wang et al. 2015). This disparity may be the result of the weak meridional gradient of water vapor in this study because of the absence of a meridional gradient of SST.

### b. Mechanism-denial experiments

We now present the results of a few mechanism-denial experiments that may help to shed some light on the importance of surface enthalpy fluxes and radiative heating for the development, maintenance, and propagation of the MJO-like disturbance. We start with the surface enthalpy fluxes. As mentioned in the previous section, although these fluxes seem to play a relatively minor role in maintaining the mature MJO-like disturbance in our simulation, they do seem to facilitate its eastward propagation.

In the first mechanism-denial experiment, HOM-SFC, the surface latent and sensible heat fluxes are zonally homogenized. Specifically, fluxes are first computed the usual way, that is, taking into account local values of surface wind and surface-to-air differences of temperature and water vapor, but then they are zonally averaged before being applied at each time step. One can see (Fig. 9b) that in this case the zonal wavenumber-1 anomaly still develops, with even larger amplitude than in the control (CTRL) case, while taking about the same time to aggregate. However, there is no coherent eastward propagation of the resulting disturbance, as there is no surface flux anomaly now, illustrating the importance of variable surface fluxes for the propagation of the disturbance. Homogenization of the surface fluxes also seems to have a noticeable effect on Kelvin waves, making them faster (Fig. 10b), and the mixed Rossby–gravity waves are slightly reduced in amplitude (Fig. 10g).

The surface flux of moist enthalpy is the product of two factors: surface wind speed and the difference or disequilibrium between the surface and air enthalpies. Homogenization of the enthalpy flux in the HOM-SFC cannot reveal the roles that these factors separately may play in the evolution and propagation of the MJO-like disturbance. For a given surface wind speed, the surface enthalpy disequilibrium should have a damping effect on MSE anomalies, as found in previous studies of self-aggregation (e.g., Andersen and Kuang 2012; Wing and Emanuel 2014; Arnold and Randall 2015) as most water vapor is concentrated in the boundary layer. This notion is supported by the next mechanism-denial experiment, NO-WISHE, in which only the surface wind speed is zonally homogenized in latent and sensible heat flux computations. The result of such a homogenization is virtually complete suppression of the zonal wavenumber-1 disturbance, which fails to develop (Fig. 9c) as the result of damping of the corresponding moisture and temperature anomalies before they have a chance to amplify. The damping effect of surface enthalpy disequilibrium is, in fact, so strong that not only does it suppress the initial development of the MJO-like disturbance, but it also suppresses the already developed disturbance in CTRL. This is demonstrated by Fig. 11, which shows the result of homogenization of the surface wind applied to the CTRL after day 280.

It is interesting that in NO-WISHE, unlike the HOM-SFC, the speed and, to a certain degree, the magnitude of Kelvin waves (Fig. 10c) do not seem to differ from the CTRL.

In the next mechanism-denial experiment, WISHE-ONLY, we eliminate the effect of thermodynamic disequilibrium on the surface flux anomalies, that is, only the water vapor and temperature at the first model layer are zonally homogenized in the surface flux calculations, while the local effect of the surface wind on the flux remains; thus, the effect of surface fluxes on the iMSE anomalies is now purely because of WISHE. The result is a strong zonal wavenumber-1 anomaly developing even faster than in CTRL (Fig. 9d), which is consistent with eliminating anomaly damping by surface disequilibrium. Unlike the HOM-SFC, though, the anomalies in WISHE-ONLY do propagate, albeit slower and not as coherently as in CTRL, mostly to the east, with some periods of hesitation and even reversal to the western direction (Fig. 9d). Thus, indeed, WISHE appears to be essential for development and propagation of the MJO-like disturbance, although it seems that some degree of added damping of the anomaly by the surface enthalpy disequilibrium is somehow also beneficial for its coherent eastward propagation.

The radiation feedback is the crucial mechanism for development and maintenance of the MJO-like disturbance in CTRL; therefore, the next logical step is to suppress the effect of radiation on the iMSE anomalies by homogenizing radiative heating, which is done in the HOM-RAD experiment. One would expect, based on dominance of radiation feedback in the development and maintenance of the MJO-like disturbance in CTRL, that the zonal wavenumber-1 mode should not then develop. Perhaps surprisingly, the mode still develops in HOM-RAD (Fig. 9e), although the variance anomaly is not as large as in the CTRL case. The corresponding frequency–wavenumber spectrum (Fig. 10e) suggests that the zonal wavenumber-1 disturbance has a somewhat faster speed (higher frequency) than the MJO-like disturbance in CTRL, with the propagation period decreasing from 45–55 to 30–35 days. Such a period puts the disturbance virtually right on the dispersion line for the Kelvin waves. Also, there seems to be no spectral gap now between the zonal wavenumber-1 disturbance and Kelvin waves, which was a feature of the CTRL. Homogenizing the radiation has virtually no effect on the asymmetric modes (Fig. 10j).

The composite precipitation anomaly in HOM_RAD (Fig. 12a) does not have as clear a swallowtail shape as the pattern seen in CTRL (cf. Fig. 7a) and that is a characteristic of the observed MJO. Also, the maximum of the surface enthalpy flux anomaly is collocated with the iMSE anomaly maximum (not shown), and, therefore, cannot be effective in moving the disturbance eastward, though it is responsible for its amplification and maintenance. As the radiative heating is horizontally homogenized, there is no other diabatic forcing to affect the iMSE anomaly and to facilitate its eastward propagation. Therefore, the zonal wavenumber-1 disturbance in RAD-HOM case propagates to the east as any Kelvin wave, that is, simply by the gravity wave mechanism.

One possible drawback of these mechanism-denial experiments is that they also change the time-mean state in ways that may affect the properties of the variability, as shown by Ma and Kuang (2016). They also showed that this effect can be mitigated by driving the model’s time-mean state back toward that of the control run. Given the expense of the current simulations, this was not deemed feasible at present but should be implemented in future work.

## 4. Linear stability analysis

We next interpret the results of our numerical simulations using linear theory. In some ways, the development here is similar to that of Fuchs and Raymond (2005) and Fuchs and Raymond (2017), who interpreted the MJO as a moisture mode whose eastward propagation is driven by WISHE, consistent with our present numerical results. We use a slightly different framework, but having many of the same key ingredients, including wind-dependent surface fluxes, small but positive effective static stability, and a representation of cloud–radiation interaction.

One of our key underlying assumptions is that deep convection always maintains a moist adiabatic lapse rate in the free troposphere, and as in previous work, this is here represented by constant saturation moist entropy *s** is the free troposphere. We also assume that the large-scale vertical velocity vanishes at the surface and at the tropopause. [This latter assumption does not work very well in the tropics, where upward radiation of wave energy is an important sink for tropospheric wave energy, particularly at higher frequencies (Yano and Emanuel 1991).] If these assumptions are satisfied, and the flow is hydrostatic, then Emanuel (1987) showed that the vertical structure of tropospheric disturbances is constrained to be the first baroclinic mode, as also demonstrated by Neelin and Yu (1994). Geopotential fluctuations at the surface are related to saturation entropy fluctuations of the free troposphere by

where is the fluctuation of the near-surface geopotential, is the mean absolute temperature at the top of the subcloud layer, and is the pressure-weighted vertical mean temperature of the free troposphere. This reduces the full horizontal momentum equations, linearized about a resting basic state, to a shallow-water form (Emanuel 1987):

where and are the fluctuations of zonal and meridional velocity at the surface, is the meridional gradient of the Coriolis parameter at the equator, and we have phrased the equations on an equatorial beta plane. This is complemented by a simple Boussinesq mass continuity equation:

where is the fluctuation of large-scale vertical velocity in the midtroposphere and *H* is a representative half depth of the troposphere.

The linearized thermodynamic equation for is written as

where and are the dry and moist adiabatic lapse rates, is the perturbation radiative heating, is the basic-state vertical gradient of dry entropy (i.e., the static stability), is the fluctuation in the convective updraft mass flux, and is the bulk precipitation efficiency. The factor converts fluctuations of dry entropy to those of saturation entropy. The term in (7) involving is just the adiabatic cooling due to large-scale vertical velocity, while the term involving is the convective heating. If the precipitation efficiency is zero, there is no net heating even if there is convection, because there is no vertically integrated latent heat release.

The convective updraft mass flux is determined through the boundary layer quasi-equilibrium hypothesis of Raymond (1995) as slightly modified by Emanuel (1995). The hypothesis holds that the import of low-entropy air into the subcloud layer by convective and clear-air downdrafts balances surface enthalpy fluxes. Its fully nonlinear form it may be expressed

where is the surface enthalpy exchange coefficient, is the magnitude of the surface winds, is the saturation entropy of the sea surface (here taken to be constant), and is a characteristic entropy of the middle troposphere. We linearize this about a mean surface easterly flow, yielding

where the overbars denote mean-state quantities. The first term on the left of (9) is the contribution of the large-scale vertical motion, the second term is the WISHE effect whereby increased perturbation easterlies yield increased surface enthalpy fluxes, which require more convection to evacuate the excess entropy, the third term is a damping of the mass flux by excess warmth of the free troposphere, and the final term is an enhancement of convection by excess midlevel moisture, which requires stronger downdrafts (and thus stronger convection) to compensate for surface fluxes.

The system is closed by a linearized equation for vertically integrated entropy, assuming that fluctuations of the latter are proportional to :

where is a form of gross moist stability,

where is the vertical structure of the vertical velocity perturbation. We will take to be positive here. We have added an artificial diffusion term at the end of (8) to account for the turbulent diffusion of entropy, with coefficient .

The system made of (4)–(7), (9), and (10) is put into nondimensional form by normalizing both the dependent and independent variables, as described in the appendix. The nondimensional equations are listed below, dropping the primes that denote the linear nondimensional variables. The cumulus updraft mass flux has been eliminated using the nondimensional form of (9), and we represent cloud–radiation interaction by taking , where *C* is a nondimensional cloud–radiation parameter. This simple representation holds that there will be more clouds in a moist atmosphere,

Here determines the degree of zonal geostrophy, is the magnitude of the cloud–radiative feedback, and measure the damping effect of boundary layer entropy perturbations on surface fluxes, governs the magnitude of the WISHE feedback, is a normalized gross moist stability, and is the normalized Fickian diffusion coefficient. The factors and account for the different ways that and have been scaled. The complete definitions of the nondimensional coefficients in (12)–(16) are provided in the appendix.

If the cloud–radiative terms and all the damping terms are dropped, as well as the terms involving and , the system reduces to the one considered by Emanuel (1987), which has unstable modes that include eastward-propagating, Kelvin-like modes but which experience zero effective stratification. Without the cloud–radiative terms, the system is similar (but not identical) to those considered by Yano and Emanuel (1991) and Emanuel (1993). It is nearly identical to the linear model of Fuchs and Raymond (2017) except that in that system, given by their (51)–(56), the Newtonian damping and WISHE terms in (15) are missing, as are the damping and cloud–radiative terms in (16). This last is because they close explicitly on moisture, whereas here we use moist static energy instead of moisture (and moist static energy is affected by radiation).

In general, normal-mode solutions to this system are given by parabolic cylinder functions multiplied by where is the zonal wavenumber and is a complex growth rate. Among these solutions is the Kelvin wave–like solution that has everywhere, and given the structure of the numerical solutions discussed in section 3 above, we focus on this solution. Eliminating between (10) and (11) for this mode gives the structure of the solutions in the meridional direction:

Thus, for the solutions to be bounded in , must be negative, corresponding to eastward-propagating disturbances. The dispersion relation for the complex growth rate is cubic, and we look for the most rapidly growing, eastward-propagating mode among the three roots.

Even with the system is controlled by eight nondimensional parameters, presenting an unwieldy phase space to explore.^{1} One can easily show, however, that none of the modes propagates when the WISHE effect is absent , consistent with the NO-WISHE simulations described in section 3. These solutions decay unless the cloud–radiative feedback parameter is sufficiently large. With relatively modest values of zonal diffusion [ in (16)], the highest growth rates of these stationary modes are at low wavenumber; in the absence of diffusion, the growth rates increase with zonal wavenumber but asymptotically approach a constant relatively quickly. Conversely, without cloud–radiation feedback (*C* = 0), WISHE modes occur consistent with those found in previous studies (e.g., Emanuel 1993), and these likewise peak at relatively low wavenumber for modest values of the entropy diffusion Their phase speeds are consistent with Kelvin waves traveling at a dimensional phase speed that represents the first baroclinic mode but with greatly reduced effective stratification. (Here is the buoyancy frequency.)

More interesting solutions are obtained when both WISHE and cloud–radiative feedback are active. Figure 13 shows (in blue) a wavenumber–frequency plot of one class of such modes, taking , , , , , , , and . The three lowest-wavenumber modes correspond to the self-aggregation instability but propagating eastward, at nearly constant frequency, owing to the WISHE effect, which also affects their growth. This frequency is somewhat lower than the extrapolation to low wavenumber of the higher-frequency modes, which behave like WISHE modes altered by cloud–radiation feedback. Their phase speeds vary from about 60% of the moist Kelvin wave speed for wavenumber 1 to 20% for wavenumber 3. We note that the similar model by Fuchs and Raymond (2017) has a short wave cutoff even though they do not have a Fickian diffusive term.

That the first three wavenumbers in Fig. 13 represent a mode distinct from the Kelvin waves can be demonstrated by examining the weak temperature gradient (WTG) approximation to the solutions to (12)–(16). These can be obtained simply by setting the perturbation of the free-tropospheric saturation entropy to zero and solving. This yields a single mode, with whose growth rate and phase speed are given by

and

For the particular values of parameters used for the full solution, only the first three zonal wavenumbers are unstable and these are plotted with red circles in Fig. 13. They are reasonably good approximations to the full solutions, but with slightly higher frequencies and slightly lower growth rates, peaking in this case at wavenumber 2. Not surprisingly, the Kelvin-like modes are entirely filtered by WTG. The form of (18) makes it clear that interactive radiation is necessary for instability and that for and sufficiently small, the growth rate will peak at some finite, nonzero wavenumber. At the same time, (17) shows that WISHE and sufficiently small gross moist stability are necessary for propagation, consistent with the results of the cloud-permitting numerical simulations and that cloud–radiation also contributes to eastward propagation.

Figure 13 bears a striking resemblance to the zonally symmetric part of Fig. 5. The structure of the wavenumber-1 normal mode is illustrated in Fig. 14, which can be compared to Figs. 7 and 8. The structure of this mode is quite similar to that of Fuchs and Raymond (2017, see their Fig. 3).

These results support the inference from the SAM numerical simulations that both WISHE and cloud–radiation feedbacks are essential to explain the MJO. WISHE alone may help explain the faster, higher-wavenumber Kelvin-like modes.

## 5. Conclusions and summary

In this study, we use a near-global aquaplanet modeling framework with uniform and constant SST and perpetual-equinox insolation to perform 260–280-day simulations using a “cloud permitting” grid spacing of 20 km over a domain as wide as Earth’s equatorial circumference and extending to “midlatitudes” at about 46°N/S, where the circulation is limited by nonslip walls. Despite the constant and uniform SST, a mean meridional Hadley cell–like circulation develops with equatorial easterlies surrounded at higher latitudes by westerlies. The simulations are sufficiently long to develop a spectrum of equatorially trapped disturbances, such as Kelvin waves and mixed Rossby–gravity waves, and a strong eastward-propagating MJO-like zonal wavenumber-1 mode with a period of 50–55 days. The MJO-like disturbance has a “swallowtail” or Gill-like pattern of precipitation and column-integrated water vapor anomalies, which closely resemble the composite structure of the observed MJO.

A budget of the spatial variance column-integrated moist static energy for the MJO-like disturbance shows a dominant role of the longwave radiative heating for its development and maintenance, with solar radiation playing a minor role. On the other hand, the surface enthalpy fluxes are found to play an equally significant role in the initiation of MJO-like disturbance, but, at the same time, a relatively minor role during its mature stage. However, it is the surface enthalpy flux anomalies, shifted a quarter-wavelength zonal distance east from the center of the MJO-like disturbance, that facilitates its eastward propagation. This aspect of the simulated MJO is at odds with several observational studies that show that horizontal moisture advection is the principle cause of eastward propagation of the MJO. The absence of background SST and associated tropospheric moisture gradients may contribute to this disparity.

We performed several mechanism-denial experiments to test the role of surface fluxes and radiation in MJO evolution. Homogenizing surface fluxes in the zonal direction does not prevent the formation of the wavenumber-1 mode, presumably by the longwave radiation feedback, but completely arrests its eastward propagation. In this respect, the development of the MJO-like disturbance supports the hypothesis that the MJO may be an example of self-aggregation on sphere, found earlier by Arnold and Randall (2015) using a superparameterized GCM, by virtue of the radiative–convective instability mechanism proposed by Emanuel et al. (2014).

Performing two additional experiments confirms the important role of WISHE in the evolution of the MJO-like disturbance. In the first experiment, zonal homogenization of the surface wind speed that goes into enthalpy surface flux computations, and which is equivalent to eliminating WISHE, completely prevents the formation of the MJO-like disturbance. Moreover, applying the same procedure to the control simulation with a mature MJO-like mode quite effectively eliminates it. In the second experiment, the surface–air enthalpy difference in the flux computations is zonally homogenized, but the wind speed feedback, WISHE, is allowed to operate. A strong MJO-like mode develops even faster than in the control case; however, the direction of its propagation is not only eastward, but can also switch to westward.

Horizontally homogenizing the radiative heating does not eliminate the eastward propagation of the zonal wavenumber-1 mode. The resultant disturbance has faster speeds than the MJO-like disturbance in the control case, corresponding to a period of 30–35 days, consistent with the moist Kelvin wave speed in this case. The disturbance core is found to be collocated with the surface flux positive anomaly, suggesting that the eastward propagation of the disturbance in this case is by the standard gravity wave mechanism of Kelvin waves. And, as shown by (18), the linear stability analysis with the WTG approximation shows no growing modes in the absence of cloud–radiation interaction. Thus, we interpret the modes with homogenized radiation as moist Kelvin waves amplified by WISHE. Homogenizing both surface fluxes and radiation heating prevents the development of any zonal wavenumber-1 disturbance.

These conclusions are bolstered by a linear stability analysis of a deep convecting atmosphere on an equatorial *β* plane, applying the approximations of an always moist adiabatic temperature profile, a convective mass flux based on boundary layer quasi equilibrium, zonal diffusion of midlevel moisture, and a simple parameterization of cloud–radiation feedback. This analysis allows pure eastward-propagating WISHE modes in the absence of cloud–radiative feedback and pure, nonpropagating self-aggregation modes in the absence of WISHE. But the two mechanisms combined lead to a slow, eastward-propagating MJO-like mode with nearly constant frequency, as well as faster, Kelvin-like modes driven mostly by WISHE. The model and the eigenmodes are in some ways similar to those found recently by Fuchs and Raymond (2017). The separation of the MJO-like modes from destabilized Kelvin waves is made clear by applying the weak temperature gradient approximation to the simple linear model. This reproduces quite well the MJO-like mode, showing it to be destabilized by radiation and driven eastward by WISHE, but completely filters all the higher-frequency modes including the Kelvin waves.

Both the numerical simulations and the linear stability analysis suggest that both WISHE and cloud–radiation feedbacks are essential to explaining the principal characteristics of the MJO.

The near-global framework (Bretherton and Khairoutdinov 2015; Narenpitak et al. 2017) applied in this study uses a nonhydrostatic, cloud-resolving model to explicitly represent deep convection in a very large equatorially centered domain that spans considerable planetary-scale distance both in longitude and latitude, but generally is smaller than used by global cloud-resolving models. A similar approach has also been used in the tropical-channel simulations with the Weather Research and Forecasting (WRF) Model (e.g., Tulich et al. 2011). Unlike the aforementioned studies, which use a cloud-resolving uniform horizontal grid spacing of just a few kilometers, the simulations in this study use a cloud-permitting grid spacing of 20 km, mostly because of the prohibitive cost of running long, multimonth, near-global simulations over such a large domain, which is almost twice as large in the zonal direction as in previous studies. We realize that use of such a coarse resolution without any cumulus parameterization may affect the conclusions of this study. However, some confidence in our results is based on the comparison of results of 4- and 20-km near-global simulations presented by Bretherton and Khairoutdinov (2015), who demonstrated that the zonal-mean distribution of cloud cover, precipitation, and radiation fluxes (their Fig. 18) are virtually independent of resolution between 4 and 20 km. Moreover, the corresponding spatial power spectra of cloud cover and precipitable water (their Fig. 17) are also found to be independent of resolution for scales larger than a few hundred kilometers, which are much smaller than the scale of the MJO-like disturbances and Kelvin waves simulated in this study.

To see if the development of the MJO-like disturbance is affected by the resolution, we conducted a relatively short 55-day simulation using a “cloud resolving” 4-km grid spacing, but otherwise using the same setup as in the CTRL case. To avoid the initial spinup, the simulation was initialized on day 20 of the CTRL, about 20 days before the development of the MJO-like mode. The latter, indeed, developed after about 20 days, with an amplitude and propagation speed similar to the CTRL. However, it would currently be computationally prohibitive to apply such a resolution to perform 280-day runs, including the mechanism-denial experiments, as each of the 4-km simulations would require at least 50 times more computational resources than the 20-km simulations.

At first blush, it is surprising that a model with 20-km grid spacing can simulate MJO-like disturbances, given that many coarse-resolution GCMs with parameterized convection struggle to do so. Certainly such processes as entrainment and detrainment are poorly simulated here, and we do not observe strong cold pools that one might take to be important. Understanding why very coarsely resolved convection succeeds where many parameterizations fail is an important issue, but one that is beyond the scope of the present study. We speculate that, given the important role of cloud–radiation interactions, realistic simulation of high-level ice clouds may be essential; this will be the subject of future work.

Another caveat of this study is the domain geometry; in particular, the existence of the walls at the northern and southern boundaries. The version of SAM used in this study is based on Cartesian geometry. We needed to limit the domain in the meridional dimension as there is no support for spherical geometry and extending the domain farther north and south would be unreasonable. At the time of writing this paper, though, a latitude–longitude version of SAM has been developed, which would allow us to extend the domain all the way to the poles to become a true aquaplanet. We have conducted a preliminary simulation with the new latitude–longitude version of SAM using a physical setup identical to the one of this study and found that the MJO-like disturbance still develops similar to the one obtained using our near-global modeling framework. Therefore, we do not believe that the existence of the walls would qualitatively affect the conclusions of the current study.

## Acknowledgments

Kerry Emanuel was supported by the NSF Grant AGS1418508 to MIT. Marat Khairoutdinov was supported by the NSF Grant AGS1418309 to Stony Brook University. All computations have been performed on Yellowstone supercomputer at the NCAR–Wyoming Supercomputer Center. Plots have been produced using the NCAR Command Language (NCL) from NCAR and MATLAB from MathWorks, Inc.

### APPENDIX

#### Definition of Nondimensional Parameters

We here introduce scalings for both the independent and the dependent variables that appear in (2)–(5), (7), and (8) in the main text. In what follows, the dimensional parameters are described in the main text, except for the mean radius of Earth, which is here denoted by *a*. We begin by defining a meridional length scale,

In the following scalings, the dimensional quantities are on the left and their nondimensional equivalents are on the right:

## REFERENCES

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

The reader interested in exploring this phase space and conversant in MATLAB may do so by contacting the second author.