The Physics of Falling Raindrops in Diverse Planetary Atmospheres

The evolution of a single raindrop falling below a cloud is governed by fluid dynamics and thermodynamics fundamentally transferable to planetary atmospheres beyond modern Earth's. Here, we show how three properties that characterize falling raindrops -- raindrop shape, terminal velocity, and evaporation rate -- can be calculated as a function of raindrop size in any planetary atmosphere. We demonstrate that these simple, interrelated characteristics tightly bound the possible size range of raindrops in a given atmosphere, independently of poorly understood growth mechanisms. Starting from the equations governing raindrop falling and evaporation, we demonstrate that raindrop ability to vertically transport latent heat and condensible mass can be well captured by a new dimensionless number. Our results have implications for precipitation efficiency, convective storm dynamics, and rainfall rates, which are properties of interest for understanding planetary radiative balance and (in the case of terrestrial planets) rainfall-driven surface erosion.


Introduction
Within a planetary condensible cycle, precipitation is the transport of the condensible species in a condensed phase (liquid or solid) through the atmosphere and, for terrestrial planets, to the surface. Extensive vertical displacement relative to the local air mass distinguishes precipitation from clouds. Because precipitating particles can fall far from the air mass where they form, they redistribute both heat and the condensible species within an atmosphere. Precipitation is a transient state, but though its effects are largely indirect, they have immense consequences for planetary climate.
On Earth, global precipitation patterns play a critical role in determining local ecology and have significant societal impacts (Margulis, 2017). On terrestrial plan-ets generally, the intensity, frequency, and spatial distribution of liquid precipitation are essential in governing surface erosion via runoff and physical weathering (e.g., Margulis, 2017) as well as chemical weathering fundamental to the carbon-silicate cycle (Walker et al., 1981;Macdonald et al., 2019;Graham & Pierrehumbert, 2020). Finally, interpreting solar system geological records shaped by fluvial erosion-e.g., ancient Mars' large-scale valley networks and crater modifications (e.g., Craddock & Howard, 2002), modern Titan's lakes and rivers (e.g., Lorenz et al., 2008), and Archean Earth's fossilized raindrops (Som et al., 2012;Kavanagh & Goldblatt, 2015)-requires an understanding of changes in precipitation events as planetary conditions vary.
Despite the importance of precipitation, little progress has been made on how precipitation physically behaves in different planetary environments (e.g., Vallis, 2020). Previous studies have tended to view precipitation behavior primarily as a function of cloud formation and evolution. Cloud physics is complicated by extreme nonlinearities, gaps in theoretical understanding of fundamental processes bridged only by empiricisms, and dependencies on spatial and temporal scales that span many orders of magnitude. In turn, planetary clouds are studied via a hierarchy of models with varying tradeoffs among complexity, robustness of included processes, and ease of interpretation (e.g., Rossow, 1978;Carlson et al., 1988;Ackerman & Marley, 2001;Lee et al., 2016;Powell et al., 2018;Gao et al., 2018). Still, climate models of terrestrial planets commonly represent clouds and precipitation using modern-Earth-tuned parameterizations with essentially ad hoc parameter sweeps (e.g., Wordsworth et al., 2013;Urata & Toon, 2013;Yang et al., 2014;Komacek & Abbot, 2019).
Given the transition from cloud to precipitation is complex and poorly modeled even on the well-observed modern Earth (e.g., Rogers & Yau, 1996;Pruppacher & Klett, 2010;Flato et al., 2014), any insights on precipitation behavior become almost impossible to extract when precipitation is considered only as an afterthought of cloud behavior. An alternative approach is to consider the behavior of individual precipitating particles independently of their formation conditions. This strategy is much more tractable in different planetary environments because individual precipitating particles are governed by thermodynamics and fluid dynamics that are both relatively well understood and fundamentally transferable to planetary regimes beyond modern Earth's.  took this approach and used key properties of individual methane-nitrogen raindrops on Titan to highlight fundamental differences between rainfall on Titan and Earth and hypothesize consequences for storm intensities.
Previous planetary science studies have also attempted to use the simplicity of raindrop physics to place constraints on paleo-air pressures on Archean Earth (Som et al., 2012;Kavanagh & Goldblatt, 2015) and early Mars (Craddock & Lorenz, 2017;Palumbo et al., 2020) via maximum raindrop sizes before breakup. This use of raindrop physics hints at the possible productivity of this approach. However, even maximum raindrop size has not been considered systematically in general planetary atmospheres before, and the conclusions reached by these studies have in some cases been inconsistent. Specifically, a recent study by Craddock and Lorenz (2017) reached opposite conclusions versus other studies (e.g., Komabayasi et al., 1964;Clift et al., 2005;Pruppacher & Klett, 2010;Som et al., 2012;Kavanagh & Goldblatt, 2015;Palumbo et al., 2020) on the dependence of maximum stable raindrop size on planetary air density, leading to entirely different conclusions about paleo-air densities.
In this paper, we build on  to establish a comprehensive and generalized picture of the "life and death" of a single raindrop over a wide range of planetary conditions. Like , we neglect the "birth" of the raindrop (i.e., the growth of a cloud particle of negligible vertical velocity to a precipitating particle below a cloud, essentially done growing). Our approach here is distinct from previous work as we present fully generalized methods, tie each component of our methodology back to fundamental physics, and focus on how the well-understood behavior of an indi-vidual raindrop can provide insight into the rest of the condensible cycle in different planetary environments. This work lays a foundation for building physically driven microphysics parameterizations for generalized mesoscale models and global circulation models (GCMs).
We limit our focus to liquid precipitating particles ("raindrops") because they have a unique shape for a given mass of condensible. The shape degeneracy of solid precipitating particles is a major challenge (e.g., Pruppacher & Klett, 2010, chapter 2.2) that we do not treat here. However, once the non-uniqueness of solid precipitating particle shapes is addressed, our methodology is applicable to solid particles as well. Water is the most familiar condensible species, but all methodology we present here is generalized for any liquid condensible, e.g., methane-nitrogen raindrops on Titan or iron raindrops on WASP-76b (Ehrenreich et al., 2020). Except when we assume the existence of a planetary surface, our methodology is also general to both terrestrial and gaseous planets.
In section 2 we present methods to calculate falling raindrop shape, terminal velocity, and evaporation rate in a generic atmosphere. We show how these characteristics can place upper and lower bounds on raindrop size in section 3. Section 4 uses the methodology developed to probe raindrop characteristics and size bounds in different atmospheres. In section 5 we discuss the implications of our results for different planetary atmospheres and microphysics parameterizations as well as possible extensions of this work to solid precipitating particles. We summarize our results in section 6.

Raindrop Characteristics
In isolation, a precipitating particle does two things: (1) it falls and (2) it evaporates. To calculate the rate at which a particle falls and evaporates requires knowledge of the relationship between particle mass and shape. Unlike solid precipitating particles, whose forms vary widely, raindrops have equilibrium shapes that can be uniquely calculated for a given mass of liquid condensible, air density, and surface gravity. A unique shape allows us, in a known external environment, to describe a raindrop with only a single size variable. Here, we use equivalent radius r eq , which is the radius a raindrop of mass m would have if it were spherical, i.e., where ρ c, is the density of the liquid condensible.

Raindrop Shape
Falling raindrops adopt a range of shapes depending on their size-though never the teardrop shape inscribed in the public imagination (Blanchard, 2004). As raindrops grow in mass, they evolve from spheres to oblate spheroids to shapes resembling the top of a hamburger bun (e.g., . (An oblate spheroid is generated from an ellipse rotated about its minor axis.) Spheres are merely a specific subset of oblate spheroid, and the more complex shapes have virtually indistinguishable dynamical properties from oblate spheroids Szakáll et al., 2010); therefore, we simplify our shape calculations by prescribing that raindrops take the shape of an oblate spheroid with semi-major axis a (oriented perpendicular to the fall direction) and semi-minor axis b (oriented parallel to the fall direction). Given this assumption, we can describe raindrop shape with only an axis ratio b/a and make use of many existing analytic expressions. Geometric properties of oblate spheroids used in this paper are given in Appendix A.
In equilibrium, the deviation of a raindrop surface from a minimum-energy-state sphere can be calculated considering the first law of thermodynamics: the change in energy from the surface's increased surface tension must be balanced by the work done to expand the surface's enclosed volume into a region of different pressure (e.g., Pruppacher & Klett, 2010, chapter 10.3.2). We follow  in accounting for the key pressures at raindrop equator and assuming an oblate spheroid shape (see Appendix B for more detail), which gives Here σ c-air is the surface tension between the liquid condensible and air, g is the local gravitational acceleration, and ρ air is the local air density. For a given r eq , equation (2) can be solved numerically for b/a to characterize the raindrop's shape.

Raindrop Velocity
The two key forces acting on a falling raindrop in a planetary atmosphere are the gravitational force F g and the aerodynamic drag force F drag . The effective gravitational force on a raindrop is after accounting for the raindrop's buoyancy within the air fluid. Drag force on a raindrop is where C D is the drag coefficient of the raindrop, A is cross sectional area of the raindrop, ρ air is the local air density, and v is the raindrop's fall speed relative to air. Raindrop cross sectional area is a function of raindrop shape and size. C D is a function of raindrop shape and flow regime. The latter can be characterized by the dimensionless Reynolds number Re where is a characteristic raindrop length scale that we take to be 2r eq and η air is the dynamic viscosity of local air. Viscosity varies with both air composition and temperature; we calculate η air in a generic air mixture following Reid et al. (1977) chapter 9.
Calculating C D theoretically for a general-shaped object in a general flow is not tractable (e.g., Stringham et al., 1969), and thus C D is typically evaluated via experimentally based parameterizations. As long as the flow regime is correctly captured by a scale analysis, defaulting to experimentally based parameterizations for generic atmospheric conditions in and of itself is not problematic-though limitations in the coverage of the parameter space used to fit expressions must be considered.
In this work, we use the drag parameterization Re 1 + 0.15Re 0.687 + 0.42 1 + 4.25 × 10 4 Re −1.16 −1 C shape (6) (Clift & Gauvin, 1970;Loth, 2008). C shape is a correction term to account for raindrop deviations from spherical shape fit by Loth (2008) for b/a ≤ 1 across a variety of falling object shapes: where f SA is the ratio of the surface area of the oblate spheroid raindrop to the surface area of a sphere of radius r eq . The formulation of C D in equation (6) is primarily based on the parameterization of drag for a sphere with Re < 3.5 × 10 5 by Clift and Gauvin (1970). C shape is most accurate when significant raindrop deviation from a sphere occurs within the Newtonian flow regime of Re between about 750 to 3.5×10 5 (Clift et al., 2005;Loth, 2008).
Compared to the other C D parameterizations we considered (Salman & Verba, 1988;Ganser, 1993;Hölzer & Sommerfeld, 2008), we found equation (6) with (7) best reproduced the velocity dependence on raindrop radius for Earth values with experimental validation over wide ranges of Re and oblate spheroid axis ratios. We note that we neglect corrections to C D for non-continuum regime effects (commonly referred to as the Cunningham or slip-flow correction factor) as we are not concerned with the behavior of very small ( 1 µm) particles (e.g., Seinfeld & Pandis, 2006, chapter 9.2). (Such corrections only become important when particle size becomes comparable to the mean free path of local air molecules.) Terminal velocity v T occurs when the raindrop is no longer accelerating, and the gravitational force F g is balanced by the aerodynamic drag force F drag . Under modern Earth atmospheric conditions, the timescale for raindrops of a fixed size to reach terminal velocity is very small compared to their lifetimes (Pruppacher & Klett, 2010, chapter 10.3.5). To test the generality of this rapid raindrop acceleration assumption, we numerically integrated water raindrop motion accounting for variable acceleration across a large range of plausible planetary conditions (air composition, surface gravity, air pressure, air temperature).
We find that the Earth-based empiricism is generally true: even the largest possible stable raindrops (described in section 3.4) with the longest acceleration timescales reach 99% of their terminal velocity after starting from rest within the first 1% of their total fall distance and 5% of their total fall time ( Figure S1). Further, when considering the effect of raindrop size changes due to evaporation on reaching terminal velocity, we find that the differences in raindrop fall speed between self-consistent treatment of raindrop acceleration and assuming terminal velocity is instantly reached are, at maximum, on the order of 10% and typically much smaller ( Figure S2). Thus, for simplicity, we henceforth make the standard assumption that raindrop falling speed relative to air v is the raindrop's terminal velocity, which can be uniquely determined for a given raindrop size and shape.
Equating F drag and F g and substituting the appropriate oblate spheroid geometry yields a terminal velocity of ρ air and g are known from planetary atmospheric properties. b/a is uniquely determined from r eq with equation (2). The nonlinear dependence of C D on v (through Re) requires we solve equation (8) numerically.
A raindrop's vertical speed relative to a planet's surface-dz/dt, its change in altitude z per unit time t-is the sum of the its velocity relative to air, here assumed to be v T , and the vertical velocity of the raindrop's local air w: We are focused on raindrop physics here and so treat w as a free parameter in the analysis that follows.

Evaporation Rate
Raindrop evaporation occurs when the atmosphere surrounding the drop is subsaturated in condensible gas. The preferred phase of the condensible molecules at the drop's surface becomes gas rather than liquid. As the condensible is transferred from the liquid phase in the raindrop to the gas phase in the air, the air closest to the raindrop surface deviates in temperature and relative humidity from the local atmospheric state; the relative humidity adjacent to the drop surface increases while the temperature drops due to the latent heat required for the liquid-to-gas phase transition. Both these effects serve to lower the thermodynamic impetus to evaporate. Thus, in addition to the environmental level of sub-saturation, the rate at which evaporation occurs is dictated by the rate at which heat and condensible gas can be transported away from the raindrop surface.
Quantitatively, the change in raindrop equivalent radius with time t can then be formulated from geometry and appropriate boundary conditions as (see Rogers & Yau, 1996, chapter 7 for a derivation). RH is the relative humidity of the local air; R is the ideal gas constant; µ c is the molar mass of the condensible in its gas phase; ρ c, is the density of the liquid condensible; D c-air is the diffusion coefficient for the condensible gas in air; p c,sat is condensible gas saturation pressure; T air is the local air temperature far from the drop's surface; T drop is the temperature at the raindrop's surface; and f V,mol is a ventilation factor that accounts for how much raindrop motion enhances condensible molecule transport relative to a stagnant drop.
Conservation of heat at the raindrop surface yields the following differential equation governing T drop : (Rogers & Yau, 1996, chapter 7). c p,c, is the specific heat at constant pressure of the liquid condensible; L c is the condensible's latent heat of vaporization at T drop ; K air is thermal conductivity of air; f V,heat is a ventilation factor that accounts for how much raindrop motion enhances heat transport relative to a stagnant drop. Without any further simplifications, equations (10) and (11) must be solved together numerically from initial conditions given their mutual dependencies. (We will discuss simplifications to calculating T drop in detail in section 3.3.) Relative humidity RH (also known as saturation) is defined as the ratio of the local condensible gas partial pressure p c to p c,sat at T air , i.e, RH ≡ p c /p c,sat (T air ). D c-air is a function of temperature, pressure, and air composition that we calculate following Reid et al. (1977) chapter 11 and Fairbanks and Wilke (1950). K air is a function of temperature as well as air composition that we calculate following the Eucken method (Reid et al., 1977, chapter 10). Note that the formulation of raindrop temperature in equation (11) only considers heat transport via conduction. For high temperature condensibles, heat transport via radiation will also need to be considered.
The ventilation factors arise from fluid dynamical effects not analytically calculable, so, as for C D , we must evaluate f V,mol and f V,heat from parameterizations based on experiments. Here, we choose to use the f V,mol parameterization of Beard and Pruppacher (1971) and Pruppacher and Rasmussen (1979): where Sc is the dimensionless Schmidt number defined as This parameterization is only experimentally validated for Re < 2600 but is hypothesized to be valid for spheres with Re < 8 × 10 4 based on theory (Pruppacher & Rasmussen, 1979 where c p,air is the specific heat at constant pressure for air. We neglect the effects of turbulence, which can act to increase ventilation. Raindrop shape impacts multiple aspects of ventilation, but in raindrop experimental data considered by Pruppacher and Rasmussen (1979), these shape effects cancel each other such that f V is independent of raindrop shape deformations at larger Re. (Note that this shape independence is only true for liquid raindrops, not solid condensibles (e.g., Pruppacher & Klett, 2010, chapter 13.3.2).) The above expressions for f V,mol and f V,heat are not definitive and could be improved by further experiments over a broader parameter space. The ventilation factors increase dr eq /dt by roughly an order of magnitude at large raindrop sizes, so they need to be accounted for despite the somewhat limited coverage of these parameterizations; as we will show in section 4, the qualitative impact of these uncertainties in ventilation for r eq (t) is often limited because large raindrops evaporate mass very gradually.
In our formulation of dr eq /dt in equation (10), we have made a number of simplifying assumptions that are valid here because we are concerned only with evaporation. (dr eq /dt can describe both evaporation and condensation depending on whether RH is less than or greater than 1.) Our focus on evaporation, rather than condensation, means we are not concerned with the behavior of very small drops, which, as we will show in Section 4, always evaporate rapidly compared to larger drops. We neglect corrections to relative humidity at the raindrop's surface due to surface tension and condensation nuclei solute effects-commonly known as Kelvin and Raoult effects, respectively (e.g., Lohmann et al., 2016, chapter 6)-that are important only at very small radii (r 1 µm). Note that in atmospheres where a gas component is soluble in the liquid condensible (e.g., N 2 in CH 4 on Titan; , solute corrections to RH cannot be neglected. (See, for example,  for how to extend what is presented here to such cases.) Consistent with our treatment of drag, we also neglect corrections to condensible vapor diffusivity and air thermal conductivity for non-continuum regime effects very close to the drop's surface, which is effectively equivalent to assuming that the mean free paths of air and condensible gas molecules are small compared to the size of the raindrop (e.g., see Lamb & Verlinde, 2011, chapter 8.2.2 for further discussion).
Beyond corrections only necessary for small drops, we also neglect corrections to the boundary conditions used to solve for dr eq /dt due to raindrop deformation from a sphere (Lamb & Verlinde, 2011, chapter 8.3). These shape corrections considerably complicate manipulating dr eq /dt and cause variations in dr eq /dt of less than 5%.

Evolution of Raindrop Size with Height below a Cloud
From the raindrop characteristics we outlined in section 2, we can calculate the change in raindrop equivalent radius with altitude z as where dr eq /dt is given by equation (10) and dz/dt by equation (9). Solving equation (10) also requires coupled evaluation of dT drop /dt from equation (11). Note that this formulation neglects the time to accelerate to a new terminal velocity as raindrop size changes (which is rapid compared to the timescale on which r eq evolves, as discussed in section 2.2).
To evaluate equation (15), we first must describe the atmospheric state variables that affect parameters required for calculating evaporation rate and terminal velocityp, T , RH-as functions of z. In this work, we prescribe planetary conditions of p, T , and RH, at a single z-either at the surface or at the cloud base depending on the calculation of interest. We also require planetary inputs of g and dry air composition, which are assumed to be fixed.
We then follow the standard assumptions for a 1D atmosphere in radiativeconvective equilibrium below saturated regions to relate atmospheric properties (e.g., Pierrehumbert, 2010;Romps, 2017): our pressure-temperature profile follows a dry adiabat, z is related to p (and thus RH and T ) assuming hydrostatic equilibrium, and RH is prescribed assuming the condensible gas is well-mixed (i.e., a constant molar concentration). Note that assuming constant T and RH from average values below the cloud base does not lead to significantly different dr eq /dz values, but such a simple profile does not allow for an internally consistent calculation of cloud base height.
We define cloud base as the "lifting condensation level" (LCL), the height at which a condensible gas reaches saturation in a parcel of air rising adiabatically: z such that p c (z LCL ) = p c,sat (T (z LCL )). The LCL errs in predicting cloud base when limited cloud condensation nuclei require supersaturation to initiate cloud particle formation. However, as we are concerned with where the raindrop starts evaporating (which requires RH < 100%), this caveat does not concern us. Equation (15) is stiff, so we integrate it using an implicit Runge-Kutta method of order 5. We define the raindrop's initial r eq at cloud base as r 0 . We calculate r eq (z) from the cloud base (z = z LCL ) to a desired z or until the raindrop fully evaporates. Here we define a "fully evaporated" raindrop as a drop of equivalent radius less than a threshold drop size ∆r.
We set ∆r = 1 µm; our results are not sensitive to this choice of ∆r as long as ∆r r 0 . ∆r must be non-zero for numerical stability, but ∆r > 0 also reflects the physical reality that the cloud drops that form raindrops very strongly thermodynamically favor condensing onto preexisting nuclei rather than forming homogeneously.

Minimum Cloud-edge Raindrop Size to Reach a Given Height
To understand the potential of a raindrop to transport condensible mass and heat within an atmosphere, we calculate the cloud-edge (where RH transitions to less than 1 and evaporation begins) size threshold where a raindrop can survive to a given height z without totally evaporating, r min (z). We define r min (z) as the r 0 such that r min (z)−∆r evaporates before reaching z, but r min (z) reaches the z, i.e., r eq (z, r 0 = r min ) ≥ ∆r. r min (z) is solved for via bisection by integrating r(z) as described in section 3.1 for initial radii between ∆r and the maximum raindrop radius described in section 3.4.
On terrestrial planets (with a surface at z surf ), clouds that can grow raindrops of r 0 ≥ r min (z surf ) can move condensible mass from the atmosphere to the surface condensible reservoir. We can place an lower bound on raindrop size from the cloud-edge size threshold where a raindrop can survive to the surface without totally evaporating: r min (z surf ), which we will henceforth abbreviate to simply r min . On gaseous planets, there is no surface, and raindrops can only evaporate, but their ability to transport mass and heat as a function of height is still important dynamically.

A Dimensionless Number Characterizing Raindrop Evaporation Regime
To better understand raindrop evaporation, we simplify equation (15) into a dimensionless number that can be more clearly interpreted-and evaluated-than a system of differential equations requiring numerical integration to solve. First, we need to simplify calculating T drop for an evaporating raindrop from the differential equation (11).
We assume T drop changes only as a function of altitude. This is justified by comparing the timescale on which atmospheric temperature changes (τ air ) to the timescale on which raindrop temperature changes (τ drop ). Assuming a dry adiabatic temperature profile, τ air ≈ (c p,air ∆T air )/(gdz/dt) where we conservatively set the characteristic change in air temperature ∆T air to 1 K. Assuming the atmosphere transfers heat to the raindrop via conduction, τ drop ≈ (r 2 eq ρ c, c p,c, )/(3K air ). Except for the largest possible raindrops, under broad planetary conditions τ drop τ air , and hence dT drop /dt=0 is a good approximation at a given altitude.
We define ∆T drop as the equilibrium temperature difference between the air and raindrop, i.e., From this definition of ∆T drop and equation (11), This transcendental equation can be solved numerically via a root-finding algorithm. It is commonly simplified to an analytic expression using Clausius Clapeyron, Taylor expansions, and series of assumptions regarding ∆T drop being small compared to T air (e.g., Rogers & Yau, 1996, chapter 7).
However, we find an analytic approximation that holds better across a broad range of planetary conditions is to evaluate equation (17) with ∆T drop values on the right hand side approximated as T drop must fall between T LCL and T air (i.e, ∆T drop ∈ [0, T air − T LCL ]) because there is no heat source for the drop once T drop = T air , and there is no heat sink for the drop once T drop = T LCL because RH=1 and evaporation ceases. Equation (18) can also be employed for a back-of-the-envelope calculation of ∆T drop . Now we can define a dimensionless number Λ to evaluate the tendency of a raindrop of radius r eq (and mass m) toward evaporation within a given vertical length scale . Λ is the ratio of evaporative mass loss during transit through to raindrop mass: Here we have made use of the chain rule, the definition of r eq , and the relation dm/dr eq = 4πρ c, r 2 eq . Expanding the terms in Λ gives Altitude-dependent values needed to calculate Λ are evaluated at the midpoint of . ∆T drop can be evaluated from equation (17) numerically (most accurate), from equation (17) and (18) algebraically, or from equation (18) (back-of-envelope). v T can be evaluated from equation (8) numerically or from a parameterized relationship for a commonly studied planet. Alternatively, v T can be estimated via Stokes law for very small drops or via v T,max ≈ 2 (σ c-air (ρ ,c − ρ air )g) 0.25 (ρ air ) −0.5 for very large drops (Clift et al., 2005, chapter 7

.C).
Λ values give the expected change in raindrop mass from evaporation relative to initial mass after falling a given distance. Λ(r eq , ) ≥ 1 indicates raindrops of size r eq will fully evaporate over distance . Therefore, the fraction of raindrop mass evaporated over can be estimated from min{Λ,1}. For a given , the r eq such that Λ = 1 approximates the minimum radius to reach that distance below the starting z without fully evaporating, r min (z start − ). For simplicity, here we only consider Λ defined when dz/dt <0, i.e., when a raindrop is falling downward. Though we do not treat raindrop formation here, we note that with some slight modifications this dimensionless number can also be employed to consider the effectiveness of cloud drop growth via condensation.

Maximum Raindrop Size Before Breakup
The final physical process we need to consider is raindrop breakup. Raindrops cannot grow to infinitely large sizes because the resistance provided by surface tension as surface area increases is limited. When surface tension ceases to be the dominant force experienced by a raindrop, the raindrop rapidly breaks apart.
A variety of approaches to estimating this maximum stable raindrop radius r max have been proposed previously; but none are expected to yield quantitatively exact values. The physics is additionally complicated in many situations (such as on presentday Earth) by the fact that the practical upper bound on raindrop size is not set from individual raindrop breakup but rather from hydrometeor collisions (e.g., Barros et al., 2010). Given the uncertainties, we are therefore primarily concerned here with how r max scales with external planetary properties. In particular, we focus on the effect of air density, which has inconsistently been claimed within the planetary literature to have no effect on r max (Som et al., 2012;Palumbo et al., 2020) and an extremely significant one (Craddock & Lorenz, 2017).
Variants of two methods have commonly been used to describe raindrop breakup; we review them here and describe their origins in more detail in Appendix C. First, we can estimate r max by considering when the base of a raindrop becomes unstable to small perturbations from a more dense fluid (liquid condensible) being on top of a less dense fluid (air)-generally referred to as Rayleigh-Taylor instability (Komabayasi et al., 1964;Grace et al., 1978;Lehrer, 1975;Clift et al., 2005;Pruppacher & Klett, 2010, chapter 10.3.4). This analysis yields a maximum length scale RT,max that can be related back to a maximum equivalent radius r max : There is not a definitive RT,max , with different authors choosing related, but often distinct, length scales.
Another common approach for estimating r max in the Earth literature is to calculate when the force of surface tension F σ is balanced by the aerodynamic drag force (e.g., Pruppacher & Klett, 2010, chapter 10.3.4). We henceforth refer to this approach as "force balance." Again, we get a relationship to be solved for r max that depends on a somewhat arbitrary length scale, here pertaining to surface tension σ,max : We evaluate equations (21) and (22) Figure 1: Raindrop altitude z versus equivalent radius r eq for equally log-spaced initial raindrop radii (r 0 ) near the minimum radius threshold for survival to surface (r min ). Gray-shaded lines are raindrops that evaporate before reaching the surface while purpleshaded lines are raindrops that successfully reach the surface. Planetary conditions are set to Earth-like as given in Table 1. or oblate spheroids. Both approaches yield similar expressions for r max with some variation in dependence on raindrop shape and constant factors depending on the choice of length scale.

Results
Having described the key physical processes that affect isolated falling raindrops in detail, we now present numerical results for a wide range of planetary conditions and circumstances pertaining to falling raindrops. We validated our shape and terminal velocity calculations against modern Earth observations, experimental results, and empirically based calculations (Figures S3-S4; Beard, 1976;. We also compared, with reasonable agreement, our results to previous planetary theoretical results on Titan's methane-nitrogen raindrops for shape, terminal velocity, and raindrop properties with altitude ( Figures S5-S8; Table S1; , using Cassini Huygens' probe data where appropriate (Fulchignoni et al., 2005;.

Raindrop Evaporation under Earth-like Conditions
Integrating equation (15), we investigated the behavior of raindrop evaporation for water raindrops falling from the cloud base to planetary surface under Earth-like conditions (Table 1). Figure 1 shows the evolution of raindrop radius as function of Earth-like a surface 300 altitude z for a number of an initial radii at cloud base until the raindrops either completely evaporate (r 0 < r min ) or reach the surface (r 0 ≥ r min ). The results imply a strong positive feedback on raindrop evaporation as raindrops grow smaller.  Figure 1, it shows the fraction of raindrop mass evaporated at the surface for a range of r 0 values. This curve approaches a step function about r min . Figure 2(b) extends 2(a) by showing fraction of raindrop mass evaporated via colormap versus surface relative humidity and r 0 . Qualitatively, Figure 2(b) shows the same sharp cut-off behavior as 2(a). Surface RH affects the quantitative value of r min because it varies RH(z), which impacts the magnitude of evaporation rate as well as the height of z LCL -a higher surface RH gives a lower z LCL . Both these effects act to make r min decrease as surface RH increases.
There is not yet an analytical method for estimating average surface RH in a generic planetary atmosphere-in part because of the poorly understood feedbacks of precipitation evaporation on average RH (Romps, 2014;Lutsko & Cronin, 2018)-so, for now, we consider surface RH a prescribed planetary parameter. In this plot, we vary surface RH from 99.9% to 25%-the former arbitrarily close to the threshold for evaporation to begin (RH < 100%) and the latter about the minimum surface RH for which the temperature at z LCL is above the freezing point of H 2 O given our other chosen planetary parameter values. Figure 3 is the same as Figure 2(b) except it probes the effect of vertical wind speed w rather than surface RH. Downdrafts (w < 0) increase the falling speed of raindrops while updrafts (w > 0) decrease the falling speed of raindrops (as long as w + v T < 0). Updrafts can transport raindrops upward once the updraft speed exceeds the magnitude of a raindrop's terminal velocity. As downdraft speed increases, there is a smoother transition through fraction mass evaporated across r 0 values. As updraft speed increases, fraction mass evaporated approaches a step function about r min until w + v T,max = 0. For updrafts speeds greater than this threshold, no r min exists as raindrops are no longer falling.
As with surface relative humidity, there is no analytic approach for estimating average w ranges in generic planetary atmospheres (though values can be probed by mesoscale models of sufficient resolution). In this plot, we bound updraft speed from where w exceeds maximum raindrop terminal velocity and then choose a symmetric downdraft speed bound. (This choice of a lower bound is arbitrary and does not represent an end-member case for downdraft speeds.) For simplicity, here we fix w as constant throughout raindrop falling and evaporation. In reality, vertical velocities vary both spatially and temporally within a given storm event (e.g., Lohmann et al., 2016), often as a result of the interaction between precipitation particles and ambient air (e.g., Rogers & Yau, 1996).

Evaluation of Dimensionless Number Characterizing Raindrop Evaporation Regime
Figure 4 compares using the dimensionless number Λ defined in equation (19) to predict raindrop evaporation behavior to using numerical integration for Earth-like atmospheric conditions (a)-(b) and across broad planetary conditions (c)-(d). Figure  4(a) shows calculations of the fraction of raindrop mass evaporated at the surface relative to initial mass at cloud base versus initial radius r 0 using both numerical integration and the minimum of Λ(r 0 , = z LCL ) and 1. Figure 4(c) is similar to (a) except it shows the difference in fraction mass evaporated at 500 m below the cloud base between these two methods at 4 r 0 for 90 different planetary conditions. The "broad" conditions in Table 1 give the ranges over which we vary T , p, and g at cloud base for background gas atmospheres of pure H 2 , N 2 , and CO 2 . Only one value among    Table 1). The horizontal light-purple line (w = 0) divides updrafts (w > 0) from downdrafts (w < 0). The purple line marks the calculated r min as a function of w. For r 0 < r min , raindrops totally evaporate before reaching the surface (hatched region). The dashed dark-purple line highlights the 10% mass evaporated contour within the more continuous shading.  (Table 1) evaluated from numerical integration (dashed gray line) and using dimensionless number Λ (purple line). (b) Altitude z versus threshold initial raindrop radius for total evaporation at z (r min (z)) for Earth-like atmospheric conditions evaluated from numerical integration (dashed gray line) and using dimensionless number Λ (purple line). (c) Difference in fraction mass evaporated between calculations using numerical integration and Λ versus four r 0 values evaluated across broad planetary conditions (Table 1)  T , p, and g is changed at a time relative to the "composition" conditions, which are used as a baseline.
Figure 4(b) shows calculations of the threshold minimum radius to reach altitude z without fully evaporating (r min (z)) using both numerical integration and Λ=1. Figure 4(d) is similar to (b) except it shows the relative error in r min (z) calculated using Λ instead of numerical integration at three values for the same 90 different planetary conditions as described for panel (c).
Unsurprisingly, we find that the accuracy of using Λ to calculate the fraction of raindrop mass evaporated at z decreases for r 0 near r min (z) and that the accuracy of using Λ to calculate r min (z) decreases as the distance between the cloud base and z increases. Nonetheless, Figure 4 demonstrates that Λ can capture the essential behavior of fraction mass evaporated and r min (z) with a small fraction of the computational cost of numerical integration ( 1%). Comparing the use of Λ and a full numerical integration for calculations at = 500 m, we found percent errors for r min (z) and percent differences in fraction raindrop mass evaporated were usually less than 10% in magnitude across a broad planetary parameter space, and even the largest errors were less than 20%. This agreement indicates that Λ is a viable way to predict and interpret raindrop evaporation regimes. In Appendix D, we use the definition of Λ to show the numerical results of Figures 1 and 2 can be understood from an analytic mathematical perspective.

Investigation into the Dependence of Maximum Raindrop Size on Air Density
Before we combine the concepts of minimum and maximum raindrop size thresholds, we return to the debate in the planetary literature over whether maximum stable raindrop size r max depends on air density. In Figure 5, we compare r max as a function of air density (or air pressure) using break up criteria as presented in section 3.4 and previous planetary literature. The two approaches we reviewed-Rayleigh-Taylor instability and force balance-are dependent on somewhat arbitrary length scales, so we plot the spread in values from different length scales. For Rayleigh-Taylor instability, we consider RT,max = 0.5πr eq , 0.5πa, 2r eq , 2a; for force balance we consider σ,max = 2πr eq , 2πa.
From the planetary literature on maximum raindrop size (Rossow, 1978;Som et al., 2012;Craddock & Lorenz, 2017;Palumbo et al., 2020), we plot in Figure 5 the approaches of , Craddock and Lorenz (2017), and Palumbo et al. (2020), which we will describe in more detail shortly. We do not plot Som et al. (2012)'s quantitative approach as they make use of empirical fits and Earth-based observations, but in practice they suggest a similar criterion to . Rossow (1978) and the Clift et al. (2005) expression cited in  use Rayleigh-Taylor instability criteria with length scales included in our range.
To calculate r max , , Craddock and Lorenz (2017), and Palumbo et al. (2020) all begin from the same criterion: the raindrop radius where the dimensionless Weber number equals 4. The Weber number We, which characterizes the ratio of drag force to the force of surface tension, is defined as This approach is equivalent to the force balance approach under the assumption that, within F drag , cross sectional area times C D is equal to πr 2 eq -an assumption that holds reasonably well under Earth surface conditions (Matthews & Mason, 1964).
To calculate r max ,  solves for the r eq satisfying We = 4 numerically (to account for the dependence of v on r eq ). As seen in Figure 5 : Maximum stable water raindrop size r max versus air pressure p air or air density ρ air , assuming fixed temperature T = 275 K and RH = 0 (for a linear relationship between p air and ρ air ), with Earth surface gravity and N 2 background gas. Different line colors correspond to different methods for calculating r max as labeled in the legend and described in the text. The Rayleigh-Taylor and force balance methods can use (arbitrarily) different length scales, so we plot the span of values from common length scale choices rather than single lines.
yields an r max that varies by about a factor of 1.5 over the range of ρ air we consider. We view this variation with density as a non-physical dependence introduced by the neglect of the dimensionless C D in representing the drag force in the formulation of We (Kolev, 2007, chapter 8). (C D nonlinearly depends on ρ air through Re.) Varying the C D parameterization within this calculation causes comparable changes in r max to varying ρ air ; thus, we consider the  method consistent with no significant dependence of r max on ρ air . To calculate r max , Palumbo et al. (2020) solves for the r eq satisfying We = 4 algebraically after assuming spherical raindrops and C D =1: Craddock and Lorenz (2017) does not present a simplified expression for calculating r max , only evaluations under different atmospheric conditions and a statement that "larger diameter raindrops are not possible at higher atmospheric pressures." When we simplify their presented equations involved in describing r max (their equations (1) and (3)) following their stated assumptions, we arrive at the same r max result as Equation (24). We are only able to reproduce the results of Craddock and Lorenz (2017)'s evaluations of r max (their Table 1) using an expression for r max inconsistent with length units: In addition to theoretical methods of estimating r max , we also plot in Figure 5 a range of claimed maximum raindrop sizes on present-day Earth, both from experiments and natural observations. This range is consistent in magnitude with all the estimates. We give a range of maximum values as the measurement is an attempt to estimate the end of the extreme tail of a stochastic process (e.g., Komabayasi et al., 1964;Grace et al., 1978;Clift et al., 2005). Single-value maxima fall between r eq of 4-5 mm (Merrington & Richardson, 1947;Beard & Pruppacher, 1969;Ryan, 1976;Hobbs & Rangno, 2004;Gatlin et al., 2015). Gatlin et al. (2015) compiled observations of 2.4 × 10 8 raindrops and found 0.4% of these raindrops had r eq ≥ 2.5 mm and only 1.9 × 10 −5 % had r eq ≥ 4 mm-statistics that suggest in practice r max is about 2.5-4 mm.
Returning to theory, in Figure 5 we see that while different assumptions about raindrop shapes in the Rayleigh-Taylor and force balance methods lead to factor of a few quantitative differences in r max , these differences are not sensitive to air density. We conclude that the effects of raindrop shape are ultimately of limited importance in estimating r max due to the ambiguity of raindrop length scales in the calculation setups. Thus the scalings of r max with planetary variables can reliably be seen analytically by assuming spherical raindrops (i.e., setting b/a = 1): The force balance, Rayleigh-Taylor instability, and Weber number methods all yield the same approximate scalings, which are effectively independent of air density. (These scalings are not novel (e.g., Clift et al., 2005); we simply present them in the context of this planetary debate on r max .) Therefore, we agree with Palumbo et al. (2020) that Craddock and Lorenz (2017)'s finding that larger raindrops become possible as time advances and Mars experiences atmospheric escape due the dependence of raindrop stability on ρ air is not justified. This conclusion is also consistent with extensive modern-Earth based literature considering r max , which does not explicitly highlight any dependence of r max on air density, which non-trivially varies from cloud to surface. Finally, we note that we  Table 1. Dark-gray lines give r eq at the onset of Rayleigh-Taylor instability (r max ); purple lines give the minimum cloud-edge raindrop radius required for the raindrop to not evaporate before reaching the surface (r min ) for three surface RH values. Purple shaded regions show cloud-edge size bounds of raindrops that can reach the surface. The gray shaded region sketches cloud condensation nuclei sizes.
are not claiming that average raindrop size is insensitive to air density, only that the instability of individual large raindrops does not have a significant dependence on air density.

Raindrop Size Bounds for Terrestrial Planets
In Figure 6, we calculated r max and r min for water raindrops across a range of planetary conditions. For clarity, we plot only one of the r max values from the methods we discussed previously (Rayleigh-Taylor with RT,max = 0.5πa) and r min values for only three representative surface relatives humidities (RH = 0.25, 0.5, 0.75) and no vertical wind (w = 0). We used as default values the Earth-like conditions of Table 1 when the x-axis planetary parameter is not varying.
The purple shading highlights how our estimates of r max and r min constrain possible raindrop sizes that can transport condensible mass from a cloud to the surface reservoir. Varying colors of purple shading correspond with size bounds set from varying surface relative humidity values. As the purple shading lightens, the bound includes smaller r min values associated with higher surface relative humidity values. The dark purple shading (for RH surf =0.25) is the strictest bound. For perspective, cloud drops grow from cloud condensation nuclei with many orders of magnitude smaller sizes as we have schematically indicated with gray shading. On Earth, typical cloud condensation nuclei are around 0.05 µm (5 × 10 −5 mm) (Lohmann et al., 2016) while, for conditions considered here, viable raindrop sizes vary by about an order of magnitude with typical values of tenths of millimeters, or about 10,000 times larger than typical cloud condensation nuclei.
In Figure 6(a), we plot raindrop size bounds as described for variable dry surface pressures. Neither of the size bounds has a strong dependence on pressure (see also section 3.4). r min depends on pressure in multiple ways that largely cancel each other out. Figure 6(b) shows the impact of surface gravity g on raindrop size bounds. As g increases, r max and r min all systematically decrease like g −0.5 . Larger raindrops are possible at lower surface gravities, and raindrops must also be larger to survive to the surface without evaporating. Figure 6(c) highlights the effect of increased evaporation rate in higher air temperatures. r min rises with T surf because evaporation rate and falling time to surface increase as long as the molar mass of the condensible gas is less than the average dry air molar mass-as considered here. Table 2 considers the effect of atmospheric composition on the time and distance from cloud base until evaporation (t evap and z evap , respectively) for H 2 O raindrops with pure H 2 , He, N 2 , O 2 , and CO 2 atmospheres. Other planetary conditions used are given under "composition" in Table 1. Atmospheric composition impacts raindrop evaporation from three main effects: (1) molar mass and heat capacity impact atmospheric structure, which governs how vertical distance maps to temperature, pressure, and relative humidity-all key parameters in calculating dr/dz; (2) molar mass impacts how a given pressure maps to a density, which impacts raindrop terminal velocity; and (3) molar mass and molecular structure impact the rate at which air can transport latent heat and condensible gas away from the raindrop.

Raindrop Evaporation with Varying Atmospheric Composition
In Table 2, we calculated t evap and z evap considering each of these effects of composition in isolation as well as all together. When only a single effect is considered, all other compositional effects are calculated with pure N 2 . Total number of molecules at cloud base is held fixed (i.e., pressure is fixed as an ideal gas assumed).
As shown Table 2, air composition acts on the integral of dr eq /dz in competing directions, so all effects must be considered in unison to understand how precipitation evaporation will vary with atmospheric composition. We find the time taken to evaporate is comparable across atmospheric conditions. The distance to evaporation is comparable for the higher molecular mass gases with the He atmosphere about 1.75 times larger and the H 2 atmosphere about 3.5 times larger. Excluding the noble gas He, the ability to transport condensible mass in units of scale heights increases as atmospheric molar mass increases.
We further probe the variations in raindrop evaporation due to atmospheric composition in Figure 7 by calculating latent heat absorbed per second (power, P evap ) as a function of (a) atmospheric pressure, (b) vertical distance from cloud base, and (c) falling time for a fixed initial condensible mass sorted into raindrops of three different radii. P evap is linearly related to the rate of condensible mass evaporation through the latent heat of vaporization, which is normalized per unit mass. Thus we choose a constant initial mass (the mass of the largest raindrop considered) to compare magnitudes of P evap for different initial raindrop sizes. We use the same planetary conditions as Table 2 (under "composition" in Table 1) with background H 2 , N 2 , and CO 2 atmospheres.
In all three vertical coordinates shown in Figure 7, the peak P evap reached is about the same for a given r 0 across all three atmospheric compositions. Maximum P evap decreases as r 0 increases (as expected in order to conserve mass with longer fall times). Figure 7 shows P evap is roughly a quadratic function of log p (a), z LCL − z (b), and t (c) between cloud base and reaching total evaporation for all atmospheric compositions and initial size values considered.
For a given composition atmosphere, increasing r 0 increases p evap , z evap , and t evap . The spread in p evap and z evap at total evaporation between different compositions increases as r 0 increases. As seen in Table 2 (given log p is essentially proportional to altitude over scale height), p evap increases with increasing dry molar mass; z evap decreases with increasing dry molar mass; and t evap is about constant across the different compositions. Note. Time to evaporate t evap and falling distance from cloud base before evaporation z evap in meters and relative to atmospheric scale height H for a water raindrop of initial size r 0 = 0.5 mm in different composition atmospheres. Planetary conditions besides dry composition are given in Table 1 under "composition." With the effect column, we consider the three main impacts of composition on t evap and z evap -atmospheric scale height H, raindrop terminal velocity v T , and transport rate of condensible gas molecules and heat away from the raindrop surface "transport"-together ("all") as well as in isolation.  scale heights versus the minimum initial radius raindrop to reach altitude z without total evaporation (r min (z)). Different colored solid lines represent ostensible atmospheric conditions (given by Table 1) for different planets as labeled. All raindrops are composed of water. Thin dashed vertical lines indicate maximum stable raindrop radius r max as estimated via Rayleigh-Taylor instability with RT,max = 0.5πa, distinguished by planet with labeled colors. r min (z) values are plotted until r min (z) = r max or z intersects with the planet's surface.

Comparison of Water Raindrops on Specific Planets
Next we move from considering raindrops in abstract conditions to studying raindrops in known (or speculated) planetary conditions. In Figure 8, we compare water raindrops on Earth; warm, wet ancient Mars; Jupiter; Saturn; and exoplanet K2-18b. Warm, wet ancient Mars is a hypothesized climate state 3-4 billion years ago where Mars was warm compared to the melting point of water and rainfall was frequent (e.g., Wordsworth, 2016).
K2-18b is an exoplanet without analog in the solar system, falling between the sizes of Earth and Neptune and receiving an Earth-like amount of stellar insolation Montet et al., 2015;Cloutier et al., 2019;Benneke et al., 2019). Though many of K2-18b's characteristics are imprecisely known, we selected this exoplanet to demonstrate the flexibility of our model as multiple teams have claimed observational detections of water vapor (Tsiaras et al., 2019;Benneke et al., 2019;Madhusudhan et al., 2020) and one team has hypothesized that observations suggest the presence of liquid water clouds (Benneke et al., 2019).
We set surface conditions on the terrestrial Earth and Mars and cloud base conditions on the gaseous Jupiter, Saturn, and K2-18b as given in Table 1. We plot altitude z from cloud base in units of (a) meters and (b) atmospheric scale heights versus r min (z), the minimum radius raindrop to reach z without totally evaporating. How r min (z) values cluster among planets varies depending on whether the distance from cloud base is measured in scale heights or meters, but for all z values tested in both unit systems we find r min (z) varies among the diverse planetary conditions considered by only a factor of few. This agreement has interesting implications for future work aimed at rigorous generalization of raindrop microphysics schemes from Earth to other planets.

Raindrops beyond H 2 O
While water is the most familiar liquid condensible, outside of the familiar Earth temperature range, many other species condense as liquids and can form raindrops. Table 3 compiles a number of liquid condensible species that are cosmochemically abundant. Detailed analysis of raindrop behavior requires specifying atmospheric conditions in addition to condensible species, so here we consider how basic condensible properties vary raindrop behavior relative to water. Again in the interest of simplicity, we do not consider mixtures of condensibles predicted by thermodynamic equilibrium in many atmospheric gas combinations, e.g., N 2 -CH 4 or NH 3 -H 2 O ( Guillot, Stevenson, et al., 2020).
We give the melting temperature T melt for each condensible to give an idea of the atmospheric temperatures where each species will be liquid. Lower temperature condensibles like CH 4 and NH 3 dominate the observable clouds of the outer solar system. The CH 4 cycle of Titan is the only active "terrestrial" condensible cycle besides Earth's we can observe in detail.
Metal and rock species like Fe and SiO 2 become condensible species at very high temperatures. Such species are predicted to be condensibles on highly irradiated exoplanets that are favored observational targets for the foreseeable future; Ehrenreich et al. (2020) have already claimed observational evidence of Fe condensing on WASP-76b. On Earth, such high temperatures can be reached during asteroid/meteoroid impacts. Geologically preserved impact spherules (e.g., Johnson & Melosh, 2012a, 2012b and micrometeorites (e.g., Tomkins et al., 2016;Payne et al., 2020) both undergo phases in their life through the atmosphere where their behaviors are described by the raindrop physics we have presented. Note. Temperature-dependent values use T = T melt . a Linstrom and Mallard (2014). b Somayajulu (1988). c Rumble et al. (2017). d Vargaftik et al. (1983). e Assael et al. (2006). f Brillo and Egry (2005). g Desai (1986). h Melosh (2007). i Bacon et al. (1960). j Kingery (1959).
Despite the very wide range of condensible species, r max values in Table 3 only vary relative to water by a factor of 0.5 to 2. As we have examined in section 3.4, the only planetary parameter beyond condensible type that affects r max is surface gravity g, which for planetary bodies varies about an order of magnitude (between about 1 and 25 m s −1 ). We therefore find that maximum stable raindrop sizes are remarkably similar across a very wide range of planetary conditions and raindrop compositions.

Discussion
This work is merely a first step toward a generalized theory of how precipitation and condensible cycles operate in planetary conditions different from modern Earth. We have considered only single raindrops, independent of their formations. Future progress will require development of theory for general planetary atmospheres on the growth of raindrops from cloud drops and extensions that include solid precipitating particles and their growth. In this context, below we discuss some future applications and extensions of this work.

Precipitation Efficiency
Precipitation efficiency measures how efficiently an atmosphere transports condensed mass from a cloud downward. Qualitatively, its distribution over time is an important metric for planetary climate as it shapes cloud coverage (both temporally and spatially), cloud radiative properties, and relative humidity profiles, which all have large consequences for radiative balance (Romps, 2014;Zhao et al., 2016;Lutsko & Cronin, 2018). On a terrestrial planet, precipitation efficiency evaluated at the surface is a particularly important quantity as it helps to set the amount of condensible mass in the atmosphere.
On short timescales, liquid precipitation efficiency is governed by the basic raindrop physics of falling and evaporation we have presented here. We have demonstrated that the dimensionless number Λ captures the essential behavior of raindrop evaporation and descent. For Λ evaluated with the length-scale from cloud base to the surface, raindrops of initial sizes with Λ > 1 (i.e., r 0 < r min ) will totally evaporate and have no condensible mass reach the surface; raindrops of initial sizes with Λ < 0.1 will experience little evaporation and have the majority of their mass reach the surface; and raindrops of initial sizes with 0.1 < Λ ≤ 1 will have both condensible mass evaporate and reach the surface. Short-scale precipitation efficiency is then fundamentally controlled by the cloud-edge condensed mass distribution among these three different size categories.
How precipitation efficiency on short timescales-governed directly by microphysicsmaps to the climatically important temporal distribution of precipitation efficiency is fundamentally influenced by both large-scale atmospheric dynamics and local-scale convection (Romps, 2014). Predicting precipitation efficiency distributions in threedimensional models requires capturing key microphysical behaviors in a parameterized representation. Our analysis suggests the key microphysical behavior we must capture to physically represent precipitation efficiency on short timescales is the sorting of condensed cloud mass into the three different size categories determined by Λ. This interpretation suggests a physical grounding for microphysics' role in controlling an important yet poorly understood climate parameter (Zhao et al., 2016;Lutsko & Cronin, 2018). Tying precipitation efficiency to the raindrop size distribution may therefore provide a framework for future improvements in generalized microphysics parameterizations.
Given this analysis, the bounds on water raindrop sizes in Figure 6 show one component of how precipitation efficiency will be shaped by different planetary conditions. A complete picture of precipitation efficiency obviously will also require better understanding of how planetary parameters influence the formation and subsequent mass distribution of raindrops. Nevertheless, the shrinking viable surface-reaching raindrop size range with rising T surf in Figure 6(c) is a striking predicted feature of warmer water cycles, which may have implications for a CO 2 -rich early Earth, for early Venus, or for exoplanets close to the runaway greenhouse threshold.

Convective Storm Dynamics
Evaporating raindrops also influence convective storm dynamics. The vertical transport of heat and condensible mass from evaporating raindrops causes local variations in air density through changes in both average molar mass and temperature. The implications of evaporating raindrops for convection depend on the ratio between an atmosphere's dry mean molecular mass µ dry and the molecular mass of its condensing species µ c as well as the removal rate of latent heat relative to the local T − p profile (Guillot, 1995;Leconte et al., 2017).
On modern-Earth, µ c /µ dry ≈ 0.6, so molar-mass-contrast effects are present but muted. Nonetheless, the interplay of rising air supplying condensible mass and sinking precipitating air plays a key role in the evolution and lifetime of a given storm (e.g., Rogers & Yau, 1996). Variations in µ c /µ dry have been hypothesized to drive storm systems completely unlike those on Earth-e.g., Saturn's giant white storms (Li & Ingersoll, 2015;Leconte et al., 2017).
As suggested by Table 2 and Figure 7, an atmosphere's background dry gas also influences the ability of raindrops to vertically re-distribute latent heat and condensible mass with respect to z and log p. This effect of background gas composition will also likely influence air density changes during storms. Future dynamical studies might explore the implications of these variations in raindrop evaporation with background gas for storm evolution and subsequent condensible gas distributions. One example solar system application of interest is better constraining how deep and how effectively ammonia raindrops (originating from melted NH 3 snowflakes or hail-like NH 3 -H 2 O "mushballs") can transport ammonia on Jupiter (Ingersoll et al., 2017;Li & Chen, 2019;Li et al., 2020;Guillot, Stevenson, et al., 2020;Guillot, Li, et al., 2020).

Rainfall Rates
On terrestrial planets, rainfall rate-the mass flux of liquid condensible hitting the surface-is a key characteristic for predicting surface erosion and flooding events (e.g., Kavanagh & Goldblatt, 2015;Craddock & Lorenz, 2017;Margulis, 2017). Rainfall rate depends on liquid condensible mass per unit air, the raindrop size distribution, and raindrop velocities as a function of size (Rogers & Yau, 1996). While we cannot yet make robust predictions for how rainfall rates should vary in different planetary conditions, we expect such a relationship will be sensitive to g and air density at cloud level because of the role of collisional kinetic energy in shaping raindrop size distributions (Low & List, 1982;Rogers & Yau, 1996;Pinsky et al., 2001;List et al., 2009). r min and Λ as a function of r eq evaluated at the surface are useful for contextualizing the size range of the in-cloud raindrop size distribution that is important for predicting rainfall rates. Future work should investigate how narrow bounds of viable surfacereaching raindrops might constrain rainfall rates across different planetary conditions in more detail.

Extension to Solid Particles
In this work, we have focused on falling liquid condensible particles because solid particle shapes are much more complicated and variable than liquid oblate spheroids. The morphology of solid condensed particles exhibits extreme variability because of the high sensitivity of crystal orientation to temperature and condensible vapor supersaturations (e.g., Libbrecht, 2017). Further, crystal structure is fixed at deposition, so because a solid condensed particle experiences variations in environmental conditions during its growth, its final shape is highly sensitive to its growth path. These considerations mean generalizing modern-Earth modeling approaches for handling ice shape degeneracies (e.g., Krueger et al., 1995) is a non-trivial exercise.
Understanding of shape is the main limiting factor for applying the methodology presented here for raindrops to solid condensible particles. Shape is a key parameter in calculating C D and f V , which impact terminal velocity and evaporation rate, respectively. Shape also needs to be accounted for more fundamentally in the derivation of evaporation rate because it controls boundary conditions, but analogous mathematics to well-investigated electrostatics means such boundary condition accounting has already been compiled (e.g., McDonald, 1963;Lamb & Verlinde, 2011, chapter 8). Extension to solid particles would also need to account for differences in saturation pressure with respect to solid and liquid condensibles and a larger assortment of latent heats because of more available phase changes, but these issues are a question of compiling thermodynamic data rather than an inherent lack of understanding.

Conclusion
We have compiled and generalized methods for calculating raindrop shape, terminal velocity, and evaporation rate in any planetary atmosphere. These properties govern raindrop behavior below a cloud, a simple but necessary component of understanding how condensible cycles operate across a wide range of planetary parameters beyond modern-Earth conditions. For terrestrial planets, raindrops sizes capable of transporting condensed mass to the surface only span about an order of magnitude, a narrow bound when compared to origins in cloud condensation nuclei many orders of magnitude smaller. We show across a wide range of condensibles and planetary parameters that maximum raindrop sizes do not significantly vary.
With more in depth calculations, we confirm the conclusion of Palumbo et al. (2020) that maximum raindrop size is only weakly dependent on air density, in contrast to the results of an earlier study (Craddock & Lorenz, 2017). By returning to the physics equations governing raindrop falling and evaporation, we demonstrate raindrop ability to vertically transport latent heat and condensible mass can be well captured by a new dimensionless number. Our analysis suggests cloud-edge, mass-weighted raindrop size distribution is a key microphysics-based control on the important climate parameter of precipitation efficiency.

Acknowledgments
Our model and all code used to generate results for this paper are available at an archived Github repository (Loftus, 2021, https://github.com/kaitlyn-loftus/ rainprops). This work was supported by NSF grant AST-1847120. K.L. thanks Miklós Szakáll for providing data on raindrop shape experiments as well as Mark Baum, Ralph Lorenz, and Jacob Seeley for helpful discussions pertaining to numerical methods, Martian raindrops, and climatic impacts of microphysics parameterizations, respectively. We thank Tristan Guillot and Ralph Lorenz for insightful reviews. ] density X air describes air property at local altitude X c describes a condensible property X drop describes a raindrop property X describes a liquid condensible property where ambiguous if "c" references liquid or gas condensible X evap describes property when raindrop finishes evaporating X LCL describes air property at lifting condensation level X surf describes air property at surface

Appendix A Oblate Spheroid Geometrical Relationships
Here, we give geometric properties of our assumed oblate-spheroid raindrop with semi-major axis a, semi-minor axis b, and equivalent radius r eq . Note that a sphere is an oblate spheroid of b/a = 1.
Raindrop semi-major axis a can be determined numerically for a given r eq and b/a from the relationship a = r eq b a −1/3 (A1) ). An oblate spheroid has cross sectional area A of (A2)  and a volume V of by definition of r eq . The ratio of the surface area of an oblate spheroid to that of a sphere f SA is (Loth, 2008).

Appendix B Raindrop Shape
The shape of a falling raindrop in air (i.e., the surface of the condensible-air boundary) can be described by the Young-Laplace equation, which governs the surface boundary between two immiscible (non-mixing) fluids. For a given point on a raindrop surface, the Young-Laplace equation can be written as Pruppacher & Klett, 2010, chapter 10.3.2 for derivation). σ c-air is surface tension between liquid condensible and air; R 1 and R 2 are the principle radii of curvature, which together describe the shape of the local surface via its curvature (De Gennes et al., 2013); and ∆p is the difference between internal and external pressures on either side of the raindrop's surface boundary.
Because we have pre-assumed an oblate spheroid geometry, we only have to evaluate equation (B1) at one point on the raindrop's surface rather than integrating over the entire surface, and the principle radii are analytic expressions rather than values which must be iteratively corrected for such integration to be self consistent (see , for shape calculation without oblate spheroid assumption). We follow  and solve equation (B1) at a point on the raindrop's equator. From geometry, for all points on an oblate spheroid's equator, and (see , Appendix B for derivation).
In a full accounting of pressures, internal pressure p int should include hydrostatic pressure, pressure from spherical surface tension, and pressure from internal circulation within the raindrop. External pressure p ext should include hydrostatic pressure, pressure from aerodynamic drag, and pressure from air turbulence. Following , we only consider internal and external hydrostatic pressures and pressure from spherical surface tension, giving ∆p = p int − p ext ≈ (ρ c − ρ air )gb + 2σ c-air r −1 eq (B4) Pressure from drag is fundamental to shaping the raindrop's shape at larger raindrop sizes (i.e., where the magnitude of F drag begins to approach the magnitude of F σ ); drag is responsible for the evolution in raindrop shape from oblate spheroid to upper hamburger bun. Drag acts to deform the raindrop so that it is no longer axially symmetric about the major axis because drag is not uniformly distributed over the raindrop's surface.
Our assumption of oblate spheroid shape is incompatible with a detailed consideration of the effects of drag on raindrop shape. We fully neglect drag because we consider the pressure balance at raindrop equator where drag pressure is minimal even when the total magnitude of the drag force on the raindrop is significant . The consequences of this neglect can be seen to be minimal from comparison of 's oblate spheroid method to computed or observed raindrops shapes where drag is included .
We neglect the effects of atmospheric turbulence on raindrop shape as we are calculating an equilibrium shape; in practice, turbulence acts to induce oscillations in raindrop shape (which are then viscously damped) rather than changing the equilibrium shape (Beard et al., 2010). We also neglect the effects of internal circulation within the raindrop as empirically its neglect has no impact on predicting Earth raindrop shapes  and theoretically internal circulation is expected to be small for a higher dynamic viscosity liquid raindrop falling through a lower dynamic viscosity air (Clift et al., 2005).

C1 Rayleigh-Taylor instability
Raindrop breakup can be studied by a linear instability analysis that incorporates capillary and gravity waves. For wavelengths above a critical wavelength λ * , total wave phase velocity on the base of the drop becomes imaginary; waves rapidly amplify in magnitude; and the raindrop becomes unstable. Assuming planar surface waves (the 3D corrections for drops near the size of r max are small (Dhir & Lienhard, 1973;Grace et al., 1978)), the critical wavelength can be written as (Komabayasi et al., 1964;Grace et al., 1978).
Physically, for a raindrop to not be disrupted by a wave, the zero mode of wavelength λ * /2 must be less than the characteristic size of the drop: Different metrics exist for the maximum raindrop length scale relative to the zero mode wave RT,max -e.g., Grace et al. (1978) uses half the equivalent radius' circumference, 0.5πr eq , while Pruppacher and Klett (2010) uses the maximum physical raindrop diameter, 2a. Different choices of length scale will result in quantitatively different results for the onset of instability by a factor of a few; this discrepancy emphasizes the role of this analysis as an estimate of when drops tend to become unstable. Regardless, any choice of length scale results in the same dependencies on physical parameters-our primary concern here. Length scales using semi-major axis a can be related to r eq via oblate spheroid geometry.
While no experimental pressure chamber data exists to explicitly test air density's effect on maximum raindrop size, numerous chemical engineering experiments have been done on the maximum-sized drops of different media falling through various other media. Such experiments consistently agree (within about 20%) with the predictions of equation (21) where "c" is replaced with the drop medium, "air" is replaced with the fall medium, and the dynamic viscosity of the drop medium is much larger than the fall medium's (like raindrops and air) (Lehrer, 1975;Grace et al., 1978;Clift et al., 2005). More complex wave analysis does not consistently produce better agreement with experiments (Grace et al., 1978).

C2 Surface Tension-Drag Force Balance
The force on a raindrop due to surface tension is where σ is a characteristic length scale of the raindrop's surface, which is conventionally taken to be 2πa or 2πr eq . Again, there is ambiguity in length scales because this setup is an estimate rather than a rigorous calculation. Drag force is given by equation (4). We note the maximum velocity relative to air of a raindrop falling in isolation is its terminal velocity (where F drag = F g ), and raindrop velocity is a monotonically increasing function with raindrop size. Thus, at the smallest r eq where F drag = F σ , the setup is analogous to solving for the raindrop size where the force of surface tension equals the gravitational force, given by equation (3).

Appendix D Analysis with Λ of the Dependence of Raindrop Evaporation on Size
We can understand the behavior of raindrop evaporation with changing raindrop size by considering the dependence of Λ on r eq . From equation (20), we can approximately represent Λ explicitly in terms of r eq as Λ ≈ Cr −(2+β) eq . (D1) C is a proportionality constant dependent on the environment across and the species of condensible; β is defined such that dz/dt ∝ r β eq . (We have neglected the dependence of ventilation factors on r eq to get a tractable expression here.) For simplicity, we consider the limiting case where w = 0, so dz/dt = v T . For very small raindrops (Re 1), the raindrop is in the Stokes regime and β = 2. For very large raindrops (as r eq approaches r max ), raindrop terminal velocity approaches a constant value (Clift et al., 2005, chapter 7.C), and β approaches 0. From the behavior of C D , β smoothly varies so that β ∈ [0, 2] (e.g., Lohmann et al., 2016, chapter 7.2.3). Therefore, as r eq decreases, Λ always exponentially increases, with the dependence on r eq moving from Λ ∝ r −2 eq to Λ ∝ r −4 eq . The exponential dependence of Λ on r eq means that the transition from minimal evaporation to full evaporation occurs over a narrow raindrop size range.
One perspective on the width of this transition regime-applicable for raindrop sizes varying by orders of magnitude-is to consider the ratio of the r eq such that Λ = 0.1 to the r eq such that Λ = 1 (i.e., the ratio of sizes between the drop that evaporates 10% of its initial mass and the drop that fully evaporates just as it reaches the end of the prescribed length scale). From equation (D1) and assuming β is constant between Λ = 1 and Λ = 0.1, this ratio is equal to 10 1/(2+β) , which evaluates to 1.78 and 3.16 for β = 2 and β = 0, respectively. Because β decreases with increasing r eq , as the r eq that evaluates to Λ = 1 increases (e.g., by increasing ), the transition width increases. Still, regardless of the exact value of β within [0,2], a change in raindrop size of a factor of 2-3 is small compared raindrop growth processes that require many order of magnitude changes in drop size. Similar analysis featuring Λ can be used to introduce the additional complexity of non-zero vertical wind speed or to isolate other variables of interest beyond r eq and consider their effects on raindrop evaporation.

Contents of this file
1. Text S1 2. Figures S1 to S8 3.  Table S1 compare our model to previous theoretical work focused on Titan raindrops with further description given in Text S1. Text S1. To compare our methods to previous theoretical work on Titan raindrops by  and , we present results (1) as reported by the previous authors, (2) from equivalent calculations done by attempting to reproduce their methods within the basic framework of our model (labeled as "reproduction"), and (3) from equivalent calculations done using the methodology we describe in this work. We note that the scatter points labeled  in Figures S7 and S8 were digitized by us from their figures and are not directly reported values. (We know there is at least one inconsistency in what is plotted in Figures S7 as our digitized starting radius is r0 = 4.80 mm rather than the reported r0 = 4.75 mm.) Where possible, we follow condensible and atmosphere properties as described or referenced in each text for the most direct comparisons; we review the most notable ones here. Under Titan's atmospheric conditions, N2 gas is expected to dissolve into liquid CH4 at significant molar concentrations (e.g., .  accounts for this mixture by prescribing a liquid raindrop density between that expected for pure CH4 and pure N2.  accounts for this mixture more rigorously: accounting for the CH4-N2 mixture thermodynamics outlined by  and evaporation of both CH4 and N2 using Raoult's law.  uses pressure, altitude, temperature, and composition data from Cassini's Huygens probe (Fulchignoni et al., 2005;. We linearly interpolate this data to be in terms of altitude to use in our numerical integration.  Figure S1. (a) Vertical distance from cloud base for raindrop to reach 99% terminal velocity (z 99% ) in units of fall distance to total evaporation (zevap) versus time for raindrop to reach 99% terminal velocity (t 99% ) in units of fall time to total evaporation (tevap) across a broad planetary parameter space. (b) Same as (a) except z 99% given in units of atmospheric scale heights and t 99% in units of seconds. For H2, N2, and CO2 background composition atmospheres (marker shape as labeled), we vary in turn TLCL, pLCL, and g (color as labeled) over ranges given in Table 1 under "broad" from baseline conditions of TLCL = 280 K, pLCL = 10 4 Pa, and Earth surface gravity. We test 90 different atmospheric conditions with H2O raindrops in all cases. Raindrop acceleration is integrated from drag and gravitational forces following Newton's second law. For maximum time to reach terminal velocity, we perform this calculation using r0 = rmax and un-physically begin integration at cloud base from dz/dt=0 m s −1 . a Properties of a CH4-N2 raindrop of initial size r 0 = 4.75 mm at 8 km after falling to Titan's surface (z = 0 km). t fall is the total fall time, and f X is the liquid molar concentration of species X in the raindrop. In the first column, we give the values reported by . In the second column, we show the results our own attempt to reproduce the methods of . In the third column, we show the results of our methods for the same calculation. Figures S7 and S8 show altitude-dependent properties of the same raindrop.
Copyright 2021 by the American Geophysical Union. 0148-0227/21/$9.00  Figure S2. (a) Median relative error in raindrop velocity v from assuming terminal velocity is instantly reached during raindrop evaporation (relative to self-consistent treatment of raindrop acceleration) versus equivalent radius at cloud base r0 across a broad planetary parameter space. (b) Same as (a) except for maximum relative error in v. For H2, N2, and CO2 background composition atmospheres (marker shape as labeled), we vary in turn TLCL, pLCL, and g (color as labeled) over ranges given in Table 1 under "broad" from baseline conditions of TLCL = 280 K, pLCL = 10 4 Pa, and Earth surface gravity. We test 90 different atmospheric conditions with H2O raindrops in all cases. Raindrop acceleration is integrated from drag and gravitational forces following Newton's second law.  Figure S3. Raindrop axis ratio b/a versus equivalent radius req for H2O raindrops and Earth surface gravity. The dark purple line shows our theoretical calculation from equation (2) (following . Dark gray points show results from more comprehensive theoretical models . Light purple points show experimental results (Beard et al., 1991;). The dashed light-gray line shows an empirical fit to additional experimental results .  Figure S4. Raindrop terminal velocity vT versus equivalent radius req for H2O raindrops. Atmospheric parameters are set to mimic Earth experimental conditions: T = 293.15 K, p = 1.01325 × 10 5 Pa, RH=0.5, Earth surface gravity, and dry air composition of 20% O2 and 80% N2. The dark purple line shows our theoretical calculation from equation (8). The light-purple scatter points represent the experimental data of  and . The dashed light-purple line shows our evaluation of the empirical relationship of Beard (1976  . Raindrop axis ratio b/a versus equivalent radius req for CH4-N2 raindrops and Titan surface gravity. The dark-purple line shows our theoretical calculation from equation (2) (following . The light-purple line shows our reproduction of the theoretical method of . (We do not give scatter points of values directly from     reproduction Earth (H 2 O), 0 km Titan (CH 4 -N 2 ), 10 km Titan (CH 4 -N 2 ), 0 km Lorenz (1993) Figure S6. Raindrop terminal velocity vT versus equivalent radius req for Earth and Titan conditions (line styles as labeled). The dark purple line shows our theoretical calculation from equation (8). The light-purple line shows our reproduction of the theoretical method of . The values of light-purple scatter points are taken from Table 1 of   this work  reproduction   Figure S7. Altitude z versus equivalent radius req for Titan conditions for a CH4-N2 raindrop of initial radius r0=4.75 mm at 8 km. The dark-purple line is calculated following the methods outlined in sections 2-3. The lightpurple line is calculated via our reproduction of the methods of . The values of light-purple scatter points are taken from Figure 2(b) of . From sensitivity tests, the majority of the difference in the calculated r(z) between our method and that of  is due to our choice to not assume the ventilation factor for heat transport is equal to the ventilation factor for molecular transport. Thus the difference is not a fundamental disagreement of models but rather the result of uncertain parameters quantifiable by fluid dynamics experiments.