A model of coupled oscillators applied to the aerosol – cloud – precipitation system

Abstract. We simulate the aerosol–cloud–precipitation system as a collection of cloud elements, each coupled through physically based interactions with adjacent clouds. The equations describing the individual clouds follow from the predator–prey model of Koren and Feingold (2011) with the addition of coupling terms that derive from the flow of air between the components resulting from surface divergence or convergence of flows associated with the life cycle of an individual cell. It is shown that some degree of coupling might stabilize clouds that would ordinarily become unstable. Varying the degree of coupling strength has significant influence on the system. For weak coupling, the clouds behave as independent oscillators with little influence on one another. As the local coupling strength increases, a point is reached at which the system becomes highly synchronized, similar to the Sakaguchi et al. (1987) model. Individual cloud oscillators in close proximity to one another can be both in-phase or out-of-phase depending on the choice of the time constant for the delay in communication between components. For the case considered, further increases in coupling strength result in reduced order and eventually unstable growth. Finally it is demonstrated that the set of coupled oscillators mimics qualitatively the spatial structure and synchronized behaviour of both closed and open-cellular cloud fields observed in satellite imagery, and produced by numerically intensive large eddy simulation.


Introduction
The study of weather and clouds can most certainly be traced back to prehistoric times; there is ample evidence that our ancestors were keen observers of clouds, rain, and other at-mospheric phenomena of immediate concern for survival.An early, written example of the human penchant for identifying patterns in the skies can be seen in Aristotle's Meteorologica, 340 BCE.The application of rigorous scientific method to weather forecast and clouds began in earnest in the late 19th century and early 20th century.Larger-scale cloud systems and associated cloud patterns often manifest at scales much larger than the field of view of a surface observer and with a few notable exceptions (e.g.Graham, 1934) it has only been in the past 70 yr that intensive study of cloud and weather systems has been made possible through computers, weather radars and, since the 1960s, satellite-borne cameras and remote sensors (Agee, 1984).
Rather than address weather patterns, this paper concerns itself with commonly observed patterns in cloud systems on scales of O(10-1000 km).The motivation of this study is to show that cloud systems often organize in ways that are common to other natural systems.Specifically, it addresses the phenomenon of synchronization within cloud systems.As is the case in other animate or inanimate systems, interactions between components of the system result in a system-wide emergence of order.When the components of the system are sufficiently well coupled, synchronization can result in particularly fascinating cloud patterns that constantly regenerate themselves (Feingold et al., 2010;Koren andFeingold, 2011, 2013) and are quite resilient to perturbation.
In the work presented below, we will apply a simple set of non-linear equations originally used to describe a cloud system (Koren and Feingold, 2011) to each of the individual clouds or cloudy elements comprising the system.The individual cloud elements grow, decay as they lose water via rain, and regenerate, much as they might in nature; each cloud element behaves as a non-harmonic oscillator.In addition, G. Feingold and I. Koren: Coupled oscillators in the cloud system the individual clouds are coupled to one another through the surface convergence (or divergence) flows associated with growing (or decaying) clouds.The system can therefore be described as a set of locally coupled oscillators.Earlier studies of such systems include Sakaguchi et al. (1987) and Blasius et al. (1999).
In Koren and Feingold (2011), parameters of the equation set were varied over a large range to identify different modes of behaviour (i) stable, non-precipitating; (ii) stable, weakly precipitating; (iii) strongly precipitating clouds with stable oscillations; and (iv) explosive growth.Because it is precipitation that drives downdrafts and surface divergence/convergence and therefore the potential mechanism for coupling, the focus here will be on the strongly precipitating state.We will demonstrate that coupling strength has significant influence on the system.Through a series of idealized one-dimensional tests we explore some of the parameters that influence system order.Results are extended to two dimensions and analogs to natural cloud systems are discussed.
2 Model description Koren and Feingold (2011), henceforth KF, presented an adaptation of the predator-prey equations, applied to a cloud system, consisting of three equations for cloud depth H , drop concentration N and rain rate R: and, where c 1 is a temperature-dependent constant, and c 2 and α are constants based on theory.H 0 is the meteorological "carrying-capacity", i.e. the H that would be reached within a few timescales (τ 1 ) in the absence of rain-related losses.Similarly, N 0 is the drop (or aerosol) concentration "carrying-capacity" that the system would reach in a few τ 2 in the absence of rain.A full description of the terms, and their values and units can be found in KF.N and R are assumed to be constant with height, which suffices for the current application (see KF).We note that R is diagnosed from the prognostic variables H and N. The delay function ensures that the current R is a function of the conditions at an earlier time; T represents the time required for cloud water to be converted to rainwater by collision and coalescence between drops (T ∼ 15 − 20 min).The approach is similar to other work that has applied simplified equation sets to help understand atmospheric systems (Wacker, 2006;Garrett, 2012), and is in the spirit of the pioneering work of Lorenz (1963).
In the current work, we apply Eqs. ( 1)-( 3) to individual cloud cells, each evolving through its individual life cycle, and connected through a coupling term (to be described below).We consider cloud cells on a grid (i, j ) so that H , N , and R in Eqs. ( 1)-(3) are replaced by H i,j , N i,j and R i,j .Thus we have and where η i,j is the coupling strength, Ḣi,j the time rate of change of H i,j for cloud oscillator i, j , and the time constant τ c,i,j is associated with the time it takes for air to move from an adjacent cell to cell i, j .System coupling is clearly the product of the η and Ḣ terms, i.e. changes in H are an essential part of the coupling.Manifestation of the interaction between the individual cloud cells might occur in nature in a number of ways: (i) a buoyant, growing cloud cell is associated with upward motion and therefore surface convergence of air within some radius of influence.Similarly a decaying cloud cell that is discharging rain is associated with downdrafts and surface divergence.Through arguments of mass continuity, convergence or divergence will influence neighbouring clouds; (ii) individual cloud cells progressing through their life cycles all interact with their environment by redistributing heat and moisture, and consuming instability (e.g.Nober and Graf, 2005).The environment thus mediates the growth of the individual cloud cells.We focus here on the first of these coupling mechanisms since Nober and Graf (2005) have to some extent addressed the first within the framework of a parameterization of atmospheric convection.
The sign of the coupling terms is such that if an adjacent cloud (e.g.H i−1,j or H i+1,j ) is decaying ( Ḣ < 0), it increases the rate of growth of H i,j , whereas if an adjacent cloud is growing ( Ḣ > 0), it suppresses the growth of H i,j ; thus, the results of practical interest are for η i,j < 0. On a spatial grid of unspecified dimension, τ c,i,j determines the length scale of the resulting features, provided one has some information on the velocity of the flow, which in cloud systems is O(10 m s −1 ) for horizontal flow and O(1 m s −1 ) for vertical flow.Such specificity seems premature and at this stage we do not delve into details of spatial scale of features.The situation is further complicated by the fact that τ c,i,j is expected to decrease as Ḣi,j increases but again we defer such exercises to a later stage.Finally, horizontal flow is a result of cloud couplings only; we assume no additional horizontal wind, or wind shear, which would break the pattern formation.
A simple one-dimensional schematic of the coupled system is shown in Fig. 1.Growing clouds are positively buoyant and associated with updrafts.Once clouds thicken sufficiently to generate rain, the evaporating rainfall evaporates, causing colder and denser, downward moving air.Downdrafts reaching the surface are forced to diverge, providing the opportunity for interaction between the cells as outflows converge at the surface.Note that some cloud systems are capped by a strong inversion, which although not impenetrable, acts much like the rigid surface boundary.This is implicitly accounted for in two ways: (i) a growing cloud, such as the one at centre top in Fig. 1, will tend to hasten the demise of adjacent cells, either by dragging surface air from those clouds, and/or by generating subsiding air from above.(We make no distinction between these mechanisms here.)(ii) The formulation to some extent includes this upper boundary via the H 0,i,j parameter: an approximately constant value implies that vertical development is capped, much as in stratocumulus clouds, unless contributions from adjacent cloud elements are particularly strong (e.g. for large |η|).For cases with highly variable H 0,i,j as in the case of a field of cumulus clouds, rising air motions can manifest as much deeper clouds, and there exists significant variability in H i,j .We defer discussion of cumulus cloud field analogs to later work.
In the subsequent calculations, it will be assumed, for simplicity that a cloud oscillator is only directly affected by its closest neighbours so that for a one-dimensional system, the summation term on the right hand side of Eq. (4) reduces to contributions from two neighbouring cells i − 1 and i + 1, and for a two-dimensional system, it includes contributions from eight adjacent points on a square grid or 6 adjacent points on an hexagonal grid.Coupling is weighted by the distance from cloud (i, j ); for a a square grid of unit length, adjacent points are weighted by unity and diagonal points are weighted by (1/ √ 2); on an hexagonal grid, all adjacent points have equal weighting.Calculations that include coupling with oscillators farther afield, and with 1/(distance) 2 weighting have also been performed, and while the qualitative nature of the results is often much the same, on occasion, this might cause a system to shift into a different mode.
For simplicity we assume τ c,i,j = τ c = constant and η i,j = η = constant for any given simulation, unless otherwise noted.Thus the temporal variation in the coupling is determined by the Ḣ terms.The discrete derivatives are calculated over t = 5 s, and updated every time step.Tests with different time intervals showed that results are not particularly sensitive, provided the interval does not span values approaching T or τ c .Lateral boundaries are assumed to be rigid, i.e. oscillators i (or j ) = 1 and i (or j )= I , only experience influ- ence from neighbouring oscillator i + 1 (or j + 1) and I − 1, respectively.
The system of Eqs. ( 4)-( 6) is represented by five primary parameters, H 0 , N 0 , τ 1 , τ 2 , T , and two coupling-related parameters η and τ c .Fortunately, many years of boundary layer cloud research provide guidance for these choices for given meteorological conditions.The examples to be presented below draw on that experience but are presented as phenomenological illustrations and for the most part do not attempt to mimic detailed simulation or observation.Instead the focus is on showing how a simple set of equations can represent previously documented emergent aspects of a coupled aerosolcloud-precipitation system.The problem is similar in principle to earlier theoretical treatments of locally coupled systems (e.g.Sakaguchi et al., 1987or Blasius et al., 1999) although the (multiple) delay times T and τ c add significant complexity and the potential for multiple synchronized solutions (e.g.Schuster and Wagner, 1989).

Results
Since parameter space for H 0 , N 0 , τ 1 , τ 2 , and T , was explored in KF, we focus here on the coupling parameters η and τ c .A base case simulation is chosen with the same parameters as those chosen in KF's Fig. 7 (H 0 = 670 m; N 0 = 515 cm −3 , τ 1 = 80 min, τ 2 = 84 min, and T = 21.5 min) and the influence of η and τ c is explored.Unless otherwise stated, these will be the parameters of choice in subsequent simulations.A constant H 0 implies constant frequency content for each of the i, j non-harmonic cloud oscillators in the absence of coupling; Gaussian distribution of H 0 is used in the two-dimensional examples as noted below.These parameters tend to create clouds with H ≈ 500 m and R ≈ 5 mm d −1 , similar to what one might find in precipitating stratocumulus clouds (e.g.Feingold et al., 2010).The weakly-or nonprecipitating regime tends to occur for smaller H 0 and/or larger N 0 and was explored fairly extensively in KF.Because precipitation drives strong changes in Ḣ , and therefore coupling strength, we prefer to focus on a more strongly precipitating case.However, as will be shown below, our systems typically start out with little or weak precipitation (based on initial conditions) and only later evolve to strongly precipitating solutions, thus affording us the opportunity to study both.
The duration of simulations is typically 1000 min (≈ 17 h), which is usually sufficient to explore robust solutions, and is also an atmospherically relevant time frame.On occasion these simulations have been extended for twice as long to persuade ourselves that the results are representative.
It is expeditious to demonstrate some basic ideas with a one-dimensional array before moving to the two-dimensional array.In the interests of brevity, results focus entirely on cloud depth H ; for a qualitative idea of R and N behaviour, we note that the R field lags that of H by ∼ T , and the N field is anticorrelated with R (see KF).

Illustrative examples
We start with results for 21 coupled oscillators all with constant H 0 distributed along a line for η = −0.51 and τ c = 5 min (∼ 1/4 × T ; Fig. 2).For atmospheric conditions this τ c is probably too small but is a useful start point.For clarity, only a sample of the 21 oscillators are included.Red lines are for i = 5 and i = 7 and black lines are for i = 10 and i = 12.
During the first part of the simulation all oscillators quickly enter a coupled state with H ≈ 500 m.It is of interest that these individual cloud elements have much smaller amplitude excursions than the KF simulation for the same conditions (their Fig. 7), but applied to the entire system (Eqs.1-3).Thus for this configuration, the coupling of the individual elements tends to provide a strong stabilizing influence on the system.At about 375 min, the odd (red line) oscillators begin to amplify their oscillation, due to the development of an instability.(This disturbance can be reproduced under multiple different conditions and is therefore regarded as meaningful.)This in turn perturbs the even (black line) oscillators.By about 440 min, the odd oscillators have much the same frequency and phase as the even oscillators, however the odd and even sets are locked in opposite phase, as one might expect from Fig. 1.Their amplitude excursion is large with clouds reaching close to 1500 m by 1000 min (i.e.significantly larger than experienced in KF's Fig. 7, where the maximum H was about 550 m).The extent to which H i can exceed H 0 is related to the magnitude of |η| (Eq.4).2, all clouds are initially synchronized with H 500 m.System undergoes a period of no phase or frequency coherence from about 400-600 min.Thereafter the system behaves much like Fig. 2 with clouds growing at the expense of their neighbours.Note the slower frequency of oscillation associated with the longer τ c .
We now increase τ c to 10.75 min (1/2 × T ) and decrease η to −0.48 (Fig. 3).The selected oscillators are 5; 7 (red) and 6; 8 (black).Figure 3 again shows almost constant H during the early part of the simulation, and as in Fig. 2, development of instability.This time, however, all four selected oscillators experience a period of time (∼ 400 to ∼ 600 min) when there is little coherence in frequency or phase.After ∼ 600 min, behaviour is similar to that in Fig. 2, albeit with a slightly lower frequency, due to the increase in τ c .While analyzing the model output we found that oscillators 10 and 12 never   and 3. Cloud depths H i,j are significantly higher than in Fig. 3 because of the larger |η|.For the case under discussion, H of up to ∼ 5 × H 0 can be sustained before the system becomes unstable.Note that the increase in |η| from −0.48 to −0.54 does not result in an acceleration of the entrainment of all the oscillators to a common frequency.First, this is because the added forcing terms associated with η have a strong influence on each individual H i , which then feeds back to N i and R i , and changes the inherent frequency of the oscillator.Second, the system has the added complexity of both the microphysical delay T and the coupling delay τ c , superimposed on the meteorological and drop concentration timescales (τ 1 and τ 2 ); this has the potential to produce a multitude of synchronized solutions (Schuster and Wagner, 1989).
Figure 4 intimates the possibility of a threshold of |η|, above which the system might exhibit stronger coherence.Figure 5 is the same as Fig. 4d but with |η| reduced to −0.53.Indeed one sees how this small change has significant influence on the individual components of the system; the even (black) oscillators (10; 12) are approximately in phase, whereas the odd (red) oscillators (5; 7) exhibit multiple modes deriving from the various terms in Eq. ( 4) with their different frequencies.The threshold behaviour manifested in the differences in results between Figs. 4 and 5 is reminiscent of the Sakaguchi et al. (1987) model, as will be discussed in Sect.4.1.
In the next simulations we increase the value of τ c to 43 min (2 × T ) and set η = −0.63.In addition to the further reduction in frequency compared to smaller τ c , the most striking feature of this simulation is that the odd (red) oscillators (5; 7) and even (black) oscillators (6; 8) all have similar phase and frequency (Fig. 6).Further examination indicates that this is generally true, regardless of the specific choice of oscillators.This response is distinctly different from the previous examples (at shorter τ c ) where adjacent cloud elements tended to be out of phase.We repeat the simulation but reduce |η| to −0.621 (Fig. 7).In this case the oscillators are incoherent from 550 to 1000 min (Fig. 7a).Extending the simulation out to 2000 min shows that there are periods of time when the adjacent cloud elements prefer to be out of  phase (Fig. 7b; 1400-1600 min), but that this is short-lived and incoherent oscillations return (Fig. 7c; 1800-2000 min).

Sensitivity to coupling strength
We now provide a broader analysis of the system for different levels of coupling strength η.For the base conditions used above (H 0 = 670 m; N 0 = 515 cm −3 , τ 1 = 80 min, τ 2 = 84 min, and T = 21.5 min), a series of many simulations covering 0.1 ≤ |η| ≤ 0.6 were performed at small increments in η to probe the system more deeply.The number of clouds was increased to 41, in recognition of the dependence of system-wide order on the number of oscillators (Sakaguchi et al., 1987).The size of the η increment was adjusted in the region of particularly interesting and/or transitional behaviour.Before delving into the broader analysis, an illustrative example upon which some of the subsequent analysis rests is shown for η = −0.53 (Fig. 8).As shown below this choice of η gives the maximum system coherence.The time evolu- tion of H i for i = 1, 5, and 21 is shown in Fig. 8a, followed by a pair of matrices of the maximum correlation between a pair of clouds i, i (Fig. 8b) and the lag time (in min) at which this maximum correlation occurs (Fig. 8c).The correlation matrix shows that large parts of the cloud field are highly correlated (r > 0.9) but with variability in the time lags (e.g.lower left quadrant of Fig. 8c).In regions of high r, the chequered pattern in the lag matrix is indicative of clouds that are in phase (high r and close to zero lag) alternating with clouds that are out of phase, as seen in earlier illustrative examples.In general, high correlation and smooth patterns in the lag matrix are indicative of coherence.Sometimes alternating "warm" and "cool" colours follow a snake-like pattern suggesting a deeper underlying structure.Parts of the domain, in the vicinity of i = 20, are characterized by relatively poor correlation (0.6 < r < 0.7), with no concomitant structure in the lag matrix.
A broader view of the system is now shown for a few select η (= −0.3, −0.475, −0.53, and −0.55062) that illustrate characteristic modes of behaviour (Fig. 9).Analysis is performed over the latter part of the simulation period (1188-2000 min).The frequency plots (first column) show the average of the periodicities of each of the 41 clouds based on the power spectra.The amplitude of these plots is another measure of the coherence of the system.They show a distinct shift from a dominant period of 80 min at η = −0.3, with a gradual shift to ∼ 60 min at the higher |η|.The strong peak in spectral density at η = −0.53reflects the maximum system coherence (see also Fig. 8).Further increases in |η| result in unstable, explosive growth of the clouds.(Note we were not able to achieve stable solutions for this set of conditions for |η| > 0.550627.)Column 2 of Fig. 9 presents another view of the cloud system in polar coordinates of correlation (radius of the circle) and phase (angle between 0 and 2π).The frequency and phase of cloud i = 1 was used as a reference oscillator, against which each of the 41 clouds were compared.Phase lags between 0 and 90 • and 0 and 270 • occur because neither the reference cloud oscillator, nor the others are fixed.With increasing |η| the system moves from one with relatively low correlation to a highly correlated and more coherent system at η = −0.53.As we approach the solutions with largest permissible |η|, correlation decreases significantly.At this point a few clouds experience very large amplitude oscillations, in contrast to the situation at lower |η| when the H i are more homogeneous (figures not shown).Column 3 of Fig. 9 presents a histogram of the correlation matrix (e.g. a summary of correlations such as those in Fig. 8a), providing further evidence of the initial migration to larger correlation as |η| increases, followed by significant decrease in correlation close to the unstable, strongly coupled state.The last column of Fig. 9 summarizes the correlation and lag matrices as two-dimensional probability distribution functions with correlation on the vertical axis and lag on the horizontal access.As the coupling strength increases there is a steady decrease in the time lags associated with maximum correlation, to some extent because of the shift to shorter periodicity (first column Fig. 9).The strong order in the system at η = −0.53 is shown in a compact fashion.

Oscillators in a two-dimensional array
Simulations are extended to two dimensions as specified by Eqs. ( 4)-( 6).Rather than explore concepts addressed in Sect.3.1, we turn our attention instead to pattern development in cloud systems experiencing significant precipitation, as studied by Feingold et al. (2010), Koren andFeingold (2011, 2013).Through a combination of large eddy simulation and satellite observation, those studies demonstrated that non-precipitating clouds achieve an approximately steady state in which losses of water to evaporation and precipitation are balanced by condensation of water.The clouds take on a closed cellular pattern that maintains a fixed structure The progression to larger |η| is initially associated with an increase in system coherence and correlation (first 3 rows) but at sufficiently large |η| (= −0.5506), there is a shift to an unstable system with poorer coherence and correlation.In this state, a few clouds have very large H i and the system cannot support further increases in |η| without exploding.
over long periods of time (∼ day).In contrast, in the presence of strong precipitation, the cloud system experiences repeated larger amplitude cycles of cloud water build-up followed by significant reduction in cloud water by rain, akin to a predator-prey cycle.The precipitating cloud pattern is an open-cell structure that constantly rearranges itself as positively buoyant cloudy regions produce rain that generates negative buoyancy.Here we attempt to reproduce these modes with the equations on a two-dimensional array.
We focus on larger values of τ c since these are likely more relevant to boundary layer cloud systems exhibiting mesoscale cellular convection.(Shorter, variable coupling times are more appropriate to cumulus cloud fields and will be explored in later work.) The first example uses the same base conditions used in Sect.3.1 (H 0 = 670 m; N 0 = 515 cm −3 , τ 1 = 80 min, τ 2 = 84 min, and T = 21.5 min) with τ c = 43 min and η = −0.20.To provide a little more realism to the simulations a small amount of spatial variability is added to H 0 (±100 m) and to η (±0.03), and maintained for the duration of the simulation.The variability in H 0 translates to a spread in the frequency distribution of the cloud oscillators, as in other theoretical studies.The relatively large number of individual os-cillators, even for a 21 × 21 grid, calls for alternate ways of displaying the model output.For the first, we use the "Hovmöller" diagram, commonly applied in meteorological analysis (Fig. 10), to display the of progression of H i,j (t).The figure shows the spatial dimension on the abscissa and time on the ordinate.To do so, some degree of averaging is performed over the j (y) dimension, in this case, over all 21 j points.Because there is no preferred symmetry in the current equation set, the main features of the field manifest when either the i or the j dimension is displayed on the abscissa, although they are not the same.
For the first 350 min of simulation, the cloud field is characterized by a spatially rigid structure.The field rises and falls in unison but the relative heights of H i,j are approximately constant with time, as borne out by the series of snapshots in Fig. 11.This reflects a frequency and phase coherence of all the oscillators, akin to the early stages of simulations in Figs. 2 and 3, or the duration of Fig. 6.After ∼ 400 min, the rigid structure is broken by the (delayed) effects of stronger rainfall; large values of H i,j are replaced by small values, both spatially and temporally, and the cloud field constantly rearranges itself with a period of ∼ 90 min.The migration in time of the peaks in H is reminiscent of the The display shows contours of H i,j as a function of i (x distance; abscissa) and time (ordinate).The linear colour scale is from 0 to 800 m (purple through red/brown) Model output has been averaged over all j -grid points.Prior to about 300 min the all components are synchronized and rise and sink in unison (see also Fig. 11).Thereafter, the spatial symmetry is broken by the onset of rain and a gridded structure merges depicting spatiotemporal oscillations in H i,j .
"travelling wave structure" discussed by Blasius et al. (1999) for coupled ecological systems.
The individual snapshots associated with the rearranging cloud field are examined more closely with a second twodimensional case that exhibits this more clearly.Conditions are the same as in Figs. 10 and 11 but |η| is increased to −0.22 ± −0.03.In addition, the domain size is increased to 41 × 41 points to allow more structure to develop.As in Figs. 10 and 11, this case starts out as a non-precipitating system with a temporally invariant structure.Strong rain develops from thick clouds and forces cellular cloud structures that rearrange with similar periodicity (∼ 2 × τ c ; Fig. 12).The cell-like structure and their rearrangement has been documented in large eddy simulation of cloud fields (Feingold et al., 2010) and in geostationary satellite imagery (Koren and Feingold, 2013), with similar periodicity (see Sect. 4).

Threshold behaviour
Threshold behaviour is prevalent in the system of equations under discussion -in both one-and two-dimensional cases.In the simulations shown here and in our experience with the many other simulations performed, small changes in coupling strength η result in a rapid transition from incoherent to synchronized oscillations.However, further increases in η lead to unconstrained growth in H i,j to unrealistically high values and instability of the system.The fact that the cloud oscillators are non-harmonic, and also include delays, makes for a rather complex system; individual oscillators have multiple intrinsic frequencies and their interaction can generate constructive or destructive coupling, with the result that the system can migrate from an ordered to a chaotic state depending on the specified parameters.In cloud systems, rain itself exhibits threshold behaviour (e.g.Kessler, 1969), which generates threshold behaviour in Ḣ .This likely influences the timing of coupling and therefore the manifestation of synchronization.The connection to the work of Sakaguchi et al. (1987), Blasius et al. (1999) and to a host of other more recent studies in a wide range of disciplines is quite clear.One study in the field of atmospheric physics of particular interest (Yair et al., 2009) showed how component clouds in a lightning storm tend to discharge in unison when the group coupling strength in the network is higher than a threshold value.In that study, the coupling occurred through a network-wide sharing of the change in the electric field strength E following a lightning strike at a given location in the system: an individual cloud element (i, j ) in the evolution of its life cycle of E, experienced a contribution E based on its distance to the cloud experiencing discharge.In their case coupling is global, as in Kuramoto (1975).In the case of electrified cloud systems, the delay time for communication of E is extremely short (seconds or fraction thereof), thus facilitating the communication across large distances (order 100 km).The analogy to our system is clear: individual clouds evolve through their life cycle of H and experience some positive or negative contribution to their growth based on inflows or outflows from adjacent clouds.The analogy is strengthened by the fact that in natural clouds, lightning is itself a function, inter alia, of the vertical development of the cloud and the formation of rain.Therefore, both systems have similar delay terms T (Eqs.4-6).However, there are distinct differences in the manner in which these systems are coupled.In the case of electric field build-up and discharge, communication is very rapid (τ c → 0) while for evolution of H , material flow, limited by the atmosphere's ability to move matter (air) from one location to another, provides the vehicle for communication, and τ c is substantial.

Cloud patterns: mesoscale cellular convection
Vast areas of the mid-latitude and subtropical oceans are covered by shallow boundary layer clouds.On the western coasts of continents where cold ocean currents prevail, the cloud system tends to exist in a closed-cellular state characterized by a high cloud fraction and high reflectance.The bright cellular cloud structures denote upward moving air and have dimensions of order 10 km.They are ringed by narrow, darker regions of downward moving air.These structures have been tracked by geostationary satellites, and after removal of the ambient wind, shown to be remarkably rigid in structure (Koren and Feingold, 2013).The sequence of cloud fields in Fig. 11 sampled prior to about 350 min (Fig. 10) also maintain a rigid structure and can be considered to be an analog to the closed-cellular state.
In stark contrast, over warmer ocean waters driven by surface heating, cloud systems tend to have an open-cellular, almost honeycomb-like structure with buoyant and reflective cloudy cells surrounding dark, negatively buoyant opencellular regions.Sometimes, open cells can be found over cold oceans that prefer the closed-cellular state.In these cases, rain has been shown to play a central role in creating and maintaining this structure.Our earlier work has identified the presence of rain as being responsible for a spatiotemporal oscillation in the cloud field: positively buoyant clouds thicken to the point that they can develop rain; rain falls and evaporates below cloud, generating negative buoyancy and downward airflow (downdrafts); at the surface, downdrafts are forced to diverge; here they encounter outflows from adjacent raining cells that result in surface convergence and the formation of positively buoyant updrafts that generate clouds, and so on.The conversion of cloud water (a proxy for H ) to rain (R) thus forces positively buoyant air to become negatively buoyant, resulting in a constant rearrangement in the spatial structure of the cloud field.Figure 10 shows such rearrangement after ∼ 400 min when rain has become more substantial.This is illustrated more clearly in Fig. 12 where a larger value of |η| creates better spatial coherence.Similar patterns can be seen in the output of cloud fields simulated by large eddy simulation (Feingold et al., 2010), or in geostationary satellite imagery (Koren and Feingold, 2013).

Summary
We have extended the predator-prey equations for the aerosol-cloud-precipitation system first presented by Koren and Feingold (2011) to individual subcomponents of the system and shown how solutions to the coupled system of equations act in ways that mimic some of the characteristics and spatial patterns developed by computationally intensive large eddy simulation and observed by satellite imagery.The individual cloud components can be considered "oscillators"; they grow in depth H based on the meteorological environment that spawns them, and decay as cloud water is converted to rain R. The atmospheric aerosol, represented here by drop concentration N acts to regulate the production of rain.Our focus here is on a system that develops significant amounts of precipitation because of the richness that it adds to the coupled system via the Ḣ terms in Eq. ( 4) (e.g.Fig. 1); The individual cloud elements (oscillators) are influenced by one another through the airflows associated with evaporating precipitation.The negatively buoyant air in a raining cloud moves toward the (rigid) lower surface where it is forced to diverge.The outflows of adjacent cloud elements that are in phase with one another will collide; this will generate upward moving air in which new clouds will form.Adjacent clouds that are out of phase will create a feedback that hastens the demise of the decaying cloud while simultaneously fuelling the growth of the growing cloud (e.g.Figs. 2 and 3).Exactly how these components interact depends on a number of different parameters associated with the model equations.The coupling strength η and delay parameters T and τ c in particular have the potential to generate complex solutions.The number of tunable parameters intimates the potential for significant further study.
For the precipitating case under consideration, we have shown that some degree of coupling can stabilize a cloud that might have been unstable in the absence of coupling.Progressive increases in coupling strength η tend to transform a disorganized system of clouds to one exhibiting a high degree of order.However, beyond this point, further increases in coupling strength result in lower coherence and eventually explosive growth (Fig. 9).Sometimes behaviour is quite complex (e.g. in Fig. 7 oscillators evolve from being incoherent to an ordered state in which adjacent oscillators are out of phase, only to migrate back to incoherence).
One of the exciting aspects of this study is that it opens up the possibility for exploring spatial cloud patterns in a relatively simple framework.The two-dimensional rigid cloud structure (early stages of Fig. 10) is similar to the rigid structures generated by large eddy simulation and observed by geostationary satellite for non-precipitating clouds; similarly the spatiotemporal oscillations of the cloud field (later stages of Figs. 10 and 12) are reminiscent of the same structures simulated by large eddy models and observed by satellites under precipitating conditions.This presents opportunity for improved understanding of the robustness of these states and the potential for transitions between them -a topic of great interest for climate studies.

Fig. 1 .
Fig. 1.Schematic showing how each cloud element is coupled to its neighbours via airflow associated with up or downward movement.For simplicity, the cloud elements are depicted as separate clouds, although for stratocumulus clouds, cloud cover might be full.Note that growing clouds are positively buoyant and are associated with upward movement.Raining clouds are negatively buoyant as a result of the evaporative cooling of precipitation.Downdrafts reaching the surface are forced to diverge.

Fig. 2 .
Fig. 2. One-dimensional solution to the equations for I = 21, τ c = 5 min (∼ 0.25× the microphysical delay T ), and η = −0.51.Red lines: i = 5 and i = 7; black lines: i = 10 and i = 12.Initially all clouds are synchronized with H 500 m.At about 380 min red lines and black lines diverge.Soon after, odd oscillators (red) are in phase, but out of phase with even oscillators (black).According to Fig. 1, this represents a cloud growing at the expense of its neighbours.

Fig. 5 .
Fig. 5.As in Fig. 4 but with |η| reduced to −0.53.The small reduction in |η| results in incoherent oscillations and suggests a threshold for synchronization.

Feingold
Fig. 7.As in Fig. 6, but for | | reduced to -0.621.(a) 0 < t < 1000 min; prior to ∼ 400 min oscillators are in phase but after that are i (b) 1400 < t < 1600 min, adjacent oscillators prefer to be out of phase; (c) 1800 < t < 2000 min, a return to incoherent oscillatio Fig. 7.As in Fig. 6, but for |η| reduced to −0.621.(a) 0 < t < 1000 min; prior to ∼ 400 min oscillators are in phase but after that are incoherent; (b) 1400 < t < 1600 min, adjacent oscillators prefer to be out of phase; (c) 1800 < t < 2000 min, a return to incoherent oscillations.

Fig. 8 .
Fig. 8. Illustrative example for a one-dimensional array of 41 clouds and η = −0.53:(a) time series for i = 1 (blue), 5 (red), 35 (black); (b) cross-correlation matrix where the position in the matrix indicates the maximum correlation r between any pair of clouds; (c) the cross-lag matrix showing the time lag at which the maximum r in (b) occurs.

Fig. 9 .
Fig. 9. Summary of the influence of η on an array of 41 clouds for η = −0.3000,−0.4750, −0.5300, and −0.5506 (top to bottom in any column).Column 1: average power spectral density in the time domain; Column 2: phase and correlation of the 41 clouds relative to the frequency and phase of cloud 1; Column 3: histogram of the correlations comprising the cross-correlation matrix (e.g.Fig. 8b); Column 4: two-dimensional probability distribution function of the correlation and lag based on the cross-correlation and cross-lag matrices.The progression to larger |η| is initially associated with an increase in system coherence and correlation (first 3 rows) but at sufficiently large

Fig. 10 .
Fig. 10.Two-dimensional simulations (21 × 21 oscillators) with τ c = 43 min (2 × T ), η = −0.20 ± 0.03 and H 0 = 670 m ± 100 m.The display shows contours of H i,j as a function of i (x distance; abscissa) and time (ordinate).The linear colour scale is from 0 to 800 m (purple through red/brown) Model output has been averaged over all j -grid points.Prior to about 300 min the all components are synchronized and rise and sink in unison (see also Fig.11).Thereafter, the spatial symmetry is broken by the onset of rain and a gridded structure merges depicting spatiotemporal oscillations in H i,j .
Fig. 11.Snapshots of model output associated with Fig. 10 at (a) 200 min, (b) 250 min, and (c) 300 min.Horizontal axes are x and y distances (or grid number i, j ).Vertical axes are height in m.Note the temporal invariance of the structure.
Fig. 12.As in Fig. 11 but for larger domain (I = 41) and η = −0.22 ± 0.03 at (a) 840 min, (b) 880 min, and (c) 920 min.Note the rearrangement of the cells as regions of high H rain out and convert to low H .