the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A simple model linking radiative–convective instability, convective aggregation and large-scale dynamics
Peter Haynes
A simple model is presented which is designed to analyse the relation between the phenomenon of convective aggregation at small scales and larger-scale variability that results from coupling between dynamics and moisture in the tropical atmosphere. The model is based on single-layer dynamical equations coupled to a moisture equation to represent the dynamical effects of latent heating and radiative heating. The moisture variable q evolves through the effect of horizontal convergence, nonlinear horizontal advection and diffusion. Following previous work, the coupling between moisture and dynamics is included in such a way that a horizontally homogeneous state may be unstable to inhomogeneous disturbances, and, as a result, localised regions evolve towards either dry or moist states, with divergence or convergence respectively in the horizontal flow. The time evolution of the spatial structure of the dry and moist regions is investigated using a combination of theory and numerical simulation. One aspect of the evolution is a spatial coarsening that, if moist regions and dry regions are interpreted as convecting and non-convecting respectively, represents a form of convective aggregation. When the weak temperature gradient (WTG) approximation (i.e. a local balance between heating and convergence) applies, and horizontal advection is neglected, the system reduces to a nonlinear reaction–diffusion equation for q, and the coarsening is a well-known aspect of such systems. When nonlinear advection of moisture is included, the large-scale flow that arises from the spatial pattern of divergence and convergence leads to a distinctly different coarsening process. When thermal and frictional damping and f-plane rotation are included in the dynamics, there is a dynamical length scale L_{dyn} that sets an upper limit for the spatial coarsening of the moist and dry regions. The f-plane results provide a basis for interpreting the behaviour of the system on an equatorial β plane, where the dynamics implies a displacement in the zonal direction of the divergence relative to q and hence to coherent equatorially confined zonally propagating disturbances, comprising separate moist and dry regions. In many cases the propagation speed and direction depend on the equatorial wave response to the moist heating, with the relative strength of the Rossby wave response to the Kelvin wave response determining whether the propagation is eastward or westward. Within this model, the key overall properties of the propagating disturbances, the spatial scale and the phase speed, depend on nonlinearity in the coupling between moisture and dynamics, and any linear theory for such disturbances therefore has limited usefulness. The model described here, in which the moisture and dynamical fields vary in two spatial dimensions and important aspects of nonlinearity are captured, provides an intermediate model between theoretical models based on linearisation and one spatial dimension and general circulation models (GCMs) or convection-resolving models.
- Article
(8419 KB) - Full-text XML
- BibTeX
- EndNote
Much theoretical and modelling work over the past few decades has focused on the coupling between dynamics and moisture in the tropical atmosphere, which has been argued should be taken into account at leading order to explain many tropical phenomena. Two topics that continue to attract significant attention are on the one hand convective aggregation, identified as a behaviour in numerical simulations in convection-representing models, and on the other hand the Madden–Julian Oscillation (hereafter MJO), identified in observations as a dominant mode of intraseasonal variability of the real tropical atmosphere. Convective aggregation (e.g. Wing et al., 2017; Muller et al., 2022) has been identified as a pattern of behaviour of a hypothetical tropical atmosphere in which there is no imposed spatial inhomogeneity but which exhibits spontaneous organisation of the circulation and convection into regions of two types, one type with active convection and large-scale ascent and the other with convection suppressed by large-scale subsidence. The relevance of convective aggregation to the behaviour of the real atmosphere remains a topic of debate, but the study of aggregation in convection-representing numerical models (hereafter CRMs) has provided a great deal of insight into the physics of the tropical atmosphere, particularly the interactions between convecting and non-convecting regions, and the way in which this physics is represented in models. It has been suggested that the same physics that is responsible for convective aggregation in numerical simulations is also part of the mechanism for the MJO (Bretherton et al., 2005; Arnold and Randall, 2015). Investigating this possibility in numerical simulations has been challenging because the numerical resolution for CRM simulations is a few kilometres or less, and the spatial structure of the MJO has scales of several thousand kilometres. A small number of papers (e.g. Arnold and Randall, 2015; Khairoutdinov and Emanuel, 2018) have described simulations that bridge this gap and have given insight into the relation between convective aggregation and the spontaneous generation of large-scale MJO-like disturbances. Other papers (e.g. Carstens and Wing, 2022, 2023) have used CRM simulations to explore more generally the relation between convective aggregation and larger-scale dynamics. Since such simulations are at the very edge of current computational capacity and the scope for thorough examination of parameter space is limited, it is desirable to find a simpler theoretical and modelling framework within which the link between processes such as aggregation on the mesoscale and larger-scale organisation of dynamics can be investigated further. The focus of this paper is the formulation and study of a simple model for such a purpose.
Several different physical processes have been proposed as being important for aggregation, often supported by the results of mechanism-denial experiments in CRMs in which the effects of particular processes have been altered or omitted altogether. As emphasised by Muller et al. (2022), there remains considerable uncertainty over which of these various descriptions of the aggregation process is most relevant to convective aggregation and, within each of them, the relative importance of different physical mechanisms. Examples of suggested mechanisms include the spatio-temporal propagation of triggering of convection through different effects, either through the formation and propagation of cold pools (e.g. Hirt et al., 2020) or through the propagation of gravity waves within the boundary layer (Yang, 2021). The study described in this paper is prompted in particular by the proposal that aggregation occurs via an instability of a spatially homogeneous radiative–convective equilibrium state resulting from feedbacks between moisture and radiation (Raymond, 2000; Emanuel et al., 2014), primarily in the free troposphere, though the boundary layer (Yang, 2018) may play a crucial role in the dynamics of these feedbacks. Instability, generally described as radiative instability or radiative–convective instability, may result if these feedbacks can overcome the dynamical stability of the convecting atmosphere. The competition between these processes is typically represented by the gross moist stability (GMS; e.g. Raymond et al., 2009), though changes in convective moisture and heat transport should be properly taken into account (Beucler et al., 2018). The growing instability is expected to ultimately saturate at finite amplitude with the result that locally the system tends towards a moist or a dry state. Such behaviour has been demonstrated using single-column radiation–convection calculations using various standard radiative and convective parametrisations and adopting the weak temperature gradient (WTG) assumption where the environmental temperature is specified, and the vertical mass flux is allowed to be non-zero (Sobel et al., 2007; Emanuel et al., 2014). Under suitable conditions the radiative–convective equilibrium (RCE) state is unstable, and the system evolves towards a moist state or a dry state, depending on the imposed initial perturbation.
A more generic and fundamental approach to the study of convective aggregation, starting with instability as the initial mechanism for the growth of moisture inhomogeneities on the homogeneous state but then seeking a description of the spatial evolution of the resulting moist and dry regions, has been undertaken by Craig and Mack (2013) and Windmiller and Craig (2019) (hereafter CMWC). The model system considered in this work is an evolution equation for a time-evolving two-dimensional moisture concentration field q, incorporating a source–sink term that is a nonlinear function of q, say G(q), with three zeros, each corresponding to possible steady states. The form of the function G(q) is such that the large-q (moist) and small-q (dry) states are stable, and the intermediate-q state is unstable. Transport of moisture is assumed to be diffusive. This system is equivalent to a reaction–diffusion equation with bistable reaction, sometimes known as the Allen–Cahn equation, which has been extensively studied using theoretical and numerical approaches. Such systems exhibit coarsening where, after the initial separation into high-q and low-q regions, typically on small scales, the scale of these regions increases monotonically with time, until it is constrained by the large-scale geometry imposed on the system. This is presented by CMWC as a mechanism for aggregation. The two essential components required for this mechanism to operate are the q dependence, and hence the bistable nature of the source–sink term, which CMWC argue results from the dependence of subsidence drying and convective moistening on free-tropospheric moisture, and the diffusive transport. Windmiller and Craig (2019) argue that diffusive transport can be justified on the basis of a simple model in which a stochastic convective cloud moistens the environment of a moist region, and they provide an estimate for the resulting diffusivity as 4×10^{2} m^{2} s^{−1}. A larger value of the diffusivity, 10^{5} m^{2} s^{−1}, is used in Craig and Mack (2013) and is envisaged as being based on the typical horizontal velocity and length scales of convective systems. The latter justifies this as an eddy diffusivity based on the typical horizontal velocity and length scales of convective motions. The appropriate value of the diffusivity therefore depends significantly on what the diffusivity is intended to represent. This reaction–diffusion model for aggregation is interesting, and it is discussed further in the Muller et al. (2022) review, but, since large-scale dynamics is omitted from the model, it cannot be used to investigate the link between convective aggregation and large-scale dynamics.
The model presented in this paper combines aspects of the CMWC model with a simple description of large-scale dynamics. The structure of this paper is as follows. In Sect. 2 we define the mathematical model to be studied, following previous approaches that use the shallow-water equations augmented by a prognostic equation for moisture, with the moisture coupling to the shallow-water dynamics. We then set out the behaviour expected of the model on the basis of previous work together with various scaling arguments, explaining the relation to CMWC. In Sect. 3 we present results from numerical simulations of the system on a doubly periodic domain. In particular, we verify that aggregation occurs and that the mechanism for aggregation can be dominated by horizonal diffusion of moisture or by horizontal advection of moisture, depending on the external parameters defining the system. In Sect. 4 we then show that rotation, thermal damping and frictional damping can each, or in combination, lead to a finite upper limit on the aggregation scale. This is with both theoretical arguments and results from numerical simulations on a doubly periodic f plane (including the zero-rotation case f=0). In Sect. 5 we exploit the f-plane results derived previously to consider the system on an equatorial β plane and show that the process of aggregation is then confined to a low-latitude region with the result of aggregation being the formation of coherently propagating disturbances. The scale and propagation speed of these disturbances depend on the external parameters, in particular, in some regimes, the relative strength of the equatorial Kelvin and Rossby wave responses to the moist heating. In Sect. 6 we discuss the results and present overall conclusions.
Whilst we have mentioned above the MJO as a key motivation for further investigation of tropical dynamics, the model as presented in this paper does not reach the stage of direct relevance to the MJO or other aspects of tropical intraseasonal variability. That requires further development to be reported in a future paper. We simply note that the model presented here can be regarded as belonging to the class of “moisture-mode” models studied by Sugiyama (2009a), Sugiyama (2009b), Sobel and Maloney (2012), Sobel and Maloney (2013), Adames and Kim (2016), and others but that this class is one of several different classes of theoretical models for the MJO that continue to be studied, as recently reviewed by Jiang et al. (2020) and Zhang et al. (2020).
2.1 Model equations
The model to be studied in the remainder of the paper follows much previous work in being based on the dynamical equations for a rotating shallow-water system, describing the evolution of horizontal velocity u and free-surface displacement h, augmented by a moisture variable q, which is transported by the horizontal velocity. Such single-layer equations can be derived from primitive equations and the corresponding moisture equation for a stratified three-dimensional atmosphere following the systematic procedure set out by Neelin and Zeng (2000) based on a vertical-mode decomposition for dynamical quantities and for moisture. Specific further simplifications, following e.g. Sugiyama (2009b), over the set of equations presented by Neelin and Zeng (2000) (their 5.1–5.4) are that the barotropic flow and the nonlinear terms in the dynamical equations are neglected.
To emphasise the overall structure of the equations, rather than the details determined by the modal decomposition, we simply write the dynamical equations in the standard form for the motion, linearised about a state of rest, of a shallow layer of fluid with undisturbed depth H:
and
where g is the gravitational acceleration, and f is the Coriolis parameter, which will be taken to be either constant f=f_{0}, corresponding to the f plane, or linearly dependent on the y coordinate, f=βy, corresponding to the equatorial β plane. k is the unit vector in the vertical. α is a linear friction coefficient and λ a thermal damping rate. The corresponding equation for the moisture variable q is
including both horizontal advective transport and horizonal diffusive transport, the latter with diffusivity κ, assumed to be constant. In these equations, u is to be interpreted as representing the lower-tropospheric horizontal velocity field, and the displacement h is a surrogate for mid-troposphere temperature, with temperature increasing as h decreases. Given typical choices for the vertical basis function for moisture (Zeng et al., 2000), q is to be interpreted as the lower-tropospheric moisture concentration. The values of the constants H, α and λ could be chosen to match the values implied by the modal decomposition and by corresponding approximations to detailed parametrisations of physical processes, as set out in Neelin and Zeng (2000) and Sugiyama (2009b), but for the purposes of the work reported here, we simply chose values typical of those chosen in previous work on simple modelling of the tropical atmosphere.
The term F_{h}(q), included on the right-hand side of Eq. (2), represents a moisture-dependent cooling term (cooling because of the relation between h and temperature), potentially including both latent heating and radiative heating. If cooling decreases with moisture, as is physically plausible, F_{h}(q) will then be a decreasing function of q. Correspondingly, the term F_{q}(q) on the right-hand side of Eq. (3) represents the combined effects of evaporation and precipitation. The part of F_{q}(q) representing precipitation (note that F_{q}(q) is a negative contribution to precipitation) is also expected, on the basis of the observed correlation between precipitation and moisture in the free troposphere (e.g. Holloway and Neelin, 2009), to be a decreasing function of q. Q is a suitably chosen constant so that q is the perturbation away from a background state where the “total” moisture variable is Q, and it is convenient to choose Q such that F_{q}(0)=0, i.e. such that in the background state there is a balance between evaporation and precipitation. Previous papers developing and exploiting the simplified set of equations, Eqs. (1), (2) and (3), such as Neelin and Zeng (2000), Zeng et al. (2000), and Sugiyama (2009b), have given detailed arguments for how these moisture-dependent terms F_{h}(q) and F_{q}(q) might be constructed from specific parametrisation schemes for radiation and convection. Here we choose simplified ad hoc forms for these terms. A specific simplification made is that both terms F_{q} and F_{h} are assumed to be independent of h; i.e. the evaporation–precipitation and the moisture-dependent part are assumed to be independent of temperature. Whilst temperature (or h) dependence of these moisture coupling terms is neglected in several previous papers (Adames and Kim, 2016; Sobel and Maloney, 2012, 2013), it is included by Sugiyama (2009a, b). We briefly examine the effect of such dependence, with results reported in Appendix E and further comments in Sect. 6.
Within the constraints of the very simple model specified above, we can identify the possibility of choosing H and Q such that $h=q=\mathrm{0}$ corresponds to a spatially homogeneous radiative–convective equilibrium (RCE) state with u=0 and F_{h} and F_{q} satisfying the conditions ${F}_{h}\left(\mathrm{0}\right)={F}_{q}\left(\mathrm{0}\right)=\mathrm{0}$. We can restrict the forms of F_{h}(q) and F_{q}(q) to those for which the system is stable to spatially homogeneous perturbations, which holds if $\partial {F}_{q}/\partial q<\mathrm{0}$, with the partial derivative being evaluated at q=0. A key question is then whether, within this restriction, the radiative–convective equilibrium (RCE) state is unstable to spatially inhomogeneous disturbances. We demonstrate below that such instability is possible for relatively simple choices of the functions F_{h}(q) and F_{q}(q), and we interpret this instability as the analogue, within this simple model, of “radiative–convective instability” that has previously been discussed in several papers and demonstrated in suitable single-column radiation–convection calculations (e.g. Sobel et al., 2007; Emanuel et al., 2014).
The key dimensional quantities that define the above system include g and H, which determine the dry gravity wave speed $c=\sqrt{gH}$; the horizonal diffusivity κ; the Coriolis parameter f; the thermal and frictional damping rates, λ and α respectively; a typical background value of moisture, say Q; and μ, an inverse timescale for the moist processes represented by F_{h} and F_{q}. It is convenient to take the dimensions of the moisture Q to be the same as those of the thickness H, and indeed this corresponds to a simple re-scaling of the parameters in F_{h}. To assess the importance of the advective term in Eq. (3), an additional dimensional quantity is needed which sets the magnitude of the spatially inhomogeneous part of the moisture field and, hence, via the leading-order balances operating in Eq. (2) and (3), the magnitude of the corresponding horizontal flow. These magnitudes are set by the nonlinear dependence of F_{q} and F_{h} on q, but it is convenient to choose the magnitude, say D, of the divergence ∇⋅u , as the relevant dimensional quantity. The reason for this becomes clear from the discussion in Sect. 2.2 below.
The fact that the magnitude of the advective nonlinearity arising in Eqs. (1), (2) and (3), as written above, depends on the choice of F_{h}(q) and F_{q}(q) potentially complicates interpretation of the behaviour of the system, and there are advantages to allowing the advective nonlinearity to be varied independently of this choice. We therefore introduce the parameter ϵ into Eq. (3) to give
The linear instability analysis of the RCE state, for example, takes ϵ=0 in the above equation. Note that in the form of the equations considered by Sugiyama (2009b), derived from the Neelin and Zeng (2000) Quasi-equilibrium Tropical Circulation Model (QTCM) equations, there is a distinct constant multiplying the nonlinear advective term. This constant is determined in principle in the derivation of the single-layer equations by the projection of a horizontal moisture advection term that varies in height onto the single basis function used to represent the moisture field, but it can also be conveniently varied as an independent parameter, and ϵ introduced here plays the same role as that parameter. We illustrate the role of advective nonlinearity by comparing ϵ=1 behaviour with ϵ=0 behaviour, but note that, for the reasons just given, ϵ=1 cannot be regarded as the only “correct” choice for including advective nonlinearity.
2.2 WTG and the relation to the CMWC reaction–diffusion system
A standard approach, particularly at low latitudes where f is small, to analysing the system defined by Eqs. (1)–(3) is to make the weak temperature gradient (WTG) approximation (e.g. Sobel et al., 2001). This neglects horizontal variation in h and can be justified provided that the horizontal length scale L satisfies $L/c\ll {T}_{q}$, i.e. that the timescale for gravity wave propagation through L is much less than the timescale T_{q} for moist processes. T_{q} could either be a timescale μ^{−1} set by an appropriate combination of F_{q} and F_{h} (see below) or be an emergent property of the system. Additionally, when damping and rotation are included, it must be the case that L≪L_{dyn}, where L_{dyn} is a dynamical length scale that is typically determined by c together with some combination of f, α and λ. We focus on the zero-damping case in this section and return to the dynamical effects of damping and rotation in Sect. 4.
Whilst under WTG h is constant in space, it may not be constant in time. Taking the spatial average of Eq. (2), the overline notation is used to denote the spatial average. It follows that
The spatially varying part of Eq. (2) then has the following form:
implying that ∇⋅u, and hence the irrotational part of the velocity field, is determined instantaneously by the moisture field q. Under the assumption $f=\mathit{\alpha}=\mathrm{0}$ in this section, the rotational part of u is constant in time. When provided with this initial rotational part of the flow, assumed to be zero for the purposes of this section, Eq. (4) becomes a self-contained equation for the evolution of the q field with the following form:
where the second equality defines the function G_{hq}. Note that whilst evaluation of $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$ requires knowledge of the q field, for the purposes of expressing the right-hand side of the equation as a function of q, $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$ is simply a parameter that appears in the definition of that function. The notation u[q] simply expresses the fact that at each instant u is determined completely, but non-locally, by the q field, through Eq. (6).
Neglecting for a moment the advection term ϵu[q]⋅∇q, this may be recognised as a reaction–diffusion equation of the type studied by CMWC. The difference is that whilst the nonlinear “reaction” term on the right-hand side of Eq. (7) was in CMWC's case entirely determined by the q dependence of precipitation and evaporation, in this case the reaction term is a combination of the moisture-driven heating/cooling F_{h}(q) and the moisture source/sink F_{q}(q). A further structural difference from the system considered by CMWC is the evolving quantity $\stackrel{\mathrm{\u203e}}{h}\left(t\right)$. The effect of this is felt by the system through the corresponding $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$ appearing in the definition of the reaction term. The reaction term is therefore not completely specified in advance as a function of q but contains the spatially constant term $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$, which also drives changes in $\stackrel{\mathrm{\u203e}}{h}$, as specified by Eq. (5). The Craig and Mack (2013) model, on the other hand, imposes an integral constraint on q, equivalent to specified domain-integrated precipitation, and then accommodates this constraint by multiplying G(q) by an unknown function of t, which is determined from the constraint. Again this means that the complete reaction term requires knowledge of the spatial distribution of q.
Simple theory of the reaction–diffusion system with the specified reaction term G(q) is that (i) homogenous steady states are possible with q equal to the constant value q_{s}, if G(q_{s})=0, and (ii) those homogeneous states are stable if ${G}^{\prime}\left({q}_{s}\right)<\mathrm{0}$ and unstable if ${G}^{\prime}\left({q}_{s}\right)>\mathrm{0}$. CMWC consider a bistable system with three possible values for q_{s}, ${q}_{-}<{q}_{\mathrm{0}}<{q}_{+}$ such that $G\left({q}_{-}\right)=G\left({q}_{\mathrm{0}}\right)=G\left({q}_{+}\right)=\mathrm{0}$ and $G\left({q}_{-}\right)<\mathrm{0}$, ${G}^{\prime}\left({q}_{\mathrm{0}}\right)>\mathrm{0}$ and ${G}^{\prime}\left({q}_{+}\right)<\mathrm{0}$. G^{′}(q_{0}) provides a useful definition of a reaction inverse timescale μ. The generic behaviour for a nonlinear reaction diffusion equation of this type is that locally q tends to one of the stable values, partitioning the domain into two regions, one with $q={q}_{+}$ and the other with $q={q}_{-}$, separated by interfaces of thickness $(\mathit{\kappa}/\mathit{\mu}{)}^{\mathrm{1}/\mathrm{2}}$. In the absence of rotation and damping, WTG will break down on length scales of order cT_{q}, so we require $c{T}_{q}\gg (\mathit{\kappa}/\mathit{\mu}{)}^{\mathrm{1}/\mathrm{2}}$; i.e. the reaction timescale μ^{−1} is such that κμ≪c^{2}. The initial geometry of these two regions is set by the initial conditions. A useful simple solution is a one-dimensional propagating reactive–diffusive wave solution with $q={q}_{+}$ on one side of the wave and $q={q}_{-}$ on the other. The speed of propagation of the wave, ${c}_{\mathrm{RD}}\sim (\mathit{\kappa}\mathit{\mu}{)}^{\mathrm{1}/\mathrm{2}}$, is determined by the form of the reaction function G(q). Defining V(q) by $\mathrm{d}V/\mathrm{d}q=G$ so that V(q) has turning points where G(q) has zeros, then if $V\left({q}_{+}\right)>V\left({q}_{-}\right)$ the region with $q={q}_{+}$ propagates into the region with $q={q}_{-}$. The corresponding result for the initial value problem, in one or more space dimensions, is that q tends everywhere to q_{+}. Similarly, if $V\left({q}_{+}\right)<V\left({q}_{-}\right)$, then q eventually tends everywhere to q_{−}. Only in the case of $V\left({q}_{+}\right)=V\left({q}_{-}\right)$, which applies in particular to the Allen–Cahn equation, do both regions $q={q}_{+}$ and $q={q}_{-}$ persist. Note that V(q) represents the area under the graph of G(q). To be precise, if the choice V(q_{0})=0 is made, then V(q_{+}) is the area under the graph of V(q) in the interval $[{q}_{\mathrm{0}},{q}_{+}]$, and V(q_{−}) is the corresponding (positive) area in the interval $[{q}_{-},{q}_{\mathrm{0}}]$.
An important effect in two dimensions is that the reaction–diffusion velocity c_{RD} becomes a local property of each point on each interface, depending not only on the form of G(q) but also on the curvature of the interface. The reaction velocity c_{RD} should be replaced by ${c}_{\mathrm{RD}}+\mathit{\kappa}/R$, where R is the (signed) radius of the curvature of the interface (such that the propagation is towards the interior of the curve). If $\left|{c}_{\mathrm{RD}}\right|<\mathit{\kappa}/\left|R\right|$, the velocity of the boundary may even change sign. Therefore, c_{RD} decreases as the curvature of the interface increases (R decreases) (Rubinstein et al., 1989). This tends to smooth out the boundary between moist and dry regions, as small-scale irregularities or indeed small-scale regions will tend to disappear. Larger moist regions can therefore expand, whilst smaller moist regions shrink. This is the standard coarsening behaviour, i.e. the geometric simplification of the geometry between the two regions through an increase in spatial scales, observed in reaction–diffusion systems (Bray et al., 2003).
The inclusion of an integral constraint on q by Craig and Mack (2013) is important because this ensures that even if $V\left({q}_{+}\right)\ne V\left({q}_{-}\right)$ (with V(q) defined as above), the system does not simply evolve to $q={q}_{-}$ or $q={q}_{+}$ everywhere. Both values of q persist as coarsening proceeds. Indeed, if this sort of constraint is not applied, then in most cases the reaction–diffusion system with a bistable reaction evolves everywhere towards one of the stable states. Indeed, Windmiller and Craig (2019) do not apply such a constraint and note that at long timescales one or the other of the q_{−} or q_{+} states occupies an increasingly large fraction of the domain. It is demonstrated below that the same property holds, with both values of q persisting as coarsening proceeds, for the model system being considered in this study, i.e. for a reaction–diffusion system with the reaction term as specified by the right-hand side of Eq. (7).
Following the arguments presented above, the stability of the spatially homogeneous RCE state will be determined by the derivative with respect to q of ${G}_{hq}(q,\mathrm{0})={F}_{q}\left(q\right)-(Q/H){F}_{h}\left(q\right)$ at q=0, with instability if the derivative is positive, provided that the domain size is large enough that diffusion does not stabilise the system through the action of the κ∇^{2}q term. (This derivative is later identified as proportional to the negative of the gross moist stability; i.e. the RCE state will be unstable if the gross moist stability is negative.) It is assumed that the derivative is indeed positive and furthermore that G_{hq}(q;0) is bistable in the sense that there are q_{+}(0) and q_{−}(0) such that ${q}_{-}\left(\mathrm{0}\right)<\mathrm{0}<{q}_{+}\left(\mathrm{0}\right)$, with ${G}_{hq}\left({q}_{-}\right(\mathrm{0});\mathrm{0})={G}_{hq}\left({q}_{+}\right(\mathrm{0});\mathrm{0})=\mathrm{0}$ and ${G}_{hq}^{\prime}\left({q}_{-}\right(\mathrm{0});\mathrm{0})<-\mathrm{0}$, ${G}_{hq}^{\prime}\left({q}_{+}\right(\mathrm{0});\mathrm{0})<\mathrm{0}$. Note that this property of G_{hq}(q;0) implies a similar property, with corresponding ${q}_{-}\left(\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}\right)$, ${q}_{\mathrm{0}}\left(\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}\right)$ and ${q}_{+}\left(\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}\right)$, for the more general right-hand side of Eq. (7), ${G}_{hq}(q;\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})$, provided that $\left|\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}\right|$ is not too large. For notational convenience the explicit dependence of e.g. ${q}_{-}\left(\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}\right)$ on $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$ is not displayed unless it is essential.
Numerical solutions below show that if G_{hq} is bistable in the sense defined, then the system indeed evolves towards two values of q and coarsening occurs. However some further insight can be obtained by assuming that after the initial adjustment the region $q={q}_{+}$ fills an area A_{+}, and the region $q={q}_{-}$ fills an area A_{−}. The area of the interfaces between the regions is assumed to be negligible. The configuration is therefore determined by the three unknowns, q_{−}, q_{+} and A_{+} (or A_{−}).
Then the above equations imply
These determine any two of the three variables q_{−}, q_{+} and A_{+} in terms of the third. For the system to be at a steady state an additional constraint, obtained by integrating Eq. (3) over the domain, might seem to be ${A}_{+}{F}_{q}\left({q}_{+}\right)+{A}_{-}{F}_{q}\left({q}_{-}\right)=\mathrm{0}$; however this can be deduced from Eqs. (8) and (9) above and provides no extra information. Therefore it has to be accepted that one piece of further information in addition to the above is required for a unique solution. In general these equations determine a state that is quasi-steady rather than exactly steady. After the initial adjustment, when the required extra information will be determined by the initial conditions (which might, for example, set A_{+}), we expect a further slow time evolution of the three variables and correspondingly of the geometry of the dry and moist regions.
Assume that the overall effect of this slow time evolution can be captured by the classical theory for one-dimensional reactive–diffusive waves, describing the propagation of the thin interfaces between regions of piecewise constant q. This suggests the (slow) time evolution equation
where L_{interface} is the length of the interface between the regions, and c_{RD} is the reaction diffusion velocity, with the convention that this is positive if the region with $q={q}_{+}$ propagates into the region with $q={q}_{-}$. L_{interface} will vary in time but is certainly positive. This equation allows a steady state when ${c}_{\mathrm{RD}}({q}_{+},{q}_{-},\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})=\mathrm{0}$.
That such a steady state exists can be deduced by considering the graphs of the relevant functions of q. As noted previously, for a given $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$ the reaction function is ${G}_{hq}(q,\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})={F}_{q}\left(q\right)-(Q/H){F}_{h}\left(q\right)+(Q/H)\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})$. Consider first the graph of ${G}_{hq}(q,\mathrm{0})={F}_{q}\left(q\right)-(Q/H){F}_{h}\left(q\right)$ as shown by the curve in Fig. 1, which intersects the q axis at q_{−}(0), q_{0}(0) and q_{+}(0). The value of ${G}_{hq}(q,\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})$ is represented by the vertical distance between this curve and the horizontal line $-(Q/H)\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$, shown in the figure for various values of $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$, varying between C_{min}<0 and C_{max}>0. For C outside this range, G_{hq}(q;C) no longer has three roots. Areas V_{+} and V_{−} are marked in the figure for a particular value of $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$. The condition ${c}_{\mathrm{RD}}({q}_{+},{q}_{-},\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})=\mathrm{0}$ is satisfied if and only if ${V}_{+}={V}_{-}$. It is clear from the figure that there is one value of $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$, say C_{s}, for which this holds, lying in the range (C_{min},C_{max}). Substituting this value into Eqs. (8)–(10) gives the corresponding values of q_{−} and q_{+}; the dry and moist values of q; and A_{+}, the fractional area occupied by the moist region.
A further question concerns the stability of this steady state. It is clear from Fig. 1 that c_{RD} is an increasing function of C (the area V_{+} increases and the area V_{−} decreases as C increases. Suppose that C>C_{s} so that ${V}_{+}\left(C\right)>{V}_{-}\left(C\right)$ and c_{RD} is positive; i.e. regions of q_{+} will propagate into regions of q_{−}. The consequence will be that the relative area occupied by q_{+} will increase, resulting in a decrease in $C=\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$, if F_{h}(q) is a decreasing function of q. Similarly, if C<C_{s}, then C will increase, indicating that the steady state C=C_{s} is stable. Note that the above arguments do not describe the process of coarsening but indicate that the two values of q persist, just as they do for the special case of the Allen–Cahn equation and for the system considered by CMWC, suggesting that coarsening is relevant.
We show examples in Sect. 3 to demonstrate that our model system, in regimes where the approximations leading to Eqs. (5) and (7) can be justified, naturally evolves to a piecewise constant configuration with both of the values of q, q_{−} and q_{+}, consistent with the relevant form of ${G}_{hq}(q,\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)})$.
The strict WTG form of the evolution equation, Eq. (7), suggests that for the system Eqs. (1), (2) and (4), locally q will tend to one of two values, q_{−} or q_{+}, as was the case for the CMWC pure reaction–diffusion system. The distinct additional feature of the system being considered here is that there is an associated pattern of convergence and divergence. The WTG balance in Eq. (2) implies that the divergence ∇⋅u will also tend to one of two values, ${D}_{-}=-{H}^{-\mathrm{1}}{F}_{h}({q}_{-},\stackrel{\mathrm{\u203e}}{h})={Q}^{-\mathrm{1}}{F}_{q}({q}_{-},\stackrel{\mathrm{\u203e}}{h})>\mathrm{0}$ or ${D}_{+}=-{H}^{-\mathrm{1}}{F}_{h}({q}_{+},\stackrel{\mathrm{\u203e}}{h})={Q}^{-\mathrm{1}}{F}_{q}({q}_{+},\stackrel{\mathrm{\u203e}}{h})<\mathrm{0}$ respectively. (Note that these estimates of convergence in moist regions and divergence in dry regions are the basis for the discussion above concerning the introduction of the parameter ϵ to control the magnitude of advective nonlinearity.) Assume that the corresponding values of ∇⋅u are $-{D}_{-}<\mathrm{0}<{D}_{+}$. Since area-integrated ∇⋅u is zero, we expect that the areas A_{−} and A_{+} filled by dry regions and moist regions respectively satisfy ${A}_{-}{D}_{-}\sim {A}_{+}{D}_{+}$.
This non-zero divergence has no consequence for the evolution of the system under the strict WTG approximation, but if this approximation is relaxed, as is likely to be required at large horizontal scales, then the coupling between moisture and divergence may lead to distinctly different behaviour from that predicted by reaction–diffusion alone. Furthermore, even if WTG can be justified in Eq. (2), the non-zero divergence will potentially be important if the nonlinear advection term ϵu[q]⋅∇q is included, since Eq. (2) implies an advecting velocity field that depends on the distribution of q. A typical velocity at length scale L will be U∼DL. The advective velocity becomes comparable with the reaction–diffusion velocity when $\mathit{\u03f5}DL\sim (\mathit{\kappa}\mathit{\mu}{)}^{\mathrm{1}/\mathrm{2}}$; i.e. $L={L}_{\mathrm{adv}}\sim (\mathit{\kappa}\mathit{\mu}/{D}^{\mathrm{2}}{\mathit{\u03f5}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$. Since moist regions are associated with convergence and dry regions with divergence, the effect of advection will be to reduce the areas of moist regions relative to those of dry regions. This suggests that in a steady state the reaction–diffusion velocity c_{RD} has to be positive rather than zero, i.e. that ${V}_{+}>{V}_{-}$ rather than ${V}_{+}={V}_{-}$. It also suggests that the values of q in both moist and dry regions are increased relative to their values without the advective nonlinearity. On scales larger than those of a single aggregated region, the combination of regions of divergence of opposite signs will generate a large-scale velocity field. Nonlinear advection will therefore cause aggregation when there is a high density of distinct convergent regions, on a timescale of $\mathit{\u03f5}/D$. The advective regime of aggregation becomes dominant for length scales L>L_{adv}. The consequence of this is illustrated in Sect. 3.2.
We can also now give a better estimate for T_{q} and hence the scale at which WTG will fail. In the linear instability phase a possible estimate is ${T}_{q}\sim {\mathit{\mu}}^{-\mathrm{1}}$. However in the nonlinear aggregation phase a potentially more relevant estimate is ${T}_{q}\sim L/{c}_{\mathrm{RD}}$; i.e. the timescale increases as the length scale increases. In this case where $\mathit{\alpha}=\mathit{\lambda}=f=\mathrm{0}$, then for WTG the first estimate would require $L\ll c{\mathit{\mu}}^{-\mathrm{1}}$. However the second estimate would require $L\ll cL/{c}_{\mathrm{RD}}$, suggesting that if c_{RD}≪c, then in the nonlinear aggregation phase the WTG description remains valid at all scales; i.e. aggregation simply proceeds until the length scale is the largest allowed by the geometry.
3.1 Model details
The model equations defined above in Eqs. (1)–(3) can be integrated numerically, and in this section we use this to confirm the previous results and further investigate model behaviour. Numerical details are given in Appendix A. Recall that the thickness variable h and the moisture variable q represent the departure from the spatially homogeneous RCE state. In all simulations reported below the initial condition is taken to be u=0, h=0 and q small with $\left|q\right|\ll \left|{q}_{\pm}\right|$. The values of q are chosen randomly at each grid point.
The theoretical discussion in the previous section does not depend on the precise form of the functions F_{h} and F_{q} appearing in Eqs. (2) and (3) respectively, requiring only that they together lead to a bistable moisture equation, Eq. (7). However, in order to investigate the behaviour numerically we need to define these forms explicitly. For illustrative purposes we use a piecewise linear construction, with
With ${\mathit{\mu}}_{\mathrm{2}}>{\mathit{\mu}}_{\mathrm{1}}>\mathrm{0}$ this represents larger effective latent heat release of precipitation near to the RCE state. Note that the key quantity $\mathit{\mu}=(\mathrm{d}/\mathrm{d}q){G}_{hq}(q;\mathrm{0}){|}_{q=\mathrm{0}}$ is equal to $-{\mathit{\mu}}_{\mathrm{1}}+{\mathit{\mu}}_{\mathrm{2}}Q/H$, which we write as $-{\mathit{\mu}}_{\mathrm{1}}(\mathrm{1}-{\mathit{\mu}}_{\mathrm{2}}Q/H{\mathit{\mu}}_{\mathrm{1}})=-{\mathit{\mu}}_{\mathrm{1}}M$, where M is a normalised gross moist stability of the RCE state. Throughout the paper we choose μ_{1} and μ_{2} such that M<0, implying that the RCE state is unstable. This simple formulation of F_{h} and F_{q} has the advantage that we can easily tune the locations of the fixed points q_{±} using the parameters q_{p} and q_{m}. Note that ${q}_{+}>{q}_{p}$ and ${q}_{-}<{q}_{m}$. Q and H are chosen such that $\mathrm{1}-Q/H>\mathrm{0}$, implying that the fixed points q_{p} and q_{m} are stable according to the analysis in Sect. 2.2. We choose $\left|{q}_{p}\right|>\left|{q}_{m}\right|$, corresponding to a more extreme value of moisture in moist regions than in dry regions. Together with the constraint of zero net heating in steady state, this implies small moist regions with strong upwelling and large dry regions with weak downwelling, as typically observed in convective aggregation (e.g. Muller and Bony, 2015).
Now that the equations are fully defined, the most basic starting point is to consider the system with no rotation or damping, $f=\mathit{\lambda}=\mathit{\alpha}=\mathrm{0}$. The computational domain is taken to be square with sides of length 10^{7} m. We initially take the other parameters in the system as $g=\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{2}}$, H=30 m, ${\mathit{\mu}}_{\mathrm{1}}^{-\mathrm{1}}=\mathrm{36}\phantom{\rule{0.125em}{0ex}}\mathrm{000}\phantom{\rule{0.125em}{0ex}}\mathrm{s}$, μ_{2}=3μ_{1}, $\mathit{\kappa}={\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, Q=15 m and ${q}_{p}/Q=\mathrm{0.1}$, ${q}_{m}/Q=-\mathrm{0.025}$. Note that the corresponding value of M for the RCE state is −0.5. This might be considered a relatively large negative value; e.g. the corresponding value in Sugiyama (2009b) is about −0.05, but the primary reason for this choice is to ensure that the RCE state is unstable, and this value has little effect on the evolution when the system has evolved away from that state. The results of a simulation with GMS of the RCE state increased to −0.05 are shown in Fig. 12. Overall, we argue that the precise values of these parameters are unimportant as here we aim to understand the general behaviour of the system, and indeed the parameters are varied throughout the paper; however the values are chosen to give similar moisture timescales to those deduced from the system studied by Sugiyama (2009b). We use the Craig and Mack (2013) value of $\mathit{\kappa}={\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ rather than $\mathit{\kappa}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{5}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ used by Sugiyama (2009b). The parameter values used in all of the two-dimensional simulations discussed in the paper are given in Table 1. See Sect. 4.3 for a brief discussion of the non-zero values of α, λ and f used later in the paper.
3.2 Simulation results
We first consider the case ϵ=0, without nonlinear advection of moisture. The evolution of the moisture distribution with time is shown in Fig. 2. The qualitative behaviour is similar to CMWC. There is an initial adjustment phase on the timescale μ^{−1} as the small-scale noise grows. During this initial phase, WTG applies only up to $L\sim {\mathit{\mu}}^{-\mathrm{1}}c$, about 10^{6} m for the parameters chosen. Hence distinct regions of enhanced and suppressed moisture on this scale or smaller form, with values corresponding to the effective stable fixed points q_{±} such that G_{hq}=0, $\partial {G}_{hq}/\partial q<\mathrm{0}$.
Once this has occurred, the coarsening process proceeds, with the scale of moist and dry regions slowly evolving, consistent with understanding of the reaction–diffusion system as discussed in the previous section. In particular, the evolution of the boundaries occurs on a slow timescale determined by the reaction–diffusion velocity c_{RD} and the smaller curvature-associated velocity $\mathit{\kappa}/R$, both of which are smaller than the gravity wave speed c. Hence the WTG approximation continues to hold as the scale of moist and dry regions increases. In this regime the proportion of the domain filled by each of the moist and dry states is changing, so there is a slow evolution of the mean heating $\stackrel{\mathrm{\u203e}}{{F}_{h}\left(q\right)}$. This causes a slow change in the locations of the stable moisture fixed points q_{±}; however this is on a longer timescale than that determined by G_{hq}, so the q distribution quickly adjusts to the new stable values. These features of the long-term evolution of the moisture distribution, the rapid adjustment to values of q close to q_{±} and the subsequent slow evolution of the values of q_{±} are shown in Fig. 3. This diffusive growth proceeds to the domain scale, when the areas of the regions are such that there is net zero heating and precipitation, and, consistent with theory (Rubinstein et al., 1989), the length of the boundary is minimised (forming either a circular or a band-shaped structure, depending on the geometry of the computational domain). At this point a steady state has been reached.
We now briefly consider the system with advective nonlinearity in the moisture equation, ϵ>0. As discussed and noted in Sect. 2.2 above, we expect this to lead to a reduction in the spatial scale of moist regions, an increase in the spatial scale of dry regions and a distinct advective mechanism for aggregation. The simulations show that the effect of the latter is that the system evolves towards a steady state, as was the case for ϵ=0. However advection changes the geometry of this final steady state, as can be seen from the steady-state distributions for different values of ϵ shown in Fig. 4. Unlike with ϵ=0, where the steady state is governed only by reaction–diffusion, the final state is now governed by a balance between reaction–diffusion and advection. The square periodic domain considered in this paper leads to the moisture forming either a band shape or a cross shape. At smaller values of ϵ (ϵ=0.5 is shown), a band forms, and at larger ϵ (ϵ=1 is shown), a cross is preferred. However at intermediate values either shape, according to details of the initial conditions, may be reached and persist. The changes in geometry associated with increasingly strong advection shown here may be relevant to realistic atmospheric flows even if the precise steady-state configurations are not. We emphasise that, as is common to a large class of nonlinear reaction- and diffusion-type problems (e.g. Rubinstein et al., 1989), these steady-state configurations are almost certainly strongly influenced by the domain geometry.
When frictional and thermal damping and rotation are included in the system, then, as noted previously, there will be an upper limit L_{dyn}, depending on the dynamical parameters α, λ and f, on the scale to which WTG balance can apply. It is expected that the coarsening to the domain scale exhibited in the previous section will be substantially modified, and perhaps halted, when the scale L_{dyn} is reached. Then in Sect. 4.2 we present semi-quantitative scaling arguments that are potentially relevant to the evolution beyond the linear instability phase. We then present results from numerical simulations in Sect. 4.3, both for ϵ=0 and for ϵ>0, including a regime diagram that summarises the overall pattern of behaviour as the dynamical parameters vary.
4.1 Key insights from the linear instability problem
Since the growth of disturbances to the RCE state is the precursor to the coarsening phase, we begin by considering an explicit solution of the linear instability problem when α, λ and f are non-zero.
There have been previous theoretical studies (e.g. Adames et al., 2019) of linear wave propagation and linear instability in systems equivalent to Eqs. (1)–(3), but it is useful to establish some of the basic properties of the particular model system that we consider in this paper. For the case discussed in Sect. 2.2 and illustrated in Sect. 3, where $\mathit{\alpha}=\mathit{\lambda}=f=\mathrm{0}$ and the WTG approximation is valid, the linear stability properties of the RCE state are very straightforward and are determined by the sign of the derivative with respect to q of G_{hq}(q,0). We consider the linear stability problem in more detail for α, λ and f non-zero values since this gives insight into the behaviour of the full nonlinear system as revealed by numerical simulation. Since the f plane is isotropic, we can assume that perturbations vary only in the x direction. Assuming small-amplitude perturbations of the form $u=\mathrm{\Re}\mathit{\left\{}\widehat{u}\mathrm{exp}\right(\mathit{\sigma}t+ikx\left)\mathit{\right\}}$, with analogous notation for other variables, the linearised forms of Eqs. (1), (2) and (3) are
where ${\mathit{\mu}}_{\mathrm{2}}=-{F}_{h}^{\prime}\left(\mathrm{0}\right)$ and ${\mathit{\mu}}_{\mathrm{1}}=-{F}_{q}^{\prime}\left(\mathrm{0}\right)$, matching the notation used in Eqs. (13) and (12) respectively. As is standard, these define an eigenvalue problem, the solution of which leads to a dispersion relation in the form of a quartic equation for σ, given explicitly in Appendix B as Eq. (B1).
An important simple case is the strict WTG limit with $\mathit{\alpha}=\mathit{\lambda}=f=\mathrm{0}$. This may be considered directly by neglecting the $\mathit{\sigma}\widehat{h}$ and $-\mathit{\lambda}\widehat{h}$ terms in Eq. (16) and then substituting for $\widehat{u}$ in Eq. (17) to deduce, neglecting the κk^{2} term,
This expression for σ motivates the previously noted definition of the normalised gross moist stability for the moist shallow-water equations,
As noted previously, in this paper parameter values are chosen such that M<0, implying σ>0 and moisture-mode instability with an inverse timescale $\mathit{\mu}={\mathit{\mu}}_{\mathrm{1}}\left|M\right|$. Note that the WTG approximation applies at small scales, with $k\gg \mathit{\mu}/c$. At large scales, with $k\ll \mathit{\mu}/c$, the moisture adjusts on a timescale shorter than that of the dynamics, and hence a steady-state balance in the moisture equation (Eq. 17), neglecting the term $\mathit{\sigma}\widehat{q}$, is appropriate. (The large-k and small-k limits in this problem correspond to the moisture-mode and gravity-mode limits identified by Adames et al., 2019.) Using the moisture equation to eliminate the moisture dependence in the height equation, Eq. (16), then gives the standard shallow-water equations with the gravity wave speed adjusted from c^{2} to
This defines the moist gravity wave speed ${c}_{m}^{\mathrm{2}}={c}^{\mathrm{2}}M$ and implies unstable growth rather than propagation if M<0. M<0 can therefore be identified as a criterion for instability whether or not the WTG approximation is valid. Re-including the effects of moisture diffusion will potentially inhibit instability on length scales comparable to or smaller than $\sqrt{\mathit{\kappa}/{\mathit{\mu}}_{\mathrm{1}}}$.
Analysis of the full dispersion relation, Eq. (B1), given in Appendix B, shows that, if M<0, there is a real positive root for σ only if
Note that the final factor represents the square of a length scale which we identify in the following subsection as the dynamical length scale L_{dyn}. However this does not give complete information on when instability is possible because there may be non-zero complex conjugate roots for σ with a positive real part. Some analytical progress can be made in describing the dependence of roots on the different parameters, but the algebra is complicated. Further details are given in Appendix B.
The behaviour found is illustrated in Fig. 5, which maps out different regions of the (α,λ) plane for a specific choice of κ and for six different choices of f, including f=0. Other parameters, μ_{1}, μ_{2}, c^{2}, Q and H, take the same values as specified in Sect. 3. The (α,λ) plane may be divided into four regions, each corresponding to a different regime of behaviour. Regime I is where ℜ(σ)>0 occurs for some k only when σ is real. Regimes IIa and IIb are where there are complex σ with ℜ(σ)>0 and with non-zero imaginary parts, and Regime III is where there is no instability, i.e. ℜ(σ)<0 for all k. The distinction between IIa and IIb is that in IIb, σ corresponding to the fastest growth over all k has a non-zero imaginary part.
For f=0 only Regimes I and III are present, and the region of instability simply corresponds to Eq. (21). For non-zero values of f that are not too large, the boundary between Regimes I and III is again described by Eq. (21), but, Regimes IIa and IIb exist in some part of the region α<f. Note that Regimes IIa and IIb are therefore confined to smaller and smaller values of α as f→0. For the largest value of f shown, again only Regimes I and III are present, with Regime I corresponding to Eq. (21). It can be shown in general that if the inequality
is satisfied, there are no Regimes IIa and IIb, and the transition between stability and instability is fully described by Eq. (21). For the parameters used to generate Fig. 5, the disappearance of regions IIa and IIb between $f=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ s^{−1} and $f={\mathrm{10}}^{-\mathrm{4}}$ s^{−1} is consistent with Eq. (22).
Whilst the parameter dependence that is found in the linear stability problem is complicated, two general rules that seem to hold are, first, increasing κ (unsurprisingly) inhibits instability, and, second, for a given κ, increasing f also tends to inhibit instability; i.e. the region of the (α,λ) plane in which there is instability reduces as these parameters increase. In particular, for any specified non-zero values of α, λ and κ, there is an f_{stab} such that f>f_{stab} implies stability. The solid curves shown in Fig. 5 bounding Regime III are therefore, for the chosen value of κ, contours of the function f_{stab}(λ,α) in the (λ,α) plane. The existence of f_{stab} is exploited in the description of the equatorial β-plane behaviour in the following section. (Note in particular that f_{stab} depends on κ as well as on α and λ, but we leave the dependence on κ unstated since in the simulations discussed later, κ is in practice kept fixed, and only α and λ are varied.)
4.2 Dynamical arguments
The description above is focused on the linear instability problem and cannot be assumed to extend to the evolution once the growing unstable disturbances have saturated, e.g. in an aggregation phase. More general insight into the evolution can be obtained by considering possible balances in the equations at horizontal scale L. Assume that the timescale of evolution is T_{q}. In the linear instability phase ${T}_{q}\sim {\mathit{\mu}}^{-\mathrm{1}}$. However after the unstable growth has saturated, T_{q} may be larger than this. For example, in the aggregation of moist and dry regions described previously, T_{q} is determined by the diffusivity κ and is large if κ is small.
If α≠0 and λ≠0, then, assuming that ${T}_{q}\gg {\mathit{\alpha}}^{-\mathrm{1}},{\mathit{\lambda}}^{-\mathrm{1}}$, a quasi-steady-state balance is possible in the dynamical equations; i.e. u is instantaneously determined by the q field according to the quasi-steady balance,
where δ is divergence, and ζ is vorticity. In this respect we can identify this case with Regime I in the previous section. Eliminating ζ gives
Substituting into Eq. (25) implies that the local h, and hence the local u, is determined by q in a surrounding region of scale:
This defines a dynamical length scale L_{dyn}. WTG balance, the local balance between divergence and heating, applies only length scales smaller than L_{dyn}. On length scales larger than L_{dyn}, the dominant balance in Eq. (25) is between λh and F_{h}(q), and hence, from Eq. (26) the divergence is proportional to ∇^{2}F(q) rather than to F(q). Another implication of the above balance, from Eq. (24), is that the flow will be dominated by the rotational component if f≫α and by the irrotational component if f≪α.
The quasi-steady balance above cannot hold when α=0, when Eq. (26) would imply δ=0. However if q is to be maintained away from the q=0 steady state, which is known to be unstable, Eq. (4), neglecting the term multiplied by ϵ, requires δ to be non-zero. This suggests that a distinct dynamical argument is required when α is small, analogous to the distinct nature of Regime II discussed in the previous section. If f≠0, a possible balance assuming a small Rossby number, i.e. that the timescale of the evolution of the q field is much larger than f^{−1}, is
implying
that is, a form of the quasi-geostrophic potential vorticity equation with the potential vorticity changing due to the effects of heating F_{h}(q) and thermal damping. Thus, in this system there are two prognostic equations, one for potential vorticity and one for moisture. In this case it is the second term in the time derivative that corresponds to divergence so that WTG applies if $L\ll (c/f)min(\mathrm{1},(\mathit{\lambda}{T}_{q}{)}^{-\mathrm{1}/\mathrm{2}})$. Here the length scale appearing is the Rossby radius $c/f$, and the flow can be expected to evolve on the timescale λ^{−1} at large timescales. Note that in this case the rotational part of the velocity field will be stronger than the divergent part.
In both the above cases it appears that the relation between moist heating and divergence becomes non-local at a sufficiently large scale, L_{dyn}, when α is large enough to bring the dynamical balance to a quasi-steady state and $c/f$ when α is smaller. Therefore the aggregation behaviour seen previously is likely to be halted, or at least strongly modified, when these scales are reached. In Sect. 4.3 below the nature of this modification is examined by numerical simulation.
The scale L_{dyn} defined above decreases as f increases, suggesting that the scale of aggregated moist and dry regions will also decrease as f increases. Furthermore, the underlying instability of the RCE state requires WTG dynamics to apply, and L_{dyn} therefore also represents an upper limit on the scale of the instability. As f increases, L_{dyn} will reduce to the diffusion scale $\sqrt{\mathit{\kappa}/\mathit{\mu}}$, and the instability of the RCE state will disappear. This provides an estimate for f_{stab},
where Γ is a non-dimensional parameter. The expression Eq. (21), provided by the linear stability calculation, is consistent with this reasoning and provides an explicit expression for Γ as equal to $\sqrt{-M}/(\sqrt{\mathrm{1}-M}-\mathrm{1})$. For the case when α is small, corresponding to Regime II, we have no expression for L_{dyn} or for the boundary of stability between Regimes II and III in Fig. 5, so a corresponding analytic result has not been determined. We do expect a qualitatively similar situation where the system is stabilised by diffusion once the maximum length scale L_{dyn} becomes sufficiently small. However, the behaviour, as f changes, of the boundary between Regime IIb and Regime III shown in Fig. 5 is geometrically complicated. This suggests that an easily interpretable expression for the entire form of the boundary will be difficult to find.
4.3 Numerical simulations
The link between the maximum length scale L_{dyn}, determined by non-zero values of f, α and λ, over which WTG is expected to apply and the spatial scale of aggregation are now illustrated using numerical simulation. The same numerical scheme as in Sect. 3 is used. The parameters, μ_{1}, μ_{2}, c^{2}, Q and H, take the same values as specified in Sect. 3, unless otherwise stated. Additionally, f, α and λ are chosen so that there is linear instability, corresponding to Regimes I, IIa and IIb in Fig. 5. Recall that details of the parameter values chosen for these simulations are given in Table 1. The values of f are chosen for illustrative purposes ($f={\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{1}}$ corresponds to about 5°). Representative values of α and λ are chosen to be similar to those in other papers on large-scale tropical dynamics; e.g. Sugiyama (2009b) takes $\mathit{\alpha}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathit{\lambda}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{7}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, and Adames and Wallace (2014) take $\mathit{\alpha}={\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathit{\lambda}={\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, Adames and Kim (2016) take $\mathit{\alpha}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$ and $\mathit{\lambda}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, but we deliberately diverge from these representative values in some simulations to investigate the effect of varying α and λ in different ways. Note that the justification of values for linear friction coefficients, or indeed the inclusion of linear friction at all, in models for tropical circulation remains a topic of active discussion (e.g. Romps, 2014).
4.3.1 ϵ=0 (advective nonlinearity excluded)
We begin with ϵ=0, i.e. excluding nonlinear advection of moisture. A selection of time series of the moisture distribution of the system for different choices of f, α and λ is shown in Fig. 6, with each simulation corresponding to a row. In all cases the moisture field q was initialised with small-scale random noise. In the first case, in row A, corresponding to Regime I in the (α,λ) plane shown in Fig. 5, the system initially evolves similarly to the case without damping or rotation, with the formation of distinct moist regions, which evolve and enlarge through aggregation. However the aggregation does not proceed to the domain scale but halts at a smaller scale, with quasi-steady circular moist regions. This is as expected from the previous dynamical discussion. The inclusion of non-zero f, α and λ implies that WTG balance can hold only up to the scale L_{dyn} and aggregation halts at this scale. Other simulations with f, α and λ corresponding to Regime I in Fig. 5 show similar evolution.
Rows B and C in Fig. 6 correspond to Regime II in the (λ,α) plane, with B corresponding to IIa and C to IIb. For case B there is again an initial segregation and then an aggregation process, leading to distinct moist and dry regions at some finite scale. However, the long-time distribution is no longer stationary but continues to evolve in time (without there being any further systematic increase in scale). In case C, whilst there is segregation, there is no clear aggregation stage, and the moist and dry regions evolve in time in a manner that is more wave-like than that seen in case B.
To establish that the division of the (α,λ) plane, originally motivated by the linear instability properties, provides a useful guide on the behaviour of the ultimate nonlinear evolution, Fig. 7 repeats the depiction of the (λ,α) plane shown previously for a single value of f, with superimposed symbols indicating whether the nonlinear evolution was aggregated and quasi-steady, as in case A above; aggregated and unsteady, as in case B; or propagating, as in case C. Regime IIa (case B) has been split into two sub-regimes, with a slow regime corresponding to transitional behaviour in which aggregated regions form, but propagation is sufficiently slow so that these remain round.
A possible interpretation of the apparent relation between the properties of the linear instability problem and the evolution observed in the numerical simulations is as follows. In Regime I, as illustrated by simulation A, the linear instability behaviour is essentially that described by the WTG approximation, with the relevant unstable mode having real σ. Therefore the system evolves through the instability to the segregated state determined by the bistability. In Regime IIa, as illustrated by simulation B, the relevant unstable mode is similar to that in A, with σ real, and the process of segregation is correspondingly similar. However the existence of slower-growing propagating (complex σ) unstable modes at larger scales is relevant to the nonlinear evolution post-segregation (even if the linear instability modes themselves do not provide a complete description of the behaviour). (A reduced mathematical model describing the evolution of the segregated state might make this relevance clearer.) In Regime IIb, illustrated by simulation C, there is a pair of the fastest-growing modes with complex conjugate σ rather than a single mode with real σ, and the mechanism for growth is therefore completely different from that described by the WTG approximation. In fact, these modes are better regarded as moisture-destabilised inertial waves and are not moisture modes since they rely on the fact that the dynamics is not slaved to the moisture field. Consequentially the nonlinear evolution is not so clearly a segregation into the two states allowed by bistability and instead is better characterised as an evolving field of nonlinear moisture-inertial waves. Note that systematic propagation in case C and the clear anisotropy of the instantaneous q distribution visible in Fig. 6i are an indication of spontaneous symmetry breaking rather than any systematic anisotropy of the system as specified.
We now focus on the behaviour of the model with parameters chosen from Regime I, where there is aggregation to a quasi-steady state. The behaviour can be usefully summarised by using the spatial autocorrelation. The autocorrelation scale, L_{auto}, is defined as the minimum radius at which the spatial autocorrelation is a factor of $\mathrm{1}/e$ less than its maximum value. The time evolution of L_{auto} for simulations with various damping and rotation rates, as well as diffusivities, is shown in Fig. 8.
Cases (a) and (b) have $\mathit{\alpha}=\mathit{\lambda}=f=\mathrm{0}$, so, as previously demonstrated in Sects. 2.2–3, aggregation is expected to eventually proceed to the domain scale. The evolution of L_{auto} for both case (a) and case (b) is consistent with this expectation. For case (b), L_{auto} reaches a limiting value within the time period shown in the figure. For case (a), a limiting value is not reached, but L_{auto} continues to increase throughout the period shown. The difference between (a) and (b) can be explained by the fact that κ for (b) is 4× larger than that for (a), therefore recalling the established theory on reaction–diffusion systems noted in Sect. 2.2 that aggregation proceeds more rapidly as κ increases.
Other cases have non-zero values of α, λ and f. All these show an approach to a finite limiting value, indicating that aggregation ceases. A candidate value for the length scale at which this occurs is ${L}_{\mathrm{dyn}}=c(\mathit{\alpha}/\mathit{\lambda}{)}^{\mathrm{1}/\mathrm{2}}({\mathit{\alpha}}^{\mathrm{2}}+{f}^{\mathrm{2}}{)}^{-\mathrm{1}/\mathrm{2}}$. The corresponding values deduced from Fig. 8 are consistent with this expression in the sense that the ordering as parameters are changed is consistent with the expression. Note in particular that for a given α and λ, the scale is smaller with f>0 than it is with f=0.
4.3.2 ϵ>0 (advective nonlinearity included)
The effects of advective nonlinearity in the presence of non-zero f, α and λ are illustrated in Fig. 9. Note that each of these cases has (λ,α) corresponding to Regime I. We have previously noted that the effect of advective nonlinearity is in the two-dimensional case to give moist regions that are more filamentary than quasi-circular. (For example, recall the steady-state moisture distributions, where aggregation has proceeded to the domain scale, for $\mathit{\alpha}=\mathit{\lambda}=f=\mathrm{0}$ shown in Fig. 4.) The effect of non-zero α, λ and f is both to limit any aggregation to a finite scale and to determine the flow pattern resulting from the moisture distribution. When f≠0 this flow has a substantial rotational component, and the advective effect of this on the filamentary moisture structures is apparent in Fig. 9a and b. The example with f=0, shown in Fig. 9c, where the advecting flow is irrotational, is distinctly different, with any curvature of the filamentary structures being weaker and resulting from deformation by a spatially structured irrotational flow. As has been noted previously, advective narrowing of moist regions means that the maximum magnitude q is affected by diffusion and not simply equal to the predicted value q_{+}.
In this section we consider the model on the equatorial β plane, i.e. with f=βy. It has been shown that on the f plane aggregation tends to be inhibited by rotation in two ways: (i) the upper limit of the scale for the underlying instability of the system is a decreasing function of f, and the lower limit is an increasing function of κ; therefore when κ is non-zero, the instability disappears altogether for $f>{f}_{\mathrm{stab}}=\mathit{\beta}{y}_{\mathrm{stab}}$, with the latter equality defining y_{stab}. (ii) The upper limit on the aggregation scale is a decreasing function of f. This suggests the possibility on the β plane of disturbances largely confined to some equatorial band with $\left|y\right|<{y}_{\mathrm{stab}}$. Such disturbances do indeed form, and we proceed to describe their behaviour.
The regime diagram shown in Fig. 5 in Sect. 4.1 provides some insight into how the dynamics might vary with latitude. At the Equator, with f=0, either Regime I or Regime III must apply, with Regime III implying that the RCE state is stable. For large f Regime III applies. Whether or not the transition from Regime I to Regime III passes through Regime II will be determined by the values of α and λ. Generally speaking this will occur when $\mathit{\alpha}/\mathit{\lambda}$ is relatively small. If Regime II is encountered then this is likely to manifest as more complicated behaviour (recall Fig. 6) as y_{stab} is approached. However we see later in this section that there are effects on the β plane that are not captured by the f-plane behaviour as described in Sect. 4.
It has been noted in the previous section that, on the f plane, whilst there was sometimes evidence of a selection of a preferred direction (recall Fig. 6c), this selection is purely random. On the β plane, however, there is a genuine east–west asymmetry, e.g. as manifested in the well-known Matsuno–Gill steady response to localised heating, which has been generalised by Wu et al. (2001) to the case where α and λ are not equal. It is of particular interest to determine whether this leads to zonal propagation of moist and dry regions and how such propagation varies with model parameters.
We begin this section by discussing an adjustment to the previous constant-f dynamical arguments and its impact on the local distribution of aggregated regions. We then describe the effects of the larger-scale equatorial circulation. Following this, we discuss the behaviour of a series of numerical experiments, in both the ϵ=0 and the ϵ>0 cases.
5.1 Implications of equatorial β-plane dynamics
Much of the scale analysis of the f-plane equations presented in Sect. 4.2 was based on a quasi-steady balance in the dynamical equations, which led to the relation (Eq. 26) between δ and h and hence an estimate ${L}_{\mathrm{dyn}}=c(\mathit{\alpha}/\mathit{\lambda}{)}^{\mathrm{1}/\mathrm{2}}({\mathit{\alpha}}^{\mathrm{2}}+{f}^{\mathrm{2}}{)}^{-\mathrm{1}/\mathrm{2}}$ for the scale on which WTG breaks down, serving as an effective upper limit on the scale of aggregation. The same approach of assuming a quasi-steady balance in the dynamical equations is now applied to the β plane. The scale L_{dyn} as defined previously is still useful but will now vary with latitude. It is convenient to use the notation L_{dyn,f} to represent the value of L_{dyn} for a particular value of f. At this stage it is also useful to note that the system being considered has no imposed inhomogeneity in x, i.e. in longitude, and there is therefore no systematic change in the character of the disturbances in x; i.e. there is representation of a “warm-pool” range of longitudes in which moisture has a stronger role in the dynamics than elsewhere. Correspondingly, any below-mentioned reference to integration over the domain refers to the entire domain and is not restricted to a particular range of x values.
The balance in Eq. (26) must be modified on the β plane because $\partial f/\partial y$ is non-zero and becomes
where $\mathit{\beta}=\partial f/\partial y$ and $\widehat{\mathit{n}}=({\mathit{\alpha}}^{\mathrm{2}}-{f}^{\mathrm{2}},\mathrm{2}f\mathit{\alpha})/({f}^{\mathrm{2}}+{\mathit{\alpha}}^{\mathrm{2}})$ is a unit vector. Substituting Eq. (33) into the thickness equation (Eq. 2) gives the steady response to a heating as
One important difference from the corresponding equations for the f plane is that there is now a preferred direction in the relation between δ and h, allowing a systematic anisotropy. A second difference is that the coefficients in the equations are now functions of y. A local analysis, treating f as constant, may therefore not always be valid. The expressions above suggest a change in the character of the system when the length scale is larger than $\mathit{\alpha}/\mathit{\beta}$. Below this scale, the first, isotropic term in Eq. (33) is dominant, and, furthermore, the coefficient appearing in this term does not vary significantly on this scale, so a local description will be valid. For length scales larger than $\mathit{\alpha}/\mathit{\beta}$, however, the second term is dominant, suggesting significant anisotropy in the dynamics. The coefficients and the vector $\widehat{\mathit{n}}$ will also vary significantly with f, so significant latitudinal variation is expected, and, more fundamentally, a local description may not be valid.
Along with the above it must be taken into account that aggregation is initiated by instability of the RCE state and that this instability is confined to $\left|y\right|<{y}_{\mathrm{stab}}$. Additionally, whilst it is expected that moisture anomalies are largely confined within this region, note that Eq. (34) implies that dynamical effects extend outside, on a length scale of ${L}_{\mathrm{dyn},f=\mathit{\beta}{y}_{\mathrm{stab}}}$. When M is large and negative, Eq. (32) suggests this scale is $\sqrt{\mathit{\kappa}/\mathit{\mu}}$: a diffusive response extending outside of the unstable equatorial region. When M is closer to zero, however, this length scale may be significantly larger.
The f-plane analysis predicts L_{dyn} as an upper limit on the scale for aggregation. This suggests ${L}_{\mathrm{dyn},f=\mathrm{0}}=c/\sqrt{\mathit{\lambda}\mathit{\alpha}}$ as the corresponding scale at the Equator on the β plane. Therefore, on the basis of the arguments above, aggregation at the Equator may be isotropic if ${L}_{\mathrm{dyn},f=\mathrm{0}}<\mathit{\alpha}/\mathit{\beta}$ and ${y}_{\mathrm{stab}}>{L}_{\mathrm{dyn},f=\mathrm{0}}$. If the second condition is not satisfied then the geometry will not allow isotropy. Furthermore, since ${L}_{\mathrm{dyn},f=\mathit{\beta}y}$ is a decreasing function of y, isotropy will extend to y=y_{stab}, but the characteristic length scales will decrease as $\left|y\right|$ increases. In other words, we expect the aggregation to be qualitatively similar to the f-plane case, except that it will be confined to a region slightly larger than $\left|y\right|<{y}_{\mathrm{stab}}$, and the characteristic length scale will vary with y. If ${L}_{\mathrm{dyn},f=\mathrm{0}}<\mathit{\alpha}/\mathit{\beta}$ and ${y}_{\mathrm{stab}}<{L}_{\mathrm{dyn},f=\mathrm{0}}$, on the other hand, then the aggregated structures in moisture will be largely confined to $\left|y\right|<{y}_{\mathrm{stab}}$ and therefore extended in x relative to their scale in y. However the dynamical signatures will extend to $\left|y\right|={L}_{\mathrm{dyn},f=\mathrm{0}}$. This suggests two distinct regimes of behaviour in a parameter space defined by ${y}_{\mathrm{stab}}\mathit{\beta}/\mathit{\alpha}$ and ${L}_{\mathrm{dyn},f=\mathrm{0}}\mathit{\beta}/\mathit{\alpha}=c\mathit{\beta}{\mathit{\alpha}}^{-\mathrm{3}/\mathrm{2}}{\mathit{\lambda}}^{-\mathrm{1}/\mathrm{2}}$. We label these Regime A (${L}_{\mathrm{dyn},f=\mathrm{0}}<\mathit{\alpha}/\mathit{\beta}$, ${L}_{\mathrm{dyn},f=\mathrm{0}}<{y}_{\mathrm{stab}}$) and Regime E (${L}_{\mathrm{dyn},f=\mathrm{0}}<\mathit{\alpha}/\mathit{\beta}$, ${y}_{\mathrm{stab}}<{L}_{\mathrm{dyn},f=\mathrm{0}}$). These two regimes are shown in the $({L}_{\mathrm{dyn},f=\mathrm{0}},{y}_{\mathrm{stab}})$ plane in Fig. 10. We add further regimes to this diagram and discuss them in more detail below.
We now consider the case where ${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$. The non-isotropic terms in Eqs. (33) and (34) must be taken into account. Furthermore, close to the Equator, a local analysis is no longer adequate. A natural y scale close to y=0, obtained by requiring a balance between the first and second terms in the middle expression in Eq. (34), is ${L}_{\mathrm{eq}}=(c/\mathit{\beta}{)}^{\mathrm{1}/\mathrm{2}}(\mathit{\alpha}/\mathit{\lambda}{)}^{\mathrm{1}/\mathrm{4}}$, as obtained by Wu et al. (2001) in their generalisation of the Matsuno–Gill problem. The corresponding x scale, however, agrees with the length scale from the local analysis at the Equator, ${L}_{\mathrm{dyn},f=\mathrm{0}}$. Note that ${L}_{\mathrm{eq}}=\sqrt{{L}_{\mathrm{dyn},f=\mathrm{0}}\mathit{\alpha}/\mathit{\beta}}$ and that L_{eq} therefore always lies between $\mathit{\alpha}/\mathit{\beta}$ and L_{dyn}, i.e. when ${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, $\mathit{\alpha}/\mathit{\beta}<{L}_{\mathrm{eq}}<{L}_{\mathrm{dyn},f=\mathrm{0}}$. The structure of Eq. (34) implies that a moisture anomaly localised within region $\left|y\right|<{L}_{\mathrm{eq}}$ will force dynamical anomalies that extend across the region $\left|y\right|\sim {L}_{\mathrm{eq}}$. This implies a distinct Regime D ( $\mathit{\alpha}/\mathit{\beta}<{L}_{\mathrm{dyn},f=\mathrm{0}}$, y_{stab}<L_{eq}) similar to Regime E, with latitudinally confined moisture anomalies driving broader dynamical anomalies.
Now consider $({L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, y_{stab}>L_{eq}). In this case, since y_{stab}>L_{eq}, the unstable region is large enough to allow multiple aggregated structures in the y direction. Since we have ${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, there is anisotropy at the Equator, but since L_{dyn} is a decreasing function of f, the anisotropic region is expected to extend only to y such that ${L}_{\mathrm{dyn},f=\mathit{\beta}y}=\mathit{\alpha}/\mathit{\beta}$; hence
which defines y_{iso}. For $\left|y\right|>{y}_{\mathrm{iso}}$ we expect isotropic aggregation. This implies two further regimes, Regime B (${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, ${y}_{\mathrm{stab}}>max\mathit{\{}{L}_{\mathrm{eq}},{y}_{\mathrm{iso}}\mathit{\}})$, in which there is anisotropic aggregation in $\left|y\right|<max\mathit{\{}{L}_{\mathrm{eq}},{y}_{\mathrm{iso}}\mathit{\}}$ and isotropic aggregation in $max\mathit{\{}{L}_{\mathrm{eq}},{y}_{\mathrm{iso}}\mathit{\}}<\left|y\right|<{y}_{\mathrm{stab}}$, and Regime C ( ${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, ${L}_{\mathrm{eq}}<{y}_{\mathrm{stab}}<{y}_{\mathrm{iso}})$, in which there is anisotropic aggregation across the entire region $\left|y\right|<{y}_{\mathrm{stab}}$.
Figure 10 shows all five of the regimes. It should additionally be noted that if y_{stab} is very small compared to dynamical length scales then instability may be inhibited. This might justify defining a further distinct regime, but since the priority has been to interpret behaviour seen in cases where there is instability, the criteria for such a regime have not been investigated in detail.
Having identified the different dynamical regimes that characterise the system on the β plane, we now consider the implications for spatial propagation of the aggregated moist and dry regions. It is useful to consider the moisture evolution equation. Using Eq. (25) to substitute for ∇⋅u, Eq. (4) becomes
This differs from the WTG form, Eq. (7), by the final term on the right-hand side, which is non-zero unless h is spatially uniform. Note that the contributions to q_{t} that arise under WTG are not expected to lead to systematic propagation since the relation between these terms and q is isotropic.
The extra term $(\mathit{\lambda}Q/H)(h-\stackrel{\mathrm{\u203e}}{h})$ can potentially cause systematic propagation if its relation to q, as expressed by Eq. (34), is anisotropic. As has been noted previously, this relation is isotropic on the f plane, implying no propagation in that case. Under circumstances where a local analysis of Eq. (34) is appropriate, it is useful to exploit the analogy between Eq. (34) and a damped advection–diffusion equation with a source term F_{h}(q). The advecting velocity is $gH\mathit{\beta}/({f}^{\mathrm{2}}+{\mathit{\alpha}}^{\mathrm{2}})\widehat{\mathit{n}}$. For a positive q anomaly, given that F_{h}(q) is negative, h will therefore be negative in the $\widehat{\mathit{n}}$ direction, and hence, according to Eq. (36), q_{t} will be negative in the $\widehat{\mathit{n}}$ direction and positive in the $-\widehat{\mathit{n}}$ direction, implying propagation of the q anomaly in the $-\widehat{\mathit{n}}$ direction. This direction is westward if f<α and eastward if f>α.
Local analysis of Eq. (34) applies for all y in Regimes A and E. These have been identified previously as isotropic at leading order, but there is weak anisotropy, and this will lead to westward propagation of aggregated moist regions close to the Equator. It follows that a sufficient condition for westward propagation at the Equator is ${L}_{\mathrm{dyn},f=\mathrm{0}}<\mathit{\alpha}/\mathit{\beta}$. In case A, the instability and hence aggregation may extend to y such that $\left|f\right|>\mathit{\alpha}$, i.e. $\left|y\right|>\mathit{\alpha}/\mathit{\beta}$, and there propagation will be eastward. The same analysis is relevant to regions such that $\left|y\right|>{y}_{\mathrm{iso}}$ in Regime B. Since this corresponds to $\left|f\right|>\mathit{\alpha}$, the propagation in these regions will be eastward.
Since the equatorial region $\left|y\right|<{L}_{\mathrm{eq}}$ in Regimes A, B and C is broader than $\mathit{\alpha}/\mathit{\beta}$, and the direction of the vector $\widehat{\mathit{n}}$ will change direction when varied across this region, the implication of the local analysis is that there will be a relative westward shift of the moisture anomalies near the Equator and a relative eastward shift elsewhere, generating a “<” shape. The overall propagation speed is likely to depend on the latitudinal extent of the equatorial moist region, i.e. on the ratio ${L}_{\mathrm{eq}}\mathit{\beta}/\mathit{\alpha}$. A wider region with a larger value of this ratio will have a larger proportion in which $-\widehat{\mathit{n}}$ is directed towards the eastward region shift and is hence likely to have greater eastward propagation.
The dependence of propagation speed on latitude suggests that when there are multiple aggregated structures in latitude, i.e. in Regimes A, B and C, there will not be a single propagation speed. If there is a dominant speed for structures at the Equator, then that will be different from and more westward than speeds for structures at higher latitudes.
An alternative way to understand the behaviour near the Equator in Regimes B, C and D is through decomposition into equatorial waves, following Wu et al. (2001). This offers a different approach to describing solutions of Eq. (34), for h given F_{h}(q), which does not require regarding the coefficients as locally constant. For reference, the method of calculation is set out in Appendix C. In this case the correction $(\mathit{\lambda}Q/H)(h-\stackrel{\mathrm{\u203e}}{h})$ appearing in Eq. (36) to the WTG convergence associated with a moisture-driven equatorial heating anomaly will be divergent responses to the east from Kelvin waves, acting to shift the moisture anomaly to the west, and to the west from Rossby waves, acting to shift the moisture anomaly to the east. There will also be a similar response, with convergence rather than divergence, associated with the negative moisture (i.e. dry) anomaly. It is the combination of these two wave responses that determines the net propagation of moist regions (and dry regions) in this system, with a stronger Rossby wave response implying eastward propagation. The fact that the Kelvin wave response decays away from the Equator on a scale L_{eq}, whereas the Rossby wave response can be excited at any latitudinal scale suggests that the Kelvin wave response will be relatively weaker, and therefore there will be eastward propagation if moist regions are significantly wider than L_{eq} (though it should be noted that it has previously been argued that L_{eq} is the natural latitudinal scale of aggregated moisture anomalies).
The arguments presented above can be used to formulate a quantification of contributions of different processes to zonal propagation that can be applied to the simulations discussed in the following section. This quantification is useful where the disturbances can be considered to be coherently propagating at a speed U, which is possible when the latitudinal structure of the disturbances is dominated by a single moisture anomaly centred on the Equator, corresponding to Regimes D and E. In a reference frame moving at constant speed U in the zonal direction, the moisture distribution will then be steady, and the moisture equation (Eq. 4) can be rewritten as
Multiplying by q_{x} and integrating over the domain imply that the propagation speed U is given by
This may be interpreted as an expression of the propagation speed as a sum of contributions from different individual terms in the moisture equation. We can exploit this, decomposing the divergence into a weak temperature gradient part, D_{WTG}, defined from Eq. (6), and its departure from WTG, D^{′}.
Since F_{q} and F_{h} are functions of q only and ${\stackrel{\mathrm{\u203e}}{F}}_{h}$ is constant, the corresponding terms in the integral are total derivatives and so vanish when integrated over a periodic domain. We can also write the diffusive term as ${q}_{x}{\mathrm{\nabla}}^{\mathrm{2}}q=\left(\frac{\mathrm{1}}{\mathrm{2}}\right[{q}_{x}^{\mathrm{2}}-{q}_{y}^{\mathrm{2}}]{)}_{x}+({q}_{x}{q}_{y}{)}_{y}$. The diffusive term therefore also does not contribute to the integral. The coherent propagation of moist regions can therefore be decomposed into parts due to the departure of the divergence from WTG and nonlinear advection:
Other terms in the moisture equation do not contribute. Note however that this does not mean that the propagation speed is independent of the value of κ, for example. Diffusivity still plays a role in setting the width of the unstable region $\left|y\right|<{y}_{\mathrm{stab}}$ and the shape of individual moist regions. The term including $-Q{D}^{\prime}$ can be further decomposed into separate contributions $-Q{D}_{\mathrm{Rossby}}^{\prime}$ and $-Q{D}_{\mathrm{Kelvin}}^{\prime}$ from Rossby and Kelvin wave parts respectively of the dynamical response to the heating implied by the moisture field. The method for calculating the Rossby and Kelvin wave responses is set out in Appendix D.
5.2 Numerical simulations
As with the f-plane case discussed earlier, much of the understanding of the behaviour on an equatorial β plane can be gleaned from numerical simulations. Numerical details are given in Appendix A, and details of the parameter values for simulations for which results are displayed are included in Table 1.
A particular focus in the analysis of the simulations is on the zonal propagation of disturbances. For each simulation a zonal phase speed U is estimated by identifying the dominant zonal wavenumber in the moisture field at the Equator and then tracking its time evolution. The estimated values U are noted in the figures showing results from the simulations. In some cases, where the latitudinal structure of the disturbances is dominated by a single moisture anomaly centred on the Equator, the disturbances can be considered to be uniformly propagating. When the latitudinal structure is more complicated, with multiple moist regions at different latitudes, the assumption of uniform propagation is not applicable, since different moist regions may propagate at different speeds. The estimated phase speed U then corresponds to the phase speed of moisture features along the Equator. Features further from the Equator tend to propagate eastward relative to those at the Equator. The expression, Eq. (40), which potentially diagnoses the mechanisms for propagation, is unlikely to be useful in cases with multiple phase speeds.
5.2.1 ϵ=0
We again start with the case where nonlinear advection is excluded. The results of a set of simulations for various diffusivities and damping rates such that L_{eq}<y_{stab} are shown in Fig. 11. In all panels the black horizontal lines denoting y=y_{stab}, calculated on the basis of the f-plane linear instability analysis, mark the boundary between the low-latitude region of instability of the RCE state and the higher-latitude region where RCE is stable. The early time evolution and upscale growth are similar to the constant-f cases, controlled by the value of κ (recall Fig. 8). The aggregation ceases at a scale determined by the dynamics, as discussed previously for the f plane in Sect. 4.2 and developed further for the β plane in Sect. 5.1. At the Equator, where f=0, we expect the scale of aggregated regions to follow $c/\sqrt{\mathit{\alpha}\mathit{\lambda}}$, up to quantisation by the domain size. This is consistent with the structure seen in Fig. 11a–c, where the predicted scale decreases by a factor of $\sqrt{\mathrm{3}}$ between each panel and the zonal scale decreases accordingly. Also, as expected, the scale of aggregated regions decreases away from the Equator as f increases. This is especially clear in panel (c) and also visible in panels (a) and (b). As in the f-plane cases, the timescale for upscale growth becomes very slow at large scales. The plots are all shown after 400 d of simulation, even though there may still be continuing systematic growth of the spatial structure at that time.
In Fig. 11a, b, d, e and f, the spatial scale of moist regions at the Equator is sufficiently large, with ${L}_{\mathrm{dyn},f=\mathrm{0}}>\mathit{\alpha}/\mathit{\beta}$, so that the anisotropic term becomes dominant in Eq. (33). According to the classification defined in Fig. 10, panels (a) and (b) reside in Regime B. Consistent with this, a region of non-isotropic aggregation forms near the Equator and extends over part of the unstable region, corresponding to $\left|y\right|<{y}_{\mathrm{iso}}$, though in panel (b) the anisotropic effect is weak. The anisotropic equatorial moist regions are locally shifted in the local direction of $-\widehat{\mathit{n}}$, in contrast to the quasi-circular structures seen at larger $\left|y\right|$ where isotropy applies. The increase in damping rates from panel (b) to panel (c), so that (c) lies in Regime A, decreases the scale of equatorial moist regions sufficiently so that the aggregation is isotropic across the entire unstable region.
Panels (d), (e) and (f) reside in Regime C, and the extension of the anisotropic region up to $\left|y\right|={y}_{\mathrm{stab}}$ is consistent with this. Panel (f) is an example of the system with increased diffusivity, and all other parameters are the same as those in panel (a). The main differences in (f) relative to (a) are the decrease in y_{stab} and the slight increase in scale due to the increased width of the diffusive boundaries between the moist and dry states. Panel (f) shows, in some sense, marginal behaviour between Regimes C and D. Whilst the aggregated regions near the Equator are confined to $\left|y\right|<{y}_{\mathrm{stab}}$, they are sufficiently wide so that no further disturbances form at larger y. The aggregation takes the form of a series of uniformly propagating regions, qualitatively similar to Regime D.
Evidence of the f-plane regimes can also be seen in the change in structure with latitude. Panels (a)–(c) and (f) all have a direct transition from Regime I to III at y=y_{stab}, and accordingly aggregated regions remain circular up to this boundary. However, panels (d) and (e) have an intermediate range of y for which the local f-plane behaviour is in Regime II. The moist and dry regions near the boundaries in these cases are no longer circular and are far more transient, similar in character to the structure seen in Fig. 6c.
A case with L_{eq}>y_{stab} is shown in Fig. 12. This has parameters chosen to lie in Regime D. The value of Q has been decreased relative to previous cases shown, to 10.5 m, and this has reduced the magnitude of the negative gross moist stability at RCE, hence reducing y_{stab} to be less than L_{eq}. The unstable region is then sufficiently narrow so that the equatorial wave response to the moist heating anomalies arising from the instability spreads into y>y_{stab}. The associated convergence drives moisture anomalies, and eventually a self-consistent balance between the dynamical and moisture fields is reached, with similar length scales for each, with the zonal scale given by $c/\sqrt{\mathit{\alpha}\mathit{\lambda}}$ and the latitudinal scale by L_{eq}.
In contrast to the more complicated cases shown in Fig. 11a–e, the structures shown in Figs. 11f and 12 can be regarded as uniformly propagating, and the decomposition expressed by Eq. (40) can be usefully applied. Given that ϵ=0, it is exclusively the Q∇⋅u term that is of interest, and, as noted previously and using a method set out in Appendix D, this can be divided into the non-local quasi-steady Rossby and Kelvin wave responses to the heating F_{h}(q). The separate Rossby and Kelvin wave contributions to the propagation speed of equatorial moist regions for these two cases, calculated on the basis of Eq. (40) with the convergence term split into the two contributions, are shown in Fig. 13. Panels (a)–(b) are the cases originally shown in Fig. 11f, and panels (c)–(d) are the cases in Fig. 12.
As expected, the propagation direction does indeed depend on the relative strengths of the Kelvin wave response, which gives divergence to the east of the moisture anomaly and hence the heating anomaly, and the Rossby wave response, which gives divergence to the west, with, in both cases, corresponding convergence associated with the dry anomaly. Note that the patterns shown in Fig. 13 are made up of both the divergence associated with the moisture anomaly and the convergence associated with the dry anomaly. Panels (a)–(b) have a stronger equatorial Kelvin wave response, and hence the moist and dry regions travel westward, whereas in panels (c)–(d) the Rossby wave response is stronger, and eastward propagation is observed. In both cases the contributions from each of the wave types are of a similar magnitude, and there is significant cancellation. This has two implications: the propagation speeds will tend to be smaller than expected from a dynamical scaling argument, and the propagation speed and direction are sensitive to the latitudinal structure of the moisture, since the details of that structure determine the relative strength of the Kelvin and Rossby wave responses. Note that the mechanism described above that associates the Rossby wave response with eastward propagation is very different to mechanisms that are relevant to a similar association in other recent papers. The papers by Yano and Tribbia (2017) and Rostami and Zeitlin (2019) describe a strongly nonlinear Rossby wave that propagates eastward as a vortex pair, relying on nonlinearity in the vorticity and hence the momentum equations, which is absent in our model equations. The paper by Hayashi and Itoh (2017) describes a kind of diabatic Rossby wave, in which the vorticity field organises convective heating such that the resulting vorticity forcing implies propagation to the east. This mechanism is allowed in principle by our equations but is not consistent with the fact that the evolution as simulated can be explained by a quasi-steady dynamical balance.
The arguments in Sect. 5.1 suggest that the parameter $c\mathit{\beta}/{\mathit{\alpha}}^{\mathrm{3}/\mathrm{2}}{\mathit{\lambda}}^{\mathrm{1}/\mathrm{2}}$ is important, with eastward propagation preferred when this is large. This is broadly consistent with the numerical simulations shown in Figs. 11 and 12, though this parameter certainly does not explain the propagation speed. For example, the system also shows a weak, complicated dependence of propagation speed and direction on diffusivity, with increasing κ both widening the boundaries between moist and dry regions and reducing y_{stab}. Regions which propagate to the east, seen in Fig. 11d and e, tend to have a spatial structure shaped like <. This structure, also noted in Sugiyama (2009b), may be explained by the same asymmetry of the Kelvin and Rossby responses that determine the direction of propagation. The Kelvin wave divergence to the east is necessarily localised near the Equator; however the Rossby wave response can be at any latitudinal scale. Hence, to the east of the moist region there is reduced moisture convergence near the Equator compared to further from the Equator; however to the west of the moist region there is reduced moisture convergence at all latitudes. There is therefore a tendency for off-equatorial (y>L_{eq}) regions of moisture to shift to the east. For a large moist region at the Equator extending into this region, an eastward tilt will be generated away from the Equator.
5.2.2 ϵ>0
A set of examples with ϵ=1 is shown in Fig. 14. The change in spatial structure of the moisture field from ϵ=0 to ϵ=1 is largely similar to that seen previously in the constant-f case (Fig. 6 versus Fig. 9). The associated convergence causes moist regions to narrow, and rotational flows at larger f lead to spirals in the moisture distribution.
Of the cases shown, (a), in Regime C, and (d), in Regime D, may be regarded as leading to a coherent steadily propagating disturbance (although the off-equatorial round regions in (a) are still unsteady). Panel (b) corresponds to Regime B and is closer to local quasi-isotropic aggregation, with features at different latitudes propagating at different speeds. This is highly unsteady, though the narrow propagation of the nonlinear advective convergence close to the Equator is preferentially in the zonal direction. Panel (c) resides near the boundary of Regimes B and C. The evolution is still very unsteady, with more evidence of rotational advection, particularly away from the Equator where f is non-zero.
For the steadily propagating cases, the same technique as before, based on Eq. (40), can be used to decompose the different contributions to propagation as shown in Fig. 15. The nonlinear moisture advection term ϵ∇⋅(qu) is now present. It is helpful to separate the contribution from any zonal mean zonal flow. This is present because the segregation–aggregation process, with the form chosen for F_{h}(q) and F_{q}(q), tends to lead to a systematic zonal mean latitudinal structure in q, with greater q at low latitudes. This leads to a corresponding structure in h and hence in u. This zonal mean structure is present when ϵ=0 but then does not have any advective effect on q. The zonal mean u usually takes the form of two off-equatorial jets, centred where the latitudinal gradient in h is the largest. This zonal flow can potentially have a strong effect on propagation, but in the simulations shown here, the generation of a strong zonal flow is avoided by choosing λ to be sufficiently large so that variations in h are small. In practice, this means $\mathit{\lambda}=\mathcal{O}\left({\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{1}}\right)$, well within the range discussed in the uniform rotation case.
For reference, the zonal mean u for the simulation in Fig. 15 is shown in Fig. 16. In this simulation the global mean u is 0.030 m s^{−1}. The damping still allows a non-zero zonal mean u; however this is significantly reduced. The domain mean u is constrained to be near zero to ensure that there is no advection by a uniform background flow.
Returning to the decomposition into different contributions to propagation shown in Fig. 15, it may be seen that the linear divergence term (which has not been separated into Rossby and Kelvin contributions in this case, but the fact that the corresponding ΔU is positive implies that the Rossby contribution dominates), the term associated with advection by the zonal mean flow and the remaining nonlinear advection term are all comparable in magnitude.
In this study we have presented a single-layer model for convective aggregation and its connection to large-scale dynamics. The linearised shallow-water equations, governing the dynamics, are augmented with a moisture equation. The moisture field affects the dynamics via a heating term, i.e. a forcing term in the thickness equation, and the dynamics affects the moisture through convergence alone (ϵ=0) or with the additional effect of horizontal advection (ϵ>0). The forms of the moisture-dependent precipitation term and the moisture-dependent heating term are such that under the WTG approximation, the spatially homogeneous state, in which both precipitation and heating are zero (the RCE state), is unstable, and the system is bistable with moist and dry stable states (interpreted as convecting and non-convecting respectively). In this regime and with ϵ=0 the behaviour is described by a nonlinear reaction–diffusion equation for q, which is very similar to that presented by Craig and Mack (2013) and Windmiller and Craig (2019) (CMWC). As discussed by CMWC, the system exhibits a well-known spatial coarsening that may be interpreted as a representation of convective aggregation.
The difference between our model system and that in the CMWC case is that in the latter the form of the reaction term is determined solely through the dependence of precipitation on q. In our case it is determined in addition by the q-dependent heating, which couples to the moisture equation through the dynamics. This allows us to include more general dynamical effects beyond the WTG approximation, including the impact of thermal and frictional damping, nonlinear advection, and rotation, including the extension to the equatorial β plane. Under WTG dynamics with ϵ=0, the system, similar to the one studied by CMWC, coarsens to the largest available scale. This remains the case when nonlinear advection (ϵ>0) is included; however the coarsening process is modified, and the end state is different. Thermal and frictional damping and f-plane rotation in combination set a dynamical scale L_{dyn}, which is an upper limit on the validity of WTG. The result is that coarsening proceeds only to this scale and then ceases. Depending on the relative values of thermal and mechanical damping and rotation, the final state may essentially be steady, or it may be unsteady, with, in some cases, a symmetry breaking leading to loss of spatial isotropy and to propagating structures. The nature of the linearly unstable modes provides some guidance on the type of the final state that is observed. If f is large enough, with the critical value f_{stab} depending on moisture diffusivity κ as well as other parameters defining the system, then the RCE state is stable.
Many aspects of the behaviour of a β plane may be interpreted in terms of the previously discussed f-plane behaviour. A latitude y_{stab} may be defined by f_{stab}=βy_{stab}. Disturbances resulting from the instability of the RCE state are largely confined within the latitudinal band $(-{y}_{\mathrm{stab}},{y}_{\mathrm{stab}})$, with some penetration outside of this band as a result of the non-locality of the dynamics, typically with a scale of L_{eq} or the local value of L_{dyn}, whichever is the largest. When thermal and mechanical damping is strong enough, specifically when $c\mathit{\beta}{\mathit{\alpha}}^{-\mathrm{3}/\mathrm{2}}{\mathit{\lambda}}^{-\mathrm{1}/\mathrm{2}}<\mathrm{1}$, the evolution of the system is similar to that observed on the f plane, with spatial modulation corresponding to the spatial variation in the value of f. The weak anisotropy introduced by the β effect leads to zonal propagation. The propagation is incoherent, with structures at different latitudes propagating at different speeds; however those on the Equator propagate to the west. When $c\mathit{\beta}{\mathit{\alpha}}^{-\mathrm{3}/\mathrm{2}}{\mathit{\lambda}}^{-\mathrm{1}/\mathrm{2}}>\mathrm{1}$, the structures that form in the q field are more strongly anisotropic and coherently propagating, and it is helpful to formulate a description of the dynamics in terms of equatorial Kelvin and Rossby waves. If the former dominates the dynamical response to the q anomalies, there is propagation to the west; if the latter dominates, there is propagation to the east. Nonlinear moisture advection enhances eastward propagation, but the enhancement is not large. Whilst the model equations allow for time-dependent dynamics, the propagating disturbances that are seen in simulations can essentially be explained by considering the quasi-steady velocity fields forced by moisture-determined heating and consequently the effect of those velocity fields through convergence/divergence and nonlinear advection on the moisture field, indicating that these disturbances are moisture modes.
The β-plane results that we present, particularly for the $c\mathit{\beta}{\mathit{\alpha}}^{-\mathrm{3}/\mathrm{2}}{\mathit{\lambda}}^{-\mathrm{1}/\mathrm{2}}>\mathrm{1}$ regime, have significant common ground with those presented by Sugiyama (2009b). The behaviour observed in both studies is similar, with small-wavenumber and slowly propagating moist regions forming at the Equator. Indeed, as noted previously, one might argue that we are re-examining aspects of the models developed by Sugiyama (2009a, b). These models include WISHE, which we do not, and have a physically derived formulation of the forcing terms F_{q} and F_{h}, which in some circumstances leads to bistability (Sugiyama, 2009a). We, on the other hand, emphasise the general implications of the bistability of the system, motivated by various lines of investigation, such as single-column radiative–convective calculations in WTG models (e.g. Sobel et al., 2007; Emanuel et al., 2014), and then choose very simple ad hoc forms of F_{q}(q) and F_{h}(q) to provide such bistability. We also deliberately trace the behaviour of the system through a sequence starting with WTG dynamics, and hence the reaction–diffusion behaviour discussed by CMWC, and finishing with the β plane with thermal and mechanical damping and nonlinear moisture advection that was the focus of Sugiyama (2009b), thereby providing a new perspective on the latter.
One difference between the conclusions of this study and those of Sugiyama (2009b) is that the latter suggests that diffusion is an important contributor to the eastward propagation of moist regions in the nonlinear advection case. We, however, conclude that the propagation is due to a combination of nonlinear advection and the displaced convergence associated with the Rossby wave response. This difference in interpretation arises in part from the diagnostic approaches used to measure the impact of different terms of the moisture tendency. We assume a uniformly propagating disturbance and then calculate a contribution from each moisture tendency term to the speed of propagation using Eq. (40), taking into account the entire spatial distribution. As noted previously, according to this approach the net contribution of the diffusive term to propagation is zero, whereas Sugiyama (2009b) compares the tendencies at y=0. Only comparing at y=0 will overestimate the contribution of the diffusive tendency to the overall propagation. In an aggregated region shaped like <, the negative curvature at the Equator to the east increases the diffusive speed of the boundary, whereas the opposite is true away from the Equator where the sign of the curvature changes. The fact that according to Eq. (15) the net contribution of diffusion to propagation is identically zero does not, of course, rule out the possibility that the form of the terms that do contribute to net propagation are affected by diffusivity. Certainly, there is evidence of diffusivity dependence of propagation speed in the results presented, for example in cases (a) and (f) in Fig. 11, where increasing diffusivity increases westward propagation speed, and cases (a) and (b) in Fig. 14, where increasing diffusivity increases eastward propagation speed. A further difference from Sugiyama (2009b) is that we have not assumed h dependence in F_{h} and F_{q}. As has been noted previously, this follows other moisture-mode models, including Sobel and Maloney (2012, 2013) and Adames and Kim (2016), but a preliminary assessment of the effect of including h dependence is given in Appendix E and concludes that the formation, evolution and propagation of moist and dry regions, as described in Sect. 5 above, remain broadly unchanged when this is included.
Recent review articles on convective aggregation, for example Wing et al. (2017) and Muller et al. (2022), include the CMWC reaction–diffusion model in their discussion of different mechanisms and models. The model we have presented above generalises CMWC by linking the moisture and large-scale dynamical equations. In the categorisation given by Wing et al. (2017), our model fits the category of a long-wave radiation feedback with an additional advective process feedback included if ϵ>0. The CMWC model, on the other hand, which relies on the form of the moisture dependence of precipitation, is categorised as a moisture feedback. Unlike the CMWC model, our model, by including dynamics, provides an upper limit on the scale of the aggregation scale, which is finite and determined by thermal and momentum damping rates when f=0 and reduces as f increases. How this relates to evidence from GCM and CRM simulations remains to be determined, although recent CRM studies show the potential for multiple aggregated regions in larger domains and more complex structures with a characteristic length scale (Yanase et al., 2022; Patrizio and Randall, 2019).
Whilst most high-resolution three-dimensional modelling of convective aggregation has focused on the non-rotating case, there has been some recent investigation of the f-plane and β-plane cases (Carstens and Wing, 2022, 2023). Some aspects of the behaviour reported in those papers are seen in our much simpler model. On the f plane, circular moist (i.e. convecting) regions form, and the scale of these regions decreases as f increases, scaling as $\mathrm{1}/f$ in both models at large rotation rates. However the physics of this behaviour is likely to be very different between the two models. In the CRM studies of Carstens and Wing (2022), the structure for larger values of f is dominated by tropical cyclones. These features have no clear analogue in our model, which neglects advective nonlinearity in the momentum equation. Carstens and Wing (2022) identify an intermediate range of f within which convective aggregation simply does not occur. One possibility is that this corresponds to our f>f_{stab}, with formation of tropical cyclones being a distinct process that occurs at a larger f value in the CRM but which is simply absent in our model. On the β plane, Carstens and Wing (2023) identify the dominant structures that arise from convective aggregation at low latitudes as convectively coupled Kelvin waves, whereas in our model the structures are clearly characteristic of moisture modes, with quasi-steady dynamical fields. In future work it would be interesting to investigate further whether there are parameter regimes in which our model also shows moisture-modified Kelvin waves as the dominant low-latitude structures.
One goal of formulating and studying the model presented in this paper was to provide a basis for understanding the MJO. But in several respects, with the model in its current form, this goal has not been met. One aspect of this is the phenomenon of aggregation, which might have provided an explanation for the large longitudinal scale of the MJO without requiring a scale-selection mechanism for an underlying instability, such as radiative–convective instability, which is consistent with such a large scale. However the rate of increase in scale due to aggregation is in our model determined by diffusivity of moisture and, with the value assumed, $\mathit{\kappa}={\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$, the time taken to reach scales of thousands of kilometres is (at least) many tens of days. This does not seem consistent with the observed MJO evolution in which convection appears to develop over a large, perhaps 10^{4} km region of the Indian Ocean on a timescale of a few days. A larger value of diffusivity would reduce the time required but might be difficult to justify – e.g. see the arguments in Biagioli and Tompkins (2023). It might be that aggregation acts alongside large-scale-selective processes such as the radiative destabilisation suggested by Adames and Kim (2016).
A second aspect is the horizontal structure and propagation characteristics of disturbances predicted by the model on an equatorial β plane. Our model in this form is similar to that considered previously by Sugiyama (2009b), and our conclusions are also similar. There is organisation into distinct moist and dry regions on the Equator, and these propagate in the zonal direction. Similar to Sugiyama (2009b), we emphasise that these propagating disturbances are fundamentally nonlinear, and their structure and propagation characteristics are not captured by a linear stability analysis. However, the direction of propagation can be westward or eastward, and the speed of propagation is consistently much less than 1 m s^{−1} when the observed MJO phase speed is around 5 m s^{−1}. Furthermore, even when the propagation is eastward-propagating disturbances, the spatial structure of eastward-propagating moist regions also necessarily forms a < shape compared to the > structure of the observed MJO (Adames and Wallace, 2014).
It seems likely that, for a model of this nonlinear moisture-mode form to produce more realistic MJO-like behaviour, further physics needs to be included (though it should not be forgotten that MJO-like behaviour can also be produced by models with a different form). Recent simple one-dimensional moisture-mode models for the MJO (Sobel and Maloney, 2013; Adames and Kim, 2016), which are based on a prescribed latitudinal structure of the flow variables, have included extra effects, such as synoptic eddy drying, boundary-layer convergence of moisture (Adames and Wallace, 2014) and latitudinal gradients in background moisture (Adames and Kim, 2016), and have shown that these effects lead to enhanced eastward propagation. However the results we have presented in Sect. 5 of this paper demonstrate that propagation speed and direction are sensitive to latitudinal structure, so there are disadvantages to prescribing this a priori. We currently include some of these effects in the two-dimensional model described in this paper with the aim of capturing the behaviour of nonlinear MJO-like disturbances in a model where the two-dimensional structure emerges rather than being prescribed. Such a model could also potentially be used to study a broader class of moist tropical variability (e.g. Wang and Sobel, 2022).
To solve Eqs. (1)–(3) numerically, we use standard methods for solving the shallow-water equations. The equations are discretised on an Arakawa C grid with the grid spacing chosen to be $d=\mathrm{4}\times {\mathrm{10}}^{\mathrm{4}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$. The system is then numerically stepped forward in time using a third-order Adams–Bashforth scheme with a time step of 112.5 s.
On the f plane, we take (doubly) periodic boundaries, whereas on the equatorial β plane, rigid north and south boundaries are used along with periodic east and west boundaries. At rigid boundaries, a sponge layer is included to avoid edge effects. This takes the form of a linear damping term included in the velocity and thickness fields. On the periodic channel we use on the β plane this takes the following form:
The magnitude decays over a scale of $\mathrm{1}/\mathrm{70}$th of the domain. Any numerical results discussed in this paper do not depend on the magnitude or decay scale of this damping term.
This appendix gives further details of the linear stability problem considered in Sect. 4.1. Considering small-amplitude perturbations of the form $u=\mathrm{\Re}\mathit{\left\{}\widehat{u}\mathrm{exp}\right(\mathit{\sigma}t+ikx\left)\mathit{\right\}}$, with analogous notation for other variables, the linearised forms of Eqs. (1), (2) and (3) lead to (14), (15), (16) and (17), with ${\mathit{\mu}}_{\mathrm{2}}=-{F}_{h}^{\prime}\left(\mathrm{0}\right)$ and ${\mathit{\mu}}_{\mathrm{1}}=-{F}_{q}^{\prime}\left(\mathrm{0}\right)$, matching the notation used in Eqs. (13) and (12) respectively. These define an eigenvalue problem for the growth rate σ, which is a root of the quartic equation
with coefficients that are either linear or quadratic in k^{2}. Simple limits of the above in the case where $f=\mathit{\alpha}=\mathit{\lambda}=\mathrm{0}$ have been noted previously, in particular in the WTG limit valid at small scales such that $k\gg \mathit{\mu}/c$ and the corresponding large-scale results for $k\ll \mathit{\mu}/c$.
To further analyse the behaviour of the full system as described by Eq. (B1), we first consider the case where α and λ are non-zero, but there is no rotation; i.e. f=0. It is helpful to note the solutions for σ in the large-scale (k small) and small-scale (k large) limits. Neglecting diffusion, it is straightforward to show that at small k the roots of Eq. (B1) are $\mathit{\sigma}\simeq -\mathit{\lambda},-\mathit{\alpha},-\mathit{\alpha}$ and −μ_{1} (i.e. the $\mathit{\sigma}\simeq -\mathit{\alpha}$ root is repeated in the limit as k→0), implying stability at small k. Correspondingly, at large k there are roots $\mathit{\sigma}\simeq \pm ick-\mathit{\lambda}-{\mathit{\mu}}_{\mathrm{2}}Q/H,-\mathit{\alpha},{\mathit{\mu}}_{\mathrm{2}}(Q/H)-{\mathit{\mu}}_{\mathrm{1}}=-{\mathit{\mu}}_{\mathrm{1}}M$, representing two gravity waves damped by moisture effects and thermal damping, a vorticity disturbance damped by friction, and a moisture mode (stable or unstable according to whether M is positive or negative) respectively. There is therefore stability at small k and, consistent with the WTG analysis, instability at large k if M>0, with the latter limit unaffected by α and λ being non-zero. Introducing diffusion by taking κ>0 will lead to all roots having negative real parts as k→∞; i.e. the system is stabilised by diffusion of moisture.
Given that all roots for σ have negative real parts for k→0, it is possible to deduce whether or not there is instability by seeking conditions, in particular a value of k, under which one of the roots has zero real parts. One possibility is σ=is, where s is real and non-zero. Substituting this form for σ into Eq. (B1), with f=0, shows that such a root is not possible. (A key simplification in the f=0 case that can be exploited here is that the fourth-order polynomial appearing in Eq. (B1) has a factor (σ+α).) The other possibility is that σ=0, implying the condition
This is a quadratic for k^{2}, which has real roots only if
where, as before, $M=\mathrm{1}-{\mathit{\mu}}_{\mathrm{2}}Q/{\mathit{\mu}}_{\mathrm{1}}H$. If this condition is satisfied, then there is instability. If it is not satisfied, then there is no instability. Furthermore, it may also be deduced that roots with positive real parts must be real and that for any value of k there can be at most one such root. Given that roots crossing the real axis must do so at σ=0, complex roots with positive real parts would be possible only if there were two such crossings, from below to above the real axis, as k increased, and those real roots then combined to give a complex conjugate pair. But this is not possible since the above shows that there are at most two crossings for any k, and, given that all roots are below the real axis for small k and (with κ>0) for large k, one crossing is from below the real axis to above, and the other is from above to below.
When f is non-zero the conditions for instability are more complicated. The large-k limits of the roots of Eq. (B1) are unaffected by the expressions given above (which were for κ=0; the same holds for κ>0). However the small-k limits are now $\mathit{\sigma}\simeq -\mathit{\lambda},\pm if-\mathit{\alpha},-{\mathit{\mu}}_{\mathrm{1}}$, implying stability at small k, but note that the previously repeated root −α now becomes the complex conjugate pair $\pm if-\mathit{\alpha}$ (i.e. frictionally damped inertial oscillations). The transition between small k and large k now depends on the values of α, λ, f, κ and other parameters. The expression (Eq. B3) for real roots of Eq. (B1) to cross the real axis generalises to (repeating Eq. 21)
however this does not give complete information on when instability is possible because with f non-zero complex conjugate roots may also cross the real axis. Some analytical progress can be made in describing the dependence of roots on the different parameters, but the algebra is complicated.
As described in Sect. 4.1, numerical investigation shows that, for a given κ, the (α,λ) plane may be divided into four regions, each corresponding to a different regime of behaviour. Regime I is where ℜ(σ)>0 occurs for some k only when σ is real. Regimes IIa and IIb are where there are complex σ with ℜ(σ)>0 and with non-zero imaginary parts, and Regime III is where there is no instability, i.e. ℜ(σ)<0 for all k. The distinction between IIa and IIb is that in IIb, σ corresponding to the fastest growth over all k has a non-zero imaginary part. (An example is given in Fig. 5.) For f=0 only Regimes I and III are present, and Regime I, the region of instability, simply corresponds to Eq. (B3).
The boundary between Regime I and Regime IIa, indicated by the dashed curves in Fig. 5, is predicted by the condition that Eq. (B1) has a double root at σ=0. This is because close to the λ=0 axis a real root crosses the real axis. As λ increases, a bifurcation point emerges on the solution branch, joining to that root, and then moves towards the real axis. A transition between the crossing of the real axis by a real root and the crossing by a complex conjugate pair occurs when the bifurcation point reaches the real axis, i.e. at the point where there is a double root at σ=0. For all f such that Regimes IIa and IIb exist, the boundary between Regime IIa and Regime I asymptotes to α≃λ as λ→0. The condition for a double root at σ=0 must be compatible with the condition given by Eq. (B4). Given that the gradient of the dashed curves is observed to reduce as λ increases, there will not be compatibility if the gradient of the curve defined by equality in Eq. (B4) is larger than 1 as λ→0, i.e. if
When this inequality is satisfied there are no Regimes IIa and IIb, and the transition between stability and instability is fully described by Eq. (B4) or, equivalently, Eq. (21). The double root at σ=0 occurs when both the constant term and the coefficient of the linear term in Eq. (B1) vanish simultaneously. Eliminating k^{2} gives a constraint on λ and α (for given values of other parameters such as f and κ), which defines the boundary between regions I and IIa in Fig. 5. For illustrative purposes, the resulting expression when κ=0 is
In the small-α limit it corresponds to λ≃α. Since this expression for λ_{*}(α) ignores diffusivity, it will provide an accurate description of the boundary between regions I and IIa only when it is approximately independent of diffusivity, which, as seen in Fig. 5, appears to be the case when λ and α are small.
The boundary between regions IIb and III corresponds to a double root of Eq. (B1) of the form σ=is, s≠0. Some information about the small-α part of this boundary can be obtained by considering the special case α=0. There is then a root that is 𝒪(k^{2}) for small k, with a real part that has the same sign as $\mathit{\lambda}{\mathit{\mu}}_{\mathrm{1}}\left({\mathit{\mu}}_{\mathrm{2}}\right(Q/H)-{\mathit{\mu}}_{\mathrm{1}})-{f}^{\mathrm{2}}(\mathit{\lambda}+{\mathit{\mu}}_{\mathrm{2}}(Q/H\left)\right)$. It follows that for ${f}^{\mathrm{2}}<\left({\mathit{\mu}}_{\mathrm{2}}\right(Q/H)-{\mathit{\mu}}_{\mathrm{1}}){\mathit{\mu}}_{\mathrm{1}}$ there is instability at small k for large λ; therefore stabilisation requires non-zero α. On the other hand, for ${f}^{\mathrm{2}}>\left({\mathit{\mu}}_{\mathrm{2}}\right(Q/H)-{\mathit{\mu}}_{\mathrm{1}}){\mathit{\mu}}_{\mathrm{1}}$ there is stability at large λ, independent of the value of α. This explains the different forms of the boundary between regions IIb and III seen in Fig. 5 for f=0 and $f={\mathrm{10}}^{-\mathrm{5}}$s^{−1} on the one hand and $f=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ s^{−1} and $f=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ s^{−1} on the other.
The condition $f={f}_{\mathrm{stab}}(\mathit{\lambda},\mathit{\alpha})$ defines a function on the (λ,α) plane as the maximum f for a given α and λ (with dependence on other parameters such as κ suppressed) for which there are unstable modes. Contours of this function are shown in Fig. 5 as the bounding curves for Regime III.
It can be seen from Fig. 5 that the variation in $f={f}_{\mathrm{stab}}(\mathit{\lambda},\mathit{\alpha})$ is non-smooth and that this continues into the region of the (λ,α) plane where α and λ are both small, implying that the limiting value of f_{stab}(α,λ) as α→0 and λ→0 cannot be straightforwardly defined. To analyse this regime further it is useful to set $\mathit{\alpha}=\mathit{\lambda}=\mathrm{0}$ in Eq. (B1), leading to
where the moist gravity wave speed is ${c}_{m}^{\mathrm{2}}={c}^{\mathrm{2}}M<\mathrm{0}$ for the potentially unstable case. The roots of this polynomial are denoted as σ_{1}, σ_{2} and σ_{3}.
First, assume that there is a purely imaginary root σ=is, where s is real. Substituting into Eq. (B7) and separating the real and imaginary parts give
Eliminating s gives
For k>0, this is true only in the dry case ${c}_{m}^{\mathrm{2}}={c}^{\mathrm{2}}$. Since we are interested in ${c}_{m}^{\mathrm{2}}<\mathrm{0}<{c}^{\mathrm{2}}$, purely imaginary roots are not possible. Therefore any change in stability happens at σ=0.
For small k, Eq. (B7) has roots $\mathit{\sigma}=-{\mathit{\mu}}_{\mathrm{1}}$ and $\mathit{\sigma}=\pm if-{k}^{\mathrm{2}}(\mathrm{2}g{\mathit{\mu}}_{\mathrm{2}}{\mathit{\mu}}_{\mathrm{1}}{f}^{\mathrm{2}}Q\pm \mathrm{2}if({f}^{\mathrm{2}}{c}^{\mathrm{2}}+{c}_{m}^{\mathrm{2}}{\mathit{\mu}}_{\mathrm{1}}^{\mathrm{2}}\left)\right)/\mathrm{4}({f}^{\mathrm{4}}+{f}^{\mathrm{2}}{\mathit{\mu}}_{\mathrm{1}}^{\mathrm{2}})$, both of which have negative real parts. For instability, a change in stability must occur at some k>0. This requires the constant coefficient of Eq. (B7) to change sign; i.e.
If the system is to be stable, this must have no real roots k^{2}>0. We non-dimensionalise with $\widehat{f}=f\sqrt{\mathit{\kappa}/{\mathit{\mu}}_{\mathrm{1}}}/c$ and $\widehat{k}=k/\sqrt{\mathit{\kappa}/{\mathit{\mu}}_{\mathrm{1}}}$, for
The transition from two roots to none happens when the discriminant of this expression (as a quadratic in k) vanishes, or
This has four roots,
corresponding to the change in global stability, and
which is a spurious solution arising from roots with k^{2}<0. (When ${\widehat{f}}^{\mathrm{2}}>\mathrm{1}-M$, the left-hand side of Eq. (B12) is strictly positive for k^{2}>0.)
With zero damping, and assuming positive f, we therefore expect instability if and only if
where the value of f_{stab}(0,0) is defined by this equality. Since only one mode may be unstable in this limiting case, the unstable mode must have a purely real growth rate. When the function f_{stab}(α,λ) is defined by Eq. (B4), the limit of this expression as α→0 and λ→0 is not well defined, since it depends on the value of $\mathit{\alpha}/\mathit{\lambda}$.
This inconsistency may be resolved as follows. There is a value f_{I} such that the region of the (λ,α) space defined by ${f}_{\mathrm{stab}}(\mathit{\lambda},\mathit{\alpha})>{f}_{I}$ contains no parameters corresponding to Regime II; i.e. the unstable roots of Eq. (B1) are real for all k. Within this region, the stability boundary is given by Eq. (B4), rewritten here as a curve in the (λ,α) plane,
Note that as f is made arbitrarily large, this curve moves closer to the α axis, but does not meet it, so there is always some range of α such that there is instability for λ=0. Within this region, the transition with varying k of a single real mode from unstable to stable corresponds to the change in sign of the constant term of Eq. (B1). When $\mathit{\alpha}=\mathit{\lambda}=\mathrm{0}$, this constant term is identically zero. We therefore have a degenerate case where the previously unstable mode is zero for all k. The limit as α→0 and λ→0 of f_{stab} in this region is therefore not expected to be well defined.
In the other section of the (λ,α) plane, where f_{stab}<f_{I}, instability both near the origin and near f=f_{stab} relies on a combination of two complex roots of Eq. (B1). If α→0 and λ→0, only one of these roots can become identically zero, so we expect the other to remain the unstable mode, and hence the limit of f_{stab} should be well defined.
Since f_{stab} is continuous (apart from at the origin), we also expect the limit along the boundary between each of the cases discussed above, f_{stab}=f_{I}, to give the correct value of f_{stab}(0,0). This is a curve of constant f_{stab}, and therefore
The curve ${f}_{\mathrm{stab}}(\mathit{\lambda},\mathit{\alpha})={f}_{\mathrm{stab}}(\mathrm{0},\mathrm{0})$ has the expected form, α=λ, when α≪f.
The full derivation of linear equatorial wave response to a heating with general thermal damping λ and friction α is derived by Wu et al. (2001). The setup and relevant results are repeated here for reference.
We look for equatorially trapped steady-state solutions to the dynamical Eqs. (1) and (3), under the long-wave approximation. These satisfy
where u, v and h→0 are $y\to \pm \mathrm{\infty}$. For simplicity, we take a domain of effectively infinite zonal length, i.e. with a length much greater than the extent of the equatorial wave response. We can also assume that F(x,y) has a mean of zero, as a non-zero spatial mean part will be balanced by a change in the mean h and not contribute to the divergence or equatorial wave response.
We must first decompose the heating in terms of parabolic cylinder functions ${\stackrel{\mathrm{\u0303}}{D}}_{n}$ as
with ${L}_{\mathrm{eq}}=\sqrt{c/\mathit{\beta}}\sqrt[\mathrm{4}]{\mathit{\alpha}/\mathit{\lambda}}$.
The relevant parts of the solutions are then the equatorial Kelvin wave, which has
and the equatorial Rossby waves, for $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots}$,
The divergence is then calculated using Eq. (C3). This gives
For our purposes it is worth noting that for a heating, F<0, the equatorial wave response has h<0, and hence the associated divergence is positive. This is the opposite sign to the WTG divergence associated with a heating.
For diagnostics of numerical simulations, we want to calculate the expressions in Appendix C numerically, including separating the divergence into weak temperature gradient and equatorial Kelvin and Rossby components. The expression for the divergence, Eq. (C7), can simply be split into parts as follows:
where
and
The first two of these are simple to calculate numerically, using the expressions in Appendix C with $F={F}_{h}-{\stackrel{\mathrm{\u203e}}{F}}_{h}$. The Rossby wave component is then calculated as the residual
This will introduce some error as the numerical fields are not exactly quasi-steady; however it also sidesteps a few issues with calculating the Rossby wave response. The first of these is that simple theoretical expressions for the Rossby wave response to a given heating assume the long-wave approximation, whereas this assumption is not needed for the Kelvin wave and may not be valid for certain simulations. The second is that the Rossby wave response requires forcing due to all terms in Eq. (C4). For a finite discretised domain, these will be both poorly resolved and cut off by boundaries. This leads to errors in the projection, since the modes are no longer orthogonal.
The impact of including a temperature dependence in the precipitation, as done by Sugiyama (2009a, b), is briefly discussed in this section. The coupling terms between the moisture and thickness equations in these papers are assumed to be functions of the precipitation P, which is assumed to depend only on the moist static energy, c_{p}T+L_{v}q in standard units, and the nonlinear combination of P with a surface flux term. The terms thus take the form ${F}_{h}(q-h/\mathit{\gamma},h)$ and ${F}_{q}(q-h/\mathit{\gamma},q)$, where γ is proportional to ${L}_{v}/{c}_{p}$. This form is motivated in Sugiyama (2009a, b) as a Betts–Miller-type representation, but the alternative interpretation as a representation of the dependence of precipitation on free-tropospheric humidity, which emerged at the time of these papers, is noted, and this interpretation would not require the h dependence.
The forcing in this paper has been defined differently: we have taken a simple linear form for all forcing terms except for a nonlinear F_{h}(q), which provides an effective limit on the magnitude of the moisture variable. Two initial questions which arise from the temperature dependence of F_{h} and F_{q} are as follows: how does the behaviour change if the nonlinear term limits $q-h/\mathit{\gamma}$ rather than q, and what is the effect of a linear h term in the moisture equation?
In the first case we can gain some insight by redefining the moisture variable. In this case the moisture terms become ${F}_{h}(q-h/\mathit{\gamma})$ and ${F}_{q}(q-h/\mathit{\gamma})$. Subtracting $\mathrm{1}/\mathit{\gamma}$ times Eq. (2) from Eq. (3) gives
Defining a new variable, ${q}^{\prime}=q-h/\mathit{\gamma}$, along with corresponding ${Q}^{\prime}=Q-H/\mathit{\gamma}$ and ${F}_{{q}^{\prime}}={F}_{q}-{F}_{h}/\mathit{\gamma}$, we can write this as
We have returned to the original form of the equations to be studied with a different interpretation of the moisture variable and some terms due to h on the right-hand side of the moisture equation. Physically, it is unclear whether the diffusive forcing should act on moisture q or on moist static energy q^{′}, so the ∇^{2}h term may be unnecessary. If this term is neglected, we end up in the situation described by the second question above. The distinction between the two setups can therefore also be interpreted as whether the diffusivity is believed to act upon the moisture or moist static energy.
The situation described by Sugiyama (2009b) is one step further, with ${F}_{h}={F}_{h}(q,h)$ and ${F}_{q}={F}_{q}(q,h)$. We now investigate this case, but for simplicity we continue to take ϵ=0. In the strict WTG limit, h=h(t), using Eq. (2) to eliminate the divergence from Eq. (3) gives
This is still a reaction–diffusion equation for q; however now the stable states depend on h, ${q}_{\pm}={q}_{\pm}\left(h\right)$. These vary slowly as h varies slowly with time.
When L_{dyn} is larger than the domain size, the arguments in Sect. 2 hold, and aggregation proceeds as expected. Otherwise, WTG cannot be valid across the entire domain. Moist regions, with their associated heating, will have decreased h, and the opposite will apply to dry regions. On a β plane, the region of reduced h associated with a heating at the Equator will extend zonally with the equatorial Rossby and Kelvin wave responses. It is therefore expected that there will be local variation in the stable values of q. Of particular interest will be whether the adjustment to the fixed points due to the equatorial wave response to heating will facilitate moistening to the east or west of an equatorial moist region.
We know that the fixed points will depend on h, but what do we expect this variation to look like? We discuss two informative special cases. For the first, we assume that the moisture and thickness coupling terms depend only on the moist static energy, ${F}_{q}(q-h/\mathit{\gamma})$ and ${F}_{h}(q-h/\mathit{\gamma})$. Note that since the spatial mean of F_{h} must be subtracted when applying the WTG approximation, an additional linear relaxation term in the thickness equation will not affect the WTG moisture equation. In this case the stable fixed points are related to the stable fixed points q_{±0} when h=0 by
Thus, the stable values of the moisture increase with h.
For the second case we assume that the effect of h in the moisture equation is weak and therefore may be represented by a linear term. Taking inspiration from the previous case but Taylor-expanding the moisture about RCE, the moisture equation, Eq. (3), is
Note that ${F}_{q}^{\prime}<\mathrm{0}$.
Now, rather than the condition that at fixed points the reaction function G_{hq}=0, we require that ${G}_{hq}={F}_{q}^{\prime}\left(\mathrm{0}\right)h/\mathit{\gamma}$. This corresponds to a translation of the curve shown in Fig. 1 upwards with increasing h. Therefore as h increases the stable moisture values will correspondingly increase.
These two situations agree on the expected effect of h-dependent forcing in the height equation on the moisture distribution; however the impact on the propagation speed of aggregated regions is unclear. Since the deviation from isotropy in the β-plane case is due to the equatorial wave response to the heating, we expect significant cancellation between the equatorial Rossby and Kelvin wave components, as seen in the divergence response discussed in Sect. 5. We investigate this numerically.
Examples of numerical simulations from both cases are shown in Fig. E1. The first case, with both coupling terms being functions of only moist static energy, leads to only very small changes from the corresponding case in Fig. 11 – the spatial distribution and wavenumber are similar, and the speed has increased by 0.1 m s^{−1}. The specific heat parameter γ=22.2 has been chosen to match the combination of geopotential and temperature in Sugiyama (2009b). The fact that the effect is small is consistent with the large value of γ. The height dependence in the second case has had a larger effect, with the wavenumber reduced from 4 to 3 and the speed reduced by a factor of 3.
The code used for numerical simulations is available upon request to the corresponding author.
No data sets were used in this article.
Videos of the time evolution of the moisture and thickness fields for a selection of figures shown in this paper are available in Davison (2024) (https://doi.org/10.17605/OSF.IO/TUEMC).
MD and PH both contributed to the conceptualisation, investigation, analysis, writing and editing.
The contact author has declared that neither of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Both authors thank the reviewers for their valuable comments, which helped improve the paper during revision, as well as George Craig, Spencer Ressel and Tomoro Yanase for their helpful discussions.
This research has been supported by the UK EPSRC via a PhD studentship to Matthew Davison (grant no. EP/T517847/1).
This paper was edited by Peter Knippertz and reviewed by three anonymous referees.
Adames, Á. F. and Kim, D.: The MJO as a dispersive, convectively coupled moisture wave: Theory and observations, J. Atmos. Sci., 73, 913–941, https://doi.org/10.1175/JAS-D-15-0170.1, 2016. a, b, c, d, e, f, g
Adames, Á. F. and Wallace, J. M.: Three-dimensional structure and evolution of the vertical velocity and divergence fields in the MJO, J. Atmos. Sci., 71, 4661–4681, https://doi.org/10.1175/JAS-D-14-0091.1, 2014. a, b, c
Adames, A. F., Kim, D., Clark, S. K., Ming, Y., and Inoue, K.: Scale analysis of moist thermodynamics in a simple model and the relationship between moisture modes and gravity waves, J. Atmos. Sci., 76, 3863–3881, https://doi.org/10.1175/JAS-D-19-0121.1, 2019. a, b
Arnold, N. P. and Randall, D. A.: Global-scale convective aggregation: Implications for the Madden-Julian Oscillation, J. Adv. Model. Earth Sy., 7, 1499–1518, https://doi.org/10.1002/2015MS000498, 2015. a, b
Beucler, T., Cronin, T., and Emanuel, K.: A Linear Response Framework for Radiative-Convective Instability, J. Adv. Model. Earth Sy., 10, 1924–1951, https://doi.org/10.1029/2018MS001280, 2018. a
Biagioli, G. and Tompkins, A. M.: A Dimensionless Parameter for Predicting Convective Self-Aggregation Onset in a Stochastic Reaction-Diffusion Model of Tropical Radiative-Convective Equilibrium, J. Adv. Model. Earth Sy., 15, e2022MS003231, https://doi.org/10.1029/2022MS003231, 2023. a
Bray, A. J., Sluckin, T. J., McLeish, T. C., Blumenfeld, R., Hinch, E. J., Magerle, R., and Ball, R. C.: Coarsening dynamics of phase-separating systems, Philos. T. R. Soc. A, 361, 781–792, https://doi.org/10.1098/rsta.2002.1164, 2003. a
Bretherton, C. S., Blossey, P. N., and Khairoutdinov, M. F.: An energy-balance analysis of deep convective self-aggregation above uniform SST, J. Atmos. Sci., 62, 4273–4292, https://doi.org/10.1175/JAS3614.1, 2005. a
Carstens, J. D. and Wing, A. A.: A Spectrum of Convective Self-Aggregation Based on Background Rotation, J. Adv. Model. Earth Sy., 14, 1–14, https://doi.org/10.1029/2021MS002860, 2022. a, b, c, d
Carstens, J. D. and Wing, A. A.: Regimes of Convective Self-Aggregation in Convection-Permitting Beta-Plane Simulations, J. Atmos. Sci., 80, 2187–2205, https://doi.org/10.1175/JAS-D-22-0222.1, 2023. a, b, c
Craig, G. C. and Mack, J. M.: A coarsening model for self-organization of tropical convection, J. Geophys. Res.-Atmos., 118, 8761–8769, https://doi.org/10.1002/jgrd.50674, 2013. a, b, c, d, e, f
Davison, M.: Video supplements for “A simple model linking convective aggregation and large-scale dynamics”, OSF Home [video], https://doi.org/10.17605/OSF.IO/TUEMC, 2024. a
Emanuel, K., Wing, A. A., and Vincent, E. M.: Radiative-convective instability, J. Adv. Model. Earth Sy., 6, 75–90, https://doi.org/10.1002/2013MS000270, 2014. a, b, c, d
Hayashi, M. and Itoh, H.: A New Mechanism of the Slow Eastward Propagation of Unstable Disturbances with Convection in the Tropics: Implications for the MJO, J. Atmos. Sci., 74, 3749–3769, https://doi.org/10.1175/JAS-D-16-0300.1, 2017. a
Hirt, M., Craig, G. C., Schäfer, S. A., Savre, J., and Heinze, R.: Cold-pool-driven convective initiation: using causal graph analysis to determine what convection-permitting models are missing, Q. J. Roy. Meteor. Soc., 146, 2205–2227, https://doi.org/10.1002/qj.3788, 2020. a
Holloway, C. E. and Neelin, J. D.: Moisture Vertical Structure, Column Water Vapor, and Tropical Deep Convection, J. Atmos. Sci., 66, 1665–1683, https://doi.org/10.1175/2008JAS2806.1, 2009. a
Jiang, X., Adames, Á. F., Kim, D., Maloney, E. D., Lin, H., Kim, H., Zhang, C., DeMott, C. A., and Klingaman, N. P.: Fifty Years of Research on the Madden-Julian Oscillation: Recent Progress, Challenges, and Perspectives, J. Geophys. Res.-Atmos., 125, 1–64, https://doi.org/10.1029/2019JD030911, 2020. a
Khairoutdinov, M. F. and Emanuel, K.: Intraseasonal variability in a cloud-permitting near-global equatorial aquaplanet model, J. Atmos. Sci., 75, 4337–4355, https://doi.org/10.1175/JAS-D-18-0152.1, 2018. a
Muller, C., Yang, D., Craig, G., Cronin, T., Fildier, B., Haerter, J. O., Hohenegger, C., Mapes, B., Randall, D., Shamekh, S., and Sherwood, S. C.: Spontaneous Aggregation of Convective Storms, Annu. Rev. Fluid Mech., 54, 133–157, https://doi.org/10.1146/annurev-fluid-022421-011319, 2022. a, b, c, d
Muller, C. J. and Bony, S.: What favors convective aggregation and why?, Geophys. Res. Lett., 42, 5626–5634, https://doi.org/10.1002/2015GL064260, 2015. a
Neelin, J. D. and Zeng, N.: A Quasi-Equilibrium Tropical Circulation Model–Formulation, J. Atmos. Sci., 57, 1741–1766, https://doi.org/10.1175/1520-0469(2000)057<1741:AQETCM>2.0.CO;2, 2000. a, b, c, d, e
Patrizio, C. R. and Randall, D. A.: Sensitivity of Convective Self-Aggregation to Domain Size, J. Adv. Model. Earth Sy., 11, 1995–2019, https://doi.org/10.1029/2019MS001672, 2019. a
Raymond, D. J.: The Hadley circulation as a radiative-convective instability, J. Atmos. Sci., 57, 1286–1297, https://doi.org/10.1175/1520-0469(2000)057<1286:THCAAR>2.0.CO;2, 2000. a
Raymond, D. J., Sessions, S., Sobel, A. H., and Fuchs, Ž.: The Mechanics of Gross Moist Stability, J. Adv. Model. Earth Sy., 1, 9, https://doi.org/10.3894/james.2009.1.9, 2009. a
Romps, D. M.: Rayleigh Damping in the Free Troposphere, J. Atmos. Sci., 71, 553–565, https://doi.org/10.1175/JAS-D-13-062.1, 2014. a
Rostami, M. and Zeitlin, V.: Eastward-moving convection-enhanced modons in shallow water in the equatorial tangent plane, Phys. Fluids, 31, 021701, https://doi.org/10.1063/1.5080415, 2019. a
Rubinstein, J., Sternberg, P., and Keller, J.: Fast Reaction, Slow Diffusion, and Curve Shortening, SIAM J. Appl. Math., 49, 116–133, https://doi.org/10.1137/0149007, 1989. a, b, c
Sobel, A. and Maloney, E.: Moisture modes and the eastward propagation of the MJO, J. Atmos. Sci., 70, 187–192, https://doi.org/10.1175/JAS-D-12-0189.1, 2013. a, b, c, d
Sobel, A. H. and Maloney, E. D.: An idealized semi-empirical framework for modeling the Madden-Julian oscillation, J. Atmos. Sci., 69, 1691–1705, https://doi.org/10.1175/JAS-D-11-0118.1, 2012. a, b, c
Sobel, A. H., Nilsson, J., and Polvani, L. M.: The weak temperature gradient approximation and balanced tropical moisture waves, J. Atmos. Sci., 58, 3650–3665, https://doi.org/10.1175/1520-0469(2001)058<3650:TWTGAA>2.0.CO;2, 2001. a
Sobel, A. H., Bellon, G., and Bacmeister, J.: Multiple equilibria in a single-column model of the tropical atmosphere, Geophys. Res. Lett., 34, L22804, https://doi.org/10.1029/2007GL031320, 2007. a, b, c
Sugiyama, M.: The moisture mode in the quasi-equilibrium tropical circulation model. Part I: Analysis based on the weak temperature gradient approximation, J. Atmos. Sci., 66, 1507–1523, https://doi.org/10.1175/2008JAS2690.1, 2009a. a, b, c, d, e, f
Sugiyama, M.: The moisture mode in the quasi-equilibrium tropical circulation model. Part II: Nonlinear behavior on an equatorial β plane, J. Atmos. Sci., 66, 1525–1542, https://doi.org/10.1175/2008JAS2691.1, 2009b. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x
Wang, S. and Sobel, A. H.: A Unified Moisture Mode Theory for the Madden–Julian Oscillation and the Boreal Summer Intraseasonal Oscillation, J. Climate, 35, 1267–1291, https://doi.org/10.1175/JCLI-D-21-0361.1, 2022. a
Windmiller, J. M. and Craig, G. C.: Universality in the spatial evolution of self-aggregation of tropical convection, J. Atmos. Sci., 76, 1677–1696, https://doi.org/10.1175/JAS-D-18-0129.1, 2019. a, b, c, d
Wing, A. A., Emanuel, K., Holloway, C. E., and Muller, C.: Convective Self-Aggregation in Numerical Simulations: A Review, Surv. Geophys., 38, 1173–1197, https://doi.org/10.1007/s10712-017-9408-4, 2017. a, b, c
Wu, Z., Sarachik, E. S., and Battisti, D. S.: Thermally driven tropical circulations under Rayleigh friction and Newtonian cooling: Analytic solutions, J. Atmos. Sci., 58, 724–741, https://doi.org/10.1175/1520-0469(2001)058<0724:TDTCUR>2.0.CO;2, 2001. a, b, c, d
Yanase, T., Nishizawa, S., Miura, H., and Tomita, H.: Characteristic Form and Distance in High-Level Hierarchical Structure of Self-Aggregated Clouds in Radiative-Convective Equilibrium, Geophys. Res. Lett., 49, e2022GL100000, https://doi.org/10.1029/2022GL100000, 2022. a
Yang, D.: Boundary Layer Diabatic Processes, the Virtual Effect, and Convective Self-Aggregation, J. Adv. Model. Earth Sy., 10, 2163–2176, https://doi.org/10.1029/2017MS001261, 2018. a
Yang, D.: A shallow-water model for convective self-aggregation, J. Atmos. Sci., 78, 571–582, https://doi.org/10.1175/JAS-D-20-0031.1, 2021. a
Yano, J.-I. and Tribbia, J. J.: Tropical Atmospheric Madden–Julian Oscillation: A Strongly Nonlinear Free Solitary Rossby Wave?, J. Atmos. Sci., 74, 3473–3489, https://doi.org/10.1175/JAS-D-16-0319.1, 2017. a
Zeng, N., Neelin, J. D., and Chou, C.: A Quasi-Equilibrium Tropical Circulation Model–Implementation and Simulation, J. Atmos. Sci., 57, 1767–1796, https://doi.org/10.1175/1520-0469(2000)057<1767:AQETCM>2.0.CO;2, 2000. a, b
Zhang, C., Adames, Á. F., Khouider, B., Wang, B., and Yang, D.: Four Theories of the Madden-Julian Oscillation, Rev. Geophys., 58, e2019RG000685, https://doi.org/10.1029/2019RG000685, 2020. a
- Abstract
- Introduction
- Model system to be studied
- Numerical simulations in the WTG regime
- Breakdown of WTG and implications for aggregation
- Equatorial β plane
- Conclusions
- Appendix A: Numerical details
- Appendix B: Details of the f-plane linear stability problem
- Appendix C: The equatorial wave response to a moist heating
- Appendix D: Decomposition of the divergence into equatorial wave components
- Appendix E: The dependence of precipitation parametrisation on temperature
- Code availability
- Data availability
- Video supplement
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Abstract
- Introduction
- Model system to be studied
- Numerical simulations in the WTG regime
- Breakdown of WTG and implications for aggregation
- Equatorial β plane
- Conclusions
- Appendix A: Numerical details
- Appendix B: Details of the f-plane linear stability problem
- Appendix C: The equatorial wave response to a moist heating
- Appendix D: Decomposition of the divergence into equatorial wave components
- Appendix E: The dependence of precipitation parametrisation on temperature
- Code availability
- Data availability
- Video supplement
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References