## Abstract

The authors introduce a simple model for the time-dependent evolution of tropical “hot spots,” or localized regions where the sea surface temperature (SST) becomes unusually high for a limited period of time. The model consists of a simple zero-dimensional atmospheric model coupled to an ocean mixed layer. For plausible parameter values, steady solutions of this model can become unstable to time-dependent oscillations, which are studied both by linear stability analysis and explicit time-dependent nonlinear simulation. For reasonable parameter values, the oscillations have periods ranging from intraseasonal to subannual. For parameter values only slightly beyond the threshold for instability, the oscillations become strongly nonlinear, and have a recharge–discharge character.

The basic mechanism for the instability and oscillations comes from cloud-radiative and wind-evaporation feedbacks, which play the same role in the dynamics and are lumped together into a single parameterization. This is possible because, under the assumption that the shortwave and longwave radiative effects of high clouds cancel at the top of the atmosphere, their net effect is only to transfer energy from the ocean to the atmosphere exactly as a surface flux does, and because the two processes are observed to be approximately in phase on intraseasonal timescales. Both feedbacks move energy from the ocean to the atmosphere in convective regions, intensifying the convection and thus destabilizing the system. The same energy transfer cools the ocean, which eventually (but not instantaneously, because of the mixed layer’s heat capacity) reduces the SST enough to render the model stable to deep convection, shutting off the convection. At that point the SST begins warming again under the resulting clear skies, starting the cycle over.

The authors also examine the forced linear response of the model, in a weakly stable regime, to an imposed atmospheric oscillation. This is meant to crudely represent forcing by an atmospheric intraseasonal oscillation. The model’s response as a function of mixed layer depth is not monotonic, but has a maximum around 10–20 m, which happens to be close to the observed value in the western Pacific warm pool.

## 1. Introduction

Observations clearly show intraseasonal variability in sea surface temperature (SST) in tropical regions of high mean SST (Krishnamurti et al. 1988; Lin and Johnson 1996; Zhang 1996; Lau and Sui 1997; Fasullo and Webster 1999; Woolnough et al. 2000; Vecchi and Harrison 2002). It is clear that intraseasonal atmospheric variability, much but not all of which is associated with the Madden–Julian Oscillation (MJO), is largely responsible for generating the SST variability (Waliser and Graham 1993; Zhang 1996; Hendon and Glick 1997; Jones et al. 1998). On the other hand, since SST is generally expected to modulate deep atmospheric convection over tropical oceans, it is natural to expect that this SST variability may influence atmospheric variability on the same intraseasonal timescale, leading to coupled dynamics. We present a simple model that aims to capture one key mechanism of coupled tropical ocean–atmosphere variability on the intraseasonal timescale. This mechanism results from local interactions between deep convection, cloud–radiative processes, surface energy fluxes, and ocean mixed layer thermodynamics.

Since the MJO is the dominant mode of intraseasonal variability in the Tropics, the present study is relevant to the following question: ‘to what degree is the MJO a coupled mode?’ General circulation model (GCM) studies focused on this question suggest that the thermodynamic response of the ocean mixed layer to the MJO has a nonnegligible effect on the atmosphere, amplifying the MJO and lowering its frequency (Flatau et al. 1997; Wang and Xie 1998; Waliser et al. 1999; Hendon 2000; Kemball-Cook et al. 2002; Maloney and Kiehl 2002; Inness and Slingo 2003), though Hendon (2000) found essentially no effect from coupling. While we are in part motivated by these studies, the model presented here is not truly a model of the MJO, because it lacks horizontal structure and propagation characteristics, which are essential components of the observed MJO (Madden and Julian 1994; Zhang and Hendon 1997).

A more direct observational motivation for our model is provided by Waliser (1996). Waliser’s study was motivated by the apparent inability of SST to exceed a limit somewhere around 30°C, except by a small margin on relatively short timescales. This had led to much debate regarding the mechanism that sets this limit (Newell 1979; Ramanathan and Collins 1991; Fu et al. 1992; Wallace 1992; Waliser and Graham 1993; Hartmann and Michelsen 1993; Waliser 1996; Lindzen et al. 2001; Hartmann and Michelsen 2002a,b; Lindzen et al. 2002). As pointed out by Waliser (1996), much of the disagreement is relevant only to the tropical mean SST. For a given tropical mean SST, the question of what prevents the local SST in the warmest regions from rising above a maximum value, presently around 30°C (a number that may perhaps vary on long timescales as a function of the mean climate), is somewhat easier to address. Waliser (1996) reasoned that whatever the limiting mechanism is, it probably acts in a time-dependent fashion to damp out “hot spots,” or localized regions where the SST exceeds the limiting value on spatial and temporal scales above some thresholds (which he chose to be a few hundred kilometers and about a month). This motivated his detailed observational study of the evolution of composite hot spots, which suggested that the reduction in surface shortwave radiation by high clouds is the dominant mechanism that limits SST, although the increase in surface latent heat flux also plays a significant role. Cubukcu and Krishnamurti (2002) analyzed the mechanisms controlling western tropical Pacific SST on a number of different timescales in a fully coupled atmosphere–ocean GCM and found that, on the intraseasonal timescale, the shortwave feedback is dominant while the surface evaporation feedback is secondary—in agreement with Waliser’s observational result—while the converse was true on the semiannual timescale.

We present an idealized model for the time-dependent behavior of SST hot spots, or localized regions of high SST. Our model represents, in a highly parameterized way, feedbacks between SST, deep convection, high cloud-radiative feedbacks, and surface fluxes. We focus on this subset of processes because observational and GCM studies suggest that the modulation of shortwave radiation by cloudiness is either the dominant mechanism, or at least a significant one, by which tropical SST is limited on these timescales (Waliser 1996; Hendon and Glick 1997; Lau and Sui 1997; Jones et al. 1998; Fasullo and Webster 1999, 2000). Our model can, under a certain interpretation, be considered to represent the evaporation feedback as well as the cloud-radiative feedbacks. However, for reasons to become clear shortly, even when seen in the most favorable light the model can in no way be construed as providing any independent evidence for the relative importance of the radiative and evaporative feedbacks in controlling SST.

Particularly because of its severe truncation of atmospheric dynamics, our model is considerably simpler than the model of Wang and Xie (1998), while at the same time being somewhat more explicitly physical than those of Waliser and Graham (1993) and Hartmann and Michelsen (1993), all of which have been used to address related questions. Since our model contains no horizontal dimension and is formally valid only in the limit of small horizontal scale (see the appendix), it is not a model of the MJO per se, although it may nonetheless be helpful in understanding the role of the ocean in modifying the MJO.

The following section introduces the model. We then present a linear stability analysis, show some nonlinear numerical solutions, discuss the mechanism of the oscillation, and then consider the forced problem of the same system driven by an imposed atmospheric oscillation. We conclude in section 4.

## 2. Model description

The atmospheric model is a zero-dimensional reduction (i.e., to a single column with a fixed vertical structure) of the Neelin–Zeng Quasi-Equilibrium Tropical Circulation Model [QTCM; Neelin and Zeng 2000; Zeng et al. 2000], with the atmospheric temperature held fixed, as has been done in a number of tropical single-column models (SCM) studies (Raymond 2000; Sobel and Bretherton 2000; Zeng 1998; Zeng and Neelin 1999; Chiang and Sobel 2002), and as is discussed in detail by Sobel and Bretherton (2000). The fixed temperature allows the large-scale vertical motion to vary, as required in order to balance diabatic heating variations generated by the model physics. This large-scale vertical motion also controls the moisture convergence and so feeds back on the physics. The resulting dynamics is rather different from that obtained under fixed large-scale vertical velocity (e.g., zero, as in radiative–convective equilibrium) but variable temperature. The SCM studies of intraseasonal radiative–convective oscillations by Hu and Randall (1994, 1995) used the latter assumption, limiting their precipitation anomalies to very small amplitudes (usually less than 0.5 mm day^{–1}). Fixing the temperature allows our precipitation anomalies to be much larger (Sobel and Bretherton 2000). This is important because cloud–radiative and surface flux feedbacks, which are parameterized here as a function of precipitation, are essential for the mechanism of the oscillations on which we focus.

The equations for the atmosphere are

Equation (1) originates as an equation for the atmospheric temperature *T,* but the rate of change *dT*/*dt* has been dropped. Here *T* is actually the heat capacity of air (*C*_{p}) times the temperature, while *q* is the latent heat of vaporization of water times the specific humidity (in practice, we then divide both *T* and *q* by *C*_{p} and express them in degrees). The *δ* is the divergence of the horizontal velocity, **v**. Here *T, **q,* and **v** are functions of *t* and implicitly multiply a fixed, specified vertical structure functions in the vertical coordinate (Neelin and Zeng 2000). To obtain the physical fields for *T* and *q,* specified reference profiles must be added to the variable components that solve the above equations. The basis functions for *T* and *q* have a single sign in the troposphere, while that for **v** changes sign once with height, the sign in the equations above representing its value in the upper troposphere. [We use the structure functions, reference profiles, and constants and coefficients derived from them, from the version 2.2 of the QTCM; see Table 1 of Neelin and Zeng (2000) or Table 1 of Bretherton and Sobel (2002). The temperature structure function is close to a moist adiabat. The humidity structure function is not exactly equivalent to constant relative humidity with height, but is almost equivalent to that in practice, since only bulk, integrated aspects of the profile are ever used in the model.] The *P* is precipitation, *R* radiative cooling (assumed positive), and *E* surface evaporation. We have neglected surface sensible heat flux since it is small compared to the evaporation over tropical oceans. Here *M*_{S} and *M*_{q} are the dry static stability and gross moisture stratification (Neelin and Held 1987), taken in Neelin and Zeng (2000) to be linear functions of *T* and *q*:

with *M*_{sr}, *M*_{sp}, *M*_{qr}, *M*_{qp} constants derived from the structure functions for *T* and *q.* Here, since *T* is constant, *M*_{S} is also. Here *b̂* is a constant coefficient that results from the projection of the original three-dimensional moisture equation on the single vertical structure function for *q.*

Physically, (1) states that the imbalance between convective heating and radiative cooling must be balanced by horizontal convergence of dry static energy. Equation (2) states that the imbalance between precipitation and evaporation can be balanced either by horizontal moisture convergence or the tendency of specific humidity, or both. Equations (3) and (4) state that the rates of convergence of dry static energy and moisture per unit mass convergence are linear functions of the atmospheric temperature and humidity, respectively. These equations are linearizations of the true dependences, which are weakly nonlinear (Neelin and Zeng 2000).

In (2), the horizontal advection term has been dropped, equivalent in the primitive equations to the approximation ∇·(**v***q*) ≈ *q*∇·**v** with ∇. the horizontal divergence and **v** the horizontal vector velocity. We make this approximation because it is impossible to model the horizontal gradient of any quantity satisfactorily in a single-column model. The primary undesirable effect on our simulations is that, as the atmosphere moistens locally, the system as written above assumes that the atmosphere in adjacent regions moistens at the same rate, so that the moisture convergence associated with mass convergence also increases, due to the dependence of *M*_{q} on *q* in (4). Under certain conditions, this feedback can cause an instability, but we view this instability as spurious. A more natural physical assumption would be that as the atmosphere moistens locally in one spot, the humidity of adjacent regions is unchanged. This implies a greater moisture gradient, in the sense that greater drying will occur as air is advected into the moist region at low levels (as it will be during times of significant precipitation), yielding a negative feedback. We attempt to mimic this crudely, and to suppress the positive feedback resulting from the q dependence in (4), by taking *M*_{q} constant. Thus, (3) and (4) are not used in the model, except to render *M*_{S} consistent with the value of *T* used in the radiation parameterization [see (5)], and to render *M*_{q} consistent with a typical mean value of *q* taken from a nonlinear simulation in which (4) was used and the other parameters were close to our control values.

The precipitation is parameterized using the simplified Betts–Miller scheme

with *τ*_{c} the convective timescale. The radiative cooling is parameterized by

where

with *τ*_{R} a radiative timescale, *T*_{e} a radiative equilibrium temperature, and *r* a nonnegative dimensionless cloud–radiative feedback parameter. Equation (7) represents clear-sky conditions, while (6) accounts for the cloudy case. Evidence justifying this simple parameterization of the greenhouse warming effects of high clouds is presented by Bretherton and Sobel (2002). Because *T* is fixed, *R*_{0} is constant once *τ*_{R} is chosen. The explicit restriction that *R* be nonnegative is an admittedly crude way of preventing the unrealistic occurrence of net radiative warming of the troposphere. The evaporation is parameterized by the bulk formula

where *q*^{*}_{s} is the saturation specific humidity at the sea surface temperature, *T*_{S}, and the surface specific humidity, *q*_{s}, is related linearly to *q* by the vertical structure function and reference profile. Here *τ*_{E} is the basic surface flux timescale. Fixing it is equivalent to taking both the surface wind speed and exchange coefficient (the latter of which might in general be a function of, e.g., near-surface stability) to be fixed. This a major limitation of the model, though one that can in part be remedied by a reinterpretation of the parameter *r* discussed shortly below.

The lower boundary condition is a dynamically passive bulk mixed layer, satisfying

where *C* is the heat capacity of the mixed layer and *S* represents the combined effects of shortwave and longwave radiation. Ocean dynamics is neglected on the basis that it has been found to exert a weak influence on SST in the warm pool regions (Seager et al. 1988; Gent 1991) where hot spots tend to develop (Waliser 1996) and the MJO is most active. The *S* will be modified by a cloud–radiative feedback which is considered to act on its shortwave component,

with *Ŝ* = constant. The longwave component of *S* is considered fixed. Together (6) and (10) imply that the net effect of high clouds (the only kind we attempt to explicitly represent in this model) vanishes in the top-of-the-atmosphere (TOA) radiation budget; high clouds warm the atmosphere and cool the ocean in equal measure, as seems to be approximately true in observations (Ramanathan et al. 1989; Harrison et al. 1990; Hartmann et al. 2001).

An interesting aspect of the TOA cancellation is that it renders the cloud–radiative feedback equivalent to a surface flux, in that its net effect is to transfer energy from the ocean to the atmosphere with no net loss or gain to the coupled system. This means that, to the extent that variations in surface latent heat flux and surface shortwave radiative flux variations are in phase for the modes of variability of interest, the cloud–radiative feedback parameterization represented by Eqs. (6) and (10) can be reinterpreted as representing the sum of cloud–radiative and surface flux feedbacks, with a suitable increase in the chosen value of *r.* Strictly viewed this way, *r* incorporates the wind-induced component of the surface flux feedback, since surface humidity variations are explicitly modeled and do affect the surface evaporation through (8). Observational studies of intraseasonal-to-subannual timescale SST and surface flux variability show variations in surface latent heat flux and shortwave radiative energy flux to be nearly in phase (Hendon and Glick 1997; Lau and Sui 1997; Jones et al. 1998; Fasullo and Webster 2000), within a week or so, which is fairly short compared to the period of interest here.

Key model parameters are shown in Table 1 [see Zeng et al. (2000) or Bretherton and Sobel (2002) for the basic QTCM parameters]. Of these the least well-constrained by observations are probably the cloud-radiative/surface flux parameter *r* and the convective timescale *τ*_{c}. Based on monthly mean observations, Bretherton and Sobel (2002) suggested *r* ≈ 0.2 for radiative processes alone, but with considerable uncertainty. They suggested that a range of 0.1–0.3 is at least plausible. J. Lin and B. Mapes (2003, personal communication) found a value of 0.1–0.15 for the infrared cloud-radiative feedback on the MJO specifically. (Lin and Mapes also found that for the MJO the shortwave forcing is considerably larger than the longwave one at TOA, in contrast to what is found at longer timescales; we do not try to incorporate this result here.) We choose *r* = 0.25; this seems reasonably conservative if it is viewed as including the surface evaporation feedback as previously discussed. The convective timescale *τ*_{c} is often taken to be a few hours, based largely on the performance of general circulation models using the Betts–Miller scheme, with some dependence on the grid spacing (Betts 1997). A timescale with an equivalent dynamical role arises in theoretical or modeling studies that relate the column-integrated precipitable water or saturation deficit to precipitation (Raymond 2000; Fuchs and Raymond 2002; Sobel and Bretherton 2003), and these studies obtain values of the order of 1 day. We take 1 day as our default value for the linear stability plots below, but examine sensitivity to this parameter among others. In a closely related study (Gildor et al. 2003) we used 0.3 days.

The values for *τ*_{R} and *τ*_{E} are chosen to give typical observed values for (clear sky) *R* and *E,* 1–2 K day^{–1} and 3–5 mm day^{–1}, respectively. Here *τ*_{E} = 12 days corresponds to surface wind speeds of a few meters per second, for typical exchange coefficients and surface air–sea humidity differences in over tropical ocean regions with frequent deep convection such as the western Pacific warm pool. Sensitivity to these parameters is not shown in the next section, but the system is less sensitive to both than to either *r* or *τ*_{c}. The linear problem is not sensitive at all to *τ*_{R}, as it does not enter the linear perturbation equations because of the assumption that atmospheric temperature perturbations are negligible.

## 3. Results

### a. Linear stability

Writing overbars for basic-state quantities and primes for perturbations, using the explicit forms of the parameterizations of (5), (6), (8), and (10) (assuming that in the basic state *P* > 0), the linearized perturbation equations are

where the constant *γ* comes from linearizing the Clausius–Clapeyron equation:

Assuming solutions proportional to *e*^{σt} gives the dispersion relation,

where

is the gross moist stability (Neelin and Held 1987; Neelin 1997; Neelin and Zeng 2000) and

is the “effective gross moist stability” (Su and Neelin 2002; Bretherton and Sobel 2002) [defined slightly differently than in Bretherton and Sobel (2002)]. Here *b*_{s} is the proportionality coefficient between the surface specific humidity and *q.*

The solutions are

All the parameters and coefficients are real, so the square root is either pure real or pure imaginary. If it is imaginary, it gives the frequency of oscillation, and the term outside the square root gives the growth rate. In that situation, the first term under the square root is proportional to the square of the growth rate, and the second term is strictly positive as long as *M* is (which we assume is the case, in general). From this we can see that if the absolute value of the growth rate is sufficiently large—which in practice means, on the order of the frequency—*σ* becomes pure real, and the linear problem has only nonoscillatory solutions, either growing or decaying.

To illustrate the behavior of the growth rate and frequency as a function of parameters, we contour them as a function of *r* on the *y* axes, and the convective timescale *τ*_{c} and mixed layer depth (*mld*) on the *x* axes in Fig. 1. The quadratic dispersion relation (15) has two solutions, and at every point in the figures the solution shown is the one with the larger growth rate, or real part of *σ,* which can require switching branches.

Over most of the space shown, and particularly near the control values of the parameters (shown by the plus symbols), the growth rate increases (becomes either more positive or less negative) with *r* and *mld* and decreases with *τ*_{c}. The dependences on *r* and *τ*_{c} are much stronger than that on *mld,* or on other parameters which are not shown such as *τ*_{E}. Only for sufficiently large *r* and sufficiently small *τ*_{c} is the system actually unstable. As shown below, the transition from the stable to the unstable regime appears to occur via a subcritical Hopf bifurcation. Therefore, if the system is close enough in parameter space to the bifurcation point, even a small finite perturbation in one of the state variables may shift it from the fixed point onto the limit cycle, with an apparently discontinuous jump in amplitude of the cycle from zero (i.e., at the fixed point) to a finite value. In Gildor et al. (2003) such transitions result from the effect of ocean biota on SST. At fixed *τ*_{c}, the stability of the system depends mainly on the strength of the cloud–radiative feedbacks, or the combination of cloud–radiative and surface evaporation feedbacks that *r* may represent. At small *τ*_{c} there is a minimum growth rate as a function of *r* around *r* = 0.17, which occurs because of branch switching as we plot the fastest-growing or slowest-decaying solution, as mentioned above. Individual solutions are monotonic in *r* at given *τ*_{c}.

Over most of the parameter space in which the linear frequency is nonzero, it corresponds to periods that are intraseasonal (0.1–0.2 corresponds to approximately 30–60 days because the period in days is 2*π*/*ω,* where *ω* is the frequency) or somewhat longer. In *r*–*τ*_{c} space, the frequency is nonzero only in a finite region roughly centered on the marginal stability curve [Re(*σ*) = 0], something that we can understand from inspection of the dispersion relation as mentioned above. To get a little more insight into how different parameters control the frequency, we can consider the special case of the marginal stability curve, on which the real part of *σ* vanishes and the frequency simplifies to

Despite that this represents only a particular slice through parameter space, the result gives a qualitative insight into what controls the frequency in nearby regions. Apparently several different parameters contribute comparably to the frequency. Of those that are most important, the gross moist stability *M* and convective timescale *τ*_{c} are the most poorly constrained by observations. The frequency depends on each only as the square root, however, so that modest changes in any of them will not make a dramatic difference. The frequency depends on *r* only as (1 + *r*)^{1/2}, a very weak dependence for *r* < 1.

The condition for neutral stability is

with instability for more negative *M*_{eff}. From this we can see the role of the finite mixed layer depth and convective timescale in stabilizing the system, compared to the fixed-SST, “strict quasi-equilibrium” (SQE; *τ*_{c} = 0) model studied by Bretherton and Sobel (2002). Their model also had a horizontal dimension, allowing horizontal moisture advection **u** · ∇*q,* but under the weak temperature gradient (WTG) approximation and SQE one has *q* = *T* = constant in convective regions, so the horizontal advection term vanishes in such regions and the dynamics becomes locally describable there as a single-column system similar to ours here (though the horizontal moisture advection term is important in determining the size of the convective region). Bretherton and Sobel (2002) found that a steady solution was not possible when *M*_{eff} [defined by them slightly differently as *M*_{eff} = (*M* – *rM*_{q})/(1 + *r*)], became negative, implying instability to time-dependent disturbances. Sobel (2003) extended essentially the same model to the coupled case and found that this criterion was no longer exactly correct. Our model corresponds locally well to Sobel’s within the convective region. Equation (18) shows that we recover the criterion of Bretherton and Sobel (2002) in the SQE limit, but that finite *τ*_{c} stabilizes the system even for fixed SST. Increasing mixed layer depth destabilizes the system for finite *τ*_{c}, but will have no effect on linear stability in the SQE limit. Even for finite *τ*_{c}, the mixed layer depth has only a weak effect on stability because for typical parameters *γ*/*C* is less than unity, while *b*_{s} = 1; for our control parameters *γ*/*C* = 0.27.

### b. Nonlinear behavior

The linear dynamics offers only limited information about the actual time-dependent behavior of the system. Comparison with nonlinear simulations shows that it predicts quite accurately the onset of instability, Re(*σ*) = 0, and gives a broadly reasonable estimate of the oscillations’ period. However, slightly past the onset of instability (whatever control parameter is used), the oscillations become strongly nonlinear. The time dependence becomes highly nonsinusoidal, and in some cases the period becomes somewhat different than the linear calculation predicts although not by anything like an order of magnitude. Figure 2 shows a nonlinear simulation in which *r* = 0.3; other parameters are as in Table 1. The system is initialized with an essentially arbitrary initial condition, to which the longtime behavior is entirely insensitive. The time “zero” on the plots in Fig. 2 is arbitrarily chosen as a point at which all memory of the initial conditions has been lost. There are no time-dependent external forcings. If all parameters are as in Table 1, the system is (very weakly) stable, and the system reaches a steady state; the time-dependent behavior shown in the figure occurs because of the choice *r* = 0.3 rather than the default value *r* = 0.25. The precipitation *P,* evaporation *E,* radiative cooling *R,* surface shortwave flux *S,* SST, and surface humidity *q*_{s}, are shown as functions of time. Note that the evaporation is that computed directly from (8); it does not explicitly include the “extra” evaporation variations that we might take to be included with the cloud–radiative feedback.

We vaguely see square-wave-type variations in precipitation and (accordingly) in *R* and *S.* Detailed features, such as the particular shape of the precipitation pulses and their duration relative to the duration of nonprecipitating periods, are dependent on parameters. (In particular, a regime with shorter precipitating pulses and longer clear periods, which may be more consistent with observations, can be obtained by reducing both *τ*_{R} and *τ*_{c}, e.g., to 25 and 0.3 days, respectively.) The SST, here and over a wide range of parameters, has more sawtoothlike variations. The oscillation overall has characteristics of a recharge–discharge oscillation. The evaporation (with the caveat above) has variability that is of the wrong sign compared to observations (Waliser 1996; Zhang 1996). Because of the constant surface wind speed, the evaporation is determined only by the difference between the saturation specific humidity at the SST and the value of *q* at the surface. This difference has a local minimum at the time of maximum precipitation. This error in *E* is not of great concern to us, because the *E* variations are small compared to those in *S* so that the latter dominate the behavior, and because as discussed above we might take the “total” evaporation to be the sum of that computed from (8) and some component of the time-variable part of *S.*

### c. Bifurcation behavior

The analysis above suggests that somewhere between *r* = 0.2 and *r* = 0.3 there is a Hopf bifurcation where the model behavior changes from a stable steady state to nonlinear oscillations. As was also shown above, similar transitions from a steady state to an unstable, oscillatory regime can result from variations of other model parameters. We now look more closely at these transitions using bifurcation diagrams. The analysis was done with the software AUTO (Doedel et al. 1997) using XPPAUT (Ermentrout 2002) as an interface. [In computing the following results, the step functions in Eqs. (5) and (6) were approximated by sharp tanh functions. The results from direct simulations using the tanh functions were compared to solutions using the step functions and no significant difference was found.]

In Fig. 3 we see how the steady-state value of SST varies as *r* is changed, with all other parameters kept fixed at their default values. As *r* increases, the steady state goes from being stable (solid line) to unstable through a Hopf bifurcation at *r* ≈ 0.255. The dramatic jump from stable solution to large amplitude oscillations suggests a subcritical Hopf bifurcation (Strogatz 1994). After the Hopf bifurcation we get a stable limit cycle, with the solid circles representing the maximum and minimum values. The minimum value decreases as *r* increases, as expected from the discussion above. In contrast, the maximum value of SST stays constant once the minimum precipitation over the cycle reaches zero, around *r* = 0.275. In the lower panel we see that the period of the oscillation for the plotted range of *r* increases slightly, from 70 to 82 days.

Qualitatively similar behavior is seen when *τ*_{c} is varied while the rest of the parameters are kept at their default values, with the general sense of the behavior change as *τ*_{c} decreases being the same as that as *r* increases, as we can understand to some degree from the linear analysis (Fig. 4). As *τ*_{c} decreases in the oscillatory regime, however, the maximum SST does decrease, while the difference between the maximum and minimum SST stays approximately constant. This presumably reflects the very small sensitivity of the period to *τ*_{c}, as shown in the figure. If the period were to change greatly, we would expect the amplitude of the SST oscillations to change, with amplitude increasing for longer period as the positive and negative changes to the net surface energy flux would have longer times in which to act.

### d. Discussion

In this section we discuss the basic instability mechanism in order to interpret some of the results in Figs. 1 and 2. For simplicity of exposition, we initially discuss the mechanism as though *r* were to represent only cloud–radiative feedbacks.

Consider an initial condition in which the SST has some value below that associated with the (unstable) steady solution, *P* = 0, and *q* is less than *T,* corresponding to the stability of the sounding to deep convection. The surface evaporation *E* adjusts *q* toward the saturation value set by the SST. If this value is greater than *T* then at some point convection will commence and *P* will become positive.

The infrared cloud–radiative feedback taken alone is positive, as we can see by considering the moist static energy equation, derived by adding (1) and (2):

In all our calculations *M* remains positive, as required if the circulation is always to be thermally direct (Neelin and Held 1987). The storage of moisture in the atmosphere is small, so upward motion (*δ* > 0), which must occur if *P* > *R,* can be sustained by negative ∂*q*/∂*t* only for a short time. Beyond that, a positive net source of moist static energy, (*E* – *R*) > 0, is required to balance the export, as in the steady case (Neelin and Held 1987). Since under our parameterization *R* decreases as *P* increases, the infrared cloud–radiative feedback amplifies the convection. In other words, as convection increases in amplitude, *R* decreases, which makes the rhs *E* – *R* increase. Since the storage term is small, the lhs can only balance the rhs by an increase in *Mδ,* which implies a stronger divergent circulation, meaning stronger convection and thus smaller *R,* etc. This mechanism, in a spatially resolved model with a fixed SST lower boundary condition, will tend to collapse the convection to smaller horizontal scales in addition to increasing its amplitude (Raymond 2000; Grabowski et al. 2000; Lee et al. 2001; Bretherton and Sobel 2002).

The shortwave cloud feedback, on the other hand, is negative. When *P* is positive, *S* decreases. This tends to reduce the SST. The surface evaporation acts to adjust *q* toward saturation at the surface, so as SST is reduced, so is *q,* and therefore also *P.* For some parameter values, this can lead to a stable steady state under conditions for which the fixed-SST system would be unstable. In a steady state, the larger *r* and *P* are, the smaller *E* will be, because in a steady state *E* = *S.* Thus the shortwave feedback damps the moist static energy source *E* – *R,* which provides the driving for *P.*

The heat capacity of the mixed layer induces a lag between the onset of precipitation and associated reduction in SST, and the reductions in *E* and *q,* which carry the negative shortwave feedback to the atmosphere. On the other hand, the reduction in *R* due to the longwave feedback happens on a timescale of the order of *τ*_{c}, which is nearly instantaneous compared to the rate of SST change. Therefore, for finite mixed layer depth, the oscillation occurs because the positive longwave feedback has some time to act, increasing *P,* before the shortwave negative feedback finally becomes significant, and eventually reduces *P* by reducing the SST.

Again, if we wish to do so, we can lump together the surface evaporation feedback together with the radiative feedbacks, without changing the substance of the argument, because both take energy out of the ocean and put it into the atmosphere. The surface evaporation, like the longwave cloud–radiative feedback, affects the precipitation on a timescale *τ*_{c}, while the SST responds equally slowly to changes in both shortwave radiative flux and evaporation.

From this explanation, we can see why larger mixed layer depth is generally associated with larger growth rate. If the mixed layer depth is small, the shortwave feedback is able to act more quickly to reduce the SST and thus suppress convection. The tendency of frequency to decrease as mixed layer depth increases is a more straightforward consequence of the fact that *C* appears directly in the position of a time constant (though it is a nondimensional quantity) in (9).

### e. Dependence on mixed layer depth and forced problem

The simplicity of our model means that its applicability to observations is mainly indirect. Its value is primarily in building intuition that can be applied to interpretations of more complex models. In particular, some of the parameters in our model have analogs in GCMs. We can easily vary these parameters and see how the behavior of our model changes. If SST variability in a GCM coupled to a mixed layer ocean changes similarly in response to varying analogous parameters, the simple model’s behavior may provide some insight into the mechanisms.

One parameter that is straightforward to vary, in a GCM or our model, is the mixed layer depth, and Watterson (2002) has recently performed a GCM study examining the sensitivity to this parameter. He found that the intensity of subannual and intraseasonal SST variability increased as the mixed layer depth was decreased from 50 to 10 m. This is in direct contrast to our model in which the growth rate increases (albeit weakly) with mixed layer depth. In agreement with our results, as well as naive intuition, Watterson did find that the spectral peak (in low-level zonal wind) moved to lower frequency as the mixed layer depth increased. However, the contradiction between our simple model and the GCM regarding the amplitude, inasmuch as our linear growth rate might perhaps be taken as an indicator of it, requires explanation.

We believe that the discrepancy is due to the stationary nature of the disturbances in our model, as opposed to the propagating MJO signal in the GCM. Clearly these affect SST variability differently, as the phase relations between the atmospheric and SST variability must be different in the two cases. A rough way to examine this issue is to consider the response of our model to forcing by an imposed time-dependent atmospheric disturbance. Clearly this is an oversimplification of the coupled interaction of an ocean mixed layer with a true atmospheric MJO, but it allows us to gain some beginning insight into the problem.

It seems most appropriate to put the forcing in the atmospheric temperature equation, so (1) becomes

with *σ*_{0} and *F*_{0} real; the factor *M*_{S} on the rhs is inserted for convenience later. This is a slightly different formulation than used by Chiang and Sobel (2002) to study the response of a (different) single-column atmosphere/mixed layer ocean model to imposed interannual temperature variations; they imposed the variations directly in *T* itself. Their model also did not include cloud-radiative feedbacks, and had more sophisticated parameterizations of convection and gaseous radiative transfer. In overall formulation, however, the present problem is similar to theirs.

We linearize the system, and assume that the solution variables *q, **T*_{s}, *δ* all oscillate with the frequency *σ*_{0},

where *q*_{0}, *T*_{s0}, *δ*_{0} are amplitudes that may be complex. The problem then is to find these amplitudes as a function of *σ*_{0} and mixed layer depth, and examine their dependence on mixed layer depth for values of *σ*_{0} in the intraseasonal–subannual range. As far as the overall magnitude of the response, measured by the absolute value of the amplitude, it is sufficient to find one of the amplitudes, as they will all be proportional to *F*_{0} and thus to each other.

With the exception of the new forcing term, the linear system is unchanged from (11) to (13), so it is sufficient to write down the solution, which is, expressed in terms of *δ*_{0},

where

In Fig. 5, we show the absolute value of the normalized amplitude | *δ*_{0}/*F*_{0} | as a function of mixed layer depth, for other parameters set to their default values. Two forcing frequencies, corresponding to periods of 31 and 63 days are shown. The response amplitude is not monotonic as a function of mixed layer depth for fixed *σ*_{0}. There is a maximum around 10–30 m, which happens to be the climatological mixed layer depth in the western Pacific warm pool (Lukas and Lindstrom 1991). We can understand this quite simply by noting that, at neutral stability, our system is just a typical undamped, unforced linear oscillator, and must have a resonance if it is forced at its natural frequency. For parameter choices that render the system weakly stable, we therefore expect a peak in the response near that frequency. Because the parameter values used in Fig. 5 put the system in a very weakly stable regime, the peak near the natural frequency at *σ*_{0} = 0.1 day^{–1} is quite strong. Further from resonance (as for the *σ*_{0} = 0.2 day^{–1} curve shown), or for a more stable system, the peaks are less pronounced.

The result of interest is that the amplitude of the forced response is nonmonotonic in mixed layer depth. There is a range over which the response does become larger as mixed layer depth is decreased, as found by Watterson (2002). However, if the mixed layer depth is made smaller than the value at peak response, around 10–30 m, the response decreases again. This is contradictory to the expectation we might have that because the surface temperature can vary more easily for smaller mixed layer depth, variability must be larger. Physically, at zero mixed layer depth, the negative feedbacks on the SST from cloud–radiative and surface flux perturbations become as instantaneous as the positive feedbacks on the atmospheric moist static energy budget from the same perturbations. This is as opposed to the situation at finite mixed layer depth, where the thermal intertia of the mixed layer depth creates an effective phase lag between the two feedbacks, allowing growth of perturbations through the atmospheric feedbacks during the time before the SST can respond appreciably.

Given the potential importance of nonlinearity, interactions with other atmospheric phenomena, the broad range of frequencies involved, etc., it is impossible to go directly from our very simple calculations to a quantitative prediction of the variance changes as a function of mixed layer depth that should be found in GCM experiments like Watterson’s. However, the smallest value used by Watterson was 10 m, near the values at which our maxima occur. We cautiously predict that if similar studies to his are taken to smaller mixed layer depths, nonmonotonic behavior may be found in the amplitude of tropical variability at intraseasonal to subannual frequencies. Variance may decrease with mixed layer depth for larger mixed layer depths, as Watterson found, but at smaller mixed layer depths variance may increase with mixed layer depth. This is a testable prediction of the present study. The results of such a test, naturally, may be model dependent, particularly as cloud–radiative feedbacks are among the aspects of GCMs that we generally expect to vary the most.

## 4. Conclusions

We have presented a simple model for localized SST hot spots. Our model consists of a single-column atmosphere, obeying the QTCM physics (single baroclinic mode, Betts–Miller-type convective adjustment), over a dynamically passive bulk mixed layer ocean. Atmospheric temperature perturbations are neglected, which is consistent for small horizontal scales. The model can be thought of as representing a region large enough that its SST anomalies can modulate atmospheric deep convection locally, but not so large that the resulting convection anomalies drive significant local atmospheric temperature anomalies on the timescales of interest (i.e., longer than the periods of gravity waves whose horizontal wavelengths are comparable to the size of the region).

A key ingredient in the model is the radiative effect of high clouds, which is taken to be zero at the top of the atmosphere, with equal longwave warming of the atmosphere and shortwave cooling of the ocean, both proportional to the precipitation rate. Although we cannot physically model the surface flux variations induced by surface wind speed variations, they can be considered to be included in a rough empirical way. This is because the surface flux variations are observed to be nearly in phase with cloud–radiative flux variations on the timescales of interest, and enter the moist static energy budget in the same way, as a transfer of energy from the ocean to the atmosphere.

The steady solutions of the model are unstable to small perturbations for plausible parameter values, particularly if we imagine the radiative feedback parameter *r* to represent the wind-induced surface flux feedback as well, implying a larger value. However, there is sufficient uncertainty in this parameter, across a range of values that straddles the stability boundary, that we cannot say that the actual ocean–atmosphere system is unstable. A similarly large uncertainty is associated with the other parameter that most strongly controls stability, the convective timescale *τ*_{c}.

The linear problem has oscillatory solutions only over a finite region of parameter space near the marginal stability boundary. Outside that region, linear solutions are purely growing or decaying. The nonlinear problem oscillates in the unstable region where the linear frequency is zero; in that region the oscillation mechanism is purely nonlinear. This is consistent with the notion that oscillations of this type play a role in limiting the maximum value of the SST, as suggested by Waliser (1996). Since the maximum is present in long-term means, and thus requires a rectified effect from any time-dependent mechanism, that mechanism cannot be a linear oscillation.

The transition from stable to unstable regime occurs via subcritical Hopf bifurcation, meaning that just after the bifurcation point the amplitude of the oscillations is quite large. In a related study (Gildor et al. 2003), we argued that such a transition might result from the effect of ocean biota on light penetration and hence on SST. Although the direct effect of the biota is small, if the atmosphere is close enough to the transition point, this finite amplitude perturbation might be strong enough to cause the atmosphere to cross the bifurcation point and thus stimulate the nonlinear oscillations. Of course, besides biological processes, many other sources of “noise” from internal atmospheric or oceanic dynamics could also force the system to cross the bifurcation point.

Whether linear or nonlinear, the oscillation frequency tends to be intraseasonal to subannual. The frequency depends on several parameters, particularly the cloud-radiative/wind-evaporation feedback parameter, the convective timescale, and the gross moist stability, but is not strongly sensitive to any of them, so that the frequency is robust at least to an order of magnitude.

The mechanism of the oscillation is easy to understand. Under clear skies, the SST increases due to high surface shortwave radiative energy flux and low surface evaporation. The atmospheric humidity increases along with the SST, in the absence of any precipitation to dry the atmosphere (though moisture divergence does partially compensate surface evaporation, reducing the increasing humidity tendency). Eventually, the system becomes unstable to deep convection. The longwave greenhouse effect of the resulting clouds amplifies the convection by increasing the net source of moist static energy to the column, requiring more convection to export that energy. At the same time, however, the surface shortwave flux is dramatically reduced, which reduces the SST, and thus humidity, until the atmospheric column again becomes stable to convection, precipitation ceases, and the skies clear, starting the cycle over. Surface flux modifications by wind-evaporation feedbacks contribute to the cycle in essentially the same way as cloud-radiative feedbacks, providing anomalous warming (in the generalized sense of moist static energy) of the atmosphere and equal cooling of the ocean during convective episodes.

We have explored the effect of varying several different parameters on the behavior of the system. Some of these parameters have analogs in GCMs, either as input parameters or as “effective” parameters that can be diagnosed from their output. Our results could in principle be used to aid in interpreting simulations of intraseasonal to subannual SST variability with GCMs coupled to ocean mixed layer models. The effect of varying mixed layer depth in particular has already been studied by Watterson (2002), who finds increased variability as the mixed layer depth decreases. In our model, the linear growth rate increases with mixed layer depth. However, in the stable regime, if the model is forced by an imposed atmospheric disturbance with intraseasonal frequency, the amplitude of the response increases as mixed layer depth decreases until a maximum around 10–20 m (near the minimum value used by Watterson), then decreases again if the mixed layer depth is decreased further. It would be interesting to see whether this behavior can be reproduced in a GCM, by varying mixed layer depth down to very small values. At least inasmuch as a GCM is felt to be several steps closer to reality than the present model, this nonmonotonic dependence on mixed layer depth is a testable prediction that emerges from our results.

Due to its lack of horizontal structure, our model cannot simulate the MJO per se. Nonetheless, it may have some relevance to the question of the degree to which the MJO is a coupled mode. If the effective parameters that best characterize the real system (*r, **τ*_{c}, etc.) were to place our model system in an unstable regime, we would expect that the coupled feedbacks in our model would play a large and destabilizing role in MJO dynamics, while we would expect a smaller role if those parameters were to place our system in a stable regime. If the real system is near the stability boundary, small changes in those parameters can be expected to yield large differences in simulated MJO behavior. It is possible that this sensitivity, in addition to the well-known difficulty of physical parameterization overall, could be part of the reason for the notorious difficulty GCMs have in simulating the MJO.

In any case, simple models such as ours may be useful in interpreting GCM results. Since the key parameters in our model relate different observable quantities to one another (e.g., precipitation and longwave cloud-radiative forcing), the effective parameters of a GCM can be estimated directly from scatterplots of those variables versus one another, with each point a space–time location under suitable averaging. We may hope that some insight may be gained by locating a GCM (or better, a number of different GCMs with different simulated intraseasonal variability) in the parameter space of a simple model whose behavior is relatively easily understood.

## Acknowledgments

We thank Chris Bretherton, Harry Hendon, Eric Maloney, Brian Mapes, David Neelin, and Duane Waliser for discussions, and Duane Waliser and an anonymous reviewer for insightful reviews leading to improvements in presentation. AHS acknowledges support for this work from the National Science Foundation under Grant DMS-01-39830 and a fellowship from the David and Lucile Packard Foundation. HG is supported by the NOAA Postdoctoral Program in Climate and Global Change, administered by the University Corporation for Atmospheric Research. Any opinions, findings, and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the sponsors.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Validity of WTG Approximation for Small Spatial Scales

We consider the linear system from section 3a, but add a temperature tendency term to (11) so that it becomes

where *â* ≈ 0.4 (Neelin and Zeng 2000). We also add a linear, zonal, inviscid momentum equation, which in the QTCM is

where *κ* is the gas constant for air divided by its heat capacity at constant pressure. With these changes, the system now supports a convectively coupled gravity wave. We assume slab symmetry so that *δ* = ∂*u*/∂*x,* and assume all variables vary as exp(*σt* – *ikx*). The system can then be written as

Using (A3) to substitute for *T*′, (A4) and (A5) become

Note that (A6) contains no *x* derivatives, and so does not need to be considered here. It is clear that for sufficiently large *k,* which means that

the first term in (A7), will become negligible relative to the second. If we do neglect that term, the remaining system becomes identical to the one considered in this study, if we again write *iku* as *δ.* That system contained no spatial derivatives and hence no dependence of *σ* on *k* whatsoever, so our neglect of atmospheric temperature variations is a consistent approximation in the limit of large *k.*

The expression (A9) physically means that the spatial scale of the WTG disturbances must be small compared to that of dry gravity waves with equal periods, since the phase speed of dry gravity waves in this system is easily shown to be *c* = . For the parameters used here *c* ≈ 50 m s^{–1}, so for *σ* = 0.2 day^{–1}, the wavelength of the disturbance must be short compared to a number on the order of 100 000 km, easily satisfied by any disturbance on earth. In reality, rotation sets a more stringent limit (Sobel et al. 2001; Majda and Klein 2003), but disturbances below the equatorial deformation scale should still be well described by WTG.

## Footnotes

*Corresponding author address:* Dr. Adam H. Sobel, Department of Applied Physics and Applied Mathematics, Columbia University, 500 West 120th St., Rm. 217, New York, NY 10027. Email: ahs129@columbia.edu