Turbulence, raindrops and the l1/2 number density law

Using a unique data set of three-dimensional drop positions and masses (the HYDROP experiment), we show that the distribution of liquid water in rain displays a sharp transition between large scales which follow a passive scalar-like Corrsin–Obukhov (k-5/3) spectrum and a small-scale statistically homogeneous white noise regime. We argue that the transition scale lc is the critical scale where the mean Stokes number (= drop inertial time/turbulent eddy time) Stl is unity. For five storms, we found lc in the range 45–75 cm with the corresponding dissipation scale Stη in the range 200–300. Since the mean interdrop distance was significantly smaller (≈ 10 cm) than lc we infer that rain consists of ‘patches’ whose mean liquid water content is determined by turbulence with each patch being statistically homogeneous. For l>lc, we have Stl<1 and due to the observed statistical homogeneity for l<lc, we argue that we can use Maxey's relations between drop and wind velocities at coarse grained resolution lc. From this, we derive equations for the number and mass densities (n and ρ) and their variance fluxes (ψ and χ). By showing that χ is dissipated at small scales (with lρ, diss≈lc) and ψ over a wide range, we conclude that ρ should indeed follow Corrsin–Obukhov k-5/3 spectra but that n should instead follow a k-2 spectrum corresponding to fluctuations scaling as Δρ∝l1/3 and Δn∝l1/2. While the Corrsin–Obukhov law has never been observed in rain before, its discovery is perhaps not surprising; in contrast the Δn≈l1/2 number density law is quite new. The key difference between the Δρ, Δn laws is the fact that the microphysics (coalescence, breakup) conserves drop mass, but not numbers of particles. This implies that the timescale for the transfer of the density variance flux χ is determined by the strongly scale-dependent turbulent velocity whereas the timescale for the transfer of the number variance flux ψ is determined by the weakly scale-dependent drop coalescence speed. We argue that the l1/2 law may also hold (although in a slightly different form) for cloud drops. Because they are consequences of symmetries, we expect the l1/3, l1/2 laws to be robust. Since the large-scale turbulence determines the n and ρ fields which are the 0th and 1st moments of the drop-size distribution, they constrain the microphysics: dimensional analysis shows that the cumulative probability distribution of nondimensional drop mass should be a universal function dependent only on scale; we confirm this empirically. The combination of number and mass density laws can be used to develop stochastic compound multifractal Poisson processes which are useful new tools for studying and modelling rain. We discuss the implications of this for the rain rate statistics including a simplified model, which can explain the observed rain rate spectra.

3 viewed as a highly intermittent, highly heterogeneous process with turbulent energy, passive scalar variance and other fluxes concentrated into a hierarchy of increasingly sparse fractal sets; over wide ranges they are multifractal (see e.g. Anselmet (2001) for a recent review). In contrast, in the precipitation literature it is still common for turbulence to be invoked as a source of homogenization, an argument used to justify the use of homogeneous (white noise) Poisson process models of rain. Incredibly, data from our main source of quantitative rain information-weather radars-are still interpreted with the help of these unrealistic drop homogeneity assumptions (essentially those of Marshall and Hitschfeld (1953) and Wallace (1953)). In the notation of this paper, this is tantamount to assuming that the patch scale l c is kilometric rather than decimetric. Even disdrometer experiments commonly assume that the turbulence leads to at most weakly heterogeneous drop statistics (e.g. Jameson andKostinski 1999, Uijlenhoet et al 1999).
Today, high-powered lidars with logarithmic amplifiers and with small (metric scale) pulse lengths produce turbulent atmospheric data sets of unparalleled space-time resolution. Analysis of such reflectivity data from aerosols has shown that if the classical turbulence laws are given multifractal extensions to account for intermittency and scaling anisotropic extensions to account for atmospheric stratification and scaling wave-like space-time extensions to account for wave behaviour, that they account remarkably well for the statistics of aerosol pollutants (Lilley et al 2004, 2007, Lovejoy et al 2008b, Radkevitch et al 2007a. While aerosols are not purely passive, they have low Stokes numbers and low sedimentation rates and were considered as reasonable passive scalar surrogates (St l ; the Stokes number at scale l, is the dimensionless ratio of the particle inertial response time to the turbulent fluctuation time). Of more direct relevance to rain is the fact that the Corrsin-Obukhov law (applied to horizontal cross-sections) also appears to hold for liquid water concentration in clouds over wide ranges of scale (Davis et al 1996, Lovejoy andSchertzer 1995a) . In addition, starting in the 1980s, a growing body of literature has demonstrated that-at least over large enough scales involving large numbers of drops-rain has nontrivial space-time scaling properties (see e.g. Lovejoy and Schertzer (1995b) for an early review).
Combining the turbulence theory with raindrop physics involves several related difficulties. Firstly, rain is particulate. This is usually handled by considering volumes of space large enough so that many particles are present and spatial averages can be taken. However, since rain is strongly coupled to the (multifractal) wind field, even these supposedly continuous average fields turn out to be discontinuous (Lovejoy et al (2003)-at least for scales larger than a critical scale l c ; see below). This means that due to the systematic and strong dependence on the scale/resolution over which the rate is estimated that the classical treatment of the rain rate R(x, t) as a mathematical space-time field (without explicit reference to its scale/resolution) is not valid. Finally, rain does not trivially fit into the classical turbulence framework of passive scalars: rain moves with speeds different from that of the ambient air (it also modifies the wind field, but this is a smaller effect).
In recent years, the treatment of particles in flows has seen great progress especially with the pioneering work of Maxey (1987), Maxey and Riley (1983) relating the particle and flow velocities. The difficulty is that strictly speaking Maxey's relations only apply to particles whose dissipation scale Stokes number St η < 1, i.e. to small (cloud, aerosol) but not large (rain) drops. For the former, they have indeed led to fairly rigorous results (see especially Bec et al 2007, Falkovich and Pumir 2004, 2007, Falkovich et al 2002, and have notably helped in explaining the role of turbulence in enhancing the initiation of rain in clouds, especially the 4 importance of small vortices in augmenting the coalescence rate. The relations have also been used to estimate the effects of turbulence on cloud drop collision efficiencies (Franklin et al 2005, 2007, Wang et al 2005a, b, 2006, Pinsky et al 1999. Because rain particles can readily have St η > 100, rigorous application of Maxey's equations are not possible and rain has typically been treated with ad hoc parameterizations. However, we empirically find that below a critical scale l c -which we argue is the scale such that St lc = 1-the observed particle distributions (number and mass) are nearly that of a white noise. Since for scales l > l c we have St l < 1, it is therefore tempting to coarse grain the field over the statistically homogeneous (white noise) range l < l c (with St l > 1) and apply Maxey's equations to the coarse grained fields at l > l c for which St l < 1. Although the results cannot be considered theoretically rigorous, we use them to justify phenomenological inertial range turbulence models and to explain two basic scaling laws the classical Corrsin-Obukhov l 1/3 law for the mass density ρ, and a new l 1/2 scaling law for the number density n. Falkovich and Pumir (2004), Falkovich et al (2006) and Bec et al (2007) have admirably attacked the full drop/turbulence interaction at its most fundamental level, but progress has been difficult. Fortunately, at a more phenomenological level, things have been easier. Schertzer and Lovejoy (1987) proposed that even if rain is not a passive scalar it nevertheless has an associated scale-by-scale conserved turbulent flux. This proposal was made in analogy with the energy and passive scalar variance fluxes based on scale invariance symmetries but it lacked a more explicit, detailed justification. It nevertheless led to a coupled turbulence/rain cascade model which predicted multifractal rain statistics over wide ranges of scale. Today, this cascade picture has obtained wide empirical support, particularly with the help of recent large scale analyses of global satellite radar reflectivities which quantitatively show that from 20 000 km down to 4.3 km, that at these scales the reflectivities (and hence presumably the rain rate) are very nearly scale-by-scale conserved fluxes (Lovejoy et al 2008a). These global satellite observations directly show that a cascade structure behaviour is followed to within ±4.6% over nearly four orders of magnitude. Other related predictions-based essentially on the use of scaling symmetry arguments-are that rain should have anisotropic (especially stratified) multifractal statistics (see Lovejoy et al 1987) and that rain should belong to three-parameter universality classes Lovejoy 1987, 1997), see the review (Lovejoy and Schertzer 1995b). However, these empirical studies have been at scales larger than drop scales and outstanding problems include the characterization of low-and zero-rain rate events and the identification of the conserved (cascaded) flux itself. In other words, up until now the connection with turbulence has remained implicit rather than explicit.
In order to bridge the drop physics-turbulence gap and to further pursue phenomenological approaches, data spanning the drop and turbulence scales were needed. Starting in the 1980s, various attempts to obtain such data have been made; this includes experiments with chemically coated blotting paper (Lovejoy and Schertzer 1990) and lidars, (Lovejoy and Schertzer 1991). The most satisfactory of these was the HYDROP experiment (Desaulniers-Soucy et al 2001) which involved stereophotography of rain drops in ≈10 m 3 volumes typically capturing the position and size of 5000-20 000 drops. The nominal rain rates were between 2 and 10 mm h −1 ; for information about HYDROP see table 1 and Lilley et al (2006), for an example see figure 1. Analyses to date (Lilley et al 2006, Lovejoy et al 2003 have shown that at scales larger than a characteristic scale determined by both the turbulence intensity and the drop-size distribution (DSD) (but typically around 40 cm, see e.g. figure 2), that the liquid water content (LWC; i.e. the mass density), number density and other statistics behave in a scaling manner as predicted Table 1. Various characteristics of the HYDROP dataset. The critical scale l c is determined from the spectral minimum estimated from the drop mass spectra. The drop mean relaxation length is estimated from the mean drop diameterL using formula (3b). The energy flux ε is estimated using equation (12). The Stokes number St, and the sedimentation number Sv are estimated from equation (13) using the dissipation scale (determined from l c estimated from the spectral minimum). In all cases, the spread ('±') is the storm-to-storm variability based on three triplets per storm (the exception being storm 207 for which there were 7 triplets). The mean diameter of the 19 triplets = 1.19 ± .17 mm (i.e. ±14%; this spread is triplet-to-triplet variation, not standard deviation for each triplet), the mean relaxation distance =1.97 ± 0.29 m (i.e. ±15%), the mean number of drops is 23 200 ± 11 800 (i.e. ±51%). The numbers of drops are larger than those in Lilley et al (2006) since all the reconstituted drops were used, not only those in the more reliable central region. The LWC statistics is for the welllit central region only, averaged at 70 cm scale (roughly the relaxation scale, l c ). The coalescence speed was estimated from the formula ϕ l = g −1/4 l 1/4 R ε 1/2 l . The mean interdrop relaxation speed difference v R,drop is the mean difference in relaxation speed averaged over all pairs of drops in the volume. v R,n and v R,ρ are calculated by averaging the relaxation speed over cubical regions 70 cm on a side, (weighted by n and ρ, respectively) and then calculating the mean differences between the neighouring 'cubes' (70 cm is roughly l c so that the differences are at the small-scale end of the density variance flux cascade).  by the cascade theories. While these results suggest that rain is strongly coupled to the turbulent wind field at scales larger than 40 cm (potentially explaining the multifractal properties of rain observed at much larger scales), the analyses did not find an explicit connection with the standard turbulence theory. See also Marshak et al (2005) for similar empirical results for cloud drops.
The key new result of this paper concerns the fluctuations in the particle number density n, which we find-both theoretically and empirically-follows a new n ≈ l 1/2 law with prefactors depending on turbulent fluxes. This explicit link between the particle number density fluctuations and turbulence is expected to hold under fairly general circumstances (including perhaps for small cloud drops) and provides the basis for constructing compound multifractal-Poisson processes. While the traditional approach to drop modelling (see e.g. Pinsky 1995, Srivastava andPassarelli 1980) is to hypothesize specific parametric forms for the DSD and then to assume spatial homogeneity in the horizontal and smooth variations in the vertical, our approach on the contrary assumes extreme turbulent-induced variability governed by the turbulent cascade processes and allows the DSD to be constrained by the turbulent fields n and ρ, i.e. by the 0th and 1st moments of the number size density (or equivalently on the 0th and 3rd moments of the DSD).
This paper is structured as follows. In section 2, we discuss the relation between the drop and wind speeds and use this to interpret the break in the liquid water spectrum at scale l c as the transition between turbulent-dominated (large) scales and drop inertia-dominated (small) scales; the scale at which the Stokes number St lc ≈ 1. This interpretation allows us to estimate the turbulent dissipation rate and the dissipation scale Stokes number St η , which we find of the This shows the 3D isotropic (angle-integrated) spectrum of the 19 stereophotographic drop reconstructions, for ρ, the particle mass density. Each of the five storms had 3-7 'scenes' (from matched stereographic triplets) with ≈5000−40 000 drops (see table 1) each taken over a 15-30 min period (orange = f207, yellow = f295, green = f229, bluegreen = f142 and cyan = f145; the numbers refer to the different storms). The data were taken from regions roughly 4.4 × 4.4 × 9.2 m 3 in extent (slight changes in the geometry were made between storms). The region was broken into 128 3 cells (3.4 × 3.4 × 7.2 cm 2 , geometric mean = 4.4 cm); we use the approximation that the extreme low wavenumber (log 10 k=0)correspondstothegeometricmean,i.e. 5.6m,theminimacorrespondtoabout40-70cm;seetable2).Thesinglelowest wavenumbers (k = 1) are not shown since the largest scales are nonuniform due to poor lighting and focus on the edges. The reference lines have slopes −5/3, +2, i.e. the theoretical values for the Corrsin-Obukhov (l 1/3 ) law and white noise, respectively.
order ≈200-300. At large scales l > l c , where the raindrops have low Stokes numbers (St l < 1), we argue that Maxey's relations can be used. Using them in section 3, we derive the basic equations for number density, mass density, the corresponding flux densities (including rain rate) and the coalescence speed. In section 4, we derive the l 1/3 and l 1/2 laws for mass and number densities respectively, and verify the results empirically on the HYDROP data; and in section 5, we conclude the paper.

Drop relaxation times, distances and velocities
Raindrops are inertial particles moving in turbulent flows under the influence of gravity. For an individual drop moving with velocity v in a wind field u, Newton's second law implies: where we have taken the downward acceleration of gravity (g) as positive, and assumed a power law drag with exponent η d . V R (M) is the 'relaxation velocity' of a drop mass M and δv is the difference in the particle and wind velocities. In still air, after an infinite time, V R is the terminal drop velocity (we use capitals for single drop properties, lower case for various averages). The corresponding relaxation times T R (M) and distances L R (M) are given by: (2) These are the typical time and length scales over which the drops adjust their velocities to the ambient air flow. According to classical dimensional analysis for the drag on rigid bodies, the low Reynolds number regime is independent of fluid density while the high Reynolds number (Re) regime is independent of fluid viscosity. This implies that in the first regime, the drag is linear (η d = 1; 'Stokes flow'), while in the second, it is quadratic (η d = 2) so that the drag coefficient (C D ) is constant. For these regimes we have: where L is the drop diameter, ρ w and ρ a are the densities of water and air, C D is the standard drag coefficient and η v is the dynamic viscosity. In the Stokes flow, we have ignored numerical factors of order unity. In order to decide which relations are more appropriate for rain drops, consider the theory and experiment reported in Le Clair et al (1972). They show that for spheres the high Re limit is obtained roughly for Re = 10 2 where C D ≈ 1 (although it decreases by a factor of 2-3 as Re increases to ≈5000). According to semi-empirical results presented in Pruppacher and Klett (1997), Re ≈ 1 is attained for drops with diameter (L) ≈ 0.1 mm, whereas Re = 10 2 is attained for drops diameter 0.5 mm. On this basis, Pruppacher proposes breaking the fall speed dynamics into three regimes with (roughly) L < 0.01 mm, 0.01 < L < 1 mm and L > 1 mm. The small L 'Stokes' regime has η d = 1, the second regime is intermediate and the third (high Re) regime has roughly η d = 2 up to a 'saturation' at L ≈ 4 mm. Although for rain, this regime is more complicated than for rigid spheres due to both drop flattening and internal drop flow dynamics, according to data and numerics reviewed in Pruppacher and Klett (1997), the basic predictions of the η d = 2 drag law (V R ∝ L 1/2 and L R = α D L) are reasonably well respected over the range 0.5-4 mm. For example, using Pruppacher's data on L R at 700 mb and a drag coefficient of C D = 0.8, the predictions of equation (3b) are verified to within ±25% over this range. Since at 20 • C, pressure = 1 atmosphere, ρ w /ρ a ≈ 1.29 × 10 3 , this implies α D = 1690, so that for drops with diameters 0.2 < L < 4 mm (the range detected by the HYDROP experiment discussed below), we expect L R to be typically of the order 40 cm to 7 m (hence, T R is in the range 0.2-0.8 s). 9

Empirical estimates of l c using drop stereophotography
The transition from rain drop behaviour dominated by turbulence to behaviour dominated by drop inertia will occur at a scale L c , which is determined by both L R and the level of turbulence. Before discussing other aspects of the turbulence theory, we give the direct estimates of L c in rain averaged over the drop sizes (denoted l c ) using the HYDROP data.
We recall that each of the five HYDROP storms had 3-7 'reconstructions' (from matched stereographic triplets) each with 5000-20 000 drops and taken over a 15-45 min period (see table 1). By 'reconstruction' we mean the determination of the position and size of nearly all (>90%) of the particles >0.2 mm in diameter. In Lilley et al (2006), the scale dependence of various densities (number, liquid water, etc) was estimated by spatial averaging over spheres of increasing diameter. It was noted that while statistical multiscaling-characteristic of turbulent cascade processes-was evident at scales greater than about 40 cm or so (depending somewhat on the reconstruction)-that at smaller scales there was no clear behaviour. Therefore, in order to get a clear statistical characterization over the full range of available scales, we turn to Fourier techniques.
To estimate the spectrum, we first discretize the liquid water density ρ on a grid (whose pixel resolution ≈4.4 cm corresponds to the accuracy with which the positions were estimated), and then we determine the discrete estimate of the Fourier transform: We introduce the power η for future use since by varying η we can conveniently consider various fields: for example, if we adopt the convention that x 0 = 1 for all x > 0 and x 0 = 0 for x = 0, and we discretize ρ η , then η = 0 corresponds to the particle number density (in roughly 8% of the cases where there is more than one drop in a 'pixel', we sum ρ η over all the particles in the pixel). Similarly, η = 1 is the usual mass density, η = 7/6 is proportional to the 'nominal rain rate' (using V R for fall speeds and equation (3b)) and η = 2 is the radar reflectivity factor. We may note that the discretization scale used here corresponds roughly to the accuracy of the position measurement; its use prevents aliasing of the spectrum (which would occur, for example, if spuriously precise positions were used).
Fromρ η k , we then determine the power spectral density: where the ' . ' indicates ensemble averaging, here estimated by averaging over the different reconstructions for each storm. Finally, the isotropic spectrum E η is estimated by angle integration in Fourier space: Figure 2 shows the 3D isotropic (angle-averaged) spectrum of the 19 stereophotographic drop reconstructions averaged over each of the five storms for η = 1 (the usual mass density). The low-wavenumber reference line indicates the theoretical Corrsin-Obukhov k −5/3 angleintegrated spectrum. Since for Gaussian white noise P η is constant, in d-dimensional space the angle-integrated spectrum varies as k d−1 (here d = 3, hence as k 2 ) and is indicated by the high-wavenumber reference line. It can be seen that the transition between the passive scalar behaviour-where the drops are highly influenced by the turbulence and the high-wavenumber For each storm, this shows the mean distance L min (in m) corresponding to the minimum of the spectrum of ρ η as a function of the order of moment η (from top to bottom on the right-hand side: the experiments (corresponding to different storms) numbered 142, 145, 295, 207, 229-see table 1 for descriptions). In the text, the value of L min for η = 1 is taken as an estimate of the critical scale l c at which St = 1. Since higher η values weight the spectra to the larger drops, so there is a slow increase of l c with increasing η. The main exception is the 142 case, which is the most dominated by small drops and shows a near doubling in l c when comparing the number density (η = 0) with the mass density (η = 1). regime, where the drops are totally chaotic (white noise) is quite sharp and occurs at critical scales l c of roughly 40-75 cm (see table 1). The overall physical interpretation of the spectrum is the following: due to the turbulent wind field, the drops are concentrated in patches of size l c (each with many drops; the mean interdrop distance here being ≈10 cm), but within each small patch, due to the transition to a white noise spectrum-the drop distribution is fairly uniform. This suggests a simple model (explored below) of the drop collision microphysics as one of statistically independent drops colliding within patches whose overall drop number and water concentration varies tremendously from patch to patch due to the turbulence. In this way, the turbulence can drive the process and constrain the microphysics.
It is significant that the transition is very well pronounced: even though the measured drop diameters varied by a factor of about ≈20, the transition systematically occurs over a range of factors ≈1.5 or so in scale. Figure 3 shows the effect of the relatively small albeit systematic changes in the positions of the spectral minima that occur by weighting the spectra more towards the small drops (small η), or the large drops (large η). Section 4 and figure 8 explore the spectra for η = 0 (corresponding to the number density) in more detail and table 1 gives some information about the various experimental conditions as well as parameter estimates obtained from the spectra.
More details on the HYDROP experiment can be found in Desaulniers-Soucy et al (2001), and for these datasets see Lilley et al (2006). The main differences between the part of the data analysed in the latter and that analysed below is that in the latter, a very conservative choice of sampling volume was used. By using only data from the region best lit and most sharply in focus, >90% of the drops with diameter >0.2 mm were identified. In this study, we sought to get a somewhat wider range of scales with as many drops as possible by exploiting the fact that Fourier techniques are very insensitive (except at the lowest wavenumbers) to slow falloffs in the sensitivity (and hence drop concentrations) near the edges of the scene (this is a kind of empirically produced spectral 'windowing' close to the numerical filters used to reduce leakage in spectral estimates). We therefore took all the data in a roughly rectangular region about 4.4 × 4.4 × 9.2 m 3 in size (geometric mean = 5.6 m), and then analysed only the spectrum for wavenumbers k > 2 (i.e. spatial scales 2.3 m). This somewhat larger scene size with respect to the previous HYDROP analyses yielded an increase in the range of scales by a factor of about 2.8 and about 2-3 times more drops.

Applying Maxey's equations to rain drops in the turbulent inertial range
In regimes of 'weak turbulence' (where Du Dt g), and power law drag (η d > 0), we may follow Maxey (1987), Maxey and Riley (1983), Falkovich and Pumir (2004) and Falkovich et al (2006) (who considered η d = 1), and obtain an expansion for the drop velocity v in terms of the wind velocity u: although he only investigated the η d = 1 case in detail, Maxey (1987) pointed out the insensitivity of this result to the exact form of the drag law (up to the first order, it is independent of η d ; note that we use the notation D/Dt = ∂/∂t + u · ∇ for the Lagrangian derivative of the wind). We should note that this equation implies that even if the wind field is incompressible, the drop velocity field is compressible: where ω is the vorticity, and s is the strain rate (Maxey 1987); this relation will be used below.
In turbulent flows in the turbulent inertial range, we may use the standard turbulence estimate of derivatives at scale l: where τ e,l = l/ u l is the eddy turnover time (i.e. eddy lifetime) for eddies of size l. Applying this in equation (7) we obtain: where St l is the scale l Stokes number evaluated at the eddy turnover time: τ e,l = l 2/3 ε −1/3 l and Sv l is the scale l 'sedimentation number' (see e.g. Grabowski and Vaillancourt 1999).
We see that at least if l is the dissipation scale η and St l < 1, the series converge. In what follows, we make the plausible but unproven assumption that equation (10) continues to be 12 approximately true for drop inertial scales l as long as the derivatives are estimated at scale l l c (see, however, Falkovich and Pumir 2007, Falkovich et al 2002, Wilkinson et al 2006. In this case, the observed statistical homogeneity (the white noise for l < l c in figure 2) of the scales with l < l c makes the assumption plausible.
In the turbulent inertial range, we obtain: where a l is the acceleration fluctuation. Since Sv l ∝ St l l 1/3 , we see that at large enough scales l, the sedimentation effect will dominate the drop inertial effects. However, if we seek to interpret the scale break in the spectra (figure 2), then it is the drop velocity fluctuations v l as functions scale l we must consider, not v l directly. From equation (10a) we see that as long as the fluctuations in the sedimentation velocity v R,l , are small compared to the fluctuations in the turbulent velocity (i.e. as long v R,l u l ) then the critical break scale l c in figure 2 is the scale at which St l = 1 so that for l < l c , the drop inertia is dominant. Although the common meteorological assumption that v R is horizontally homogeneous (and hence v R,l ≈ 0) is unjustified (cf figure 5), it seems at least plausible that it is smaller than u l . The fluctuation v l can be estimated either by first averaging equation (10a) over a scale l and then taking differences between neighbouring l sized patches, or by taking the average difference between v R for drops separated by distance l or less. The results of both of these definitions are given in table 1; we see that at l = 70 cm (≈ l c ) in table 1; the former is about a factor 2 smaller than the latter, and we see below that in turn this is about a factor 2 smaller than u l . To compare this with u l , we need an estimate of the turbulent energy flux ε. As indicated below, this can be obtained by neglecting v R,l in comparison with u l using the condition St l,c = 1. When this is done (table 1) we find that u l c is indeed 2-6 times larger than v R,l . This is an ex post facto justification of the assumption that at l c , v R,l is negligible with respect to u l .
If this interpretation is correct, l c is the critical decoupling scale at which the velocity and acceleration terms in equation (7) are of equal magnitude. It is also the scale at which (according to theoretical considerations in Wang and Maxey (1993) and experimental results in Fessler et al (1994), we expect the maximum effect in causing preferential concentration of particles). With this assumption, using equation (11) with the drop averaged l R , we see that the critical scale l c at which St l = 1 is: (t R is the relaxation time averaged over the drops). If we now consider the typical velocity difference v l between two drops separated by a distance l: where a l is the fluctuation (difference) in acceleration of the wind at the scale and we have ignored higher order velocity derivatives in the l l c case. As discussed above, the l l c result follows if we assume that the spatially averaged fluctuations in v R at scale l tend to either decrease as l increases or at least increase more slowly with scale than the turbulent u l which increases as l 1/3 . Empirically it implies that the drop velocity spectrum is the same 13 as the Kolmogorov wind spectrum, and thus at least compatible with the observed Corrsin-Obukhov passive scalar spectrum, figure 2. We have assumed that the sedimentation term gives a small contribution (i.e. v R ≈ 0 even if the drop averaged mean relaxation velocity v R is not negligible; see section 3 for more discussion of this). Also, encouraged by the relatively sharply defined spectral minimum in figure 2, we have used a mean scale l c rather than a more precise individual drop diameter (equivalently mass)-dependent relation. For scales l > l c since the drop velocity fluctuation statistics are essentially the same as those for the wind (i.e. they are both k −5/3 ), it is plausible that equation (13) can explain the observed rain drop mass density spectrum ( figure 2). Similarly, at scales l < l c , where the differences between particle velocities v depend on the highly variable wind acceleration gradients, it is perhaps not so surprising that the spectrum follows Gaussian white noise: the drop inertia dominated regime is 'chaotic'. It is therefore natural to associate the scale at which the spectral minimum occurs with the transition equation (13) from viscous forces to drop inertial forces at the high wavenumbers.
The interpretation of the break in figure 2 as the critical drop inertial (St l = 1) scale is sufficiently important that it is worth mentioning that an alternative transition between viscous and gravitational forces (i.e. sedimentation) is sometimes invoked (for much smaller particles, cloud drops, see e.g. Grabowski and Vaillancourt (1999), and see Lilley et al (2006) for a similar argument). In this case, the critical sedimentation scale l cs is obtained by balancing the drop size-averaged relaxation ('terminal') velocity v R with the turbulent velocity gradient ε 1/3 l cs l 1/3 cs . This leads to the following estimate of the critical sedimentation scale l cs : To distinguish the two explanations for the spectral transition (i.e. at l c or at l cs ?), we can use the overall mean values l R ≈ 2.0 m, t R = 0.45 s, v R = 4.27 m s −1 (these are the means of l R , v R,n ,see table 1; v 2 R /l R is not exactly equal to g since v R and l R are averages of different moments of the drop volumes). These values combined with the use of the spectral minima to estimate l c and l cs , and with equations (12) and (14) can be used to estimate the intensity of the turbulence (ε) implied by the two different scenarios (i.e. St lc = 1 or Sv lcs = 1). Assuming the critical spectral scale L min is l c (equation (12)), the results for individual storms are shown in table 1; they are all not too far from ε ≈ 4 m 2 s −3 , which is a large but still plausible value (the value ε = 10 −4 m 2 s −3 is often considered a typical atmospheric value; however, the intermittency is huge so that during storms values several orders of magnitude larger may indeed be realistic: Grabowski and Vaillancourt (1999) suggest that ε = 0.1 m 2 s −3 is a large but plausible value in clouds but they do not mention the scale at which the value should apply: due to intermittency, large values are more common at smaller scales and values in rain are plausibly larger than that in clouds). In comparison, using equation (14) and identifying the spectral minimum instead with l cs rather than with l c , we obtain much larger and probably unrealistic values: the overall mean being ε ≈ 160 m 2 s −3 . This makes it unlikely that equation (14) could explain the observations. We therefore conclude that equation (12) is indeed valid. With this assumption, we can estimate Sv lc , the sedimentation number at the decoupling scale l c (table 1); the overall mean is ≈3.7. This indicates that the sedimentation scale l cs is larger than l c (see table 1); the overall mean is ≈1.2 m, i.e. 2-3 times larger than l c .
If we accept the position of the spectral minima as an estimate l c , then we can also estimate the values St η and Sv η , i.e. the values at the dissipation scale l η = ν 3/4 ε −1/4 . First with kinematic viscosity η ≈ 1.5 × 10 −5 m 2 s −3 (roughly the value for air at 20 • C, pressure 1 atmosphere), we obtain an overall average l η ≈ 170 µm (see table 1). There is surprisingly little storm-tostorm spread so that the overall mean values St η ≈ 220 and Sv η ≈ 60 give a reasonable idea of the relative magnitudes (table 1). We can now compare these values with those of Bec et al (2007), who performed large scale numerical simulations of monodisperse inertial particles in 3D hydrodynamic turbulence (see also the experimental results of Aliseda et al (2002)). The main differences were: (i) the Stokes numbers in the simulations were much smaller; 0.16 < St η < 3.5, (ii) there was no gravity (so Sv η = 0), (iii) that the particles were small enough so as to be in the Stokes flow (η d = 1) range (in order to simulate cloud drops). Perhaps the most important result was the finding that particles cluster right through the turbulent inertial range. Characterizing the clustering by a fractal correlation dimension (D 2 ), Bec et al (2007) found that D 2 is dependent on St η (hence on drop size). Indeed, the correlation codimension (= 3 − D 2 ) varied from about 0.7 for St η = 1 to 0 for St η ≈ 3.5 suggesting that with the full mix of particle sizes, the number density measure would be multifractal as found in the HYDROP experiment (Lilley et al 2006).

The drop mass density function N
Having argued that St l < 1 at scales l > l c , we now systematically use equation (7) to relate the turbulent wind field to the drop velocities. To do this, we first introduce N = N (M, x, t), which is the drop-mass density function, i.e. the number of drops with mass between M and M + dM per unit volume of space at location x , time t (the corresponding function in terms of drop radius the DSD). Ignoring condensation and drop break up, the usual coalescence (Smoluchowski) equation (used, for example, in cloud and rain modelling (Srivastava and Passarelli 1980)) can be written: The right-hand side term is the coalescence operator in compact notation (see Lovejoy et al 2004), N 1 |H |N 2 : If we now assume that the drop-drop collision mechanism is space-time independent, then the full coalescence kernel H can be been factored H M , M, x, t = ϕ x, t h (M , M) so that space-time variations in the coalescence rate are accounted for by ϕ, which is the coalescence speed and h(M , M) is a time-independent kernel characterizing the drop interaction mechanism; we return to this in section 3.4. The budget equation for N is thus: The validity of this turbulence-drop coalescence equation can be verified by integrating over a volume; the left-hand side is the change in the volume of particles between M and 15 M + dM, the first right-hand side term is the flux of such particles across the bounding surface, the second is the change in the number within the volume due to coalescence. In this equation, v(M, x, t) is the velocity of a particle mass M. While this equation was invoked in Falkovich and Pumir (2004), only the static v = 0 case was studied. Equation (15) is also used in conventional DSD modelling but typically with the additional assumption of horizontal homogeneity (only smooth vertical variability is considered). Although a great deal of effort has gone into studying various possible kernels, few attempts have been made to study the constraints placed on the microphysics by the larger scale turbulence dynamics-our goal here.

The mass and number densities and fluxes
We now consider the first two moments of the turbulence-drop coalescence equation (15); the number density (n) and drop mass density (ρ), as well as their fluxes: the mass flux r and the number flux N : Note that the vertical component of r is the usual rain rate; R = r z (see section 3.4). By respectively integrating equation (17) with respect to M and by multiplying equation (17) by M and then by integrating with respect to M, we can obtain the budget equations for the particle number and mass densities: where we have used the fact that the coalescence operator conserves drop mass: We see that because of equation (23), coalescence is not directly relevant for the ρ equation (22) whereas it is relevant for the n equation (21). Also for the same reasons, taking into account drop break up will add a new term to equation (21) but will not affect equation (22). We now introduce the average relaxation velocities v R,n and v R,ρ weighted by number and mass densities, respectively: Substituting v from equation (7) into equations (19) and (20) and using definitions (24) and (25), we obtain: so that α n and α ρ are the effective velocities of the drops. Note that here and below, the derivatives are understood to be estimated at the scale l = l c . It is of interest to consider the spatial variability of v R,ρ and v R,n and also how they are related to each other. Figure 4 shows the scatter plots obtained from the HYDROP data by using the high Re theoretical relaxation formula (equation (3)b). Table 1 shows the storm-bystorm mean as well as the spread between scenes within the same storm, averaging over a 70 cm resolution, i.e. at about the relaxation scale l c . From the table, we can see that the mean over each storm is fairly stable, whereas from figure 4 we see that the variability is very large. Although there is clearly no one-to-one relation between v R,ρ and v R,n , the following relation is roughly satisfied: v R,n = av b R,ρ with a = 0.95 ± 0.07, b = 0.70 ± 0.10 (the standard deviations being the storm-to-storm spread in values, for v in units of m s −1 ). This shows that v R,ρ and v R,n cannot be taken to be equal. To get an idea of the spatial variability, we also calculated the mean horizontal spectra of v R,ρ and v R,n ( figure 5(a)). Although the range of scales is small, the figure suggests that the spectrum is dominated by the low rather than high wavenumbers (of interest below); possibly even with a roughly Corrsin-Obukhov spectrum at low wavenumbers. While this result follows if ρ has a Corrsin-Obukhov spectrum (since the v R s are the powers of ρ; see section 4.4, figure 8(c)), note that it is contrary to the usual meteorological assumption that v R is horizontally homogeneous (being determined by a spatially homogeneous DSD).  To investigate the variability in the vertical gradients, we determined the spectra of the vertical gradients of v R,ρ and v R,n (using finite difference estimates of ∂v R,n /∂z and ∂v R,ρ /∂z), these are needed in section 3.3 and are shown in figures 6(a) and (b). We see that in both cases the spectrum shows a slight tendency to rise at higher wavenumbers; if the spectrum is E(k z ) ≈ k −β z , then roughly β ≈ −0.3. Recall that whenever β < 1, that the variance is dominated by the small scales/large wavenumbers (an 'ultraviolet catastrophe'). Since the spectrum of ∂v R /∂z is k 2 z times the spectrum of v R , this is consistent with a near-Kolmogorov value β ≈ 1.7 in the vertical gradients.

The mass and number fluxes
To put the mass flux vector into a more useful form, we can appeal to the dynamical fluid equations: where we have used the incompressible continuity equation for u adequate for our purposes, see Pruppacher and Klett (1997) and F r is the reaction force of the rain on the wind and p is the pressure. The reaction force per volume of the rain on the air is: plus higher order corrections. We therefore obtain: We can now note that the critical decoupling scale l c is typically much greater than the turbulence dissipation scale so that at the scale l c , we can neglect the viscous term. In addition, since we consider scales l > l c and St l < 1, we neglect the D 2 u Dt 2 term. We therefore obtain the approximation: and hence using this in equations (26) and (27), we obtain: Finally, if we make the hydrostatic approximation (∇ p = g(ρ + ρ a )ẑ), we see that this reduces to a simpler form: i.e. the effective velocities are equal to the wind velocities plus the appropriate mean relaxation speeds. To judge the accuracy of the hydrostatic approximation, we can consider the dimensionless difference 1 − 1 gρ a ∂ p ∂z , which in nonprecipitating air at scales above the viscous scale, should be equal to the dimensionless vertical acceleration: ≈ 1 g Dw Dt . Although we did not have measurements during HYDROP, this difference has been evaluated empirically using stateof-the-art drop sondes with roughly 5-10 m resolution in the vertical during the 'NOAA Winter Storms 04' experiment (North Pacific Ocean). It was found that the typical mean value at the surface was about 0.02 with fluctuations of the order ±0.005 so that, even at 10 m scales, we may use the hydrostatic approximation to reasonable accuracy. In addition, at higher wavenumbers, the vertical spectrum of 1 − 1 gρ a ∂ p ∂z ≈ 1 g Dw Dt (figure 7) increases nearly linearly with wavenumber showing its strong dependence on the small scales. This suggests that the (notoriously difficult to measure) vertical velocity (the integral of the acceleration) has a spectrum k 2 times smaller, i.e. roughly k −1 (see also the discussion on the vertical velocity in Lovejoy et al (2008b), where similar conclusions are reached).
Using the hydrostatic approximation (for α n and α ρ , but not for their divergences, where we only neglect the Du/Dt · ∇v R term), we therefore obtain the following equations for the number and mass densities: where we note the additional coalescence term in the equation for n.

The rain rate
In a similar manner, we can obtain the z-component of the mass flux r , which is the rain rate R: where for the vertical wind we have used the notation w = (u) z . Although the full ramifications of this equation for rain will be developed elsewhere, it should be noted that this approximation will break down in regions where the particle number density n (strongly correlated with ρ) is so small that there is low probability of finding a drop (see the discussion of the necessary compoundcascade/Poissonmodelinsections4.3and4.4andequation (63)).Inthisway,zerorain rate regions simply correspond to regions with very small n.
To understand how this may affect the rain statistics, consider the following simple model. Since n and ρ are highly correlated, we may set ρ to zero whenever it is below some threshold ρ t . If for the moment we ignore the term v R,ρ the model would be: We now take ρ to have Corrsin-Obukhov statistics ( ρ ≈ l H with H = 1/3 corresponding to a spectrum k −β ; ignoring intermittency corrections, β = 1 + 2H = 5/3 as found in the HYDROP experiment), and w to have H ≈ 0 statistics (β ≈ 1; as inferred indirectly from the drop sonde data, figure 7). We seek the spectrum of R. First consider the effect of the threshold. At large scales it imposes a fractal support, which flattens the spectrum; at high wavenumbers it remains smoother than w so that the product ρ t,w will have high wavenumber statistics dominated by the vertical wind ≈k −1 , whereas at low wavenumbers, it will be much shallower (depending somewhat on the threshold and the value of H ). In addition, if we confine our rain analyses to regions with no zeroes, then the break disappears and we expect roughly k −1 statistics.
If we now consider the (neglected) ρv R,ρ in equation (38), we see that it will have much smoother, nearly Corrsin-Obukhov (β ≈ 5/3) statistics; its contribution to the spectrum will thus be dominated by the vertical wind term ρw. These conclusions have been substantiated by numerical simulations and will be the subject of a future paper. This simple model may well be sufficient to explain observations of rain from high-(temporal) resolution raingauges: several studies have indeed found corresponding ω −0.5 spectra at long times but ω −1 spectra at short times (see e.g. Fraedrich and Larnder 1993) with transitions occurring at scales of 2-3 h (see de Montera et al 2007).

Energy flux, number and mass variance fluxes
In figure 1, we noted that at scales l > l c , ρ accurately obeyed the Corrsin-Obhukov law for passive scalars: E ρ (k) ≈ k −5/3 . The reason is that the variance flux χ = −(∂ρ 2 /∂t) is 'conserved' by the nonlinear terms and can therefore only be dissipated at small scales. To see this, multiply equation (22) by −2ρ and use equation (27) to obtain: Integrating over a volume and using the divergence theorem, we see that only the ρ 2 ∇ · α ρ represents a dissipation of the variance flux. We therefore obtain: (with the hydrostatic approximation). These terms are only important at small scales. This is true of the vorticity/shear term since E ω (k) = E s (k) = k 2 E u (k) (the standard result in isotropic turbulence, which is a consequence of the fact that ω and s are the derivatives of velocity). In Kolmogorov turbulence, therefore both ω 2 and s 2 have roughly k 1/3 spectrum; figure 5(a) shows that v R,ρ varies more smoothly so that the far right term in equation (41) is dominated by high wavenumbers. Furthermore, we have seen from figure 6 that the spectrum of ∂v R /∂z is also dominated by the small scales and also has a near k 1/3 spectrum. If we estimate ω 2 /2 − s 2 ≈ ( u l /l) 2 ≈ ε 2/3 l l −4/3 , then from table 1 we see that at l c , ∂v R /∂z can be neglected with respect to (v R,ρ /g)(ω 2 /2 − s 2 ). Since the dissipation depends only on t R = v R /g and 21 on ε l , using dimensional analysis we obtain l ρ,diss = ε 1/2 l t 3/2 R = l c (cf equation (12); alternatively, we obtain the same result by estimating ρ in equation (41) from ρ in the scaling regime discussed later, equation (55). Since the above arguments are only valid for l > l c , this at least demonstrates the self-consistency of the model.
Similarly, for the dissipative part of the number variance flux (from equation (21)): so that we have: (again using the hydrostatic approximation). The dissipation of number density flux thus has a first term of the same form as the mass density flux, but an additional coalescence term, which does not depend on spatial gradients. We therefore expect it to contribute to the dissipation of ψ over a wide range of scales with intensity strongly dependent on the local drop concentration. Indeed, since n is related to ρ via the DSD, it will be the scale l ρ,diss = l c which will be the inner scale for the n field, roughly as observed (cf figures 3 and 8). However, n does directly define an important scale: the local inter-drop distance l inter = n −1/3 , where the field description breaks down; this is discussed below.

Coalescence
From the empirical spectrum and supported by the theoretical considerations of section 2, we have developed a picture of rain being composed of l c -sized statistically homogeneous (white noise) 'patches'. While each patch has a highly variable liquid water content (determined by the turbulent dynamics), within each patch (corresponding to scales with St l > 1), the variation is white noise suggesting that within a patch the particles are statistically independent. If the mean inter-drop distance l inter is smaller than l c , so that each patch contains many drops, then the key coalescence processes are greatly simplified. In this section, we examine the consequences. We now consider the coalescence term in more detail: where we have factored the general kernel H(M, M , x, t) into a purely mass-dependent h(M, M ) and a time-space varying speed term ϕ. The h(M, M ) kernel describes the basic collision mechanism, whereas the ϕ determines the coalescence speed. Fairly generally, we may write: where v is the difference in velocities of the colliding drops and E takes into account the geometry and collision efficiency (typically taken as a power law of the mass ratios, see e.g. Beard and Ochs (1995)), for much recent work on these factors in cloud drops, especially the impact of turbulence, see Franklin et al (2005Franklin et al ( , 2007, Pinsky et al (1999,) Wang et al (2005a, b, 2006. Following the empirical spectrum, which indicates a rapid transition from drops following the flow at scales larger than l c to drops with apparently white noise statistics at smaller scales, it is natural to consider the simplified model that the particles follow the flow (but with an extra sedimentation velocity v R ) down to l c (roughly independent (a) (b) (c) Figure 8. (a) Same as previous but for n, the particle number density (calculated using an indicator function on a 128 3 grid. The reference lines have slopes −2 and +2, the theoretical values for the l 1/2 law and white noise, respectively. (b) This shows the ratio of the ensemble spectra E ρ (k)/E n (k); for each of the 5 storms, and the overall ensemble (purple), with the theoretical reference line slope −1/3. This is a sensitive test of the prediction of equations (55) and (59). (c) These are mean β values for the different ηs with the bars indicating storm-to-storm differences. The fits were for wavenumbers <70 cm −1 (the turbulent regime; l > l c , St l < 1). Note that for η = 0 (number density) and η = 1 (mass density), the theoretical predictions (2 and 5/3; ignoring intermittency corrections) are well verified since β = 2.02 ± 0.11 and 1.74 ± 0.12, respectively. See figure 3 for the corresponding spectral minima.
of particle mass) and then that they are independent of the flow at smaller scales. Using this approximation, when two particles collide, we may thus assume that since the moment that they decouple from the flow a distance l c away they have followed roughly linear trajectories and according to equation (13) their typical velocity difference upon collision will be: where t R is the mean relaxation time. Note that as discussed in the derivation of equation (13), we have neglected the gradient of the relaxation velocities which (empirically, table 1) are somewhat smaller than the turbulent velocity gradients. Using the turbulent inertial range estimate of the acceleration ( a l = v l /τ e,l ) and t R = (l R g) 1/2 , we obtain: a very weak dependence on l R , which gives some ex post facto justification to the lumping of all the drops together (independent of their diameters) and using mean values. Using this velocity difference in the interaction kernel (equation (45)), we obtain: or: i.e. the coalescence speed ϕ is equal to the turbulent velocity gradient at the mean relaxation scale (the subscript on the integral denotes spatial averaging at the relaxation scale), the dimensions of h are (length) 2 . Although ε(x, t) and, hence, the coalescence speed ϕ(x, t) is highly variable, it defines characteristic speeds for the coalescence process: ε 1/2 l c l 1/4 R g −1/4 . We note that this approximation implies that the inter-drop collision rate is determined primarily by the turbulence and not directly by the mass-dependent relaxation velocities (as is usually assumed). As a final comment, we could consider the possibility of applying this to cloud drop dynamics. Clearly, we would have to include terms representing condensation, however, as far as the coalescence processes are concerned, a key question is whether for cloud drops l inter < l c , the condition, which allows us to exploit the 'white noise' regime (corresponding to scales with St l > 1).

26
HYDROP experiment; we see that the values are little smaller than the turbulence-based speeds ϕ l . Using v R,l,drop as an estimate of the coalescence speed, one obtains: so that as long as the scale dependence of v R,l,drop is weak, we again have an l 1/2 law. It would be interesting to attempt to empirically check this l 1/2 law in clouds using aircraft drop data (for example, using the FSSP probe). Since table 1 shows that at l c , ϕ l = g −1/4 l 1/4 R ε 1/2 l is about double the value of v R,l,drop , we see that in the HYDROP experiment this sedimentation-based coalescence speed will be slower than the estimated turbulence-induced speed.
If we put the two models together-one for the dominant coalescence process in rain, and the other for the dominant process in clouds, we obtain the following picture for the production of rain from clouds. At first, the particles are small and numerous enough so that l inter < l c (see discussion in section 3.6) and the coalescence rate is primarily the gravity-dominated sedimentation mechanism. However, as the drops grow in size, l inter can increase more rapidly than l c . In regions where l inter > l c it is the enhanced turbulence-based coalescence speed which is dominant and this leads to the accelerated production of larger (rain-sized) drops.
It is worth briefly considering the effect of drop break up, which will be important for the larger drops and that we have ignored until now. Since breakup conserves the mass, it will not affect the equations for ρ. However, we must add another term to equation (15), which will be a linear (rather than quadratic) integral operator to account for breakup, this will add a new term to the equation for the number density (equation (21)) and a new contribution to the number variance flux dissipation (equation (43)); it will modify the effect of the microphysics on the flux dissipation (ψ diss ) but without changing its fundamental character. We will therefore still expect the number variance flux to exhibit a cascade from the large to small scales. However, breakup might make a significant modification to the timescale for the transfer τ l (equation (58)), which is currently estimated by the coalescence speed. At the moment, the coalescence speed is assumed to depend on the turbulence velocity differences across structures of size l c . Breakup mechanisms would presumably not be too sensitive to the turbulence intensity; they might add a contribution to the coalescence speed, which is independent of the turbulence. However, as long as the result has only weak scale (l) dependence, this could modify the detailed prediction equation (59) but not the basic l 1/2 law. Once again the key to the l 1/2 law is the existence of a (roughly) scale-independent number-variance flux transfer speed.

Coupled χ and ψ cascades and the l 1/2 law
If we eliminate the energy flux, we can express the mass density fluctuations ρ in terms of the number density fluctuations n via the simple relation: which has no explicit dependence on the energy flux ε or scale l. Although there is a weak scale dependence on the (intermittently varying) flux ratio χ 1/2 l /ψ 1/3 l and a weak DSD dependence of the mean relaxation length l R , equation (61) implies that the mass and number fluctuations are closely related. Equation (61) raises the question of the statistical couplings between ε, χ and ψ cascades. While the links between ε and χ, and between ε and ψ are not obvious, χ and ψ are indirectly linked through the number size density N , so that they cannot be statistically independent of each other. Indeed, the ratio χ /ψ has dimensions of mass 2 and it seems clear that at least at large enough scales the two should be related by the (ensemble, i.e. climatological) drop mass variance M 2 L ext : where L ext is the external scale of the cascades; M 2 L ext indicates the mass variance averaged over the drop distribution at the largest scale L ext . A 'strongly coupled' model would take χ l ≈ ψ l M 2 l , i.e. it would assume this relation to hold at all scales l, but this is likely to be too strong an assumption to be realistic. In section 4.5, we discuss the constraints imposed by 0th-and 1st-order moments n and ρ.
In order to exploit these statistical relations so as to make stochastic processes with the corresponding statistics, we may use the standard multifractal simulation techniques to simulate ε, χ and ψ and from them (by fractional integration to obtain the extra l 1/3 , l 1/2 scalings), the u, ρ and n fields. However, the n field determines the probability per unit volume of finding a drop; in other words, it should control a compound Poisson-multifractal process. In a future paper, we show how to make such processes which produce stochastic realizations respecting all the above statistics, and which determine implicitly the DSDs (see also Lovejoy and Schertzer (2006) for an early proposal). For the moment, we note that the particulate nature of rain in fact imposes a length scale equal to the mean interdrop distance: (n is the actual number density, not its fluctuation n). Therefore, wherever l inter > l c , l c must be replaced by l inter . This will occur at low rain rates and will modify the low rain rate statistics. This behaviour, in fact, provides a natural 'cutoff' mechanism for the transition from rain to no rain; as the rain rate becomes lower and lower, the inner scale of the cascade rapidly becomes larger and larger. We can study this using the compound cascade Poisson processes.

The empirical number density law
In order to empirically test the l 1/2 law for the number density, we used the same grid as for the mass density but with η = 0 (equation (4)), producing the spectra shown in figure 8(a). We see that the convergence to the low wavenumber theoretical k −2 behaviour (straight line) occurs at slightly smaller scales than for the k −5/3 behaviour of ρ since n is less variable (smoother) than ρ (see figure 3). In figure 8(b), we show the ratio of the ensemble spectra (all 19 triplets) for E ρ (k), E n (k); this shows that the number density field really is smoother (by about k 1/3 ) than the corresponding spectrum for the mass density. Because we have taken the ratio of the spectra, the y-axis is 'blown up' with respect to the previous; this is quite a sensitive indicator that the basic theory is roughly correct. Finally, in figure 8(c), we show the exponents of spectra of ρ η ; η = 0, 1 are the number and mass density spectra, whereas η = 1/6 and 7/6 are important for v R,n and v R,ρ , respectively, and η = 2 for the radar reflectivity factor.

The implications for the DSD
The model that has emerged combines three large-scale turbulence based fluxes-the energy, liquid water and number variance fluxes which are determined by the large-scale state of the atmosphere; the drop microphysics primarily determines the inner scales of ρ and n. Together, these three fluxes determine the highly intermittent u, n and ρ fields. Since n and ρ are respectively the 0th and 1st moments of the mass number density N (M), we see that the large scales constrain the microphysics via these two moments. In order to succinctly express the constraint, we can appeal to work done in the last few years exploring different ways of applying dimensional analysis to the analysis of DSDs (Lee et al 2004, Sempere-Torre et al 2000, Uijlenhoft et al 1999. The basic conclusion of these papers is that on purely empirical grounds a single-dimensional quantity is not enough to account for the diversity of DSDs so that at least two such quantities are needed. In Lee et al (2004) the use of two different moments to nondimensionalize the distribution is explored empirically. It is surprising that the authors never consider the mass moments order j = 0 and 1, instead they concentrate their attention on combinations of the moments of orders j/3 and j/2 or j/3 and ( j + 1)/3, where j is an integer 1.
In any case, nondimensionalization using the scale l fields n l , ρ l is straightforward-it can be done by inspection. First, the mass can be nondimensionalized by ρ l /n l = M l , which is the mean drop mass at scale l, and the cumulative drop number distribution N l (M/M l > s) of the nondimensional mass M/M l (i.e. the number density of drops with mass ratio M/M l exceeding the threshold s) can be nondimensionalized by n l which yields Pr l (M/M l > s), i.e. the cumulative probability distribution (accumulated from s to infinity rather than the more usual 0 to s). Combining both nondimensionalizations, the large-scale determination of n l , ρ l therefore implies that the nondimensional Pr l (M/M l > s) = p l (s) should be independent of n and ρ, so that p l (s) should be a universal (microphysics determined) probability distribution: only a function of the scale l and nondimensional drop mass s. This implies that the drop size density is N l (M) = n l (d p l (Mn/ρ)/dM) = (n 2 l /ρ l ) p l (s) with s = Mn l /ρ l , the consequences of which will be explored elsewhere.
In order to test this prediction, we present figure 9, which shows Pr l (M/M l > s) for the 5 storms with scales corresponding to the entire experimental volume (which was nearly the same for all the storms). From the figure, we see that the curves are so close as to be nearly indistinguishable; for Pr > 10 −3 , the standard deviation of the (log 10 ) residue of each curve compared to the overall mean is ±0.038 so that as predicted, there seems to be a universal nondimensional DSD. We also see that the extremes are perhaps not far from this algebraic form (itself a generic consequence of cascades).

Conclusions
Cloud and rain physics have been divorced from turbulence theory for too long. To date, the attempts to unify them at the most fundamental level have been arduous, and have proceeded from the small dissipation scales up, taking advantage of the general and elegant low Stokes number relation between the particle and wind velocities (Maxey 1987). In contrast, phenomenological turbulence cascades-which model the fluxes from large scales down to the dissipation scales-provide an attractive alternative. This is because they have proved to be not only highly accurate but even indispensable for understanding and modelling both rain and turbulence; they thus provide the natural framework for explicitly marrying the two.
In order to make progress, we relied on a new and unique source of data spanning the drop/turbulence scales: the HYDROP experiment. The latter provides the positions and masses of rain drops in a roughly 10 m 3 volume; the empirical HYDROP mass spectra show quite directly that rain drop (mass and number) densities are organized into statistically homogeneous (white noise) 'patches' whose overall liquid water content varies according to the classical turbulence (Corrsin-Obukhov) laws, while within each patch the distribution is homogeneous (the spectrum is white noise). In contrast, the classical theory of precipitation still used in interpreting radar echoes assumes that l c is the size of a radar pulse volume, which is typically 1 km or more.
Although rain drops have high Stokes numbers so that Maxey's relations are not valid, we argued that the patches were precisely of critical size l c where St l c = 1, so that it was plausible that Maxey's relations might hold for drop velocities coarse-grained over the drops within the homogeneous patches size l c (since for scales l > l c St l < 1). By assuming St l c = 1, we were able to empirically estimate the dissipation scales, the Stokes and the sedimentation numbers and the energy fluxes. By applying the coarse-grained Maxey relations, we formulated dynamical equations for both liquid water mass and particle number densities (ρ and n). With the exception of a mean vertical (sedimentation) velocity, the liquid water mass equation was identical to that of standard passive scalar advection and therefore predicted that rain should follow the Corrsin-Obukhov (l 1/3 and k −5/3 laws) of passive scalar advection. This was because both the microphysical processes and the sedimentation process conserve overall mass density. However, the equation for the particle number density was quite new: of the form l 1/2 (i.e. k −2 spectrum) the difference arising from the fact that while coalescence processes conserve total mass, they do not conserve total particle numbers. Although any microphysical process that defines a nearly scale-invariant coalescence speed ϕ leads to an l 1/2 number density law, we argued that for rain, the speed was proportional to ε 1/2 l , whereas for cloud drops the speed would depend rather on the differences in drop relaxation speeds. While the latter mechanism is presumably also present in rain, it could be dominated by the former energy flux-determined rate, whenever the patch size contained many particles (i.e. l c l inter ). Either mechanism allows the (weakly scale-dependent) coalescence speed to determine the rate of scale-by-scale transfer of number variance flux, whereas the speed of the corresponding mass-density variance flux depends on the strongly scale-dependent turbulence velocity. In both cases, the HYDROP data were used to verify that both the l 1/3 and l 1/2 laws were obeyed with reasonable accuracy. The weakest part of the argument was the use of Maxey's equations on the coarse-grained drop velocities. However, we only used certain symmetries of the resulting equations (notably the scale invariance and conservation of the variance fluxes by the nonlinear terms) and the conservation of mass by the coalescence processes, so that our conclusions about the l 1/3 and l 1/2 laws are likely to be robust.
Although in this paper, we only touch on it, an important consequence concerns the rain rate field R, which is the vertical component of the liquid water mass flux vector. Within the hydrostatic approximation (which is of acceptable accuracy in this context), we showed that R is the product of the mass density and the sum of the vertical velocity and mean relaxation speed. However, when the number density (n) field (which is correlated with the ρ field) is sufficiently close to zero, it will imply zero-rain rates since the probability of finding a drop becomes very low (a full study of this effect must be done using compound Poisson multifractal processes discussed below). In this way, a natural rain/no-rain mechanism emerges and we argue that it can quantitatively explain observations of rain rate spectra.
The conventional methods of modelling the evolution of raindrops (e.g. Khain and Pinsky 1995, Khairoutdinav and Kogan 1999, Srivastava and Passarelli 1980 give turbulence at most a minor (highly 'parameterized') role: the atmosphere is considered homogeneous and the spatial variability of the DSD arises primarily due to complex drop interactions. Our approach is in many ways the opposite: down to a small drop inertial scale (below which a kind of white noise 'drop chaos' exists), the drops are shown to be strongly coupled with the highly variable (multifractal) wind field. The coalescence speed is then determined by the turbulence velocity gradient at this decoupling scale. By considering the 0th and 1st moments of the DSDs (n and ρ), and showing how they are related to scale-by-scale conserved turbulent fluxes, we obtain strong (but implicit) constraints on the DSD, which we expect to be highly variable simply because it is controlled by turbulent cascades. Since these 0th-and 1st-order moments constrain the microphysics, we expect-on dimensional grounds (and there is literature on this type of DSD argument, e.g. Lee et al (2004)) that the dimensionless cumulative DSD as a function of the dimensionless drop mass should be a universal function of dimensionless mass. We empirically verified this on the HYDROP data.
The implications of this model can be taken further: if we are given a turbulent number density field, it can be used to control a compound Poisson process, which determines the locations of the drops. Once the positions are determined, the corresponding ρ field can then be used to attribute masses to the particles. In this way, we obtain a 3D model of particle positions and masses, which implicitly determines the highly variable DSD in the vicinity of each point and which respects the turbulence laws and is compatible with the microphysics (see Lovejoy and Schertzer (2006) for a preliminary model). Such models can be used to solve longstanding observers problems in the radar meteorology, satellite rain algorithms, rain-amount estimation and cloud physics.