Exploring the nonlinear cloud and rain equation Open Access

Marine stratocumulus cloud decks forming over dark, subtropical oceans are regarded as the reflectors of the atmosphere.¹ The decks of low clouds 1000s of km in scale reflect back to space a significant portion of the direct solar radiation and therefore dramatically increase the local albedo of areas otherwise characterized by dark oceans below.^2,3 This cloud system has been shown to have two stable states: open and closed cells. Closed cell cloud systems have high cloud fraction and are usually shallower, while open cells have low cloud fraction and form thicker clouds mostly over the convective cell walls and therefore have a smaller domain average albedo.^4–6 Closed cells tend to be associated with the eastern part of the subtropical oceans, forming over cold water (upwelling areas) and within a low, stable atmospheric marine boundary layer (MBL), while open cells tend to form over warmer water with a deeper MBL. Nevertheless, both states can coexist for a wide range of environmental conditions.^5,7 Aerosols, liquid or solid particles suspended in the atmosphere, serve as Cloud Condensation Nuclei (CCN) and therefore affect the concentration of activated cloud droplets.⁸ Changes in droplet concentration affect key cloud properties such as the time it takes for the onset of significant collision and coalescence between droplets, a process critical for rain formation. The onset of significant collision-coalescence process can thus be represented by a delay factor.^9,10

The emergent behavior of the coevolution of cloud and rain has been shown to be captured by a set of dynamical equations with a delayed sink term.¹¹ Numerical analysis of these equations yields bifurcation points that separate different dynamical regimes. The first point marks a shift from a steady-state (stable fixed point) in which the rain consumes the cloud at the exact rate of cloud replenishment, to oscillations (limit cycle) of stronger rain that depletes the cloud that created it and therefore dissipates until the cloud is thick enough to reform rain. These results were shown to provide insights into naturally occurring closed and open cells in marine stratocumulus cloud systems,¹² and the processes underlying transitions between these states. The oscillating branch of the solutions represents open cell clouds that typically produce stronger rain,^1,13,14 and as the rain depletes the cloud water and suppresses the updraft by evaporating below cloud base, the average cloud tends to last for shorter durations.⁶ The closed cells tend to produce very little drizzle and their morphology remains stable for more than 10 h despite the fact that their characteristic scale suggests a theoretical lifetime of only ∼1 h.^12,15,16 The one-dimensional, time-delay equations for cloud thickness (H) and for droplet concentration (N) were later coupled by a spatial dynamical feedback to a set of oscillators that produce spatial patterns of cellular convection similar to the ones produced in detailed cloud resolving models and seen in nature.¹⁷ 

The nature of the transitions from open to closed cellular convection has been studied using both large eddy simulation (LES) and the cloud and rain equations.¹⁸ This study showed that the transition between closed and open cellular cloud states shows hysteresis as function of the aerosol loading. Such behavior is expected in a Delay-Differential-Equation (DDE) as the solution depends on the history of the delayed element.

Here, we explore the cloud and rain equation response to small perturbations around the steady state, which allows analytical exploration of the nature of the damped oscillations toward the fixed points, and of the first bifurcation point that marks the transition from a fixed point to a limit cycle. For completeness, we further explore the transition toward chaotic behavior. This transition occurs in a nonphysical regime in the current simple model, yet it is possible that subsequent studies will find such a transition in a physical regime of somewhat more detailed model. The linear stability analysis allows us to better understand transitions between states and how they depend on changes in aerosol and the environmental variables, which provides new insights into the complex interactions between clouds, aerosol, and rain. As aerosol concentration strongly affects the cloud droplet concentration (N), changes in the parameter N imply here changes in the aerosol loading.

Time delay differential equations (DDE) are used in many dynamical systems,¹⁹ including population dynamics,²⁰ neural networks,²¹ El Nino,²² and more. As even a simple linear DDE requires an infinite number of initial conditions to initialize the delay term, it is formally equivalent to a partial differential equation or to an infinite system of ordinary differential equations in terms of its number of degrees of freedom. DDEs can display complex behavior, some of them chaotic, and therefore often cannot be solved analytically.¹⁹ 

Using the notation $H(t−D)$ for the $H$ value in D time units before the time t, the cloud and rain equation can be coupled to an aerosol equation as follows:¹¹ 

where the time dependent variables are the cloud depth (H), and the droplet concentration (N). R is the rain-rate, which is a function of H and N, and $α≈100 [day−1 m−2.5]$ is a scaling constant that links cloud depth, droplet concentration, and rain rate. The $α$ value was determined both theoretically and from measurements.^23–25 The environmental and droplet (aerosol) conditions are represented by H₀ and N₀ as the cloud-depth and droplet concentration carrying capacities, τ and τ₂ are characteristic times for reaching the carrying capacity values under no sink conditions, D is the delays that represents past states that control the current sinks. Note that following theoretical and modeling studies the delayed rain sink term depends on the inverse of the square root of N.²⁶ 

Here, we consider the basic cloud and rain equation in which N is a free parameter

Such a representation of the problem assumes that changes in the aerosol concentration are relatively small. It represents cases for which the source of aerosols is steady or when the aerosol consumption by drizzle is relatively small. This reduces the problem to a first order nonlinear DDE controlled by four parameters: (i) H₀—the cloud carrying capacity parameter that represents the systems maximal potential for cloud depth. (ii) τ—the characteristic cloud recovery time. (iii) N, which controls the strength of the sink term (rain), and (iv) D—the time delay that represents the time it takes to convert cloud water into rain by stochastic microphysical collection processes.¹¹ 

Equation (1) can be nondimensionalized by replacing H with the normalized height $h=H/H0$ and t with $t*=t/τ$⁠. When translating the equation for $dH/dt$ to $dh/dt*$⁠, model parameters N, H₀, and τ are replaced by a single parameter

and Eq. (1) is therefore transformed to a simpler nondimensional form

where $D*=D/τ$ is the nondimensional delay. Throughout this paper, we will refer to the nondimensional version of the problem (Eq. (3)) as the Cloud and Rain (C&R) equation. When solving Eq. (3) for the steady state case (⁠$hsts$⁠), for which the derivative (⁠$dhdt*$⁠) vanishes, $h=hsts$ is constant in time and the delay does not play a role in determining $hsts$⁠. In such a case, Eq. (3) reduces to a simple one-parameter polynomial whose physical solution (allowing only positive h values) is

The fact that N^1/2 and the cloud depth carrying capacity appear only via their ratio in a single nondimensional parameter (⁠$μ$⁠) has an important implication that will be discussed later. We note the $hsts$ is a fixed point that can be either stable (to which the C&R equation converges in the steady-state) or non-stable. For all cases, if the equation's initial conditions are equal to $hsts$ throughout the delay period, the solution will remain $hsts$ for all later times. Next, we will use linearized stability analysis to explore the system's response to small perturbations from $hsts$⁠.

When the delay is set to zero the C&R equation becomes a first order, nonlinear ordinary differential equation with a quadratic term (a Riccati equation)²⁷ that can be transformed to a second order linear ordinary differential equation and has an analytical solution of the form

where $k1$ is an integration constant and $ς=4μ+1$ is a positive nondimensional number. Prescribing the initial conditions $h(t=0)=0$⁠, to show how the cloud develops in time, yields $k1=2ςtanh−1(−1/ς)$⁠. The hyperbolic tangent function approaches 1 for positive argument, or −1 for negative argument, with a sharp transition near zero. For large enough t*, the hyperbolic tangent component can therefore be replaced with $−1$ such that the solution for h converges to the steady state solution of Eq. (4). The transition in time from: $h=0$ to $h=hsts$ is smooth and monotonic with no oscillations, as for an over-damped oscillator. This behavior will be linked to the stability analysis of the fuller equation with delay in Section IV.

As stated earlier, $hsts$ is a fixed point for all parameter values. We can therefore perform a stability analysis of the C&R equation to investigate the response to small perturbations around the fixed point for different values of the model parameters.

Let δ be a small perturbation around the fixed point $hsts$ such that $h=hsts+δ(t)$⁠, thus $dδdt*=dhdt*$⁠. Linearizing around $hsts$ and neglecting terms nonlinear in δ yields

Expressing δ as an exponent $δ=eβt*$ allows separation of the contribution of the delay, and transforms Eq. (6) to

thus yielding a transcendental equation for $β$

The exponent $β$ may, in general, be complex. Its real part, $Re{β}$⁠, determines the stability of the perturbations, and its imaginary part, $Im{β}$⁠, the frequency as the solution converges or diverges from the fixed point. For $Re{β}<0and Im{β}=0$ the system is in a state similar to an overdamped oscillator converging to the fixed point with no oscillations (a continuation of the case of no delay). The case $Re{β}<0 and Im{β}>0$ describes decaying oscillations toward the steady state. $Re{β}>0$ generates an unstable fixed point where the perturbation is amplified, leading to steady oscillations. Points in the parameter space for which $Re{β}=0$ describe the first bifurcation in which the system transforms from having a stable to a non-stable fixed point.²⁸ 

The transcendental equation for $β$ (Eq. (8)) has a closed form solution based on the Lambert W function that solves $W(z)eW(z)=z$⁠, which is often used for DDE analysis²⁹ 

Inserting the expression for $hsts$ from Eq. (4), into the argument of the Lambert W function in Eq. (9) (⁠$ξ)$ yields

For the physical parameter range, the argument $ξ$ is always negative and is composed of two factors, one a function of μ and one of $D*$⁠. $ξ$ grows monotonically with μ (becomes less negative and smaller in absolute value) and decreases with $D*$⁠. The main branch of the Lambert W function is a real and negative number as long as its argument is real, negative number equal or larger than $−e−1$ (where it reaches its global minimum $W(−e−1)=−1$⁠). Therefore, since $ξ$ is always negative, Eq. (9) reveals that as long as $ξ≥−e−1$⁠, $β$ will have a real and negative value. For such cases, small perturbation $δ=eβt*$ decays exponentially to the fixed point $hsts$ with no oscillations (an overdamped regime). Inserting the $ξ$ expression (Eq. (10)) and solving this condition explicitly shows that as long as

$β$ is real and negative. For a given nonzero value of μ, as $D*$ approaches zero, the left side of Eq. (11) increases rapidly such that the condition is always fulfilled and the solution is similar to the no-delay case of $D*=0$⁠, which we showed to be overdamped (Eq. (5)). We define a critical delay (⁠$Dc*$⁠) which satisfies the equality $1+D*eD*+1D*eD*+1=μ+4μ$⁠, for which $ξ=−e−1$ and therefore where $Re{β}$ reaches its minimal value of

and where the solution's convergence to the fixed point is fastest.

For $ξ<−e−1$⁠, $Im{β}$ is not zero. The oscillations decay toward the stable fixed point as long as $Re{β}<0$⁠, namely, based on Eq. (9), as long as $Re{W(ξ)}<D*$⁠, while the first bifurcation point, i.e., the transition to a limit cycle state, occurs when $Re{β}=0$⁠, namely, $Re{W(ξ)}=D*$⁠.

Writing the explicit equation for the value of the model parameters at the bifurcation point

reveals an equation that resembles one of the definitions of the Lambert W function (i.e., $W(zez)=z$⁠) with an additional negative scalar $ψ=(1−1+4μ)$⁠. The argument of the Lambert W function in Eq. (13) is therefore a negative real number. If the solution for $W(ψD*eD*)=z,$ and $z=x+iy$⁠, solving and for $x$ and $y$ while requiring that $Im{zez}=0$ yields $x=D*$ and $ψ=1 cos y$⁠, which implies that $|ψ|≥1$ (for a more detailed proof see the lemma in the  Appendix). Therefore Eq. (13) can be satisfied only when $|1−1+4μ|≥1$⁠, namely, when μ $≤$ 4/3. We define this limit as $μlim=4/3$⁠. The bifurcation points exist only when $μ<μlim$⁠. Therefore, any point for which $μ≥μlim$ is a stable fixed point.

Figure 1 shows β as a function of the two controlling parameters: μ and D*. Five distinct features can be seen in this parameter-space (the μ vs. D* space):

• The area in which $Im{β}=0$⁠, and therefore where the system is overdamped, is shaded light blue in panel (c), revealing that for lower delay values the perturbation decays toward the fixed point $(hsts)$ with no oscillations.

• The line along which $Re{β}$ has a local minimum is shown in panel (a) (narrow blue “valley”) and marked as a blue contour in panel (c), marking the edge of the overdamped regime. Points along this line mark the fastest decay of the perturbation toward $hsts$ for a given value of μ. For a given μ, the corresponding delay values along this line are defined as the critical delays (⁠$Dc*$⁠), which can be calculated by Eq. (12).

• The area in which $Re{β}<0 and Im{β}>0$ is between the critical delay line and the black contour shown on the left panel (white area in panel (c)). Values of $D*$ and μ in this regime will lead to damped oscillations. The shading of the period as function of the two model parameters (Fig. 1, panel (b)) reveals a decrease of the period near the critical delay line up to a local minimum followed by a monotonic increase (reduction of the frequency) as a function of the delay.

• Solutions corresponding to points in the parameter domain to the right of the vertical $μlim$ line (i.e., $μ>μlim$⁠, purple vertical line in panel (c)) can only be of stable fixed point type.

• The first bifurcation line is marked by the black contour line, for which $Re{β}=0$⁠, and therefore $Re{W(ξ)}=D*$⁠. This line shows the transition at which $hsts$ shifts from a stable to a non-stable fixed point. For small D* values, the bifurcation point described in Eq. (13) occurs for small μ values that monotonically grow as D* increases approaching the μ = 4/3 limit.

• The area in which $Re{β}>0$ is left of the black contour (panel (c)). In this part of the parameter space, $hsts$ becomes an unstable fixed point and therefore perturbations from the fixed point will shift the solution to a limit cycle or, as will be shown later, to a period doubling route to chaos.²⁸ 

To further explore the nature of the transition from stable fixed point to non-stable, we run a set of numerical simulations around the first bifurcation point. Because $hsts$ is independent of $D*$⁠, it is convenient to change it while holding μ constant such that $hsts$ is unchanged (see Eq. (4)). Figure 2 shows that the equation undergoes a supercritical Hopf bifurcation, showing that as the delay is increased, the amplitude of the oscillation indeed increases as $∼D*−D0*$⁠, away from the bifurcation point $(D0*)$⁠. For completeness, we also present the evolution for the non-physical regime for which the lower h values are negative, and find that the solution undergoes a period-doubling route to chaos, and the ratios of the distance between the values of D* at the period-doubling points converge to the Feigenbaum constant.^30,31

Recall that we solved here the nondimensional form of the C&R equation for which 3 out of the 4 model parameters were replaced by μ, reducing the dimensionality to only 2 nondimensional parameters. However, the interpretation of the results in the original dimensional parameter space is not always straightforward, because changing a dimensional parameter corresponds to moving along some curve in the 2d μ vs. D* nondimensional parameter space. In particular, we note that both nondimensional parameters D* and μ are normalized by $τ−1$ while μ is also normalized by $N/H0$⁠. Varying $N/H0$ while keeping other parameters fixed corresponds to a straight horizontal line parallel to the μ axis in the μ vs. D* space, showing that as $N/H0$ increases (i.e., increase in μ) the system's state will shift from an unstable to a stable fixed point. Changes in the dimensional delay D correspond to a straight vertical line in the μ vs. D* space, parallel to the D* axis (e.g., the dotted orange vertical line on Fig. 1, panel (c)) for which μ is held constant. An increase in D leads to a transition from stable to non-stable fixed points. As shown in Eq. (13), the D* values for which the Hopf bifurcation point occurs depend on μ, and for μ > μ_lim, the equation will have a stable fixed point for any D*. Varying τ while holding all other variables fixed is described by an inclined line in the μ vs. D* space, approaching the origin (when $τ→∞)$ for which the slope depends on $N/H0$⁠. As illustrated in the dotted magenta line in Fig. 1 (panel (c)), some of these lines can cross the Hopf bifurcation contour twice, showing a transition from stable fixed point to non-stable and then to stable again, as τ increases.

To illustrate how β affects the solution, Fig. 3 shows simulations for parameters corresponding to two lines in the μ vs. D* parameter space: first (dotted orange line in Fig. 1 panel (c)), when D is the free variable and all other parameters are fixed (similarly to the runs in Fig. 2 for which H₀ = 1000 m, N = 16 cm⁻³, τ = 20 min, therefore μ = 0.29, yielding a constant fixed point), and second, when τ is the free variable (H₀ = 1000 m, N = 16 cm⁻³, and D = 15 min, marked by the dotted magenta line in Fig. 1, panel (c)). The two numerical simulations were run for five values of β each (marked in Fig. 3, upper panels) controlled by changes in D (left column) and τ (right column), demonstrating qualitatively different solutions for β and therefore different stability regimes of the dynamics around the steady state. Each simulation was initiated with the fixed-point cloud height, $H=Hsts$⁠, and the thickness therefore remained constant in time, until a small perturbation was introduced. As shown in Eq. (4), $hsts$ is a function of μ only; hence for the cases in which the delay is the free parameter, $hsts$ remains the same for all simulations (Fig. 3 left column). For the first 3 delay varying simulations (Fig. 3, panels (b)–(d)), $Re{β}<0$⁠, and therefore, the perturbations decay back to the $Hsts$ which is a stable fixed point. In the case of panel (b), $Im{β}=0$⁠, and therefore, the system is in the overdamped regime and solutions decay to $Hsts$ with no oscillations. As the values of $β$ are closer to the bifurcation point, the decay to the fixed point is slower. On the same note when $Re{β}>0$ (panels (e) and (f)), the closer the $β$ values are to the bifurcation point, the slower is the shift to a steady limit cycle state. The τ varying simulations (Fig. 3 right column) show that the system can enter and exit the non-stable fixed point regime crossing the Hopf bifurcation contour twice, and therefore will have a stable fixed point for relatively small τ values (panels (h) and (i)) and for relatively high τ values (panel (l)) and between them be in the limit cycle state (panels (j) and (k)).

In this paper, we consider warm marine stratocumulus clouds for which all hydrometeors are liquid (cloud droplets and raindrops). Cloud droplets nucleate on aerosol particles that serve as Cloud Condensation Nuclei (CCN). Raindrops form when larger droplets collide and collect the smaller ones. The efficiency of transferring cloud droplets to rain correlates positively with the first two moments of the droplet size distribution, i.e., the average size and the variance. An increase in the aerosol concentration increases the concentration of activated cloud droplets. This implies that more droplets are competing for the available vapor and therefore their average size will decrease.^8,32 Such changes incur a suppression of drop collection, so that variance decreases. The result is a delay in the onset of rain formation, and for marine stratocumulus, a decrease in the rain-rate up to a point of complete rain suppression.^9,33

We have analyzed the cloud and rain (C&R) equation for which the droplet (aerosol) concentration is assumed fixed and is prescribed. Such an assumption can be justified when the source of aerosols is steady (either local or due to long range transport) and when the aerosol consumption by drizzle is relatively small. In such cases, fixed N (often also assumed in cloud resolving models) allows one to study the coevolution of cloud and rain in a model of reduced complexity.

The C&R question has four parameters (τ, D, H₀, and N). When nondimensionalizing the equation, the environmental parameters τ and H₀ and the droplet concentration N are replaced by a single parameter $μ=N/(ατH0)$⁠. Apart from reducing the complexity, this offers important physical insights into the interplay between the environmental and aerosol properties with respect to the system's stability.

We performed linear stability analysis which yielded an analytical expression for the stability as a function of the two parameters (μ and D*). The C&R equation forms a very rich solution space, ranging from a stable over-damped convergence through fixed point, to a limit cycle, to chaotic behavior. Each solution type occupies a distinct regime in the μ vs. D* parameter space.

The nondimensional delay D* and μ have opposite effects on the stability. Larger μ values (larger aerosol concentration, all other variables being equal, hence smaller sink term) imply a more stable solution, while an increase in D* leads to an instability of the steady solution. In the limit we have shown that any point for which μ > $μlim$ yields solutions that are characterized by a stable fixed point.

For a small D* and a large μ, perturbations decay back to the steady solution with no oscillations. For larger D*, the system shifts, through a critical delay $Dc*$⁠, to a damped oscillatory regime, and later, crossing the first bifurcation point, to a limit cycle regime, and finally to a chaotic regime.

As expected, an increase in the aerosol and therefore in the droplet concentration N decreases the rain-rate. Reduction in the rain-rate can shift the system state from a limit cycle to coexistence of steady (weak) rain and fixed cloud thickness (stable fixed point solution). The μ parameter is linear in $N/H0$ and does not depend on D*. This suggests that clouds forming in atmospheric conditions that allow thicker clouds will shift from a limit cycle state to a steady-state under higher aerosol concentrations, and therefore that the shift from open to closed cells would occur at higher aerosol concentrations for a thicker marine boundary layer.

Figure 4 illustrates how $μ$ values depend on N and H₀ for a given characteristic time for cloud recovery time (τ = 20 min). μ_lim is marked by the magenta contour. The $μ<μlim$ subspace, to the right of the contour, represents $μ$ values where $Re{β}$ values can have positive values and so the system can have unstable fixed points (depending on the D values). For larger τ values, the μ_lim contour shifts to the left, expanding the relative area over which the system can be unstable.

However, within the possibly unstable subspace in which μ < $μlim$ the system's response to changes in τ were shown to exhibit non-monotonic behavior, namely, an increase in τ will shift the solution from a stable to an unstable, and back to a stable fixed point state. Both D* and μ are linear in $1/τ$⁠, and therefore, varying τ corresponds to moving along a straight, inclined line via the origin in the μ vs. D* parameter space (dotted magenta line in Fig. 1(c)). The slope of such lines depends on the value of $N/H0$⁠. $$ Smaller $N/H0$ values imply a steeper slope and therefore a longer path in the non-stable fixed point regime, again showing that smaller aerosol concentrations reduce the size of the parameter regime in which the system is in a steady state.

We acknowledge that such a simplified description of the interplay between cloud and rain cannot capture the full complexity of these interactions in natural systems. Moreover, some of the physical processes that are important in observed clouds are not included here (radiation, for example). Nevertheless, in the spirit of simple dynamical systems analogues to complex systems, the analytical and numerical analyses shown here offer new perspectives on the role and interplay of the main physical parameters. The work demonstrates how changes in these parameters can shift the system between different regimes within the solution space. The existence of such regimes could be explored in the future with more detailed numerical models.

The research leading to these results received funding from the European Research Council (ERC) under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 306965 (CAPRI). E.T. is funded by the NSF Physical Oceanography program, Grant No. OCE-1535800, and thanks the Weizmann Institute for its hospitality during parts of this work.

Lemma 1: An extended view of the solution form of the Lambert W function for any real number. The Lambert W function maps $zez→$ $z$⁠. If

then

Expressing $z=x+iy$ and requiring that $Im(zez)=0$ implies the following link

Inserting Eqs. (A3) in (A2) yields

Therefore, Eq. (A1) becomes

Eq. (A4) suggests an extended view on the Lambert W function adding a factor of $1 cos y$ to the standard $xex$ kernel of the function. Moreover, $|1 cos y|≥1, ∀y$⁠. For the cases that $η≥−1e$⁠, a real solution for $W(η)$ exist for which y = 0 and $1 cos y=1$⁠, such that the solution collapse to the standard form of $W(η)=x or η=xex$⁠. In the general case, however, the fact that $|1 cos y|≥1$ implies that for $η>1e$⁠, $W(η)≥Re{W(−η)}$⁠. Such representation for the cases of real arguments and the derived results may have implications in other cases where solutions of the W function are involved. We note that the solution presented in Eq. (5) is not unique and that there are other W function's branches that yields x,y pairs that can satisfy the equation. The above theorem is applicable to all of these x,y pairs.